# -*- coding: cp1252 -*-
# -*- coding: utf-8 -*-

################################################
#   Módulo: appl_bij. Aplicaciones biyectivas. #
#   Autor: dserrano, bcurtu.                   #
#   Versión 1.0. Enero de 2009.                #
################################################

"""
Aplicaciones biyectivas.

Resumen:

	La clase ApplBij le permitirá extraer aleatoriamente aplicaciones
	biyectivas a partir de una matriz (booleana) de correspondencias.

Abstract:

	ApplBij class allows you to randomly extract bijective applications
	from a (boolean) matrix of correspondences.
"""

import os
import random
import time

# Comprobación de si el argumento T es una matriz (lista anidada).
def IsMatrix(T):

	"Comprobación de si el argumento T es una matriz (lista anidada)."

	if type(T)!=type([]) or len(T)==0:
		return False
	for i in range(len(T)):
		if type(T[i])!=type([]) or len(T[i])!=len(T[0]):
			return False
	return True

# Comprobación de si una matriz T es cuadrada.
def IsSquareMatrix(T):

	"Comprobación de si una matriz T es cuadrada."

	if len(T)!=len(T[0]):
		return False
	return True

# Comprobación de si una matriz T es booleana.
def IsBooleanMatrix(T):

	"Comprobación de si una matriz T es booleana."

	for i in range(len(T)):
		for j in range(len(T[0])):
			if not T[i][j] in (0, False, 1, True):
				return False
	return True

# Creación de una submatriz resultado de eliminar la fila m y la columna n en la matriz T.
def Submatrix(T, m, n):

	"Creación de una submatriz resultado de eliminar la fila m y la columna n en la matriz T."

	return [[T[i][j] for j in range(len(T[0])) if j!=n] for i in range(len(T)) if i!=m]

# Aplicaciones biyectivas.
class ApplBij:

	"""
	Clase ApplBij. Aplicaciones biyectivas.
	---------------------------------------

	Dada una matriz booleana y cuadrada A, llamada de correspondencias (o de in-
	clusiones / exclusiones), crearemos un objeto ApplBij de la forma:

		import appl_bij
		...
		apli=appl_bij.ApplBij(A, mode)	# Objeto ApplBij denominado apli.

	El parámetro mode puede valer 'SAFE' o 'UNSAFE' y es el modo de inicialización
	de la clase.

	.- Cuando un objeto ApplBij se inicializa en modo 'SAFE' se computan las apli-
	caciones y aplicaciones biyectivas existentes.
	.- Cuando un objeto ApplBij se inicializa en modo 'UNSAFE' no se computan las
	aplicaciones biyectivas.

	Para obtener una aplicación biyectiva al azar, simplemente escribiremos:

		loopSize=apli.GetLoopSize([alfa])	# (Esta sentencia sólo una vez).
		result=apli.GetRandomApplBij(loopSize)

	En el ejemplo anterior result es una lista donde se guardará la aplicación bi-
	yectiva hallada y alfa es la probabilidad (0.0<alfa<1.0) de éxito deseado en
	la búsqueda de dicha aplicación biyectiva (alfa=0.999 por defecto).
	El algoritmo de cálculo fija automáticamente el tamaño del bucle de búsqueda
	en función del valor alfa introducido. Si el modo de inicialización es 'UNSA-
	FE', no existe dependencia con alfa, por lo que es más claro llamar al método
	GetLoopSize() sin argumentos.

	Es posible que, aun existiendo aplicaciones biyectivas, no se encuentren den-
	tro del tamaño del bucle de búsqueda. En dicho caso se devolverá None.


	Un ejemplo de uso de la clase ApplBij: juego del Amigo Invisible.
	-----------------------------------------------------------------

	Supongamos que tenemos la siguiente matriz de correspondencias:

		#	Donantes de regalos
		#	 0  1  2  3  4  5
		Corr=[
			[0, 0, 1, 1, 1, 1],	# 0
			[0, 0, 1, 1, 1, 1],	# 1
			[1, 1, 0, 0, 1, 1],	# 2	Receptores de regalos
			[1, 1, 0, 0, 1, 1],	# 3
			[1, 1, 1, 1, 0, 0],	# 4
			[1, 1, 1, 1, 0, 0]	# 5
			]

	Si miramos la primera columna del ejemplo anterior vemos que el amigo 0 puede ha-
	cer un regalo a los amigos: 2, 3, 4 y 5, pero no a sí mismo ni a 1.
	...
	Si miramos la última columna del ejemplo anterior vemos que el amigo 5 puede ha-
	cer un regalo a los amigos: 0, 1, 2 y 3, pero no a sí mismo ni a 4.

	Escribamos ahora (por ejemplo):

		apli=appl_bij.ApplBij(Corr, 'SAFE')
		loopSize=apli.GetLoopSize()
		result=apli.GetRandomApplBij(loopSize)
		if result:
			print result

	Supongamos que el resultado de la sentencia print es:

		[4, 5, 0, 1, 2, 3]

	Lo anterior significa que:

	.- El amigo 0 deberá hacer un regalo al amigo 4.
	.- El amigo 1 deberá hacer un regalo al amigo 5.
	.- El amigo 2 deberá hacer un regalo al amigo 0.
	.- El amigo 3 deberá hacer un regalo al amigo 1.
	.- El amigo 4 deberá hacer un regalo al amigo 2.
	.- El amigo 5 deberá hacer un regalo al amigo 3.

	Es posible que la clase ApplBij no encuentre ninguna solución compatible con
	la matriz de correspondencias. Esto se debe a:

	.- En modo 'SAFE': porque no existe tal solución (con la probabilidad alfa
	especificada).
	.- En modo 'UNSAFE': porque no existe tal solución o porque, existiendo, la
	clase ApplBij no la ha encontrado con el tamaño especificado del bucle de
	búsqueda (loopSize).

	En dicho caso el método GetRandomApplBij() devolverá None.

	Por esta razón es siempre necesaria la comparación:

		if result:

	antes de usar la lista almacenada en result.


	Detalles técnicos sobre la clase ApplBij.
	-----------------------------------------

	La clase ApplBij sigue una estrategia de búsqueda de aplicaciones biyectivas
	que es función de la densidad de la matriz de correspondencias. Se define la
	densidad de dicha matriz como el cociente entre el número de aplicaciones y
	el número de permutaciones existentes, siendo:

	.- Número de aplicaciones: es el producto, extendido a todas las columnas, de
	la suma del número de Trues en cada columna.
	.- Número de permutaciones: para una matriz de tamaño n es Factorial(n) o n!.

	De una forma menos rigurosa, la densidad puede entenderse aproximadamente co-
	mo la relación del número de Trues entre el número de Falses contenidos en la
	matriz de correspondencias. Así, densidades mayores de 1.0 indican una preva-
	lencia de Trues sobre Falses. La prevalencia se invierte para valores de den-
	sidad inferiores a 1.0.

	La clase ApplBij es muy eficiente tanto en el caso de matrices muy densas (
	como generalmente se da en el juego del Amigo Invisible) como en el caso de
	matrices poco densas. Sin embargo es poco eficiente cuando la densidad se a-
	proxima a 1.0. En este caso:

	.- El número de aplicaciones biyectivas (soluciones) es muy pequeño en rela-
	ción al número tanto de aplicaciones como de permutaciones, por lo que, si
	la matriz de correspondencias es suficientemente grande, deberemos elegir un
	tamaño del bucle de búsqueda (loopSize) grandísimo si queremos asegurarnos
	de encontrar alguna solución. Por ello:
	.- Si el modo de inicialización es 'SAFE', la llamada al método GetLoopSize()
	devolverá un tamaño del bucle de búsqueda muy grande. En caso de existir so-
	lución, ésta se hallará, pero a costa de un tiempo de cálculo que puede lle-
	gar a ser inaceptablemente largo.
	.- Si el modo de inicialización es 'UNSAFE', la llamada al método GetLoopSi-
	ze() devolverá un tamaño del bucle de búsqueda muchas veces demasiado peque-
	ño, por lo que será probabilísticamente muy difícil encontrar una solución.
	Es claro que el usuario siempre tiene la libertad de usar el valor loopSize
	que le plazca (sin necesidad de llamar al método GetLoopSize()), pero eso no
	acortará el tiempo de cálculo necesario para encontrar una solución.

	Un ejemplo de matriz con densidad 1.0 es aquélla en que sus valores son True
	salvo debajo de la diagonal principal, donde todos son False. Esta matriz tie-
	ne las siguientes propiedades:
		.- El número de aplicaciones coincide con el número de permutaciones y
		es: Factorial(n), siendo n el tamaño de la matriz.
		.- Sólo existe una aplicación biyectiva. Ésta es: [0, 1, 2, ... , n-1].
	Observe el lector que si por ejemplo n=10, se ha de encontrar al azar la úni-
	ca solución existente entre 10!=3628800 aplicaciones. Para encontrarla con
	una probabilidad 0.999 necesitaremos un bucle de búsqueda de tamaño al menos
	¡25066860!


	Modos de inicialización de la clase ApplBij.
	--------------------------------------------

	Cuando inicializamos la clase ApplBij en modo 'SAFE', lo primero que hace és-
	ta es calcular el número de aplicaciones biyectivas (soluciones) por recur-
	sión. Esto nos brinda la seguridad de que, si existe, obtendremos una solución
	al llamar al método GetRandomApplBij() (con la probabilidad alfa especificada).
	Tiene el grave inconveniente de que si el tamaño de la matriz es grande (a
	partir de 8 o así) el tiempo de cálculo se volverá inaceptablemente largo.
	Se desaconseja por tanto operar con matrices de tamaño superior a 8 con este
	modo de inicialización.

	En el modo de inicialización 'UNSAFE' no se computan las aplicaciones biyecti-
	vas existentes, por lo que su funcionamiento es mucho más rápido. Tiene no obs-
	tante dos inconvenientes, ambos derivados de que no sabemos cuántas aplicacio-
	nes biyectivas (soluciones) hay:
		.- Puede que no exista ninguna solución y la estemos buscando en balde.
		.- El método GetLoopSize() devuelve un tamaño del bucle de búsqueda cons-
		tante (10000). Puede que este valor sea insuficiente para matrices grandes
		de densidad cercana a 1.0. En dicho caso es posible que no se encuentre
		ninguna solución (aun habiéndola).
	"""

	# Comprobación de si una aplicación es biyectiva.
	def __IsBijective(self, x):

		for i in range(len(x)):
			if x[i]==self.NullAppl[i]:
				return False

			for j in range(i+1, len(x)):
				if x[i]==x[j]:
					return False

		return True

	# Comprobación de si una permutación es compatible con la matriz de correspondencias.
	def __IsCompatible(self, x):

		for j in range(len(self.Table[0])):
			if not self.Table[x[j]][j]:
				return False
		return True

	# Verificación de la existencia de aplicaciones biyectivas.
	# (Si el método devuelve False se garantiza la no existencia de aplicaciones biyectivas).
	# (Si el método devuelve True no se garantiza la existencia de aplicaciones biyectivas).
	def __AreThereApplBij(self):

		for row in self.Table:
			if not (1 in row or True in row):
				return False
		return True

	# Cálculo del número de aplicaciones biyectivas.
	@staticmethod
	def __ComputeNumApplBij(T, numApplBij):

		if len(T)==1:
			if T[0][0]:
				numApplBij[0]+=1
		else:
			for i in range(len(T)):
				if T[i][0]:
					ApplBij.__ComputeNumApplBij(Submatrix(T, i, 0), numApplBij)

	# Inicialización del objeto ApplBij.
	def __Initialize(self):

		for j in range(len(self.Table[0])):
			for i in range(len(self.Table)):
				if self.Table[i][j]:
					self.NumCorr[j]+=1
			self.NumAppl*=self.NumCorr[j]
			self.NumPerm*=j+1
			self.density*=float(self.NumCorr[j])/(j+1)

	# Obtención aleatoria de una aplicación.
	def __GetRandomAppl(self):

		Appl=[]

		for j in range(len(self.Table[0])):
			count=0

			if self.NumCorr[j]>0:
				draw=random.randint(1, self.NumCorr[j])
			else:
				draw=1

			for i in range(len(self.Table)):
				if self.Table[i][j]:
					count+=1
				if count==draw:
					Appl.append(i)
					break
			else:
				Appl.append(i+1)

		return Appl

	# Constructor.
	def __init__(self, T, mode):

		"""
		Constructor de la clase ApplBij.

		T:	matriz de correspondencias. Es una matriz (en forma de lista anidada)
			cuadrada y booleana.
		mode:	modo de inicialización de la clase. Puede ser 'SAFE' o 'UNSAFE'.

		Se lanzará una excepción ValueError si:
		.- T no es una matriz, no es cuadrada o no es booleana.
		.- mode no es ni 'SAFE' ni 'UNSAFE'.
		"""

		if not (IsMatrix(T) and IsSquareMatrix(T) and IsBooleanMatrix(T)) \
		or not mode in ('SAFE', 'UNSAFE'):
			raise ValueError, 'Base class: ApplBij. Invalid argument(s).'

		self.Table=T
		self.mode=mode
		self.NullAppl=[len(self.Table) for i in range(len(self.Table[0]))]
		self.NumCorr=[0 for i in range(len(self.Table[0]))]
		self.NumAppl=self.NumPerm=1
		self.NumApplBij=[0]
		self.density=1.0

		self.__Initialize()

		if self.__AreThereApplBij():
			if self.mode=='SAFE':
				ApplBij.__ComputeNumApplBij(self.Table, self.NumApplBij)
			else:
				self.NumApplBij=[1]

		random.seed()

	# Obtención del número máximo de aplicaciones al azar que se han de extraer
	# para obtener, con probabilidad de al menos alfa, una aplicación biyectiva.
	def GetLoopSize(self, alfa=0.999):

		"""
		Obtención del número máximo de aplicaciones al azar que se han de extraer
		para obtener, con probabilidad de al menos alfa, una aplicación biyectiva.
		(Esta probabilidad no se puede garantizar si la clase ApplBij se ha ini-
		cializado en modo 'UNSAFE').

		Se lanzará una excepción ValueError si alfa<=0.0 o alfa>=1.0.
		"""

		if alfa<=0.0 or alfa>=1.0:
			raise ValueError, 'Base class: ApplBij. Invalid argument(s).'

		if self.NumApplBij[0]<=0 or self.density<=0.0:
			return 0

		if self.mode=='UNSAFE':
			return 10000

		if self.density>1.0:
			p=float(self.NumApplBij[0])/self.NumPerm
		else:
			p=float(self.NumApplBij[0])/self.NumAppl

		i=0
		s=0.0

		while s<alfa:
			i+=1
			s+=((1-p)**(i-1))*p

		return i

	# Obtención aleatoria de una aplicación biyectiva.
	# (Se usará un bucle de búsqueda de tamaño especificado LoopSize).
	def GetRandomApplBij(self, LoopSize):

		"""
		Obtención aleatoria de una aplicación biyectiva (en forma de lista).
		Se usará un bucle de búsqueda de tamaño especificado LoopSize.
		Si no se encuentra ninguna aplicación biyectiva se devolverá None.
		"""

		if self.density>1.0:
			Perm=[i for i in range(len(self.Table))]
			for i in range(LoopSize):
				random.shuffle(Perm)
				if self.__IsCompatible(Perm):
					self.stat=(len(self.Table), self.mode, self.density, LoopSize, i+1, True)
					return Perm
		else:
			for i in range(LoopSize):
				Appl=self.__GetRandomAppl()
				if self.__IsBijective(Appl):
					self.stat=(len(self.Table), self.mode, self.density, LoopSize, i+1, True)
					return Appl

		self.stat=(len(self.Table), self.mode, self.density, LoopSize, LoopSize, False)
		return None

	# Obtención de la estadística del sorteo (en forma de tupla).
	def GetDrawStat(self):

		"""
		Obtención de la estadística del sorteo (en forma de tupla).
		
		La tupla contiene la siguiente información:
		
		.- Tamaño de la matriz de correspondencias.
		.- Modo de inicialización de la clase ('SAFE' / 'UNSAFE').
		.- Densidad de la matriz de correspondencias.
		.- Tamaño del bucle de búsqueda usado en el sorteo.
		.- Número real de iteraciones habidas en el sorteo.
		.- ¿Sorteo exitoso? (True / False).
		
		Para obtener la estadística de un sorteo, use este método
		después de usar el método GetRandomApplBij().
		"""

		return self.stat

# Test de la clase ApplBij.
def Test():

	"Test de la clase ApplBij."

	# Creación de la matriz de correspondencias para un sorteo del Amigo Invisible
	# entre parejas. (Uno no puede regalarse ni a sí mismo ni a su pareja).
	def CreateCorr(size):

		Corr=[[True for j in range(2*size)] for i in range(2*size)]

		for i in range(len(Corr)/2):
			for j in range(2*i, 2*i+2):
				Corr[2*i+1][j]=Corr[2*i][j]=False

		return Corr

	# Creamos la matriz A donde alojamos una correspondencia.
	A=CreateCorr(6)

	# Imprimimos la matriz A.
	print 'Sorteos del Amigo invisible para %d parejas.' % (len(A)/2)
	print
	for i in range(len(A)):
		for j in range(len(A[0])):
			print '%-5s' % str(A[i][j]),
		print

	# Creamos una aplicación biyectiva a partir de la matriz A.
	# (La inicializamos en modo 'UNSAFE').
	appl=ApplBij(A, 'UNSAFE')

	# Determinamos el tamaño del bucle de búsqueda.
	loopSize=appl.GetLoopSize()

	# Buscamos e imprimimos los sorteos exitosos.
	while(True):
		result=appl.GetRandomApplBij(loopSize)
		print 'Resultado del sorteo:'
		if result:
			print result
		print 'Estadística del sorteo:'
		print appl.GetDrawStat()
		time.sleep(1)

# Ejecución del test de la clase ApplBij.
if __name__=='__main__':
	help(ApplBij)
	Test()