Surface growth simulation - Cellular Automata Domany-kinzel

#!/usr/bin/python
# -*- coding: utf-8 -*-
"""
Surface growth simulation using Cellular Automata Domany-kinzel
Autores:
  THIAGO LINHARES BRANT REIS
  VANDERCI FERNANDES ARRUDA
  Wagner Cipriano
  Email: wagnao-gmail
Disciplina: Princípios de Modelagem Matemática. CefetMG.
Prof: Thiago Mattos

Options:
 -L #  Qtde de Sítios           Def: 10**4
 -T #  Unidades de tempo        Def: 10**5
"""

#Imports::
from os import system
from time import time, sleep, strftime
import getopt, sys
from random import randint, random
import math
import numpy as np
from scipy.stats import linregress

#Global vars
global gFileName, ExitFileNameFR, ExitFileNameEC, Funcs
gFileName = 'Exit'
Funcs = 'FR+EC'
#FR: Flutuacao Rugosidade
#EC: Expoente de crescimento

def __Pot2(i, f):
    """Retorna as potencias de 2 em um intervalo [i, f]"""
    v = i
    idx = 1
    Result = []
    while(v < f):
        Result.append(v)
        v = 2**idx
        idx+=1
    Result.append(f)
    return Result
#

def __NormalizeDimList(List, Tam, Value=None):
    """
      Insere elemento fixo na lista ate que completar qtde de elementos necessária(Tam)
      O lemento inserido e o param Value ou o ultimo elemento da lista caso Value = Null
    """
    if(not List):
        raise Exception('List invalid')
    elif(Value is None):
        Value = List[-1]
    elif(len(List) >= Tam):
        print('List len Ok. List: %s, Tam: %s' %(len(List), Tam))
    while (len(List) < Tam):
        List.append(Value)
#

def __GetNewCellValueByNeighbors(ln, rn, P1, P2):
    """
      Defini o valor de um sítio baseado na sua vizinhaça no instante de tempo anterior
      e nos valores de P1 e P2
      ln:  left neighbor
      lr:  rigth neighbor
      P1:  Probability P1
      P2:  Probability P2
    """
    RNumber = random() #Pseudo random numbers!
    NeigValue = '%s%s' %(ln, rn)
    if(NeigValue == '00'):                #00
        return 0
    elif(NeigValue == '11'):              #11
        return [0,1] [RNumber <= P2]
    else:                                 #'01 or 10'
        return [0,1] [RNumber <= P1]
#

def __RMS(Height, L, mean=None):
    """
      Calcula a rugosidade da interface
      RMS (Root Mean Square)
      #Roughness of the interface at time instant t
      Height: Altura da interface
    """
    if(mean is None): mean = Height.mean()
    WLt = math.sqrt( (((Height - mean) ** 2).sum()) / L )
    return WLt
#

def __GetInitialCondition(L, dorandom=True, key=3):
    """
      Gera a condição inicial
      L: Qtde de Sítios
    """
    if(dorandom):
        first_row = [randint(0,1) for i in range(L)]
    else:
        first_row = [[0,1] [i%key ==0] for i in range(L)]
    return np.array(first_row)
#

def __SaveToFile(Text, Var='', lFileName=None, Ext='log', Mode='a+'):
    """
     Save result to txt file
    """
    global gFileName
    if(lFileName): FileName = lFileName
    else:          FileName = gFileName
    with open('./%s.%s' %(FileName, Ext), 'a+', encoding='utf-8') as f:
        f.write('%s %s\n' %(Text, str(Var)))
#

def RegLinGetEq(x, y):
    #Regressao linear: Encontrar a equação da reta a partir de um conjunto de pontos: vetor x, y
    #http://webpages.fc.ul.pt/~aeferreira/python/aulas/11_scipy_regression.html
    #y = mx + b
    x = [math.log(i,10) for i in x]
    y = [math.log(j,10) for j in y]
    m, b, R, p, SEm = linregress(x, y)
    #    m = sum(a*b for (a,b) in zip(x,y)) - sum(x) * sum(y) / len(x)
    #    m /= sum(a**2 for a in x) - (sum(x)**2)/len(x)
    #    b = sum(y)/len(y) - m * sum(x)/len(x)
    print('y = %sx + %s' %(m, b))
    return m
#

def SurfaceGrowthSim(InitCond, L, T, P1, P2, PtsTime):
    """
      Simulate Surface Growth with ACDK
      InitCond: array of inital condition
      T:        Unidades de tempo
      P1:       Probability P1
      P2:       Probability P2
    """
    MeanHeightByTime = []  #hbar
    RMSByTime = []  #Root Mean Square
    Height = InitCond.copy()
    LastOne = InitCond

    #Mean height initial condition
    MeanHeightByTime.append(Height.sum() / L)
    RMSByTime.append(__RMS(Height, L))

    # Determine each new row (t) based on the last old data (t-1).
    ti = time()
    for t in range(T-1): #Passos de tempo
        NewOne = [__GetNewCellValueByNeighbors(LastOne[ (j - 1) %L ],  LastOne[ (j + 1) %L ], P1, P2)
                   for j in range(L)]
        NewOne = np.array(NewOne)
        ActivSitGen = NewOne.sum()
        #Add height
        Height = Height + NewOne

        #Mean height
        h_bar = Height.sum() / L             #(Pag 20, Form 2.1)
        #Roughness of the interface at time instant t   --> RMS (Root Mean Square)
        w = __RMS(Height, L, h_bar)

        #Armazenar dados - filtro: somente armazena alguns pontos (PtsTime)
        if((not PtsTime or (t in PtsTime)) or (ActivSitGen <= 0)):
            MeanHeightByTime.append(h_bar)
            RMSByTime.append(w)

        #update last one to use in the next generation aa
        LastOne = NewOne.copy()
        print('t: %s   h_bar: %s   w: %s   AS: %s' %(str(t).ljust(6), str(round(h_bar, 3)).ljust(12), str(round(w, 3)).ljust(12), str(ActivSitGen).ljust(6)))

        #Avaliar se a minha geração já chegou no regime de saturação
        if(ActivSitGen <= 0):
            print('  Regime de saturação alcançado na geração %s' %(t+1+1));
            #system('speaker-test -t sine -f 500 -l 1 -p 4');
            break;
        #
    print('P1: %s   Time i: %s' %(str(P1).ljust(7), time() - ti))
    return Height, MeanHeightByTime, RMSByTime
#

def SurfaceGrowthSimWithProbParams(L, T, P1List, P2List, InitCond, PtsTime):
    """
      Simulate Surface Growth using ACDK
      Results:  array of results
      T:        Unidades de tempo
      P1:       Probability P1
      P2:       Probability P2
    """
    Results = []
    for P1 in P1List:
        for P2 in P2List:
            print('\n\n\nP2: %s, P1: %s' %(P2, P1))
            Height, MeanHeightByTime, RMSByTime = SurfaceGrowthSim(InitCond, L, T, P1, P2, PtsTime)
            Results.append({'T':T, 'L':L, 'P1':P1, 'P2':P2, 'Height': Height, 'MeanHeightByTime':MeanHeightByTime, 'RMSByTime': RMSByTime})

            #Save Result Data
            __SaveToFile(Text='\n#P1: %s, P2: %s     %s' %(P1, P2, strftime('%Y-%m-%d %H:%M:%S')))
            __SaveToFile(Text='  Final Height:       ', Var=','.join(str(i) for i in Height))
            __SaveToFile(Text='  MeanHeightByTime:   ', Var=MeanHeightByTime)
            __SaveToFile(Text='  RMSByTime:          ', Var=RMSByTime)
            __SaveToFile(Text='  Generations len:    ', Var=len(RMSByTime))

    return Results
#

def GetGraphicFlutRugos(L, T, InitCond):
    """Simulando a Evolução da flutuação da rugosidade"""
    #CONSTs
    P1 = [0.2, 0.5, 0.59, 0.595, 0.6, 0.8]
    P2 = [0.95]

    #Itens Legenda
    ListLeg = []
    for P1i in P1:
        ListLeg.append('P1 = %s' %(P1i))

    #Passos de tempo a armazenar os pontos
    PtsTime = __Pot2(i=1, f=T)
    #Save to file
    __SaveToFile(Text='\n*Evolução da flutuação da rugosidade:\n  L:       %s\n  T:       %s\n  P1:      %s\n  P2:      %s\n  PtsTime: %s' %(L, T, P1, P2, str(PtsTime)))
    #simulate surface growth
    ListItems = SurfaceGrowthSimWithProbParams(L, T, P1, P2, InitCond, PtsTime)
    RespItems = []
    for Item in ListItems:
        __NormalizeDimList(List=Item['RMSByTime'], Tam=len(PtsTime))
        RespItems.append( Item['RMSByTime'] )

    #Armazenar dados geração grafico
    Ext='py';
    __SaveToFile(Text='# -*- coding: utf-8 -*-\n#DateTime: %s\n\nX = ' %(strftime('%Y-%m-%d %H:%M:%S')), Var=PtsTime, lFileName=ExitFileNameFR, Ext=Ext)
    __SaveToFile(Text='Y = ', Var=RespItems, lFileName=ExitFileNameFR, Ext=Ext)
    __SaveToFile(Text='ListLeg = ', Var=ListLeg, lFileName=ExitFileNameFR, Ext=Ext)
    __SaveToFile(Text="\nxlabel = 'time step'\nylabel = r'Rugosidade $\delta \omega(L, t)$'\ntitle = 'Evolução da flutuação da rugosidade   L=%s  T=%s'\nScale = 'loglog';\nlinewidth = 0.5;\nLegLoc = 2;\n" %(L, T), lFileName=ExitFileNameFR, Ext=Ext)

    #Validações dos dados
    print ('\n\nPtsTime: ', len(PtsTime));
    for i in range(len( RespItems)):
        print (i, len(RespItems[i]))
        if(len(RespItems[i]) != len(PtsTime)):
            raise Exception('Erro validacao dados')
#

def GetGraphicExpCresc(L, T, InitCond):
    #Evolução do expoente de crescimento

    #CONSTs
    Ext = 'py';
    P2 = [0.1, 0.3, 0.5, 0.7, 0.9]
    #P1i = np.linspace(0.61, 0.85, 20)                  #20pts   0.61 - 0.85
    #P1 = np.append(P1i, np.linspace(0.86, 1, 10))      #7pts   0.86 - 1.0
    #P1List = list(P1); del P1i, P1;
    ##P1 = [0.50, 0.55, 0.60, 0.65, 0.70, 0.75, 0.80, 0.85, 0.90, 0.95, 1.0]
    ##P1 = [0.61, 0.62, 0.63, 0.64, 0.65, 0.69, 0.70, 0.71, 0.79, 0.8, 0.81, 0.89, 0.9, 0.91, 0.99, 1.0]
    #novo modelo para gerar os valores de P1, conforme grafico da referencia (Allbens)
    P1ListDc = { '0.1': [0.780, 0.785, 0.790, 0.795,   0.800,   0.805, 0.810, 0.815, 0.820,   0.830, 0.840,  0.850,   0.870, 0.890, 0.910, 0.930],
                 '0.3': [0.760, 0.765, 0.770, 0.775,   0.780,   0.785, 0.790, 0.795, 0.800,   0.810, 0.820,  0.830,   0.850, 0.870, 0.890, 0.910],
                 '0.5': [0.730, 0.735, 0.740, 0.745,   0.750,   0.755, 0.760, 0.765, 0.770,   0.780, 0.790,  0.800,   0.820, 0.840, 0.860, 0.880],
                 '0.7': [0.685, 0.690, 0.695, 0.700,   0.705,   0.710, 0.715, 0.720, 0.725,   0.735, 0.745,  0.755,   0.775, 0.795, 0.815, 0.835],
                 '0.9': [0.610, 0.615, 0.620, 0.625,   0.630,   0.635, 0.640, 0.645, 0.650,   0.660, 0.670,  0.680,   0.700, 0.720, 0.740, 0.760],
               }
    __SaveToFile(Text='# -*- coding: utf-8 -*-\n#DateTime: %s\n\nP1List = %s\n' %(strftime('%Y-%m-%d %H:%M:%S'), P1ListDc), lFileName=ExitFileNameEC, Ext=Ext)

    #Passos de tempo a armazenar os pontos
    PtsTime = None
    PtsTime = __Pot2(i=1, f=T)

    #Save to file
    __SaveToFile(Text='\n\n*Evolução do expoente de crescimento:\n  L:       %s\n  T:       %s\n  P1:      %s\n  P2:      %s\n  PtsTime: %s' %(L, T, P1ListDc, P2, str(PtsTime)))

    #Begin simulation
    ListLeg = []
    RespItens = []
    for ItemP2 in P2:
        BetaList = []
        P1List = P1ListDc[str(ItemP2)]
        for ItemP1 in P1List:
            #simulate surface growth
            ListItems = SurfaceGrowthSimWithProbParams(L, T, [ItemP1], [ItemP2], InitCond, PtsTime)
            w = ListItems[0]['RMSByTime']

            #Encontrar o coeficiente angular de cada um dos conjuntos de dados
            vx = [PtsTime, range(L)] [not PtsTime]
            Beta = RegLinGetEq(vx[:len(w)], w)
            BetaList.append(Beta)
            print ('\n\n\nP2=%s P1=%s\nt: %s\nw: %s\nB: %s\n\n' %(str(ItemP2).ljust(5), str(ItemP1).ljust(5), vx, w, Beta));

        #Save data to file (para conferencia)
        LogTxt = 'P2=%s\nBetaList = %s;\n' %(ItemP2, str(BetaList))
        __SaveToFile(Text=LogTxt, lFileName=ExitFileNameEC, Ext=Ext)

        #Armazenando os dados: ItemP2, P1List, BetaList
        ListLeg.append('P2 = %s' %(ItemP2))
        RespItens.append(BetaList)

    #Armazenar dados geração grafico
    __SaveToFile(Text='\n\n\nX =', Var=P1List, lFileName=ExitFileNameEC, Ext=Ext)
    __SaveToFile(Text='Y =', Var=RespItens, lFileName=ExitFileNameEC, Ext=Ext)
    __SaveToFile(Text='ListLeg = ', Var=ListLeg, lFileName=ExitFileNameEC, Ext=Ext)
    __SaveToFile(Text="\nxlabel = 'P1';\nylabel = r' $ \\beta $';\ntitle = 'Evolução expoente de crescimento   L=%s  T=%s'\nScale = 'linear';\nlinewidth = 0.1;\nLegLoc = 1;\n" %(L, T), lFileName=ExitFileNameEC, Ext=Ext)
    ##__SaveToFile(Text='ListItems = ', Var=ListItems, lFileName=lFileName, Ext=Ext)
#


def main():
    Help = """
      Cellular Automata Domany-kinzel
        Options:
          -L   Qtde de Sítios           Def: 10**4
          -T   Unidades de tempo        Def: 10**5
          -F   Funcs                    Def: FR+EC
                                             FR: Flutuacao Rugosidade
                                             EC: Expoente de crescimento
    """
    global Funcs
    #PARAMETROS
    L = 10**4 #Qtde de Sítios
    T = 10**5 #Unidades de tempo
    IniCondRand = True

    #Get params (comand line options)
    opts,args = getopt.getopt(sys.argv[1:],'L:T:F:h')
    for key,val in opts:
        if(key == '-L'):    L = int(val)
        elif(key == '-T'):  T = int(val)
        elif(key == '-F'):  Funcs = str(val).upper()
        elif(key == '-h'):
            print(Help);
            sys.exit(1);
    print('L:     %s\nT:     %s\nICR: %s\nFuncs: %s' %(L, T, IniCondRand, Funcs)); sleep(1);

    #Condição inicial é a mesma para todas as simulações da sessão ativa
    InitCond = __GetInitialCondition(L, IniCondRand)

    #Gerar grafico 01 - Fig 3.3
    if(Funcs.find('FR') != -1):
        GetGraphicFlutRugos(L, T, InitCond)

    #Gerar grafico 02 - Fig 3.4
    if(Funcs.find('EC') != -1):
        GetGraphicExpCresc(L, T, InitCond)
#

if(__name__ == '__main__'):
    gFileName = 'Exit-%s' %(strftime('%Y-%m-%d-%Hh'))
    ExitFileNameFR = 'DadosFlutRugos_%s' %(strftime('%Y%m%d%H%M%S'))
    ExitFileNameEC = 'DadosExpCresc_%s'  %(strftime('%Y%m%d%H%M%S'))

    ini = time()
    __SaveToFile(Text='--- '*20 + '\nbegin: %s ' %(strftime('%Y-%m-%d %H:%M:%S')) )

    main();

    __SaveToFile(Text='\n\nend: %s\nDur: %s min\n\n\n\n\n' %(strftime('%Y-%m-%d %H:%M:%S'), (time()-ini)/60  ))

    #Gerar os graficos
    if(Funcs.find('FR') != -1):
        system('python3 Plots.py -f %s' %(ExitFileNameFR));
    if(Funcs.find('EC') != -1):
        system('python3 Plots.py -f %s' %(ExitFileNameEC));