Simula automata neural de Kinouchi-Copelli

kc.py

import networkx as nx
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import rc

def kc( ativos0, num_estados, rede, tempo):    
    firstN = set();
    ativos = set();
    ativos.add(ativos0);
    firstN.add(ativos0);
    deterministicos = set();    
    N = len(rede);
    for instante in tempo:
#        print " Tempo "+str(instante)+"\n";
#        print " Ativos"+str(ativos)+"\n";       
        deterministicos.update(ativos);
#        print " Transicao deterministica "+str(deterministicos)+"\n";
        ptrans = primeirosVizinhosAtivos( ativos, rede);
#        print " Rede "+str(rede.edge)+"\n";
        lt = len(ptrans);
        transicoes = np.random.rand(lt);
        c = 0;
#        print " sorteados "+str(transicoes)+"\n"+" 1-p "+str(ptrans)+"\n";
        
        for t in ptrans:         
            if transicoes[c] < (1.0-ptrans[t]):
                rede[t]['estado'] = rede[t]['estado'] + 1;
                ativos.add(t);
                c = c + 1;
            else:
                c = c + 1;
        for b in deterministicos.copy():                
            rede[b]['estado'] = np.mod( rede[b]['estado'] + 1, num_estados);
            if rede[b]['estado'] == 0:
                deterministicos.remove(b);
            if rede[b]['estado'] != 1:
                ativos.remove(b);
#        print " Rede "+str(rede.edge)+"\n";

        Num = len(ativos);
        if Num == 0:
            break;
        firstN.update(ativos);
    return (len(firstN),instante);

## Selecao de probabilidade, primeira tentativa produtorio_ativos(1-p). probabilidade de nao chegar nenhum sinal aos vizinhos, a probabilidade de chegar pelo menos um sinal eh 1-p_vizinhos ##
def primeirosVizinhosAtivos(ativos, net):
    p_vizinhos = {};
    for i in ativos:
        for j in net[i]:
            if j in p_vizinhos:
                p_vizinhos[j] = p_vizinhos[j]*(1.0 - net[i][j]['P']);
            else:
                if j != 'estado' and j != 'input':
                    if net[j]['estado'] == 0:
                        p_vizinhos[j] = net[i]['input']*(1-net[i][j]['P']); 
#    for i in p_vizinhos:
#        l = len(p_vizinhos[i]);
#        prod = 1.0;
#        print "\n Neighb_num "+str(p_vizinhos)+" \n"
#        for j in range(l):
#            prod = prod*(1.0-p_vizinhos[i][j]);        
#        p_vizinhos[i] = 1 - net[i]['input']*prod;
#    print "\n Quiescent transition P "+str(p_vizinhos)+" \n";    
    return p_vizinhos;
    
def shrink(a):
    a = list(a);
    a.sort();
    for el in a:
        while a.count(el) > 1:
            a.remove(el);
    return np.array(a);

def allIndex(rede, a, diverse = False):
    N = len(rede);
    indexes = [];
    if diverse:
        for i in range(N):
            if rede[i]['estado'] != a:
                indexes.append(i); 
        return indexes;
    else:
        for i in range(N):
            if rede[i]['estado'] == a:
                indexes.append(i);              
        return indexes;

def associarEstadosRede(rede, estados):
    for i in range(len(rede)):
        rede[i]['estado'] = estados[i];

def associarEstimuloExterno(rede, individuos, taxa):
    for i in range(len(rede)):
        rede[i]['input'] = 1.0;
    for member in individuos:
        rede[individuos.pop()]['input'] = np.exp(-taxa);

def associarProbabilidades(rede0, pmax):
    for edge in rede0.edge:
        for a in rede0.edge[ edge ].items():
            if a[0] == (edge-1):
                rede0.edge[ edge ][ a[0] ][ 'P' ] = rede0.edge[ edge -1][ edge ][ 'P' ];
            else:           
                rede0.edge[ edge ][ a[0] ][ 'P' ] = np.random.uniform(0, pmax);

def histogram(data,bins, minimum = -7.5, maximum = 7.5):
    binsize = (maximum - minimum)/float(bins);
    binxpos = [minimum + binsize*j for j in np.arange(bins)];
    binypos = [];
    data.sort();
    for i in range(bins):
        s = data.searchsorted([minimum + i*binsize,minimum + (i+1)*binsize]);
        binguys = data[s[0]:s[1]];
        binypos.append(binguys.size);
        norm = float(binsize*data.size) #float(sum(binypos));
        binypos = [j/norm for j in binypos];
    plt.bar(binxpos,binypos,binsize,color='red');

numero_estados = 2;
N = 100000;
k_mean = 10.0;
# Calcula <k> = np.mean(rede0.degree(rede0).values());  
tempot = np.linspace(0, 9999, 10000);
pmax = 2.*0.8/k_mean;
e0 = np.zeros(N);
estimulo = 0.0;
## Condicao inicial aleatoria
# e0 = np.random.randint(0, numero_estados, N);


## Para o caso de haver individuo com entradas externas do tipo processo de Poisson
#inp = np.random.randint(0, N, N/10);
#inp = shrink(inp);
S_ize = [];
T_ime = [];
for i in range(100000):
    rede0 = nx.generators.gnm_random_graph(N, N*k_mean/2);

    associarProbabilidades(rede0, pmax);
    associarEstadosRede(rede0, e0);
    associarEstimuloExterno(rede0, set(), 2.0);
    inicial = np.random.randint(0,N);
    rede0[inicial]['estado'] = 1.0;

    (a,b) = kc(inicial ,numero_estados, rede0, tempot);
    S_ize.append(a);
    T_ime.append(b);
    print "\n";
    print i;


S_ize.sort()
T_ime.sort()

l1 = S_ize.count(1);
l2 = T_ime.count(0);

for i in range(l1):
    S_ize.pop(0);
for i in range(l2):
    T_ime.pop(0);