#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Mon Nov 13 13:19:11 2017

@author: laro
"""

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Thu Nov  9 21:35:40 2017

@author: laro
"""

"""
Created on Wed Nov  8 13:36:59 2017

@author: laro
"""

#!/usr/bin/python -O 
# -*- encoding: utf-8 -*-
import numpy as np
from numpy import array
from scipy import * 
from random import *
#from pylab import *
import time
from numpy import cumsum
#import matplotlib.pyplot as plt

from sklearn import datasets, linear_model
from sklearn.metrics import mean_squared_error, r2_score
# script general para hacer una corrida  a un set de parámetros,
# beta, tamaño de la red,

#beta = 1/T

def En(S):
    energia = 0
    for i in range(0,Lx):
        iright = irightl[i]#mod(i+1,Lx) 
        ileft = ileftl[i]#mod(i+Lx-1,Lx) 
        for j in range (0, Ly):
            jup = jupl[j] #mod(j+1,Ly) 
            jdown = jdownl[j] #mod(j+Ly-1,Ly) 
            energia += - S[i,j]*(S[i,jup] + S[i,jdown] + S[iright,j] + S[ileft,j])

    return energia/2.

def ising2Dpaso(S,beta):

    DEc = 0
    DM = 0

    for i in range (0,Lx):
        iright = irightl[i];#mod(i+1,Lx)
        ileft = ileftl[i];#mod(i+Lx-1,Lx)
        for j in range (0, Ly):
            jup = jupl[j];#mod(j+1,Ly)
            jdown = jdownl[j];#mod(j+Ly-1,Ly) 
            spin_old = S[i,j]
            spin_new = -S[i,j]

            DE = 2.0*spin_old*(S[iright,j] + S[ileft,j] + S[i,jup] + S[i,jdown] )
            if( DE < 0):
        	    
                S[i,j] = spin_new
                DEc += DE
                DM += (spin_new - spin_old)
            else:

                p = random()
                expbetaE = exp(-beta*DE)
                if (p < expbetaE):
                    S[i,j] = spin_new
                    DEc += DE
                    DM += (spin_new - spin_old)
                else:
                    g=1
                		  
    return DEc, DM, S


# Tamaño en x,
Lx=32
# Tamaño en y,
Ly=32
irightl=roll(arange(0,Lx),1)
ileftl=roll(arange(0,Lx),-1)
jupl=roll(arange(0,Ly),1)
jdownl=roll(arange(0,Ly),-1)

T=2
beta=1.0/T





def isin(T,Sij,npre,npasos):      
    #pretermalizo
    for n in range(1,npre):
        DE,DM,Sij=ising2Dpaso(Sij,beta)
    energia=zeros(npasos)
    magnet=zeros(npasos)            
    energia[0] = En(Sij)/(Lx*Ly)
    magnet[0] = sum(Sij)/(Lx*Ly)
    for n in range(1,npasos):
        DE,dM,Sij = ising2Dpaso(Sij,beta)
        energia[n] = energia[n-1] + DE/(Lx*Ly)
        magnet[n] = magnet[n-1] + dM/(Lx*Ly)
    return energia,magnet,Sij


#busco transicion de fase

meanE=[]
meanE2=[]
meanM2=[]
meanM=[]
stdE=[]
stdM=[]
Tvect=[]
npre=2000

npasos=8000

# propongo un estado inicial al azar
# Sij es una matriz de 1 y -1 indicando las dos proyecciones de
# espin
# AL estar afuera del loop aleatorizo solo una vez... como arranco desde temp bajas propongo un estado inicial magnetizado
m=0.95
Sij = zeros((Lx,Ly))
for i in range(0,Lx):
    for j in range (0,Ly):
        if(random() > m):
            Sij[i,j] = 1 
        else:
            Sij[i,j] = -1
            
for i in range(1,51):
    T=i/10
    Tvect.append(T)
    print("T={}".format(T))
    beta=1.0/T
    for n in range(1,npre):
        DE,DM,Sij=ising2Dpaso(Sij,beta)
       
    energia,magnet,Sij=isin(T,Sij,npre,npasos)
    meanE.append(mean(energia))
    meanE2.append(mean(array(energia)*array(energia)))
    meanM.append(mean(abs(magnet)))
    meanM2.append(mean(array(abs(magnet))*array(abs(magnet))))
    stdE.append(std(energia))
    stdM.append(std(magnet))
    
np.save("mE.npy",meanE)
np.save("mE2.npy",meanE2)
np.save("mM.npy",meanM)
np.save("mM2.npy",meanM2)
np.save("stdE.npy",stdE)
np.save("stdM.npy",stdM)
np.save("Tvect.npy",Tvect)