TransferMatrix2DIsing

import numpy
import itertools
import scipy.sparse.linalg
import math
import heapq
import numpy as np
from irlb import *
import time


# ALGORITHM TO PERMUTE WITHOUT REPETIONS 
class unique_element:
    def __init__(self,value,occurrences):
        self.value = value
        self.occurrences = occurrences
        
def perm_unique(elements):
    eset=set(elements)
    listunique = [unique_element(i,elements.count(i)) for i in eset]
    u=len(elements)
    return perm_unique_helper(listunique,[0]*u,u-1)
        
def perm_unique_helper(listunique,result_list,d):
    if d < 0:
        yield tuple(result_list)
    else:
        for i in listunique:
            if i.occurrences > 0:
                result_list[d]=i.value
                i.occurrences-=1
                for g in  perm_unique_helper(listunique,result_list,d-1):
                    yield g
                i.occurrences+=1


class TransferMatrix:
######### working with the Transfer matrix 
    def __init__(self,N,J):       
        self.N = N
        self.J = J
    def compute(self,N,J):
        self.N = N
        self.J = J
        lista = list() # defining a list
        lista2 = list() # defining a list
        for i in range(1,self.N+1):
            lista.append(1) 
        
        for j in range(1,self.N):
            lista[self.N-j] = -1
            a = list(perm_unique(lista))# making the permutations with the class (permutation without repetition)
            lista2.append(a)
                       
        # removing extra brackets    
        lista2 = list(itertools.chain(*lista2))  
        #lista2.reverse()
        ################## just adding the missing states #################
        add1 = []
        add2 = []
        
        for i in range(1,self.N+1):
            add1.append(1)
            add2.append(-1)
        
        add1 = tuple(add1)
        add2 = tuple(add2)
        lista2.insert(0,tuple(add1))
        lista2.insert(len(lista2),tuple(add2))
        
        tst = []
        for i in range(0,len(lista2)):
            tst.append(lista2[i])
        
        lista2 = list(itertools.chain(tst)) # just removing the extra unneccessery brackets
        #################### 
        
        Nu = 2**(self.N)
        p = numpy.zeros(Nu*Nu).reshape((Nu,Nu))
        
        for i in range(0,Nu):
            for j in range(0,Nu):
                
                
                sum1 = 0                                
                
                for k in range(0,self.N):                        
                    sum1 = sum1 + self.J*(lista2[i][k]*lista2[j][k])
                
                for l in range(0,N):
                    if l==N-1:
                        k = 0
                    else:
                    	k = l + 1
                    sum1 = sum1 + self.J*(lista2[i][l]*lista2[i][k])/2
                    sum1 = sum1 + self.J*(lista2[j][l]*lista2[j][k])/2
		


                p[i][j] = np.float64(math.exp(sum1))
        return p

   
"""
N = 10
J = 1
T = TransferMatrix(N,J).compute(N,J)
vals = numpy.linalg.eigvals(T)
#print vals
#print T
#print "----- -------- ------ ----"
T = T.reshape(2**(N),2**(N))
#print T
nu = 2

tic1 = time.clock()
S = numpy.linalg.eigvals(T)
toc1 = time.clock()

tic2 = time.clock()
X = irlb(T,nu)
toc2 = time.clock()


print 'vreme za prvi:',toc1-tic1,'s'
print 'vreme za drugi:',toc2-tic2,'s'
print(X[1])
print(max(S))
print max(X[1])
"""

nu=2
for k in range(0,50):
    print k


    s = 0.2 + k*0.01
        
    c = ''.join([str('koji'),str(s),str('.txt')])
    f = open(c,'w')
            
    print "N/","/Largest EigenValues","/Free Energy", "-----J = ",s
    for i in range(0,13):
        N = i+1
        print 'gde sam:',N
        J = s
	tic1 = time.clock()
        T = TransferMatrix(N,J).compute(N,J)
	toc1 = time.clock()
	print 'vreme za transfer matricu:',toc1-tic1,'s'
        #vals, vecs = scipy.sparse.linalg.eigsh(T,k=2**(N)-1)
	#vals = numpy.linalg.eigvals(T)
	T = T.reshape(2**(N),2**(N))
	tic2 = time.clock()
	X = irlb(T,nu)
	toc2 = time.clock()
	print 'vreme za dijagonalizaciju:',toc2-tic2,'s'
        print(X[1]) 
	print 'samo jedna vrednost:',max(X[1])    
        eigenvalues = np.float64(math.log(max(X[1]))/N)	
        free_en = np.float64(-eigenvalues)
        print N,eigenvalues,free_en
        ww = ' '.join([str(N),str(free_en)])
        f.write(str(ww))
        f.write('\n')
    f.close()

###########################################