From 6513f4e2b48f4b5526744f257bafe69b06a67804 Mon Sep 17 00:00:00 2001 From: timofej Date: Sat, 30 Sep 2023 19:53:06 +0200 Subject: [PATCH] more evaluations --- evaluation/1Layer Activity.py | 75 +++ evaluation/1Layer PHI.py | 91 ++++ evaluation/1Layer Suscept.py | 183 ++++++++ evaluation/1LayerPHI.py | 118 +++++ evaluation/2Layer Activity.py | 81 ++++ evaluation/2Layer Suscept.py | 155 +++++++ evaluation/2Layer fft.py | 433 ++++++++++++++++++ evaluation/2Layer phiplot.py | 98 ++++ evaluation/2axplot.py | 60 +++ evaluation/4Layer Activity.py | 77 ++++ evaluation/4Layer Bootstrap Resultant.py | 286 ++++++++++++ evaluation/4Layer Causal 4way resultant.py | 369 +++++++++++++++ evaluation/4Layer Causal Significance.py | 191 ++++++++ evaluation/4Layer Causal resultant.py | 405 ++++++++++++++++ evaluation/4Layer Res+EI.py | 76 ++- .../4Layer Resultant 4way Phase Evolution.py | 277 +++++++++++ .../4Layer Resultant Phase Evolution.py | 86 +++- evaluation/4Layer Resultant.py | 106 +++++ evaluation/4Layer Significance.py | 25 + evaluation/4Layer phiplot.py | 98 ++++ evaluation/4laysign.py | 173 +++++++ evaluation/phi.py | 173 +++++++ evaluation/plot.py | 58 ++- evaluation/plots/bootstrap_plot.py | 29 ++ evaluation/plots/plot.py | 169 +++++++ evaluation/sign.py | 374 +++++++++++++++ evaluation/simpleplot.py | 28 ++ evaluation/updated_bib_file.bib | 353 ++++++++++++++ evaluation/updated_sync.bib | 353 ++++++++++++++ examples/Simulate1Layer.py | 22 +- examples/Simulate2Layers.py | 22 +- examples/Simulate4Layers.py | 38 +- neuropercolation/automaton.py | 42 +- neuropercolation/display.py | 24 +- 34 files changed, 5032 insertions(+), 116 deletions(-) create mode 100644 evaluation/1Layer Activity.py create mode 100644 evaluation/1Layer PHI.py create mode 100644 evaluation/1Layer Suscept.py create mode 100644 evaluation/1LayerPHI.py create mode 100644 evaluation/2Layer Activity.py create mode 100644 evaluation/2Layer Suscept.py create mode 100644 evaluation/2Layer fft.py create mode 100644 evaluation/2Layer phiplot.py create mode 100644 evaluation/2axplot.py create mode 100644 evaluation/4Layer Activity.py create mode 100644 evaluation/4Layer Bootstrap Resultant.py create mode 100644 evaluation/4Layer Causal 4way resultant.py create mode 100644 evaluation/4Layer Causal Significance.py create mode 100644 evaluation/4Layer Causal resultant.py create mode 100644 evaluation/4Layer Resultant 4way Phase Evolution.py create mode 100644 evaluation/4Layer Resultant.py create mode 100644 evaluation/4Layer Significance.py create mode 100644 evaluation/4Layer phiplot.py create mode 100644 evaluation/4laysign.py create mode 100644 evaluation/phi.py create mode 100644 evaluation/plots/bootstrap_plot.py create mode 100644 evaluation/plots/plot.py create mode 100644 evaluation/sign.py create mode 100644 evaluation/simpleplot.py create mode 100644 evaluation/updated_bib_file.bib create mode 100644 evaluation/updated_sync.bib diff --git a/evaluation/1Layer Activity.py b/evaluation/1Layer Activity.py new file mode 100644 index 0000000..e7314d8 --- /dev/null +++ b/evaluation/1Layer Activity.py @@ -0,0 +1,75 @@ +# -*- coding: utf-8 -*- +""" +Created on Tue Aug 30 14:25:12 2022 + +@author: timof +""" + +import sys +import os +import json + +from plot import qtplot + + +import math as m +import numpy as np + + +vect = np.vectorize + +@vect +def log2(x): + try: + return m.log2(x) + except ValueError: + if x==0: + return float(0) + else: + raise + +def new_folder(path): + if not os.path.exists(path): + os.makedirs(path) + return path + +path = '/cloud/Public/_data/neuropercolation/1lay/steps=1000100/' +suffix = '' + +chi = chr(967) +vareps = chr(949) + +vals = [[],[]] + +runsteps = 1000100 + +eps_space = np.linspace(0.005, 0.5, 100) +eps_space = eps_space[1::2] + +dims = list(range(3,10))#+[16,49] + +mode='density' +ma=[] +s=[] +k=[] +mk=[] +lastkurt=None +for dim in dims[-1:]: + dimpath = new_folder(path + f'dim={dim:02}/') + for epsilon in eps_space: + with open(dimpath+f"eps={round(epsilon,3):.3f}_activation.txt", 'r', encoding='utf-8') as f: + activation = np.array(json.load(f)[:500]) + + qtplot(f"Activity time series for eps={round(epsilon,3):.3f}", + [list(range(500))], + [activation], + x_tag = 'time step', + y_tag = f'activity density', + y_range = (-1,1), + export=True, + path=dimpath+"evolution/", + filename=f'eps={round(epsilon,3):.3f}_evolution.png', + close=True) + +mode = 'density' +#%% diff --git a/evaluation/1Layer PHI.py b/evaluation/1Layer PHI.py new file mode 100644 index 0000000..dcb498a --- /dev/null +++ b/evaluation/1Layer PHI.py @@ -0,0 +1,91 @@ +# -*- coding: utf-8 -*- +""" +Created on Tue Aug 30 14:25:12 2022 + +@author: timof +""" + +import sys +import os +import json + +from plot import qtplot +from phi import MIP,calc_mips,smartMIP + +import math as m +import numpy as np +from datetime import datetime + +vect = np.vectorize + +@vect +def log2(x): + try: + return m.log2(x) + except ValueError: + if x==0: + return float(0) + else: + raise + +def new_folder(path): + if not os.path.exists(path): + os.makedirs(path) + return path + +path = '/cloud/Public/_data/neuropercolation/1lay/steps=1000100/' +suffix = '' + +chi = chr(967) +vareps = chr(949) +varphi = chr(981) + +vals = [[],[]] + +runsteps = 1000100 + +eps_space = np.linspace(0.005, 0.5, 100) +eps_space = eps_space[1::2] + +dims = list(range(3,10))#+[16,49] + +mode='density' +ma=[] +s=[] +k=[] +mk=[] +lastkurt=None +phinum=10000 +PHIS=[] +for dim in dims[:1]: + dimpath = new_folder(path + f'dim={dim:02}/') + for epsilon in eps_space: + try: + with open(dimpath+f"eps={round(epsilon,3):.3f}_phi_{phinum}.txt", 'r', encoding='utf-8') as f: + phis=json.load(f) + PHIS.append(np.average(phis)) + print(f'Loaded and Done eps={epsilon:.3f} with dim={dim} at {datetime.now()}') + except: + with open(dimpath+f"eps={round(epsilon,3):.3f}_states.txt", 'r', encoding='utf-8') as f: + states = np.array(json.load(f)[100:phinum+100]) + + phis = [smartMIP(dim,state.translate(str.maketrans('','','.-=')),epsilon)[1] for state in states] + with open(dimpath+f"eps={round(epsilon,3):.3f}_phi_{phinum}.txt", 'w', encoding='utf-8') as f: + json.dump(phis, f, indent=1) + PHIS.append(np.average(phis)) + print(f'Generated and Done eps={epsilon:.3f} with dim={dim} at {datetime.now()}') + #%% + qtplot(f"Average phi for different noise parameters", + [eps_space], + [PHIS], + [''], + colors='w', + x_tag = f'noise level {vareps}', + y_tag = f'effect integration {varphi}', + export=False, + path=dimpath+"evolution/", + filename=f'eps={round(epsilon,3):.3f}_evolution.png', + close=False) + +mode = 'density' +#%% diff --git a/evaluation/1Layer Suscept.py b/evaluation/1Layer Suscept.py new file mode 100644 index 0000000..14204b8 --- /dev/null +++ b/evaluation/1Layer Suscept.py @@ -0,0 +1,183 @@ +# -*- coding: utf-8 -*- +""" +Created on Tue Aug 30 14:25:12 2022 + +@author: timof +""" + +import sys +import os +import json + +from plot import qtplot + + +import math as m +import numpy as np + + +vect = np.vectorize + +@vect +def log2(x): + try: + return m.log2(x) + except ValueError: + if x==0: + return float(0) + else: + raise + +path = '/cloud/Public/_data/neuropercolation/1lay/steps=1000100/' +suffix = '' + +chi = chr(967) +vareps = chr(949) + +vals = [[],[]] + +runsteps = 1000100 + +eps_space = np.linspace(0.005, 0.5, 100) +eps_space = eps_space[1::2] + +dims = list(range(3,10))#+[16,49] + +def __calc_chi(mode='density'): + ma=[] + s=[] + k=[] + mk=[] + lastkurt=None + for dim in dims: + try: + with open(path+f"magnets/magnetisation_{mode}_dim={dim:02}.txt", 'r', encoding='utf-8') as f: + magnet = json.load(f) + with open(path+f"suscepts/susceptibility_{mode}_dim={dim:02}.txt", 'r', encoding='utf-8') as f: + suscept = json.load(f) + with open(path+f"kurts/kurtosis_{mode}_dim={dim:02}.txt", 'r', encoding='utf-8') as f: + kurt = json.load(f) + except: + magnet = [] + suscept = [] + kurt = [] + jumped=False + print('magnets or suscept or kurt file not found') + if not os.path.exists(path+"magnets/"): + os.makedirs(path+"magnets/") + if not os.path.exists(path+"suscepts/"): + os.makedirs(path+"suscepts/") + if not os.path.exists(path+"kurts/"): + os.makedirs(path+"kurts/") + for epsilon in eps_space: + try: + with open(path+f"dim={dim:02}/eps={round(epsilon,3):.3f}_activation.txt", 'r', encoding='utf-8') as f: + activation = np.array(json.load(f)[100:]) + + if mode=='absolute': + dmag = np.sum(activation*dim**2)/len(activation) + mag = np.sum(np.abs(activation*dim**2))/len(activation) + mag2 = np.sum(activation**2*dim**4)/len(activation) + mag4 = np.sum(activation**4*dim**8)/len(activation) + elif mode=='density': + dmag = np.sum(activation)/len(activation) + mag = np.sum(np.abs(activation))/len(activation) + mag2 = np.sum(activation**2)/len(activation) + mag4 = np.sum(activation**4)/len(activation) + else: + raise NotImplementedError + + mag = round(abs(mag),6) + fluct = round(mag2-mag**2,6) + tail = round(mag4/mag2**2,6) + + magnet.append(mag) + suscept.append(fluct) + kurt.append(tail) + print(f"Done dim={dim:02} eps={round(epsilon,3):.3f}: mag={mag}") + print(f"Done dim={dim:02} eps={round(epsilon,3):.3f}: fluct={fluct}") + print(f"Done dim={dim:02} eps={round(epsilon,3):.3f}: tail={tail}") + jumped=True + except: + if not jumped: + magnet.append(1) + elif jumped: + magnet.append(0) + suscept.append(0) + if not jumped: + kurt.append(1) + elif jumped: + kurt.append(3) + print(f"Missing dim={dim:02} eps={round(epsilon,3):.3f}") + with open(path+f"magnets/magnetisation_{mode}_dim={dim:02}.txt", 'w', encoding='utf-8') as f: + json.dump(magnet, f, ensure_ascii=False, indent=1) + with open(path+f"suscepts/susceptibility_{mode}_dim={dim:02}.txt", 'w', encoding='utf-8') as f: + json.dump(suscept, f, ensure_ascii=False, indent=1) + with open(path+f"kurts/kurtosis_{mode}_dim={dim:02}.txt", 'w', encoding='utf-8') as f: + json.dump(kurt, f, ensure_ascii=False, indent=1) + + if lastkurt is not None: + pos=0 + while kurt[pos]2: + val=1 + elif zm>2: + val=0 + elif zp==zm: + val=0.5 + elif zm==2: + val=0.5**(3-zp) + elif zp==2: + val=1-0.5**(3-zm) + elif zm==0 and zp==1: + val=9/16 + elif zp==0 and zm==1: + val=7/16 + else: + raise NotImplementedError(zp,zm) + return val + +path = '/cloud/Public/_data/neuropercolation/1lay/steps=100000/' + +def phi(dim,statestr,partstr,eps): + length = dim**2 + eta = 1-eps + + state = np.array([int(q) for q in statestr]) + state = list(state.reshape((dim,dim))) + state = [list([int(cell) for cell in row]) for row in state] + + part = np.array([int(p) for p in partstr]) + part = list(part.reshape((dim,dim))) + part = [list([int(cell) for cell in row]) for row in part] + + inp = [[q+sum([state[(i+1)%dim][j], + state[(i-1)%dim][j], + state[i][(j+1)%dim], + state[i][(j-1)%dim] + ]) for j,q in enumerate(row)] for i,row in enumerate(state)] + + beps = [[int(inp[i][j]>2)*eta+int(inp[i][j]<3)*eps for j,q in enumerate(row)] for i,row in enumerate(state)] + + zplus = [[q+sum([state[(i+1)%dim][j]*(part[i][j]==part[(i+1)%dim][j]), + state[(i-1)%dim][j]*(part[i][j]==part[(i-1)%dim][j]), + state[i][(j+1)%dim]*(part[i][j]==part[i][(j+1)%dim]), + state[i][(j-1)%dim]*(part[i][j]==part[i][(j-1)%dim]) + ]) for j,q in enumerate(row)] for i,row in enumerate(state)] + zminus = [[sum([(1-state[(i+1)%dim][j])*(part[i][j]==part[(i+1)%dim][j]), + (1-state[(i-1)%dim][j])*(part[i][j]==part[(i-1)%dim][j]), + (1-state[i][(j+1)%dim])*(part[i][j]==part[i][(j+1)%dim]), + (1-state[i][(j-1)%dim])*(part[i][j]==part[i][(j-1)%dim]) + ]) for j,q in enumerate(row)] for i,row in enumerate(state)] + + kplus = [[kcomb(zplus[i][j],zminus[i][j]) for j,q in enumerate(row)] for i,row in enumerate(state)] + + pi = [[eps*(1-kplus[i][j]) + eta*kplus[i][j] for j,q in enumerate(row)] for i,row in enumerate(state)] + + crossent = [[-beps[i][j]*m.log2(pi[i][j])-(1-beps[i][j])*m.log2(1-pi[i][j]) for j,q in enumerate(row)] for i,row in enumerate(state)] + + return np.sum(crossent) - length*H2(eps) + +def MIP(dim,statestr,eps): + lophi=np.inf + mip = [] + for parti in range(1,2**(dim**2-1)): + partstr = bin(parti)[2:].zfill(dim**2) + print(partstr) + curphi = phi(dim,statestr,partstr,eps) + + if curphi=halfh])/sum([specc[pos] for pos in range(len(specc)) if specc[pos]=sigpower/2: + median = pos + stepdiff = speccnorm[pos] + halfdiff = sum(speccnorm[:pos+1]) - sigpower/2 + + freqstep = freqqs[pos] - freqqs[pos-1] + freqmax = freqqs[pos] - freqstep * halfdiff/stepdiff + + break + + for pos in range(median, len(speccnorm)): + if cumspecc[pos+1]>=sigpower*3/4: + stepdiff = speccnorm[pos] + halfdiff = sum(speccnorm[:pos+1]) - sigpower*3/4 + + + freqstep = freqqs[pos] - freqqs[pos-1] + freqquart3 = freqqs[pos] - freqstep * halfdiff/stepdiff + + break + for pos in range(median+1): + if cumspecc[pos+1]>=sigpower/4: + stepdiff = speccnorm[pos] + halfdiff = sum(speccnorm[:pos+1]) - sigpower/4 + + freqstep = freqqs[pos] - freqqs[pos-1] + freqquart1 = freqqs[pos] - freqstep * halfdiff/stepdiff + + break + + + + freqfwhm = freqquart3-freqquart1 + gamma = freqfwhm/2 + + freqmax_parts.append(freqmax) + freqfwhm_parts.append(freqfwhm) + + print(gamma**2/median**2) + + lorentz = sigpower/m.pi * accum/partlen * gamma/((freqqs-freqmax)**2+gamma**2) + noisefloor + lorentz_list.append(lorentz) + + speccpdf = speccnorm/sigpower + deviation = np.sqrt(sum((specc-lorentz)**2)/len(specc))/sigpower + + maxpos = 0 + for i in range(len(freqqs)): + if abs(freqqs[i]-freqmax)lorentz[maxpos]/2: + pos-=1 + lhm=freqqs[pos] + except: + lhm=freqqs[0] + pos=maxpos + try: + while lorentz[pos]>lorentz[maxpos]/2: + pos+=1 + uhm=freqqs[pos] + except: + uhm=freqqs[-1] + + SPECC = np.average(specc_list,axis=0) + LORENTZ = np.average(lorentz_list,axis=0) + DEV = np.sqrt(sum((SPECC-LORENTZ)**2)/len(SPECC))/sum(SPECC) + + TOTP.append(np.average(totalpower_parts)) + SIGP.append(np.average(sigpower_parts)) + SN.append(np.average(sn_parts)) + STON.append(sigtonoi) + + FMAX.append(np.average(freqmax_parts)) + FWHM.append(np.average(freqfwhm_parts)) + LHM.append(lhm) + UHM.append(uhm) + Q.append(np.average(freqmax_parts)/np.average(freqfwhm_parts)) + SOS.append(DEV) + #%% + with open(path+f"signal/SN_{parts}_accum={accum}_dim={dim}_runsteps={runsteps}.txt", 'w', encoding='utf-8') as f: + json.dump(SN, f, ensure_ascii=False, indent=4) + with open(path+f"signal/FMAX_{parts}_accum={accum}_dim={dim}_runsteps={runsteps}.txt", 'w', encoding='utf-8') as f: + json.dump(FMAX, f, ensure_ascii=False, indent=4) + with open(path+f"signal/FWHM_{parts}_accum={accum}_dim={dim}_runsteps={runsteps}.txt", 'w', encoding='utf-8') as f: + json.dump(FWHM, f, ensure_ascii=False, indent=4) + with open(path+f"signal/LHM_{parts}_accum={accum}_dim={dim}_runsteps={runsteps}.txt", 'w', encoding='utf-8') as f: + json.dump(LHM, f, ensure_ascii=False, indent=4) + with open(path+f"signal/UHM_{parts}_accum={accum}_dim={dim}_runsteps={runsteps}.txt", 'w', encoding='utf-8') as f: + json.dump(UHM, f, ensure_ascii=False, indent=4) + with open(path+f"signal/TOTP_{parts}_accum={accum}_dim={dim}_runsteps={runsteps}.txt", 'w', encoding='utf-8') as f: + json.dump(TOTP, f, ensure_ascii=False, indent=4) + with open(path+f"signal/SIGP_{parts}_accum={accum}_dim={dim}_runsteps={runsteps}.txt", 'w', encoding='utf-8') as f: + json.dump(SIGP, f, ensure_ascii=False, indent=4) + with open(path+f"signal/Q_{parts}_accum={accum}_dim={dim}_runsteps={runsteps}.txt", 'w', encoding='utf-8') as f: + json.dump(Q, f, ensure_ascii=False, indent=4) + with open(path+f"signal/STON_{parts}_accum={accum}_dim={dim}_runsteps={runsteps}.txt", 'w', encoding='utf-8') as f: + json.dump(STON, f, ensure_ascii=False, indent=4) + + #%% + plotpath = path + 'spectra/' + # qtplot(f'Frequency spectrum for dim={dim} eps={eps:.3f}', + # [freqqs]*2, + # [LORENTZ,SPECC], + # [f'lorentz regression dev={DEV}','spectrum'], + # x_tag = 'frequency f', + # y_tag = 'spectral power p', + # y_log = False, + # export=True, + # lw=1, + # path=plotpath, + # filename=f'Lorentz and Cauchy fit to spectrum {parts} parts accum={accum} eps={round(eps,3):.3f}.png', + # close=True) + + + #%% + + print("Power {}, Signal {}".format(round(totalpower), round(sigpower))) + + print("Done {:.3f}: DEV {}, FMAX {}, FWHM {}".format(eps, round(deviation,3), round(freqmax,5), round(FWHM[-1],5))) + + LORENTZ_eps.append(LORENTZ) + SPECC_eps.append(SPECC) + + #%% + # vareps=chr(949) + # qtplot(f'Frequency spectra to noise level', + # [freqqs[50::100]]*len(eps_space), + # [SPECC[50::100] for SPECC in SPECC_eps], + # [f'{vareps}={eps:.3f}' for eps in eps_space], + # x_tag = 'frequency f', + # y_tag = 'energy spectral density', + # y_log = True, + # export=False, + # lw=2, + # path=plotpath, + # filename=f'Lorentz and Cauchy fit to spectrum log {parts} parts accum={accum} eps={round(eps,3):.3f}.png', + # close=False) + #%% + TOTP_eps.append(TOTP) + SIGP_eps.append(SIGP) + SN_eps.append(SN) + STON_eps.append(STON) + FMAX_eps.append(FMAX) + FWHM_eps.append(FWHM) + Q_eps.append(Q) + LHM_eps.append(LHM) + UHM_eps.append(UHM) + SOS_eps.append(SOS) + + #%% +sigpath = basepath + 'signal/' +if not os.path.exists(sigpath): + os.makedirs(sigpath) + +# qtplot(f'Signal to noise for dim={dim} eps={eps:.3f}', +# [eps_space], +# [SN], +# ['signalpower to noisefloor'], +# x_tag = 'epsilon', +# y_tag = 'ratio', +# y_log = False, +# export=True, +# path=sigpath, +# filename=f'SN_plot.png', +# close=True) + +qtplot(f'Q facor for odd automaton sizes', + [eps_space]*3, + Q_eps, + [f'dim={dim}x{dim}' for dim in dims], + y_tag = 'purity', + y_log = False, + export=True, + path=sigpath, + filename=f'Q_plot_{parts}parts_accum={accum}.png', + close=False) +#%% + +# with open(path+f"signal/FMAX_1_accum={accum}_dim={dim}_runsteps={runsteps}.txt", 'r', encoding='utf-8') as f: +# FMAX=json.load(f) +# with open(path+f"signal/FWHM_1_accum={accum}_dim={dim}_runsteps={runsteps}.txt", 'r', encoding='utf-8') as f: +# FWHM=json.load(f) + +qtplot(f'Dominant frequency', + [[0]+list(eps_space[:-1])]*3, + [[0]+FMAX[-1],[0]+LHM[-1],[0]+UHM[-1]], + # [f'dim={dim}x{dim}' for dim in dims], + ['dominant frequency', 'lesser half maximum', 'upper half maximum'], + #colors=['g','r','r'], + y_tag = 'frequency f', + y_log = True, + export=True, + path=sigpath, + filename=f'FMAX+FWHMto50_plot_{parts}parts_accum={accum}.png', + close=False) +#%% +qtplot(f'Dominant frequency', + [[0]+list(eps_space[:-1])]*3, + [[0]+FMAX[:-1] for FMAX in FMAX_eps],#, [[0]+LHM[i][:50] for LHM in LHM_eps], [[0]+UHM[i][:50] for UHM in UHM_eps], + [f'dim={dim}x{dim}' for dim in dims], + #[['dominant frequency', 'lesser half maximum', 'upper half maximum'], + #colors=['g','r','r'], + y_tag = 'frequency f', + y_log = True, + export=True, + path=sigpath, + filename=f'FMAX+FWHMto50_plot_{parts}parts_accum={accum}.png', + close=False) +#%% +qtplot(f'Total and signal energy to noise level', + [eps_space]*6, + list(np.array(SIGP_eps)/runsteps*accum)+list(np.array(TOTP_eps)/runsteps*accum), + [f'signal energy for dim={dim:02d}' for dim in dims]+[f'total energy for dim={dim:02d}' for dim in dims], + colors=['b','g','r','c','y','m'], + y_tag = f'energy', + y_log = True, + export=True, + path=sigpath, + filename=f'POWER_plot_{parts}parts_accum={accum}.png', + close=False) +#%% +qtplot(f'Signal-to-noise ratio to noise level', + [eps_space]*len(dims), + STON_eps, + [f'dim={dim:02d}x{dim:02d}' for dim in dims], + y_tag = f'ratio S/N', + y_log = True, + export=True, + path=sigpath, + filename=f'SN_plot_{parts}parts_accum={accum}.png', + close=False) +#%% +qtplot(f'Goodness of fit to Lorentz curve', + [eps_space]*len(dims), + SOS_eps, + [f'dim={dim:02d}x{dim:02d}' for dim in dims], + y_tag = f'mean deviation', + y_log = False, + export=True, + path=sigpath, + filename=f'DEV_plot_{parts}parts_accum={accum}.png', + close=False) +#%% +with open(basepath+f"signal/SN_{parts}_accum={accum}_dim={dim}_runsteps={runsteps}.txt", 'w', encoding='utf-8') as f: + json.dump(SN_eps, f, ensure_ascii=False, indent=4) +with open(basepath+f"signal/FMAX_{parts}_accum={accum}_dim={dim}_runsteps={runsteps}.txt", 'w', encoding='utf-8') as f: + json.dump(FMAX_eps, f, ensure_ascii=False, indent=4) +with open(basepath+f"signal/FWHM_{parts}_accum={accum}_dim={dim}_runsteps={runsteps}.txt", 'w', encoding='utf-8') as f: + json.dump(FWHM_eps, f, ensure_ascii=False, indent=4) +with open(basepath+f"signal/LHM_{parts}_accum={accum}_dim={dim}_runsteps={runsteps}.txt", 'w', encoding='utf-8') as f: + json.dump(LHM_eps, f, ensure_ascii=False, indent=4) +with open(basepath+f"signal/UHM_{parts}_accum={accum}_dim={dim}_runsteps={runsteps}.txt", 'w', encoding='utf-8') as f: + json.dump(UHM_eps, f, ensure_ascii=False, indent=4) +with open(basepath+f"signal/TOTP_{parts}_accum={accum}_dim={dim}_runsteps={runsteps}.txt", 'w', encoding='utf-8') as f: + json.dump(TOTP_eps, f, ensure_ascii=False, indent=4) +with open(basepath+f"signal/SIGP_{parts}_accum={accum}_dim={dim}_runsteps={runsteps}.txt", 'w', encoding='utf-8') as f: + json.dump(SIGP_eps, f, ensure_ascii=False, indent=4) +with open(basepath+f"signal/Q_{parts}_accum={accum}_dim={dim}_runsteps={runsteps}.txt", 'w', encoding='utf-8') as f: + json.dump(Q_eps, f, ensure_ascii=False, indent=4) +with open(basepath+f"signal/STON_{parts}_accum={accum}_dim={dim}_runsteps={runsteps}.txt", 'w', encoding='utf-8') as f: + json.dump(STON_eps, f, ensure_ascii=False, indent=4) +with open(basepath+f"signal/Q_{parts}_accum={accum}_dim={dim}_runsteps={runsteps}.txt", 'w', encoding='utf-8') as f: + json.dump(Q_eps, f, ensure_ascii=False, indent=4) +with open(basepath+f"signal/STON_{parts}_accum={accum}_dim={dim}_runsteps={runsteps}.txt", 'w', encoding='utf-8') as f: + json.dump(STON_eps, f, ensure_ascii=False, indent=4) + +def plot_execute(): + if sys.flags.interactive != 1 or not hasattr(qt.QtCore, 'PYQT_VERSION'): + qt.QtGui.QApplication.exec_() + diff --git a/evaluation/2Layer phiplot.py b/evaluation/2Layer phiplot.py new file mode 100644 index 0000000..a338e06 --- /dev/null +++ b/evaluation/2Layer phiplot.py @@ -0,0 +1,98 @@ +# -*- coding: utf-8 -*- +""" +Created on Tue Aug 30 14:25:12 2022 + +@author: timof +""" + +import sys +import os +import json + +from plot import qtplot + + +import math as m +import numpy as np + + +vect = np.vectorize + +@vect +def log2(x): + try: + return m.log2(x) + except ValueError: + if x==0: + return float(0) + else: + raise + +def new_folder(path): + if not os.path.exists(path): + os.makedirs(path) + return path + +phase = np.vectorize(lambda x,y: (m.atan2(y,x)+m.pi)%(2*m.pi)-m.pi) +diff = np.vectorize(lambda x,y: (y-x+m.pi)%(2*m.pi)-m.pi) +H2 = lambda x: -x*m.log2(x)-(1-x)*m.log2(1-x) + + +path='/cloud/Public/_data/neuropercolation/2lay/steps=100100/' +suffix = '' + +chi = chr(967) +vareps = chr(949) +varphi = chr(981) + +vals = [[],[]] + +runsteps = 1000100 + +eps_space = np.linspace(0.005, 0.5, 100) +eps_space = eps_space[1::2] + +dims = list(range(3,10,2))#+[16,49] + +mode='density' +ma=[] +s=[] +k=[] +mk=[] +PHI=[] +lastkurt=None +for dim in dims: + phis=[] + con_gap = 3 + cons = [(n,(2*n+m)%dim) for n in range(dim) for m in range(0,dim-2,con_gap)] + dimpath = new_folder(path + f'dim={dim:02}/') + for epsilon in eps_space: + try: + with open(dimpath+f"eps={round(epsilon,3):.3f}_ei.txt", 'r', encoding='utf-8') as f: + ei = np.array(json.load(f)[100:]) + except: + print(f'Calcing phi for eps={epsilon}') + with open(dimpath+f"eps={round(epsilon,3):.3f}_channels.txt", 'r', encoding='utf-8') as f: + channels = np.array(json.load(f)[100:]) + + ei = np.sum(channels,axis=0)*(1-H2(epsilon)) + ei=list(ei) + with open(dimpath+f"eps={round(epsilon,3):.3f}_ei.txt", 'w', encoding='utf-8') as f: + json.dump(list(ei),f, ensure_ascii=False,indent=1) + + phi=np.average(ei) + phis.append(phi) + PHI.append(phis) + #%% +qtplot(f"Mean effect integration over noise level", + [[0]+list(eps_space)]*len(dims), + [[0]+phi for phi in PHI[::-1]], + [f'dim={dim:02d}x{dim:02d}' for dim in dims[::-1]], + y_tag = f'effect integration {varphi}', + export=True, + path=dimpath+"", + filename=f'eps={round(epsilon,3):.3f}_evolution.png', + close=False) + +mode = 'density' +#%% diff --git a/evaluation/2axplot.py b/evaluation/2axplot.py new file mode 100644 index 0000000..a40f7be --- /dev/null +++ b/evaluation/2axplot.py @@ -0,0 +1,60 @@ +""" +Demonstrates a way to put multiple axes around a single plot. + +(This will eventually become a built-in feature of PlotItem) +""" + + +import pyqtgraph as pg + +pg.mkQApp() + +pw = pg.PlotWidget() +pw.show() +pw.setWindowTitle('pyqtgraph example: MultiplePlotAxes') +p1 = pw.plotItem +p1.setLabels(left='axis 1') + +## create a new ViewBox, link the right axis to its coordinate system +p2 = pg.ViewBox() +p1.showAxis('right') +p1.scene().addItem(p2) +p1.getAxis('right').linkToView(p2) +p2.setXLink(p1) +p1.getAxis('right').setLabel('axis2', color='#0000ff') + +## create third ViewBox. +## this time we need to create a new axis as well. +p3 = pg.ViewBox() +ax3 = pg.AxisItem('right') +p1.layout.addItem(ax3, 2, 3) +p1.scene().addItem(p3) +ax3.linkToView(p3) +p3.setXLink(p1) +ax3.setZValue(-10000) +ax3.setLabel('axis 3', color='#ff0000') + + +## Handle view resizing +def updateViews(): + ## view has resized; update auxiliary views to match + global p1, p2, p3 + p2.setGeometry(p1.vb.sceneBoundingRect()) + p3.setGeometry(p1.vb.sceneBoundingRect()) + + ## need to re-update linked axes since this was called + ## incorrectly while views had different shapes. + ## (probably this should be handled in ViewBox.resizeEvent) + p2.linkedViewChanged(p1.vb, p2.XAxis) + p3.linkedViewChanged(p1.vb, p3.XAxis) + +updateViews() +p1.vb.sigResized.connect(updateViews) + + +p1.plot([1,2,4,8,16,32]) +p2.addItem(pg.PlotCurveItem([10,20,40,80,40,20], pen='b')) +p3.addItem(pg.PlotCurveItem([3200,1600,800,400,200,100], pen='r')) + +if __name__ == '__main__': + pg.exec() \ No newline at end of file diff --git a/evaluation/4Layer Activity.py b/evaluation/4Layer Activity.py new file mode 100644 index 0000000..704ce38 --- /dev/null +++ b/evaluation/4Layer Activity.py @@ -0,0 +1,77 @@ +# -*- coding: utf-8 -*- +""" +Created on Tue Aug 30 14:25:12 2022 + +@author: timof +""" + +import sys +import os +import json + +from plot import qtplot + + +import math as m +import numpy as np + + +vect = np.vectorize + +@vect +def log2(x): + try: + return m.log2(x) + except ValueError: + if x==0: + return float(0) + else: + raise + +def new_folder(path): + if not os.path.exists(path): + os.makedirs(path) + return path + +path = '/cloud/Public/_data/neuropercolation/4lay/cons=27-knight_steps=100100_causal/dim=09/batch=0/' +suffix = '' + +chi = chr(967) +vareps = chr(949) + +vals = [[],[]] + +runsteps = 1000100 + +eps_space = np.linspace(0.005, 0.5, 100) +eps_space = eps_space[1::2] + +dims = list(range(3,10))#+[16,49] + +mode='density' +ma=[] +s=[] +k=[] +mk=[] +lastkurt=None +for dim in dims[-1:]: + dimpath = new_folder(path + f'dim={dim:02}/') + for epsilon in eps_space[:]: + with open(path+f"eps={round(epsilon,3):.3f}_phase_diff.txt", 'r', encoding='utf-8') as f: + phase_diff = np.array(json.load(f)[:500]) + with open(path+f"eps={round(epsilon,3):.3f}_ei.txt", 'r', encoding='utf-8') as f: + phase_diff = np.array(json.load(f)[:500]) + + qtplot(f"Phase relation time series for eps={round(epsilon,3):.3f}", + [list(range(500))]*2, + [phase_diff], + x_tag = 'time step', + y_tag = f'phase diffe´rence', + y_range = (-m.pi,m.pi), + export=True, + path=dimpath+"evolution/", + filename=f'eps={round(epsilon,3):.3f}_evolution.png', + close=False) + +mode = 'density' +#%% diff --git a/evaluation/4Layer Bootstrap Resultant.py b/evaluation/4Layer Bootstrap Resultant.py new file mode 100644 index 0000000..87ce790 --- /dev/null +++ b/evaluation/4Layer Bootstrap Resultant.py @@ -0,0 +1,286 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Mon Aug 21 14:59:22 2023 + +@author: astral +""" + +import os +import json +import math as m +import numpy as np +from numpy.linalg import norm +from datetime import datetime +from random import sample as choose +from random import random + +from plot import qtplot + +from neuropercolation import Simulate4Layers + + +def resultant(sample): + phase_x = [m.cos(ind) for ind in sample] + phase_y = [m.sin(ind) for ind in sample] + + return (np.average(phase_x), np.average(phase_y)) + +def coherence(sample): + phase_x = [m.cos(ind) for ind in sample] + + return np.average(phase_x) + + +def bootstrap_stat(whole, sub, strength=10000, estimator='resultant'): + k = len(sub) + + if estimator=='resultant': + whole_est = norm(resultant(whole)) + sub_est = norm(resultant(sub)) + + boot_dist = [] + for i in range(strength): + boot_dist.append(norm(resultant(choose(whole,k)))) + if i%1000==999: + print(f'Done {i:0{len(str(strength))}d} bootstraps') + + confidence = len([val for val in boot_dist if (val-whole_est)sub_est])/len(boot_dist) + else: + raise NotImplemented + return confidence + +def bootstrap_res_parts(wholes, subs, strength=10000): + ks = [len(s) for s in subs] + + whole = [w for whole in wholes for w in whole] + whole_est = norm(resultant(whole)) + + resultants = [] + for i in range(strength): + sample = [pha for j in range(len(wholes)) for pha in choose(wholes[j],ks[j])] + resultants.append(norm(resultant(sample))) + if i%1000==999: + print(f'Done {i:0{len(str(strength))}d} bootstraps') + + return resultants + +def new_folder(path): + if not os.path.exists(path): + os.makedirs(path) + return path + +phase = np.vectorize(lambda x,y: (m.atan2(y,x)+m.pi)%(2*m.pi)-m.pi) +diff = np.vectorize(lambda x,y: (y-x+m.pi)%(2*m.pi)-m.pi) +H2 = lambda x: -x*m.log2(x)-(1-x)*m.log2(1-x) + +extremes = None +maxdt = 250 + +stp = 1000100 +batch = 0 + +eps_space = list(np.linspace(0.01,0.5,50)) + +print(f'Started at {datetime.now()}') + +for dim in [9]: + for eps in eps_space[:1]:#list(eps_space[3:8])+list(eps_space[:3])+list(eps_space[8:]): + eps = round(eps,3) + path='/cloud/Public/_data/neuropercolation/4lay/cons=dimtimesdimby3_steps=100100/dim=09_cons=27/batch=0/' + + try: + with open(path+f"eps={round(eps,3):.3f}_phase_diff.txt", 'r', encoding='utf-8') as f: + phase_diff = json.load(f) + except: + with open(path+f"eps={round(eps,3):.3f}_activation.txt", 'r', encoding='utf-8') as f: + activation = json.load(f)[100:] + + osc = list(zip(*activation)) + phase_abs = np.array([[np.arctan2(*act[::-1]) for act in osc[i]] for i in range(2)]) + phase_diff = diff(phase_abs[0],phase_abs[1]) + phase_diff = [round(pha,6) for pha in phase_diff] + + with open(path+f"eps={round(eps,3):.3f}_phase_diff.txt", 'w', encoding='utf-8') as f: + json.dump(list(phase_diff), f, indent=1) + + + all_res = norm(resultant(phase_diff)) + av_diff = np.arccos(all_res) + + try: + with open(path+f"eps={round(eps,3):.3f}_ei.txt", 'r', encoding='utf-8') as f: + ei = json.load(f) + except: + with open(path+f"eps={round(eps,3):.3f}_channels.txt", 'r', encoding='utf-8') as f: + channels = json.load(f)[100:] + + ei = [round(np.sum(cha)*(1-H2(eps)),6) for cha in channels] + + with open(path+f"eps={round(eps,3):.3f}_ei.txt", 'w', encoding='utf-8') as f: + json.dump(ei, f, indent=1) + + strength = 100000 + extremes = 10000 #[l//2 for l in lens] + ext_rat = extremes/(stp-100) + + circparts = 100 + pha_dev = m.pi/circparts + pha_max = np.max(np.abs(phase_diff)) + + phase_in_part = lambda ph, i: abs(ph)<=pha_max/circparts if i==0 else i*pha_max/circparts0 else [] for i in range(circparts)] + + pha_parts = [[phase_diff[i] for i in part] for part in dev_parts] + syn_parts = [[phase_diff[i] for i in part] for part in top_parts] + ran_sam = sample_parts(dev_parts, top_parts) + top_sam = sample_parts(top_parts, top_parts) + + dev = [] + for part in dev_parts: + dev.extend(part) + top = [] + for part in top_parts: + top.extend(part) + + assert sorted(top)==sorted(top_sam) + + ran_res = [] + top_res = [] + tot_res = [] + + bot_ei = [] + top_ei = [] + + maxdt = 250 + for dt in range(maxdt): + tot_pha = phase_diff[dt:dt-maxdt] + ran_pha = [(phase_diff[i+dt]) for i in ran_sam] + top_pha = [(phase_diff[i+dt]) for i in top_sam] + + tot_res.append( norm(resultant(tot_pha)) ) + ran_res.append( norm(resultant(ran_pha)) ) + top_res.append( norm(resultant(top_pha)) ) + + + sampling = 'samedist_varmaxpha' + plotpath = new_folder(path+f'{sampling}_causal_roll{pha_dev:.3f}/')#'bootstrap={strength}/' + savepath=new_folder(plotpath+f'extremes={extremes}_bootstrength={strength}/') + + newpath = new_folder(path+f'extremes={extremes}_bootstrength={strength}/') + resultant_stat = bootstrap_res(phase_diff,top,strength=strength) + + confs = [] + confdt = 250 + for dt in range(confdt+1): + pha_evs = [[phase_diff[i+dt] for i in part] for part in dev_parts] + syn_evs = [[phase_diff[i+dt] for i in part] for part in top_parts] + phas = [phase_diff[i+dt] for i in top] + res_pha = norm(resultant(phas)) + + conf = len([val for val in resultant_stat if val0 else None + to_sync = lambda i: True if abs(phase_diff[i])<0.08*m.pi else False if 0.42*m.pi0] + + print(f'{len(ei_ind)} states with positive EI') + + samples = choose(ei_ind, extremes) + sampling = 'allpos_ei' + + with open(path+f"eps={round(eps,3):.3f}_states.txt", 'r', encoding='utf-8') as f: + states = json.load(f)[100:] + with open(path+f"eps={round(eps,3):.3f}_coupling.txt", 'r', encoding='utf-8') as f: + coupling = json.load(f) + coupling = [tuple(edge) for edge in coupling] + + phase_pairs = [] + ei_pairs = [] + for num,i in enumerate(samples): + causal_init = states[i].translate(str.maketrans('', '', '.-=')) + + sim = Simulate4Layers( dim, + eps, + coupling=coupling, + init=causal_init, + noeffect=0, + steps=1, + draw=None, + ) + + activation = sim._activations() + channel = sim._channels() + + osc = list(zip(*activation)) + phase_abs = np.array([[np.arctan2(*act[::-1]) for act in osc[i]] for i in range(2)]) + phasediff_c = np.round(diff(phase_abs[0],phase_abs[1]),6) + ei_c = [round(np.sum(cha)*(1-H2(eps)),6) for cha in channel] + max_ei_c = max([np.sum(cha) for cha in channel]) + + sim = Simulate4Layers( dim, + eps, + coupling=coupling, + init=causal_init, + noeffect=-1, + steps=1, + draw=None, + ) + + activation = sim._activations() + channel = sim._channels() + + osc = list(zip(*activation)) + phase_abs = np.array([[np.arctan2(*act[::-1]) for act in osc[i]] for i in range(2)]) + phasediff_i = np.round(diff(phase_abs[0],phase_abs[1]),6) + ei_i = [round(np.sum(cha)*(1-H2(eps)),6) for cha in channel] + max_ei_i = max([np.sum(cha) for cha in channel]) + + phase_pairs.append((phasediff_i[-1], phasediff_c[-1])) + ei_pairs.append((ei_i[-1], ei_c[-1])) + + savepath = path + sampling + '/' + new_folder(savepath) + + if num%100==99: + print(f'Done {num:0{len(str(extremes))}d}') + + with open(savepath+f"eps={round(eps,3):.3f}_phase_pairs.txt", 'w', encoding='utf-8') as f: + json.dump(phase_pairs, f, indent=1) + with open(savepath+f"eps={round(eps,3):.3f}_ei_pairs.txt", 'w', encoding='utf-8') as f: + json.dump(ei_pairs, f, indent=1) + + phases_i, phases_c = zip(*phase_pairs) + ei_i, ei_c = zip(*ei_pairs) + + phase_space = np.linspace(0,m.pi,100+1) + ei_space = np.linspace(0,np.max([ei_i,ei_c]),100+1) + + phase_dist_i = [len([ph for ph in phases_i if low<=ph0 else [] for i in range(circparts)] + bot_parts = [dev_parts[i][:ext_parts[i]] if ext_parts[i]>0 else [] for i in range(circparts)] + + top = [] + for part in top_parts: + top.extend(part) + bot = [] + for part in bot_parts: + bot.extend(part) + + print(len(top), len(bot), extremes, 'equal?') + sampling = 'samedist_varmaxpha' + + # pha_center = av_diff + # pha_dev = m.pi/8 + + # from_sync = lambda i: True if abs(phase_diff[i])<0.08*m.pi else False if 0.42*m.pi0 else None + # to_sync = lambda i: True if abs(phase_diff[i])<0.08*m.pi else False if 0.42*m.pi=maxdt: + path_c = path+f'causal_maxdt={causal_maxdt}/eps={round(eps,3):.3f}/{i:0{len(str(stp))}}/' + else: + path_c = path+f'causal_maxdt={maxdt}/eps={round(eps,3):.3f}/{i:0{len(str(stp))}}/' + + causal_init = states[i].translate(str.maketrans('', '', '.-=')) + + try: + with open(path_c+f"eps={round(eps,3):.3f}_phase_diff.txt", 'r', encoding='utf-8') as f: + phasediff = json.load(f) + except: + sim=Simulate4Layers(dim, + eps, + coupling=coupling, + init=causal_init, + noeffect=0, + steps=maxdt, + draw=None, + save='simple', + path=path_c, + ) + + with open(path_c+f"eps={round(eps,3):.3f}_activation.txt", 'r', encoding='utf-8') as f: + activation = json.load(f) + + osc = list(zip(*activation)) + phase_abs = np.array([[np.arctan2(*act[::-1]) for act in osc[i]] for i in range(2)]) + phasediff = diff(phase_abs[0],phase_abs[1]) + + with open(path_c+f"eps={round(eps,3):.3f}_phase_diff.txt", 'w', encoding='utf-8') as f: + json.dump(list(phasediff), f, indent=1) + + for i in []:#top: + original_maxdt=0 + for file in os.listdir(path): + f = os.path.join(path, file) + if not os.path.isfile(f): + o_maxdt = file.replace('original_maxdt=','') + if c_maxdt != file: + original_maxdt = int(o_maxdt) + + if original_maxdt>=maxdt: + path_c = path+f'original_maxdt={original_maxdt}/eps={round(eps,3):.3f}/{i:0{len(str(stp))}}/' + else: + path_c = path+f'original_maxdt={maxdt}/eps={round(eps,3):.3f}/{i:0{len(str(stp))}}/' + + causal_init = states[i].translate(str.maketrans('', '', '.-=')) + + try: + with open(path_c+f"eps={round(eps,3):.3f}_phase_diff.txt", 'r', encoding='utf-8') as f: + phasediff = json.load(f) + except: + sim=Simulate4Layers(dim, + eps, + coupling=coupling, + init=causal_init, + noeffect=-1, + steps=maxdt, + draw=None, + save='simple', + path=path_c, + ) + + with open(path_c+f"eps={round(eps,3):.3f}_activation.txt", 'r', encoding='utf-8') as f: + activation = json.load(f) + + osc = list(zip(*activation)) + phase_abs = np.array([[np.arctan2(*act[::-1]) for act in osc[i]] for i in range(2)]) + phasediff = diff(phase_abs[0],phase_abs[1]) + phasediff = [round(pha,6) for pha in phasediff] + + with open(path_c+f"eps={round(eps,3):.3f}_phase_diff.txt", 'w', encoding='utf-8') as f: + json.dump(list(phasediff), f, indent=1) + + + + bot_res = [] + top_res = [] + dis_res = [] + tot_res = [] + + bot_ph = [] + top_ph = [] + dis_ph = [] + tot_ph = [] + + bot_ei = [] + top_ei = [] + dis_ei = [] + tot_ei = [] + + + present_dt=0 + for file in os.listdir(path): + f = os.path.join(path, file) + if not os.path.isfile(f): + p_maxdt = file.replace(f'{sampling}_causal_roll{pha_dev:.3f}_maxdt=','') + if p_maxdt != file: + present_dt = int(p_maxdt) + + if present_dt>0: + try: + datapath = path + f"{sampling}_causal_roll{pha_dev:.3f}_maxdt={present_dt}/data/" + with open(datapath+f"eps={round(eps,3):.3f}_bot_dia_res.txt", 'r', encoding='utf-8') as f: + bot_res = json.load(f) + with open(datapath+f"eps={round(eps,3):.3f}_top_dia_res.txt", 'r', encoding='utf-8') as f: + top_res = json.load(f) + with open(datapath+f"eps={round(eps,3):.3f}_dis_dia_res.txt", 'r', encoding='utf-8') as f: + dis_res = json.load(f) + with open(datapath+f"eps={round(eps,3):.3f}_bot_dia_ei.txt", 'r', encoding='utf-8') as f: + bot_ei = json.load(f) + with open(datapath+f"eps={round(eps,3):.3f}_top_dia_ei.txt", 'r', encoding='utf-8') as f: + top_ei = json.load(f) + with open(datapath+f"eps={round(eps,3):.3f}_dis_dia_ei.txt", 'r', encoding='utf-8') as f: + dis_ei = json.load(f) + except: + present_dt=0 + bot_res = [] + top_res = [] + dis_res = [] + + for dt in range(present_dt,maxdt): + top_pha = [(phase_diff[i+dt]) for i in top] + bot_pha = [(phase_diff[i+dt]) for i in bot] + + dis = [] + for i in top: + causal_maxdt=0 + for file in os.listdir(path): + f = os.path.join(path, file) + if not os.path.isfile(f): + c_maxdt = file.replace('causal_maxdt=','') + if c_maxdt != file: + causal_maxdt = int(c_maxdt) + + if causal_maxdt>=maxdt: + path_c = path+f'causal_maxdt={causal_maxdt}/eps={round(eps,3):.3f}/{i:0{len(str(stp))}}/' + else: + path_c = path+f'causal_maxdt={maxdt}/eps={round(eps,3):.3f}/{i:0{len(str(stp))}}/' + + with open(path_c+f'eps={round(eps,3):.3f}_phase_diff.txt', 'r', encoding='utf-8') as f: + dis.append((json.load(f)[dt])) + + dei = [] + for i in top: + causal_maxdt=0 + for file in os.listdir(path): + f = os.path.join(path, file) + if not os.path.isfile(f): + c_maxdt = file.replace('causal_maxdt=','') + if c_maxdt != file: + causal_maxdt = int(c_maxdt) + + if causal_maxdt>=maxdt: + path_c = path+f'causal_maxdt={causal_maxdt}/eps={round(eps,3):.3f}/{i:0{len(str(stp))}}/' + else: + path_c = path+f'causal_maxdt={maxdt}/eps={round(eps,3):.3f}/{i:0{len(str(stp))}}/' + + with open(path_c+f'eps={round(eps,3):.3f}_channels.txt', 'r', encoding='utf-8') as f: + dei.append(np.sum(json.load(f)[dt])*(1-H2(eps))) + + dis_pha = dis + tot_pha = (phase_diff[dt:dt-maxdt]) + + bot_res.append( [norm(resultant(bot_pha))] ) + top_res.append( [norm(resultant(top_pha))] ) + dis_res.append( [norm(resultant(dis_pha))] ) + tot_res.append( norm(resultant(tot_pha)) ) + + # bot_ph.append( [phase(*resultant(bot_pha[i])) for i in range(1)] ) + # top_ph.append( [phase(*resultant(top_pha[i])) for i in range(1)] ) + # dis_ph.append( [phase(*resultant(dis_pha[i])) for i in range(1)] ) + # tot_ph.append( phase(*resultant(tot_pha)) ) + + bot_ei.append( [np.average([ei[i+dt] for i in bot])] ) + top_ei.append( [np.average([ei[i+dt] for i in top])] ) + dis_ei.append( [np.average(dei)] ) + tot_ei.append( np.average(ei[dt:dt-maxdt]) ) + + if dt%10==9: + print(f'Done dt={dt:{len(str(maxdt))}d}') + + plotpath = path+f'{sampling}_causal_roll{pha_dev:.3f}_maxdt={maxdt}/' + new_folder(plotpath+'data/') + + with open(plotpath+f"data/eps={round(eps,3):.3f}_bot_dia_res.txt", 'w', encoding='utf-8') as f: + json.dump(list(bot_res), f, indent=1) + with open(plotpath+f"data/eps={round(eps,3):.3f}_top_dia_res.txt", 'w', encoding='utf-8') as f: + json.dump(list(top_res), f, indent=1) + with open(plotpath+f"data/eps={round(eps,3):.3f}_dis_dia_res.txt", 'w', encoding='utf-8') as f: + json.dump(list(dis_res), f, indent=1) + with open(plotpath+f"data/eps={round(eps,3):.3f}_bot_dia_ei.txt", 'w', encoding='utf-8') as f: + json.dump(list(bot_ei), f, indent=1) + with open(plotpath+f"data/eps={round(eps,3):.3f}_top_dia_ei.txt", 'w', encoding='utf-8') as f: + json.dump(list(top_ei), f, indent=1) + with open(plotpath+f"data/eps={round(eps,3):.3f}_dis_dia_ei.txt", 'w', encoding='utf-8') as f: + json.dump(list(dis_ei), f, indent=1) + + bot_res = list(zip(*bot_res)) + top_res = list(zip(*top_res)) + dis_res = list(zip(*dis_res)) + + # bot_ph = list(zip(*bot_ph)) + # top_ph = list(zip(*top_ph)) + # dis_ph = list(zip(*dis_ph)) + + bot_ei = list(zip(*bot_ei)) + top_ei = list(zip(*top_ei)) + dis_ei = list(zip(*dis_ei)) + + + + qtplot(f'Diachronic resultant for dim={dim} eps={eps:.3f} with 4 layers', + [np.array(range(maxdt))]*4, + [tot_res, bot_res[0], dis_res[0], top_res[0]], + ['Average Resultant', f'bottom {extremes} ei' + f'dis {extremes} ei', f'top {extremes} ei'], + x_tag = 'dt', + y_tag = 'concentration', + export=True, + path=plotpath, + filename=f'All Diachronic Resultant Norm eps={round(eps,3):.3f} dim={dim} extremes={extremes}.png', + close=True) + + # qtplot(f'Diachronic resultant phase for dim={dim} eps={eps:.3f} with 4 layers', + # [np.array(range(maxdt))]*3, + # [top_ph[0], dis_ph[0], tot_ph], + # ['top {extremes} ei', 'dis {extremes} ei', + # 'Average'], + # x_tag = 'dt', + # y_tag = 'phase', + # export=True, + # path=plotpath, + # filename=f'All Diachronic Resultant Phase eps={round(eps,3):.3f} dim={dim} extremes={extremes} roll{pha_dev:.3f}.png', + # close=True) + + # qtplot(f'Diachronic ei for dim={dim} with 4 layers', + # [np.array(range(maxdt))]*4, + # [bot_ei, top_ei, dev_ei, tot_ei], + # ['EI ev of bottom {extremes} ei', 'EI ev of top {extremes} ei', + # 'EI ev of phase filtered ei', 'Average EI'], + # x_tag = 'dt', + # y_tag = 'average ei', + # export=True, + # path=path+'plots/', + # filename=f'Diachronic EI balanced for eps={round(eps,3):.3f} dim={dim} extremes={extremes} roll{pha_dev:.3f}.png', + # close=True) + + print(f'Done eps={eps:.3f} with dim={dim} at {datetime.now()}') + + # qtplot(f'Resultant and EI evolution for dim={dim} with 4 layers', + # [[0]+eps_space]*2, + # [max(av_ei)*diff_res, av_ei], + # ['Resultant', 'avEI'], + # export=True, + # path=path, + # filename=f'Resultant and EI for dim={dim}.png', + # close=True) + + \ No newline at end of file diff --git a/evaluation/4Layer Res+EI.py b/evaluation/4Layer Res+EI.py index 1f143f9..5c2c76e 100644 --- a/evaluation/4Layer Res+EI.py +++ b/evaluation/4Layer Res+EI.py @@ -14,6 +14,7 @@ from datetime import datetime from plot import qtplot eps_space = list(np.linspace(0.005,0.5,100)) +eps_space= eps_space[1::2] def resultant(sample): phase_x = [m.cos(ind) for ind in sample] @@ -24,37 +25,68 @@ def resultant(sample): def H2(x): return -x*m.log2(x)-(1-x)*m.log2(1-x) -for dim in [7]: - diff_res = [0] +dims=list(range(3,10)) + +R=[] +EI=[] +for dim in dims: + cons = [(n,(2*n+m)%dim) for n in range(dim) for m in range(0,dim-2,3)] + path=f'/cloud/Public/_data/neuropercolation/4lay/cons=dimtimesdimby3_steps=100100/dim={dim:02d}_cons={len(cons)}/batch=0/' + + diff_res = [1] av_ei = [0] for eps in eps_space: - path=f'/cloud/Public/_data/neuropercolation/4lay/steps=100000/dim={dim:02}/' + try: + with open(path+f"eps={round(eps,3):.3f}_phase_diff.txt", 'r', encoding='utf-8') as f: + phase_diff = json.load(f)[100:] + with open(path+f"eps={round(eps,3):.3f}_ei.txt", 'r', encoding='utf-8') as f: + ei = json.load(f)[100:] + except: + with open(path+f"eps={round(eps,3):.3f}_activation.txt", 'r', encoding='utf-8') as f: + activation = json.load(f)[100:] + with open(path+f"eps={round(eps,3):.3f}_channels.txt", 'r', encoding='utf-8') as f: + channels = json.load(f)[100:] - with open(path+f"eps={round(eps,3):.3f}_activation.txt", 'r', encoding='utf-8') as f: - activation = json.load(f)[100:] - - osc = list(zip(*activation)) - phase = np.array([[np.arctan2(*act[::-1]) for act in osc[i]] for i in range(2)]) - phase_diff = (phase[1]-phase[0]+m.pi)%(2*m.pi)-m.pi + osc = list(zip(*activation)) + phase = np.array([[np.arctan2(*act[::-1]) for act in osc[i]] for i in range(2)]) + phase_diff = (phase[1]-phase[0]+m.pi)%(2*m.pi)-m.pi + ei = [np.sum(cha)*(1-H2(eps)) for cha in channels] res = np.linalg.norm(resultant(phase_diff)) diff_res.append(res) - with open(path+f"eps={round(eps,3):.3f}_channels.txt", 'r', encoding='utf-8') as f: - channels = json.load(f)[100:] - - ei = [np.sum(cha)*(1-H2(eps)) for cha in channels] av_ei.append(np.average(ei)) + with open(path+f"eps={round(eps,3):.3f}_phase_diff.txt", 'w', encoding='utf-8') as f: + json.dump(list(phase_diff), f, indent=1) + + with open(path+f"eps={round(eps,3):.3f}_ei.txt", 'w', encoding='utf-8') as f: + json.dump(ei, f, indent=1) print(f'Done eps={eps:.3f} with dim={dim} at {datetime.now()}') - - qtplot(f'Resultant and EI evolution for dim={dim} with 4 layers', - [[0]+eps_space]*2, - [diff_res, av_ei], - ['Resultant', 'avEI'], - export=True, - path=path, - filename=f'Resultant and EI for dim={dim}.png', - close=True) + R.append(diff_res) + EI.append(av_ei) + + #%% +savepath=f'/cloud/Public/_data/neuropercolation/4lay/cons=dimtimesdimby3_steps=100100/plots/' + +qtplot(f'Total resultant of phase-difference', + [[0]+eps_space]*len(dims), + R[::-1], + [f'{dim}x{dim}' for dim in dims[::-1]], + y_tag = 'concentration', + export=False, + path=savepath, + filename=f'Resultant_ev.png', + close=False) + +qtplot(f'Average Effect Integrated Information', + [[0]+eps_space]*len(dims), + EI[::-1], + [f'{dim}x{dim}' for dim in dims[::-1]], + y_tag='integrated information', + export=False, + path=savepath, + filename=f'EI_ev.png', + close=False) \ No newline at end of file diff --git a/evaluation/4Layer Resultant 4way Phase Evolution.py b/evaluation/4Layer Resultant 4way Phase Evolution.py new file mode 100644 index 0000000..876a03a --- /dev/null +++ b/evaluation/4Layer Resultant 4way Phase Evolution.py @@ -0,0 +1,277 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Mon Aug 21 14:59:22 2023 + +@author: astral +""" + +import json +import math as m +import numpy as np +from numpy.linalg import norm +from datetime import datetime +from random import sample as choose + +from plot import qtplot + +eps_space = list(np.linspace(0.01,0.2,20)) + +def resultant(sample): + phase_x = [m.cos(ind) for ind in sample] + phase_y = [m.sin(ind) for ind in sample] + + return (np.average(phase_x), np.average(phase_y)) + +phase = np.vectorize(lambda x,y: (m.atan2(y,x)+m.pi)%(2*m.pi)-m.pi) +diff = np.vectorize(lambda x,y: (y-x+m.pi)%(2*m.pi)-m.pi) +H2 = lambda x: -x*m.log2(x)-(1-x)*m.log2(1-x) + +extremes = None +maxdt = 300 + +for dim in [9]: + for eps in eps_space: + path=f'/cloud/Public/_data/neuropercolation/4lay/steps=100100/dim=09/' + + try: + with open(path+f"eps={round(eps,3):.3f}_phase_diff.txt", 'r', encoding='utf-8') as f: + phase_diff = json.load(f) + except: + with open(path+f"eps={round(eps,3):.3f}_activation.txt", 'r', encoding='utf-8') as f: + activation = json.load(f)[100:] + + osc = list(zip(*activation)) + phase_abs = np.array([[np.arctan2(*act[::-1]) for act in osc[i]] for i in range(2)]) + phase_diff = diff(phase_abs[0],phase_abs[1]) + + with open(path+f"eps={round(eps,3):.3f}_phase_diff.txt", 'w', encoding='utf-8') as f: + json.dump(list(phase_diff), f, indent=1) + + + all_res = norm(resultant(phase_diff)) + av_diff = np.arccos(all_res) + + try: + with open(path+f"eps={round(eps,3):.3f}_ei.txt", 'r', encoding='utf-8') as f: + ei = json.load(f) + except: + with open(path+f"eps={round(eps,3):.3f}_channels.txt", 'r', encoding='utf-8') as f: + channels = json.load(f)[100:] + + ei = [np.sum(cha)*(1-H2(eps)) for cha in channels] + + with open(path+f"eps={round(eps,3):.3f}_ei.txt", 'w', encoding='utf-8') as f: + json.dump(ei, f, indent=1) + + pha_center = av_diff + pha_dev = m.pi/32 + + from_sync = lambda i: True if abs(phase_diff[i])<0.08*m.pi else False if 0.42*m.pi0 else None + to_sync = lambda i: True if abs(phase_diff[i])<0.08*m.pi else False if 0.42*m.pi0 else None + to_sync = lambda i: True if abs(phase_diff[i])<0.08*m.pi else False if 0.42*m.pi0] @@ -108,19 +111,22 @@ for dim in [7]: dev_ei = [] tot_ei = [] for dt in range(maxdt): - bot_pha = [abs(phase_diff[i+dt]) for i in bot_ind] - top_pha = [abs(phase_diff[i+dt]) for i in top_ind] - dev_pha = [abs(phase_diff[i+dt]) for i in dev_ind] + bot_pha = [[abs(phase_diff[i+dt]) for i in bot_in_ind], + [abs(phase_diff[i+dt]) for i in bot_out_ind]] + top_pha = [[abs(phase_diff[i+dt]) for i in top_in_ind], + [abs(phase_diff[i+dt]) for i in top_out_ind]] + dev_pha = [[abs(phase_diff[i+dt]) for i in dev_in_ind], + [abs(phase_diff[i+dt]) for i in dev_out_ind]] tot_pha = np.abs(phase_diff[dt:dt-maxdt]) - bot_res.append( norm(resultant(bot_pha)) ) - top_res.append( norm(resultant(top_pha)) ) - dev_res.append( norm(resultant(dev_pha)) ) + bot_res.append( [norm(resultant(bot_pha[i])) for i in [0,1]] ) + top_res.append( [norm(resultant(top_pha[i])) for i in [0,1]] ) + dev_res.append( [norm(resultant(dev_pha[i])) for i in [0,1]] ) tot_res.append( norm(resultant(tot_pha)) ) - bot_ph.append( phase(*resultant(bot_pha)) ) - top_ph.append( phase(*resultant(top_pha)) ) - dev_ph.append( phase(*resultant(dev_pha)) ) + bot_ph.append( [phase(*resultant(bot_pha[i])) for i in [0,1]] ) + top_ph.append( [phase(*resultant(top_pha[i])) for i in [0,1]] ) + dev_ph.append( [phase(*resultant(dev_pha[i])) for i in [0,1]] ) tot_ph.append( phase(*resultant(tot_pha)) ) bot_ei.append( np.average([ei[i+dt] for i in bot_ind]) ) @@ -130,29 +136,59 @@ for dim in [7]: if dt%100==99: print(f'Done dt={dt}') - + + bot_res = list(zip(*bot_res)) + top_res = list(zip(*top_res)) + dev_res = list(zip(*dev_res)) + + bot_ph = list(zip(*bot_ph)) + top_ph = list(zip(*top_ph)) + dev_ph = list(zip(*dev_ph)) + qtplot(f'Diachronic resultant for dim={dim} with 4 layers', [np.array(range(maxdt))]*4, - [bot_res, top_res, dev_res, tot_res], - ['Resultant ev of bottom {extremes} ei', 'Resultant ev of top {extremes} ei', - 'Resultant ev of phase filtered ei', 'Average Resultant'], + [bot_res[0], top_res[0], dev_res[0], tot_res], + ['in bottom {extremes} ei', 'in top {extremes} ei', + 'in all filtered {len(dev_in_ind)} ei', 'Average Resultant'], x_tag = 'dt', y_tag = 'concentration', export=True, - path=path+'plots/', - filename=f'Diachronic Resultant Norm balanced eps={round(eps,3):.3f} dim={dim} extremes={extremes} roll{pha_dev:.3f}.png', + path=path+'inout/', + filename=f'Diachronic Resultant Norm in eps={round(eps,3):.3f} dim={dim} extremes={extremes} roll{pha_dev:.3f}.png', + close=True) + qtplot(f'Diachronic resultant for dim={dim} with 4 layers', + [np.array(range(maxdt))]*4, + [bot_res[1], top_res[1], dev_res[1], tot_res], + ['out bottom {extremes} ei', 'out top {extremes} ei', + 'out all filtered {len(dev_out_ind)} ei', 'Average Resultant'], + x_tag = 'dt', + y_tag = 'concentration', + export=True, + path=path+'inout/', + filename=f'Diachronic Resultant Norm out eps={round(eps,3):.3f} dim={dim} extremes={extremes} roll{pha_dev:.3f}.png', close=True) qtplot(f'Diachronic resultant phase for dim={dim} with 4 layers', [np.array(range(maxdt))]*4, - [bot_ph, top_ph, dev_ph, tot_ph], - ['bottom {extremes} ei', 'top {extremes} ei', - 'all filtered ei', 'Average'], + [bot_ph[0], top_ph[0], dev_ph[0], tot_ph], + ['in bottom {extremes} ei', 'in top {extremes} ei', + 'in all filtered {len(dev_in_ind)} ei', 'Average'], x_tag = 'dt', - y_tag = 'concentration', + y_tag = 'phase', export=True, - path=path+'plots/', - filename=f'Diachronic Resultant Phase balanced eps={round(eps,3):.3f} dim={dim} extremes={extremes} roll{pha_dev:.3f}.png', + path=path+'inout/', + filename=f'Diachronic Resultant Phase in eps={round(eps,3):.3f} dim={dim} extremes={extremes} roll{pha_dev:.3f}.png', + close=True) + qtplot(f'Diachronic resultant phase for dim={dim} with 4 layers', + [np.array(range(maxdt))]*4, + [bot_ph[1], top_ph[1], dev_ph[1], tot_ph], + ['out bottom {extremes} ei', 'out top {extremes} ei', + 'out all filtered {len(dev_out_ind)} ei', 'Average'], + x_tag = 'dt', + y_tag = 'phase', + export=True, + path=path+'inout/', + filename=f'Diachronic Resultant Phase out eps={round(eps,3):.3f} dim={dim} extremes={extremes} roll{pha_dev:.3f}.png', close=True) # qtplot(f'Diachronic ei for dim={dim} with 4 layers', diff --git a/evaluation/4Layer Resultant.py b/evaluation/4Layer Resultant.py new file mode 100644 index 0000000..3dbcb42 --- /dev/null +++ b/evaluation/4Layer Resultant.py @@ -0,0 +1,106 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Mon Aug 21 14:59:22 2023 + +@author: astral +""" + +import os +import json +import math as m +import numpy as np +from numpy.linalg import norm +from datetime import datetime +from random import sample as choose + +from plot import qtplot + +eps_space = list(np.linspace(0.005,0.5,100)) + +def resultant(sample): + phase_x = [m.cos(ind) for ind in sample] + phase_y = [m.sin(ind) for ind in sample] + + return (np.average(phase_x), np.average(phase_y)) + +def H2(x): + return -x*m.log2(x)-(1-x)*m.log2(1-x) +def new_folder(p): + if not os.path.exists(p): + os.makedirs(p) + return p + +extremes = None +maxdt = 300 + +path=f'/cloud/Public/_data/neuropercolation/4lay/steps=100000_diares_manydim/' + + +dims = list(range(7,9)) + +resus = [] +phalocks = [] +for dim in dims: + resu = [] + phalock = [] + for eps in eps_space: + dimpath=f'/cloud/Public/_data/neuropercolation/4lay/steps=100000_diares_manydim/dim={dim:02d}/' + + try: + with open(dimpath+f"eps={round(eps,3):.3f}_phase_diff.txt", 'r', encoding='utf-8') as f: + phase_diff = json.load(f) + except: + with open(dimpath+f"eps={round(eps,3):.3f}_activation.txt", 'r', encoding='utf-8') as f: + activation = json.load(f)[100:] + + osc = list(zip(*activation)) + phase = np.array([[np.arctan2(*act[::-1]) for act in osc[i]] for i in range(2)]) + phase_diff = (phase[1]-phase[0]+m.pi)%(2*m.pi)-m.pi + + with open(dimpath+f"eps={round(eps,3):.3f}_phase_diff.txt", 'w', encoding='utf-8') as f: + json.dump(list(phase_diff), f, indent=1) + + + all_res = norm(resultant(phase_diff)) + av_diff = np.arccos(all_res) + all_pha = np.arctan2(*resultant(phase_diff)[::-1]) + + resu.append(all_res) + phalock.append(all_pha) + + plotpath = new_folder(dimpath + 'resplot/') + with open(plotpath+f"resultants.txt", 'w', encoding='utf-8') as f: + json.dump(list(resu), f, indent=1) + with open(plotpath+f"mainphase.txt", 'w', encoding='utf-8') as f: + json.dump(list(phalock), f, indent=1) + + resus.append(resu) + phalocks.append(phalock) + +plotspath = new_folder(path + 'resplot/') + +#%% +qtplot(f'Resultant evolution with 4 layers', + [eps_space]*len(dims), + resus, + [f'Resultant for size {dim}x{dim}' for dim in dims], + x_tag = 'epsilon', + y_tag = 'resultant', + lw=3, + export=True, + path=plotspath, + filename=f'resultant.png', + close=False) + +qtplot(f'Dominant phase evolution with 4 layers', + [eps_space]*len(dims), + phalocks, + [f'dominant phase for size {dim}x{dim}' for dim in dims], + x_tag = 'epsilon', + y_tag = 'resultant', + lw=3, + export=True, + path=plotspath, + filename=f'mainphase.png', + close=True) diff --git a/evaluation/4Layer Significance.py b/evaluation/4Layer Significance.py new file mode 100644 index 0000000..62371bd --- /dev/null +++ b/evaluation/4Layer Significance.py @@ -0,0 +1,25 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Sun Sep 10 07:07:28 2023 + +@author: astral +""" + +import sys as _sys +from sign import significance as _sign +import numpy as _np + +import gc as _gc + +_dim = 9 +_epsilons = _np.linspace(0.01,0.20,20) + +for _eps in _epsilons[_sys.argv[1]]: + _sign(_dim,_eps) + + for name in dir(): + if not name.startswith('_'): + del globals()[name] + + _gc.collect() \ No newline at end of file diff --git a/evaluation/4Layer phiplot.py b/evaluation/4Layer phiplot.py new file mode 100644 index 0000000..46d5345 --- /dev/null +++ b/evaluation/4Layer phiplot.py @@ -0,0 +1,98 @@ +# -*- coding: utf-8 -*- +""" +Created on Tue Aug 30 14:25:12 2022 + +@author: timof +""" + +import sys +import os +import json + +from plot import qtplot + + +import math as m +import numpy as np + + +vect = np.vectorize + +@vect +def log2(x): + try: + return m.log2(x) + except ValueError: + if x==0: + return float(0) + else: + raise + +def new_folder(path): + if not os.path.exists(path): + os.makedirs(path) + return path + +phase = np.vectorize(lambda x,y: (m.atan2(y,x)+m.pi)%(2*m.pi)-m.pi) +diff = np.vectorize(lambda x,y: (y-x+m.pi)%(2*m.pi)-m.pi) +H2 = lambda x: -x*m.log2(x)-(1-x)*m.log2(1-x) + + +path = '/cloud/Public/_data/neuropercolation/4lay/cons=dimtimesdimby3_steps=100100/' +suffix = '' + +chi = chr(967) +vareps = chr(949) +varphi = chr(981) + +vals = [[],[]] + +runsteps = 1000100 + +eps_space = np.linspace(0.005, 0.5, 100) +eps_space = eps_space[1::2] + +dims = list(range(3,10))#+[16,49] + +mode='density' +ma=[] +s=[] +k=[] +mk=[] +PHI=[] +lastkurt=None +for dim in dims: + phis=[] + con_gap = 3 + cons = [(n,(2*n+m)%dim) for n in range(dim) for m in range(0,dim-2,con_gap)] + dimpath = new_folder(path + f'dim={dim:02}_cons={len(cons)}/batch=0/') + for epsilon in eps_space: + try: + with open(dimpath+f"eps={round(epsilon,3):.3f}_ei.txt", 'r', encoding='utf-8') as f: + ei = np.array(json.load(f)[100:]) + except: + print('Calcing phi') + with open(dimpath+f"eps={round(epsilon,3):.3f}_channels.txt", 'r', encoding='utf-8') as f: + channels = np.array(json.load(f)[100:]) + + ei = channels*(1-H2(epsilon)) + + with open(dimpath+f"eps={round(epsilon,3):.3f}_ei.txt", 'r', encoding='utf-8') as f: + json.dump(ei,f,indent=1) + + phi=np.average(ei) + phis.append(phi) + PHI.append(phis) + #%% +qtplot(f"Mean effect integration over noise level", + [eps_space]*len(dims), + PHI[::-1], + [f'dim={dim:02d}x{dim:02d}' for dim in dims[::-1]], + y_tag = f'effect integration {varphi}', + export=True, + path=dimpath+"", + filename=f'eps={round(epsilon,3):.3f}_evolution.png', + close=False) + +mode = 'density' +#%% diff --git a/evaluation/4laysign.py b/evaluation/4laysign.py new file mode 100644 index 0000000..a09801c --- /dev/null +++ b/evaluation/4laysign.py @@ -0,0 +1,173 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Sun Sep 10 07:07:28 2023 + +@author: astral +""" + +import sys as _sys +from sign import sampling as _samp +from sign import plotting as _plot +from sign import full_stats as _stats +import numpy as _np +import json +from plot import qtplot + +import gc as _gc + +_dim = 9 +_epsilons = _np.linspace(0.005,0.20,40) +_noeffects = list(range(1,41,5)) +try: + _exeps = int(_sys.argv[1]) + _epses=_epsilons[_exeps-1:_exeps] + _nf = int(_sys.argv[2]) + _nfs = _noeffects[_nf-1:_nf] +except: + _epses=_epsilons + _nfs = list(range(1,41,5)) +_samplemethod = 'allpos_ei' +_ext = 5000 +_dt=40 + +path='/cloud/Public/_data/neuropercolation/4lay/cons=27-knight_steps=1000100/dim=09/batch=0/' +plotpath = path + _samplemethod + f'_samples={_ext}/' +savepath = path + _samplemethod + f'_samples={_ext}/dt={_dt}/' + +meanlist = [] +stdlist = [] +integrallist = [] +cohendlist = [] +norm_p_list = [] +ttest_p_list = [] +sign_p_list = [] + + +resultant_i = [] +resultant_c = [] +for _eps in _epses[11:21]: + res_evol_i = [] + res_evol_c = [] + for nf in _nfs: + #_samp(_dim,_eps,_samplemethod,_ext,_dt,noeff=nf) + resultant_i = [] + resultant_c = [] + for t in range(1,_dt): + res_i, res_c = _stats(_dim,_eps,_samplemethod,_ext,t,noeff=nf,ret='phases') + resultant_i.append(res_i) + resultant_c.append(res_c) + + res_evol_i.append(resultant_i) + res_evol_c.append(resultant_c) + # mean,std,integral,cohend,norm_p,ttest_p,sign_p = _stats(_dim,_eps,_samplemethod,_ext,1) + + # meanlist.append(mean) + # stdlist.append(std) + # integrallist.append(integral) + # cohendlist.append(-cohend) + # norm_p_list.append(norm_p) + # ttest_p_list.append(ttest_p) + # sign_p_list.append(sign_p) + + + # with open(savepath+f"phasediff_means.txt", 'w', encoding='utf-8') as f: + # json.dump(meanlist, f, indent=1) + # with open(savepath+f"phasediff_stds.txt", 'w', encoding='utf-8') as f: + # json.dump(stdlist, f, indent=1) + # with open(savepath+f"phasediff_integrals.txt", 'w', encoding='utf-8') as f: + # json.dump(integrallist, f, indent=1) + # with open(savepath+f"phasediff_cohends.txt", 'w', encoding='utf-8') as f: + # json.dump(cohendlist, f, indent=1) + # with open(savepath+f"phasediff_norm_ps.txt", 'w', encoding='utf-8') as f: + # json.dump(norm_p_list, f, indent=1) + # with open(savepath+f"phasediff_ttest_ps.txt", 'w', encoding='utf-8') as f: + # json.dump(ttest_p_list, f, indent=1) + # with open(savepath+f"phasediff_sign_ps.txt", 'w', encoding='utf-8') as f: + # json.dump(sign_p_list, f, indent=1) + +# stdlowlist = [meanlist[eps] - stdlist[eps] for eps in range(len(meanlist))] +# stdhighlist = [meanlist[eps] + stdlist[eps] for eps in range(len(meanlist))] + +# norm_conf = [1-p for p in norm_p_list] +# ttest_conf = [1-p for p in ttest_p_list] +# sign_conf = [1-p for p in sign_p_list] + res_tot_i = _np.average(res_evol_i,axis=0) + res_evol_c.append(res_tot_i) + qtplot(f'Resultant reduction disconnect eps={_eps} dim={_dim} with 4 layers', + [list(range(_dt-1))]*(len(_nfs)+1), + res_evol_c, + [f'{nf} disconnected steps' for nf in _nfs]+['connected steps'], + x_tag = 'dt', + y_tag = 'resultant', + export=True, + path=plotpath, + filename=f'Resultant reduction disconnect eps={_eps} dim={_dim} extremes={_ext}.png', + close=True) + +# qtplot(f'Mean causal phase reduction for dt={_dt} dim={_dim} with 4 layers', +# [_epsilons]*3, +# [stdlowlist, stdhighlist, meanlist], +# ['Low standard deviation', +# 'High standard deviation', +# 'Mean'], +# colors=['r','r','g'], +# x_tag = f'noise level {chr(949)}', +# y_tag = 'abs phase reduction', +# export=True, +# path=savepath, +# filename=f'Mean phase reduction dim={_dim} extremes={_ext}.png', +# close=True) +# qtplot(f'Phase reduction probability for dt={_dt} dim={_dim} with 4 layers', +# [_epsilons], +# [integrallist], +# ['probability of phase reduction'], +# x_tag = f'noise level {chr(949)}', +# y_tag = 'abs phase reduction', +# export=True, +# path=savepath, +# filename=f'Probability of phase reduction dim={_dim} extremes={_ext}.png', +# close=True) +# qtplot(f'Effect size phase reduction for dt={_dt} dim={_dim} with 4 layers', +# [_epsilons], +# [cohendlist], +# ["Absolute Cohen's d (negative)"], +# x_tag = f'noise level {chr(949)}', +# y_tag = 'effect size', +# x_range = (0,0.2), +# # y_range = (0.05,1), +# y_log = True, +# export=True, +# path=savepath, +# filename=f'Effect size phase reduction dim={_dim} extremes={_ext}.png', +# close=False) +# qtplot(f'Test p-values for dt={_dt} dim={_dim} with 4 layers', +# [_epsilons]*5, +# [[0.001 for eps in _epsilons],[0.0001 for eps in _epsilons], +# norm_p_list,ttest_p_list,sign_p_list], +# ['0.001 limit','0.0001 limit', +# 'p-value for normality', 'p-value for mean difference', 'p-value for median'], +# x_tag = f'noise level {chr(949)}', +# y_tag = 'p-value', +# y_range = (0,0.01), +# lw=2, +# export=True, +# path=savepath, +# filename=f'P-values phase reduction dim={_dim} extremes={_ext}.png', +# close=True) +# qtplot(f'Test confidences for dt={_dt} dim={_dim} with 4 layers', +# [_epsilons]*5, +# [[0.999 for eps in _epsilons],[0.9999 for eps in _epsilons], +# norm_conf,ttest_conf,sign_conf], +# ['0.999 limit', '0.999 limit', +# 'confidence for normality', 'confidence for mean difference', 'confidence for median'], +# x_tag = f'noise level {chr(949)}', +# y_tag = 'confidence´', +# y_range = (0.99,1), +# lw=2, +# export=True, +# path=savepath, +# filename=f'Confidences phase reduction dim={_dim} extremes={_ext}.png', +# close=True) + + diff --git a/evaluation/phi.py b/evaluation/phi.py new file mode 100644 index 0000000..bcce633 --- /dev/null +++ b/evaluation/phi.py @@ -0,0 +1,173 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Wed Sep 27 04:39:54 2023 + +@author: astral +""" +import os +import json +import math as m +import numpy as np +from numpy.linalg import norm +from datetime import datetime +from random import sample as choose +from random import random +from numba import jit, njit, prange + +from plot import qtplot + +from neuropercolation import Simulate4Layers + +eps_space = list(np.linspace(0.01,0.2,20)) +def new_folder(path): + if not os.path.exists(path): + os.makedirs(path) + return path + +phase = np.vectorize(lambda x,y: (m.atan2(y,x)+m.pi)%(2*m.pi)-m.pi) +diff = np.vectorize(lambda x,y: (y-x+m.pi)%(2*m.pi)-m.pi) +H2 = lambda x: -x*m.log2(x)-(1-x)*m.log2(1-x) + +@njit +def neighbor(digit0, digit1, lenght): + layer = int(lenght) + dim = int(np.sqrt(layer)) + digit0, digit1 = np.array([digit0%dim, digit0//dim]), np.array([digit1%dim, digit1//dim]) + #print(digit0,digit1) + coord_dif = list(map(abs,digit1 - digit0)) + layer_nbor = 0 in coord_dif and len(set([1,dim-1]).intersection(set(coord_dif))) != 0 + #print(coord_dif, set([1,dim-1]).intersection(set(coord_dif))) + if layer_nbor: + return True + else: + return False + +@njit +def kcomb(zp,zm): + if zp>2: + val=1 + elif zm>2: + val=0 + elif zp==zm: + val=0.5 + elif zm==2: + val=0.5**(3-zp) + elif zp==2: + val=1-0.5**(3-zm) + elif zm==0 and zp==1: + val=9/16 + elif zp==0 and zm==1: + val=7/16 + else: + raise NotImplementedError(zp,zm) + return val + +path = new_folder('/cloud/Public/_data/neuropercolation/1lay/mips/') + + +def phi(dim,statestr,partstr,eps): + length = dim**2 + eta = 1-eps + # statestr=statestr.translate(str.maketrans('','','.-=')) + + state = np.array([int(q) for q in statestr]) + state = list(state.reshape((dim,dim))) + state = [list([int(cell) for cell in row]) for row in state] + + part = np.array([int(p) for p in partstr]) + part = list(part.reshape((dim,dim))) + part = [list([int(cell) for cell in row]) for row in part] + + inp = [[q+sum([state[(i+1)%dim][j], + state[(i-1)%dim][j], + state[i][(j+1)%dim], + state[i][(j-1)%dim] + ]) for j,q in enumerate(row)] for i,row in enumerate(state)] + + beps = [[int(inp[i][j]>2)*eta+int(inp[i][j]<3)*eps for j,q in enumerate(row)] for i,row in enumerate(state)] + + zplus = [[q+sum([state[(i+1)%dim][j]*(part[i][j]==part[(i+1)%dim][j]), + state[(i-1)%dim][j]*(part[i][j]==part[(i-1)%dim][j]), + state[i][(j+1)%dim]*(part[i][j]==part[i][(j+1)%dim]), + state[i][(j-1)%dim]*(part[i][j]==part[i][(j-1)%dim]) + ]) for j,q in enumerate(row)] for i,row in enumerate(state)] + zminus = [[(1-q)+sum([(1-state[(i+1)%dim][j])*(part[i][j]==part[(i+1)%dim][j]), + (1-state[(i-1)%dim][j])*(part[i][j]==part[(i-1)%dim][j]), + (1-state[i][(j+1)%dim])*(part[i][j]==part[i][(j+1)%dim]), + (1-state[i][(j-1)%dim])*(part[i][j]==part[i][(j-1)%dim]) + ]) for j,q in enumerate(row)] for i,row in enumerate(state)] + + kplus = [[kcomb(zplus[i][j],zminus[i][j]) for j,q in enumerate(row)] for i,row in enumerate(state)] + + pi = [[eps*(1-kplus[i][j]) + eta*kplus[i][j] for j,q in enumerate(row)] for i,row in enumerate(state)] + + crossent = [[-beps[i][j]*m.log2(pi[i][j])-(1-beps[i][j])*m.log2(1-pi[i][j]) for j,q in enumerate(row)] for i,row in enumerate(state)] + + return np.sum(crossent) - length*H2(eps) + + +def MIP(dim,statestr,eps): + lophi=np.inf + mip = [] + # statestr=statestr.translate(str.maketrans('','','.-=')) + for parti in range(1,2**(dim**2-1)): + partstr = bin(parti)[2:].zfill(dim**2) + curphi = phi(dim,statestr,partstr,eps) + + if curphi= 380 and wavelength <= 440: + attenuation = 0.3 + 0.7 * (wavelength - 380) / (440 - 380) + R = ((-(wavelength - 440) / (440 - 380)) * attenuation) ** gamma + G = 0.0 + B = (1.0 * attenuation) ** gamma + elif wavelength >= 440 and wavelength <= 490: + R = 0.0 + G = ((wavelength - 440) / (490 - 440)) ** gamma + B = 1.0 + elif wavelength >= 490 and wavelength <= 510: + R = 0.0 + G = 1.0 + B = (-(wavelength - 510) / (510 - 490)) ** gamma + elif wavelength >= 510 and wavelength <= 580: + R = ((wavelength - 510) / (580 - 510)) ** gamma + G = 1.0 + B = 0.0 + elif wavelength >= 580 and wavelength <= 645: + R = 1.0 + G = (-(wavelength - 645) / (645 - 580)) ** gamma + B = 0.0 + elif wavelength >= 645 and wavelength <= 750: + attenuation = 0.3 + 0.7 * (750 - wavelength) / (750 - 645) + R = (1.0 * attenuation) ** gamma + G = 0.0 + B = 0.0 + else: + R = 0.0 + G = 0.0 + B = 0.0 + R *= 255 + G *= 255 + B *= 255 + return (int(R), int(G), int(B)) + +def plot_execute(): + if sys.flags.interactive != 1 or not hasattr(qt.QtCore, 'PYQT_VERSION'): + qt.QtGui.QApplication.exec_() + +def qtplot(titletext, spaces, vals, names=None, + x_tag=f'noise level {chr(949)}', y_tag=None, x_log=False, y_log=False, x_range=None, y_range=None, + colors=None, + export=False, path=None, filename=None, lw=3, close=False): + linewidth = lw + + if not close: + app = qt.mkQApp() + + ph = qt.plot() + ph.showGrid(x=True,y=True) + if x_range is not None: + ph.setXRange(*x_range) + else: + ph.setXRange(min([min(sp) for sp in spaces]), max([max(sp) for sp in spaces])) + if y_range is not None: + ph.setYRange(*y_range) + else: + ph.setYRange(min([min(v) for v in vals]), max([max(v) for v in vals])) + ph.setLogMode(x=x_log, y=y_log) + + ph.setTitle(title=titletext, offset=(1000,500))#.format(dim)) + ph.setLabel('left', y_tag) + ph.setLabel('bottom', x_tag) + #ph.setXRange(0, 0.15) + ph.setFlag(ph.ItemIgnoresTransformations, False) # the relevant line + + if names is not None: + ph.addLegend(offset=(400, 30),size=20) + + + #s = ph.plot(np.linspace(0.01,0.32,32), eps_max_freq0, title=sp_Title, pen='w') + #s = ph.plot(np.linspace(0.01,0.32,32), eps_max_freq1, title=sp_Title, pen='w') + + + if colors=='rgb': + colors=[__get_color(fac/(len(vals)-1)) for fac in range(len(vals))] + elif colors is None: + colors=['r', 'g', 'b', 'y', 'c', 'm', 'w', + (100,100,0), (0,100,100), (100,0,100), + (100,200,0), (0,100,200), (200,0,100), + (200,100,0), (0,200,100), (100,0,200)] + else: + colors = colors[::-1] + + + for i in range(len(vals)-1,-1,-1): + try: + s = ph.plot(spaces[i], + vals[i], + name=names[i] if names is not None else None, + pen=qt.mkPen(colors[i], + width=linewidth)) + except: + print('Failed plotting '+names[i]) + raise + + #title_item = qt.TextItem(text=titletext, anchor=(0.5, 7), color='grey') + #ph.addItem(title_item) + #title_item.setPos(ph.getViewBox().viewRect().center()) + #font = QFont() + #font.setPointSize(14) # Adjust the font size as needed + #title_item.setFont(font) + + if export: + new_folder(path) + + exporter = qt.exporters.ImageExporter(ph.plotItem) + + # set export parameters if needed + exporter.parameters()['width'] = 1200 # (note this also affects height parameter) + + # save to file + exporter.export(path+filename) + + print(f'Saving to {path+filename}') + + def handleAppExit(): + # Add any cleanup tasks here + print("closing") + + if close: + ph.close() + else: + app.aboutToQuit.connect(handleAppExit) + app.exec_() + \ No newline at end of file diff --git a/evaluation/sign.py b/evaluation/sign.py new file mode 100644 index 0000000..349c304 --- /dev/null +++ b/evaluation/sign.py @@ -0,0 +1,374 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Sat Sep 9 22:23:30 2023 + +@author: astral +""" + + + +def sampling(dim,eps,samplemethod,extremes,dtmax,noeff=1): + print(f'Starting causal simulation with dim={dim}, eps={eps:.3f}, {samplemethod}, {extremes} extremes for dt={dtmax} and noeff={noeff:03d}') + + import os + import json + import math as m + import numpy as np + from numpy.linalg import norm + from datetime import datetime + from random import sample as choose + + from plot import qtplot + + from neuropercolation import Simulate4Layers + + eps_space = list(np.linspace(0.01,0.2,20)) + + def resultant(sample): + phase_x = [m.cos(ind) for ind in sample] + phase_y = [m.sin(ind) for ind in sample] + + return (np.average(phase_x), np.average(phase_y)) + + def new_folder(path): + if not os.path.exists(path): + os.makedirs(path) + + phase = np.vectorize(lambda x,y: (m.atan2(y,x)+m.pi)%(2*m.pi)-m.pi) + diff = np.vectorize(lambda x,y: (y-x+m.pi)%(2*m.pi)-m.pi) + H2 = lambda x: -x*m.log2(x)-(1-x)*m.log2(1-x) + + #print(f'Started at {datetime.now()} with eps={eps:.3f}') + + eps = round(eps,3) + path='/cloud/Public/_data/neuropercolation/4lay/cons=27-knight_steps=1000100/dim=09/batch=0/' + + try: + with open(path+f"eps={round(eps,3):.3f}_phase_diff.txt", 'r', encoding='utf-8') as f: + phase_diff = json.load(f) + except: + with open(path+f"eps={round(eps,3):.3f}_activation.txt", 'r', encoding='utf-8') as f: + activation = json.load(f)[100:] + + osc = list(zip(*activation)) + phase_abs = np.array([[np.arctan2(*act[::-1]) for act in osc[i]] for i in range(2)]) + phase_diff = diff(phase_abs[0],phase_abs[1]) + phase_diff = [round(pha,6) for pha in phase_diff] + + with open(path+f"eps={round(eps,3):.3f}_phase_diff.txt", 'w', encoding='utf-8') as f: + json.dump(list(phase_diff), f, indent=1) + + + all_res = norm(resultant(phase_diff)) + av_diff = np.arccos(all_res) + + try: + with open(path+f"eps={round(eps,3):.3f}_ei.txt", 'r', encoding='utf-8') as f: + ei = json.load(f) + except: + with open(path+f"eps={round(eps,3):.3f}_channels.txt", 'r', encoding='utf-8') as f: + channels = json.load(f)[100:] + + ei = [round(np.sum(cha)*(1-H2(eps)),6) for cha in channels] + + with open(path+f"eps={round(eps,3):.3f}_ei.txt", 'w', encoding='utf-8') as f: + json.dump(ei, f, indent=1) + + ei_ind = [i for i,val in enumerate(ei) if val>0] + + print(f'{len(ei_ind)} states with positive EI') + + samples = choose(ei_ind, extremes) + + with open(path+f"eps={round(eps,3):.3f}_states.txt", 'r', encoding='utf-8') as f: + states = json.load(f)[100:] + with open(path+f"eps={round(eps,3):.3f}_coupling.txt", 'r', encoding='utf-8') as f: + coupling = json.load(f) + coupling = [tuple(edge) for edge in coupling] + + phase_pairs = [[] for dt in range(dtmax+1)] + ei_pairs = [[] for dt in range(dtmax+1)] + for num,i in enumerate(samples): + causal_init = states[i].translate(str.maketrans('', '', '.-=')) + + sim = Simulate4Layers( dim, + eps, + coupling=coupling, + init=causal_init, + noeffect=noeff-1, + steps=dtmax, + draw=None, + ) + + activation = sim._activations() + channel = sim._channels() + + osc = list(zip(*activation)) + phase_abs = np.array([[np.arctan2(*act[::-1]) for act in osc[i]] for i in range(2)]) + phasediff_c = diff(phase_abs[0],phase_abs[1]) + ei_c = [round(np.sum(cha)*(1-H2(eps)),6) for cha in channel] + max_ei_c = max([np.sum(cha) for cha in channel]) + + sim = Simulate4Layers( dim, + eps, + coupling=coupling, + init=causal_init, + noeffect=-1, + steps=dtmax, + draw=None, + ) + + activation = sim._activations() + channel = sim._channels() + + osc = list(zip(*activation)) + phase_abs = np.array([[np.arctan2(*act[::-1]) for act in osc[i]] for i in range(2)]) + phasediff_i = diff(phase_abs[0],phase_abs[1]) + ei_i = [round(np.sum(cha)*(1-H2(eps)),6) for cha in channel] + max_ei_i = max([np.sum(cha) for cha in channel]) + + for dt in range(1,dtmax+1): + phase_pairs[dt].append((i, phasediff_i[dt], phasediff_c[dt])) + ei_pairs[dt].append((i, ei_i[dt], ei_c[dt])) + + if num%1000==999: + print(f'Done {num:0{len(str(extremes))}d}') + + + for dt in range(1,dtmax+1): + savepath = path + samplemethod + f'_samples={extremes}/dt={dt}/noeff={noeff:03d}/' + new_folder(savepath) + with open(savepath+f"eps={round(eps,3):.3f}_phase_pairs.txt", 'w', encoding='utf-8') as f: + json.dump(phase_pairs[dt], f, indent=1) + with open(savepath+f"eps={round(eps,3):.3f}_ei_pairs.txt", 'w', encoding='utf-8') as f: + json.dump(ei_pairs[dt], f, indent=1) + +def plotting(dim,eps,samplemethod,extremes,dt): + import os + import json + import math as m + import numpy as np + from numpy.linalg import norm + from datetime import datetime + from random import sample as choose + + from plot import qtplot + + from neuropercolation import Simulate4Layers + + eps_space = list(np.linspace(0.01,0.2,20)) + + def resultant(sample): + phase_x = [m.cos(ind) for ind in sample] + phase_y = [m.sin(ind) for ind in sample] + + return (np.average(phase_x), np.average(phase_y)) + + def new_folder(path): + if not os.path.exists(path): + os.makedirs(path) + + phase = np.vectorize(lambda x,y: (m.atan2(y,x)+m.pi)%(2*m.pi)-m.pi) + diff = np.vectorize(lambda x,y: (y-x+m.pi)%(2*m.pi)-m.pi) + H2 = lambda x: -x*m.log2(x)-(1-x)*m.log2(1-x) + + + print(f'Started at {datetime.now()} with eps={eps:.3f}') + + eps = round(eps,3) + path='/cloud/Public/_data/neuropercolation/4lay/cons=27-knight_steps=1000100/dim=09/batch=0/' + + savepath = path + samplemethod + f'_samples={extremes}/dt={dt}/' + + try: + with open(savepath+f"eps={round(eps,3):.3f}_phase_pairs.txt", 'r', encoding='utf-8') as f: + phase_pairs = json.load(f) + with open(savepath+f"eps={round(eps,3):.3f}_ei_pairs.txt", 'r', encoding='utf-8') as f: + ei_pairs = json.load(f) + except: + sampling(dim,eps,samplemethod,extremes,dt) + with open(savepath+f"eps={round(eps,3):.3f}_phase_pairs.txt", 'r', encoding='utf-8') as f: + phase_pairs = json.load(f) + with open(savepath+f"eps={round(eps,3):.3f}_ei_pairs.txt", 'r', encoding='utf-8') as f: + ei_pairs = json.load(f) + + t, phases_i, phases_c = zip(*phase_pairs) + t, ei_i, ei_c = zip(*ei_pairs) + + extremes = len(t) + + phases_cdiff = [abs(phases_i[i])-abs(phases_c[i]) for i in range(len(t))] + + phase_space = np.linspace(-m.pi,m.pi,101) + absph_space = np.linspace(0,m.pi,50+1) + cdiff_space = np.linspace(min(phases_cdiff),max(phases_cdiff),51) + + phase_dist_i = [len([ph for ph in phases_i if (phase_space[j]+phase_space[j-1])/2<=ph<(phase_space[j]+phase_space[j+1])/2])/len(t) + for j in range(100)] + phase_dist_c = [len([ph for ph in phases_c if (phase_space[j]+phase_space[j-1])/2<=ph<(phase_space[j]+phase_space[j+1])/2])/len(t) + for j in range(100)] + + absph_dist_i = [len([ph for ph in phases_i if low<=abs(ph) max_bat: + # max_bat = int(bat)+1 + + for batch in range(max_bat, max_bat+batches): + savepath = dimpath + f'batch={batch}' + new_folder(savepath) + for eps in [0.3]:#list(eps_space[61::2]): eps = round(eps,3) - cons = [(n,(n+m)%dim) for n in range(dim) for m in [0,int(dim/2)]] initstate = [[0,0],[0,0]] sim = Simulate4Layers(dim, eps, coupling=cons, init=initstate, steps=stp, + noeffect=-1, + #fps=20, draw=None, - res=2, - save='simple', - path=f'/cloud/Public/_data/neuropercolation/4lay/cons={len(cons)}-2diag_steps={stp}/dim={dim:02}/batch={batch}/', + res=6, + save='all', + path=path, ) print(f'Done eps={eps:.3f} with dim={dim} at {datetime.now()}') \ No newline at end of file diff --git a/neuropercolation/automaton.py b/neuropercolation/automaton.py index c684d7f..b2183da 100644 --- a/neuropercolation/automaton.py +++ b/neuropercolation/automaton.py @@ -92,7 +92,7 @@ class Neurolattice(CellularAutomatonCreator, abc.ABC): self.epsilon = eps def init_cell_state(self, coord): # pylint: disable=no-self-use - return DEAD + return [self.init] def get_cells(self): return self._current_state @@ -195,7 +195,8 @@ class Neuropercolation(CellularAutomatonCreator, abc.ABC): for coord, old, new in zip(this_state.keys(), this_state.values(), next_state.values()): coord_c = tuple([*coord[:2],int(1-coord[2])]) old_c = this_state[coord_c] - new_state = evolution_rule(old.state.copy(), old_c.state.copy(), [n.state[0] for n in old.neighbors], coord[2], coord_c[2]) #inverse the inhibitory layer's action + new_state = evolution_rule(old.state.copy(), old_c.state.copy(), [n.state[0] for n in old.neighbors], coord[2], coord_c[2]) + #inverse the inhibitory layer's action evolve_cell(old, new, new_state) @@ -233,16 +234,25 @@ class Neuropercolation(CellularAutomatonCreator, abc.ABC): class NeuropercolationCoupled(CellularAutomatonCreator, abc.ABC): - def __init__(self, dim, eps, coupling=[], init=[[0,0],[0,0]], *args, **kwargs): + def __init__(self, dim, eps, coupling=[], init=[[0,0],[0,0]], noeffect=-1, *args, **kwargs): super().__init__(dimension=[dim, dim, 2, 2], init=init, neighborhood=VonNeumannNeighborhood(EdgeRule.FIRST_AND_LAST_CELL_OF_DIMENSION_ARE_NEIGHBORS)) self._evolution_step = 0 self.epsilon = eps self.coupling = coupling + self.noeffect = noeffect def init_cell_state(self, coord): # pylint: disable=no-self-use - return [self.init[coord[3]][coord[2]]] + if type(self.init) is list: + return [self.init[coord[3]][coord[2]]] + elif type(self.init) is str: + self.init = self.init.translate(str.maketrans('', '', '.-=')) + dim = self._dimension[0] + assert len(self.init)==4*dim**2 + coordval = [1, dim, dim**2, 2*dim**2] + digit = sum([coord[i]*coordval[i] for i in range(4)]) + return [int(self.init[digit])] def get_cells(self): return self._current_state @@ -277,11 +287,14 @@ class NeuropercolationCoupled(CellularAutomatonCreator, abc.ABC): if coord[:2] in self.coupling: coord_c = tuple([*coord[:2],coord[2],int(1-coord[3])]) old_c = this_state[coord_c] - new_state = evolution_rule(old.state.copy(), old_c.state.copy(), [n.state[0] for n in old.neighbors], 0) + new_state = evolution_rule(old.state.copy(), old_c.state.copy(), [n.state[0] for n in old.neighbors], 0) \ + if self.noeffect=self.epsilon) - if alive_neighbours > 2: + if alive_neighbours > 2.5: if CASE: new_cell_state = ALIVE else: new_cell_state = DEAD + elif alive_neighbours < 2.5: + if CASE: + new_cell_state = DEAD + else: + new_cell_state = ALIVE else: - if CASE: - new_cell_state = DEAD - else: - new_cell_state = ALIVE + new_cell_state = [int(random.random()>=0.5)] + return new_cell_state def count_channels(self): diff --git a/neuropercolation/display.py b/neuropercolation/display.py index c1965ab..d222de6 100644 --- a/neuropercolation/display.py +++ b/neuropercolation/display.py @@ -148,7 +148,7 @@ class Simulate1Layer: self._append_activation() def _append_all(self): - automaton_state = [[0 for n in range(self.__dimension)] for m in range(self.__dimension)] + automaton_state = [[0 for n in range(self.__dimension)] for l in range(self.__dimension)] activation = 0 for coord, cell in self._cellular_automaton._current_state.items(): x,y = coord @@ -321,7 +321,7 @@ class Simulate2Layers: self.__channel_list.append(self._cellular_automaton.count_channels()) def _append_all(self): - automaton_state = [[[0 for n in range(self.__dimension)] for m in range(self.__dimension)] for l in range(2)] + automaton_state = [[[0 for n in range(self.__dimension)] for l in range(self.__dimension)] for k in range(2)] activation = [0]*2 for coord, cell in self._cellular_automaton._current_state.items(): x,y,l = coord @@ -424,6 +424,7 @@ class Simulate4Layers: eps, coupling=[], init=[[0,0],[0,0]], + noeffect=-1, steps=100, draw='pygame', res=4, @@ -442,11 +443,12 @@ class Simulate4Layers: :param state_to_color_cb: A callback to define the draw color of CA states (default: red for states != 0) """ super().__init__(*args, **kwargs) - self._cellular_automaton = NeuropercolationCoupled(dim,eps,coupling,init) + self._cellular_automaton = NeuropercolationCoupled(dim,eps,coupling,init,noeffect) self.__active = True self.__cell_size = [res,res] self.__dimension = dim self.__epsilon = eps + self.__coupling = coupling self.__gridside = dim*res self.__samegap = 10 self.__othergap = 20 @@ -499,18 +501,22 @@ class Simulate4Layers: except: print('Failed to quit pygame') + def _activations(self): + return self.__activation_list - + def _channels(self): + return self.__channel_list + def _track(self): if self.__save == 'all': self._append_all() - elif self.__save == 'simple': + else: self._append_activation() self.__channel_list.append(self._cellular_automaton.count_channels()) def _append_all(self): - automaton_state = [[[[0 for n in range(self.__dimension)] for m in range(self.__dimension)] for l in range(2)] for k in range(2)] + automaton_state = [[[[0 for n in range(self.__dimension)] for l in range(self.__dimension)] for k in range(2)] for j in range(2)] activation = [[0,0],[0,0]] for coord, cell in self._cellular_automaton._current_state.items(): x,y,l,o = coord @@ -545,7 +551,9 @@ class Simulate4Layers: def _save_all(self): if not os.path.exists(self.__path): os.makedirs(self.__path) - + + with open(self.__path+f"eps={round(self.__epsilon,3):.3f}_coupling.txt", 'w', encoding='utf-8') as f: + json.dump(self.__coupling, f, indent=1) with open(self.__path+f"eps={round(self.__epsilon,3):.3f}_states.txt", 'w', encoding='utf-8') as f: json.dump(self.__state_list, f, indent=1) with open(self.__path+f"eps={round(self.__epsilon,3):.3f}_activation.txt", 'w', encoding='utf-8') as f: @@ -557,6 +565,8 @@ class Simulate4Layers: if not os.path.exists(self.__path): os.makedirs(self.__path) + with open(self.__path+f"eps={round(self.__epsilon,3):.3f}_coupling.txt", 'w', encoding='utf-8') as f: + json.dump(self.__coupling, f, indent=1) with open(self.__path+f"eps={round(self.__epsilon,3):.3f}_activation.txt", 'w', encoding='utf-8') as f: json.dump(self.__activation_list, f, indent=1) with open(self.__path+f"eps={round(self.__epsilon,3):.3f}_channels.txt", 'w', encoding='utf-8') as f: