Trying to understand Resultant evolution

evaluation/4Layer Bootstrap.py

@@ -23,15 +23,16 @@ def resultant(sample):
     
     return (np.average(phase_x), np.average(phase_y))
 
-def H2(x):
+phase = lambda x,y: (m.atan2(y,x)+m.pi)%(2*m.pi)-m.pi
-    return -x*m.log2(x)-(1-x)*m.log2(1-x)
+
+H2 = lambda x: -x*m.log2(x)-(1-x)*m.log2(1-x)
 
 extremes = None
 maxdt = 300
 
-for dim in [9]:
+for dim in []:
     for eps in eps_space:
-        path=f'/cloud/Public/_data/neuropercolation/4lay/cons=27-3diag_steps=1000100/dim=09/batch=0/'
+        path=f'/cloud/Public/_data/neuropercolation/4lay/cons=7-knight_steps=1000100/dim=07/batch=0/'
         
         try:
             with open(path+f"eps={round(eps,3):.3f}_phase_diff.txt", 'r', encoding='utf-8') as f:
@@ -64,74 +65,106 @@ for dim in [9]:
                 json.dump(ei, f, indent=1)
         
         pha_center = av_diff
-        pha_dev = m.pi/8
+        pha_dev = m.pi/32
-        orth_ind = sorted([i for i,val in enumerate(ei[:-maxdt]) if (pha_center-pha_dev)<=abs(phase_diff[i])<=(pha_center+pha_dev)], key = lambda i: ei[i])
+        
+        from_sync = lambda i: True if abs(phase_diff[i])<0.08*m.pi else False if 0.42*m.pi<abs(phase_diff[i])<0.58*m.pi else from_sync(i-1)
+        
+        dev_ind = sorted([i for i,val in enumerate(ei[:-maxdt]) if (pha_center-pha_dev)<=abs(phase_diff[i])<=(pha_center+pha_dev)], key = lambda i: ei[i])
+        dev_in_ind = [i for i in dev_ind if not from_sync(i)]
+        dev_out_ind = [i for i in dev_ind if from_sync(i)]
+        
+        devm_in_ind = [i for i in dev_in_ind if phase_diff[i]<0]
+        devp_in_ind = [i for i in dev_in_ind if phase_diff[i]>0]
+        devm_out_ind = [i for i in dev_out_ind if phase_diff[i]<0]
+        devp_out_ind = [i for i in dev_out_ind if phase_diff[i]>0]
         
         if extremes is None:
-            extremes = len(orth_ind)//4000*1000
+            extremes = min(len(dev_in_ind),len(dev_out_ind))//1000*100
 
-        print(len(orth_ind))
+        print(len(dev_out_ind), len(dev_in_ind))
         #print(all_res, av_diff)
         
         
-        bot_ind = orth_ind[         :extremes]
+        bot_in_ind = dev_in_ind[         :extremes]
-        top_ind = orth_ind[-extremes:        ]
+        top_in_ind = dev_in_ind[-extremes:        ]
+        bot_out_ind = dev_out_ind[         :extremes]
+        top_out_ind = dev_out_ind[-extremes:        ]
         
-        # bot_ei_pool = [ei for ei in enumerate(ei[:-maxdt]) if ei[1] <= bot_ei[1]]
+        bot_ind = bot_in_ind + bot_out_ind
-        # top_ei_pool = [ei for ei in enumerate(ei[:-maxdt]) if ei[1] >= top_ei[1]]
+        top_ind = top_in_ind + top_out_ind
-        
-        # bot_eis = choose(bot_ei_pool, extremes)
-        # top_eis = choose(top_ei_pool, extremes)
-        
-        # bot_ind = [enum[0] for enum in bot_ei]
-        # top_ind = [enum[0] for enum in top_ei]
         
         bot_res = []
         top_res = []
-        orth_res = []
+        dev_res = []
         tot_res = []
         
+        bot_ph = []
+        top_ph = []
+        dev_ph = []
+        tot_ph = []
+        
         bot_ei = []
         top_ei = []
-        orth_ei = []
+        dev_ei = []
         tot_ei = []
         for dt in range(maxdt):
-            bot_res.append( norm(resultant([phase_diff[i+dt] for i in bot_ind])) )
+            bot_pha = [phase_diff[i+dt] for i in bot_ind]
-            top_res.append( norm( resultant([phase_diff[i+dt] for i in top_ind])) )
+            top_pha = [phase_diff[i+dt] for i in top_ind]
-            orth_res.append( norm( resultant([phase_diff[i+dt] for i in orth_ind])) )
+            dev_pha = [phase_diff[i+dt] for i in dev_ind]
-            tot_res.append( norm( resultant(phase_diff[dt:dt-maxdt])) )
+            tot_pha = 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)) )
+            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)) )
+            tot_ph.append( phase(resultant(tot_pha)) )
             
             bot_ei.append( np.average([ei[i+dt] for i in bot_ind]) )
             top_ei.append( np.average([ei[i+dt] for i in top_ind]) )
-            orth_ei.append( np.average([ei[i+dt] for i in orth_ind]) )
+            dev_ei.append( np.average([ei[i+dt] for i in dev_ind]) )
             tot_ei.append( np.average(ei[dt:dt-maxdt]) )
             
             if dt%100==99:
                 print(f'Done dt={dt}')
             
-            
         qtplot(f'Diachronic resultant for dim={dim} with 4 layers',
                 [np.array(range(maxdt))]*4,
-                [bot_res, top_res, orth_res, tot_res],
+                [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'],
+                  'Resultant ev of phase filtered ei', 'Average Resultant'],
                 x_tag = 'dt',
                 y_tag = 'concentration',
                 export=True,
-                path=path,
+                path=path+'plots/',
-                filename=f'Diachronic Resultant for eps={round(eps,3):.3f} dim={dim} extremes={extremes} roll{pha_dev:.3f}.png',
+                filename=f'Diachronic Resultant balanced 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_ph, top_ph, dev_ph, tot_ph],
+                ['Resultant ev of bottom {extremes} ei', 'Resultant ev of top {extremes} ei',
+                  'Resultant ev of phase filtered ei', 'Average Resultant'],
+                x_tag = 'dt',
+                y_tag = 'concentration',
+                export=True,
+                path=path+'plots/',
+                filename=f'Diachronic Resultant balanced 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, orth_ei, tot_ei],
+                [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,
+                path=path+'plots/',
-                filename=f'Diachronic EI for eps={round(eps,3):.3f} dim={dim} extremes={extremes} roll{pha_dev:.3f}.png',
+                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()}')

evaluation/4Layer Percolation small eps.py

@@ -0,0 +1,153 @@
+#!/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))
+
+def H2(x):
+    return -x*m.log2(x)-(1-x)*m.log2(1-x)
+
+extremes = None
+maxdt = 300
+
+for dim in [9]:
+    for eps in eps_space[:10]:
+        path=f'/cloud/Public/_data/neuropercolation/4lay/cons=7-knight_steps=1000100/dim=07/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 = 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(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
+        orth_ind = sorted([i for i,val in enumerate(ei[:-maxdt]) if (pha_center-pha_dev)<=abs(phase_diff[i])<=(pha_center+pha_dev)], key = lambda i: ei[i])
+        
+        thres = 0.04*m.pi
+        last_sync = lambda i: i if abs(phase_diff[i])<0.08*m.pi else last_sync(i-1) 
+        last_orth = lambda i: i if 0.42*m.pi<abs(phase_diff[i])<0.58*m.pi else last_orth(i-1) 
+        
+        trav_out = [(i,pha) for i,pha in enumerate(phase_diff[:-maxdt]) if abs(pha)<thres and abs(phase_diff[i+1])>thres]
+        diff00_ind = [last_sync(i) for i,pha in trav_out]
+        diff00m_ind = [last_sync(i) for i,pha in trav_out if pha<0]
+        diff00p_ind = [last_sync(i) for i,pha in trav_out if pha>0]
+        
+        trav_in = [(i,pha) for i,pha in enumerate(phase_diff[:-maxdt]) if m.pi/2-thres<abs(pha)<m.pi/2+thres and abs(phase_diff[i+1])<m.pi/2-thres]
+        diff90_ind = [last_orth(i) for i,pha in trav_in]
+        diff90m_ind = [last_orth(i) for i,pha in trav_in if pha<0]
+        diff90p_ind = [last_orth(i) for i,pha in trav_in if pha>0]
+        
+        print(f'{len(diff00_ind)}/{len(trav_out)} synced states and {len(diff90_ind)}/{len(trav_in)} orthogonal states')
+        #print(all_res, av_diff)
+        
+        in_ei = []
+        out_ei = []
+        
+        inm_res = []
+        inp_res = []
+        outm_res = []
+        outp_res = []
+        for dt in range(maxdt):
+            inm_res.append( resultant([phase_diff[i+dt] for i in diff90m_ind]) )
+            outm_res.append( resultant([phase_diff[i+dt] for i in diff00m_ind]) )
+            
+            inp_res.append( resultant([phase_diff[i+dt] for i in diff90p_ind]) )
+            outp_res.append( resultant([phase_diff[i+dt] for i in diff00p_ind]) )
+            
+            in_ei.append( np.average([ei[i+dt] for i in diff90_ind]) )
+            out_ei.append( np.average([ei[i+dt] for i in diff00_ind]) )
+            
+            if dt%100==99:
+                print(f'Done dt={dt}')
+            
+        im_res = list(zip(*inm_res))
+        om_res = list(zip(*outm_res))
+        
+        ip_res = list(zip(*inp_res))
+        op_res = list(zip(*outp_res))
+        
+        space = np.linspace(0,2*m.pi,101)
+        
+        av_ei = np.average(ei)
+        
+        qtplot(f'Transitions resultants for eps={round(eps,3):.3f} with dim={dim} with 4 layers',
+                [im_res[0], om_res[0], ip_res[0], op_res[0], np.cos(space)],
+                [im_res[1], om_res[1], ip_res[1], op_res[1], np.sin(space)],
+                ['Resultant ev inward -', 'Resultant ev outward -',
+                 'Resultant ev inward +', 'Resultant ev outward +',
+                 'Unit Circle'],
+                x_tag = 'x',
+                y_tag = 'y',
+                lw=1,
+                export=True,
+                path=path,
+                filename=f'Transitions resultant pm for eps={round(eps,3):.3f} dim={dim}.png',
+                close=False)
+        
+        qtplot(f'Transitions ei for eps={round(eps,3):.3f} with dim={dim} with 4 layers',
+                [np.array(range(maxdt))]*3,
+                [in_ei, out_ei, [av_ei]*maxdt],
+                ['EI ev inward', 'EI ev outward', 'Average EI'],
+                x_tag = 'dt',
+                y_tag = 'average ei',
+                export=True,
+                path=path,
+                filename=f'Transitions EI for eps={round(eps,3):.3f} dim={dim} .png',
+                close=False)
+        
+        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)
+        
+        

evaluation/4Layer Resultant Evolution.py

@@ -0,0 +1,266 @@
+#!/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))
+
+def H2(x):
+    return -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/cons=7-knight_steps=1000100/dim=07/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 = 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(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
+        
+        last_sync = lambda i: i if abs(phase_diff[i])<0.08*m.pi else last_sync(i-1) 
+        last_orth = lambda i: i if 0.42*m.pi<abs(phase_diff[i])<0.58*m.pi else last_orth(i-1) 
+        
+        dev_ind = sorted([i for i,val in enumerate(ei[:-maxdt]) if (pha_center-pha_dev)<=abs(phase_diff[i])<=(pha_center+pha_dev)], key = lambda i: ei[i])
+        dev_in_ind = [i for i in dev_ind if last_sync(i)<last_orth(i)]
+        dev_out_ind = [i for i in dev_ind if last_sync(i)>last_orth(i)]
+        
+        devm_in_ind = [i for i in dev_in_ind if phase_diff[i]<0]
+        devp_in_ind = [i for i in dev_in_ind if phase_diff[i]>0]
+        devm_out_ind = [i for i in dev_out_ind if phase_diff[i]<0]
+        devp_out_ind = [i for i in dev_out_ind if phase_diff[i]>0]
+        
+        if extremes is None:
+            extremes = len(dev_in_ind)//4000*1000
+
+        print(len(dev_out_ind), len(dev_in_ind))
+        #print(all_res, av_diff)
+        
+        bot_in_ind = dev_in_ind[         :extremes]
+        top_in_ind = dev_in_ind[-extremes:        ]
+        bot_out_ind = dev_out_ind[         :extremes]
+        top_out_ind = dev_out_ind[-extremes:        ]
+        
+        bot_ind = bot_in_ind + bot_out_ind
+        top_ind = top_in_ind + top_out_ind
+        
+        botm_in_ind = [i for i in bot_in_ind if phase_diff[i]<0]
+        botp_in_ind = [i for i in bot_in_ind if phase_diff[i]>0]
+        topm_in_ind = [i for i in top_in_ind if phase_diff[i]<0]
+        topp_in_ind = [i for i in top_in_ind if phase_diff[i]>0]
+        
+        botm_out_ind = [i for i in bot_out_ind if phase_diff[i]<0]
+        botp_out_ind = [i for i in bot_out_ind if phase_diff[i]>0]
+        topm_out_ind = [i for i in top_out_ind if phase_diff[i]<0]
+        topp_out_ind = [i for i in top_out_ind if phase_diff[i]>0]
+        
+        # bot_ei_pool = [ei for ei in enumerate(ei[:-maxdt]) if ei[1] <= bot_ei[1]]
+        # top_ei_pool = [ei for ei in enumerate(ei[:-maxdt]) if ei[1] >= top_ei[1]]
+        
+        # bot_eis = choose(bot_ei_pool, extremes)
+        # top_eis = choose(top_ei_pool, extremes)
+        
+        # bot_ind = [enum[0] for enum in bot_ei]
+        # top_ind = [enum[0] for enum in top_ei]
+        
+        botm_in_res = []
+        topm_in_res = []
+        devm_in_res = []
+        
+        botp_in_res = []
+        topp_in_res = []
+        devp_in_res = []
+        
+        botm_out_res = []
+        topm_out_res = []
+        devm_out_res = []
+        
+        botp_out_res = []
+        topp_out_res = []
+        devp_out_res = []
+        
+        bot_ei = []
+        top_ei = []
+        dev_ei = []
+        
+        tot_res = []
+        tot_ei = []
+        for dt in range(maxdt):
+            botm_in_res.append( resultant([phase_diff[i+dt] for i in botm_in_ind]) )
+            topm_in_res.append( resultant([phase_diff[i+dt] for i in topm_in_ind]) )
+            devm_in_res.append( resultant([phase_diff[i+dt] for i in devm_in_ind]) )
+            
+            botp_in_res.append( resultant([phase_diff[i+dt] for i in botp_in_ind]) )
+            topp_in_res.append( resultant([phase_diff[i+dt] for i in topp_in_ind]) )
+            devp_in_res.append( resultant([phase_diff[i+dt] for i in devp_in_ind]) )
+            
+            botm_out_res.append( resultant([phase_diff[i+dt] for i in botm_out_ind]) )
+            topm_out_res.append( resultant([phase_diff[i+dt] for i in topm_out_ind]) )
+            devm_out_res.append( resultant([phase_diff[i+dt] for i in devm_out_ind]) )
+            
+            botp_out_res.append( resultant([phase_diff[i+dt] for i in botp_out_ind]) )
+            topp_out_res.append( resultant([phase_diff[i+dt] for i in topp_out_ind]) )
+            devp_out_res.append( resultant([phase_diff[i+dt] for i in devp_out_ind]) )
+            
+            bot_ei.append( np.average([ei[i+dt] for i in bot_ind]) )
+            top_ei.append( np.average([ei[i+dt] for i in top_ind]) )
+            dev_ei.append( np.average([ei[i+dt] for i in dev_ind]) )
+            
+            tot_res.append( resultant(phase_diff[dt:dt-maxdt]) )
+            tot_ei.append( np.average(ei[dt:dt-maxdt]) )
+            
+            
+            if dt%100==99:
+                print(f'Done dt={dt}')
+            
+        bm_in_res = list(zip(*botm_in_res))
+        tm_in_res = list(zip(*topm_in_res))
+        dm_in_res = list(zip(*devm_in_res))
+        
+        bp_in_res = list(zip(*botp_in_res))
+        tp_in_res = list(zip(*topp_in_res))
+        dp_in_res = list(zip(*devp_in_res))
+        
+        bm_out_res = list(zip(*botm_out_res))
+        tm_out_res = list(zip(*topm_out_res))
+        dm_out_res = list(zip(*devm_out_res))
+        
+        bp_out_res = list(zip(*botp_out_res))
+        tp_out_res = list(zip(*topp_out_res))
+        dp_out_res = list(zip(*devp_out_res))
+        
+        t_res = list(zip(*tot_res))
+        
+        space = np.linspace(0,2*m.pi,101)
+        
+        qtplot(f'Diachronic resultant +/- for eps={round(eps,3):.3f} with dim={dim} with 4 layers',
+                [bm_in_res[0], tm_in_res[0], dm_in_res[0], bp_in_res[0], tp_in_res[0], dp_in_res[0],
+                 bm_out_res[0], tm_out_res[0], dm_out_res[0], bp_out_res[0], tp_out_res[0], dp_out_res[0],
+                 t_res[0], np.cos(space)],
+                [bm_in_res[1], tm_in_res[1], dm_in_res[1], bp_in_res[1], tp_in_res[1], dp_in_res[1],
+                 bm_out_res[1], tm_out_res[1], dm_out_res[1], bp_out_res[1], tp_out_res[1], dp_out_res[1],
+                 t_res[1], np.sin(space)],
+                ['Resultant ev in of bottom {extremes} ei -', 'Resultant ev in of top {extremes} ei -',
+                 'Resultant ev in of phase filtered ei -',
+                 'Resultant ev in of bottom {extremes} ei +', 'Resultant ev in of top {extremes} ei +',
+                 'Resultant ev in of phase filtered ei +',
+                 'Resultant ev out of bottom {extremes} ei -', 'Resultant ev out of top {extremes} ei -',
+                 'Resultant ev out of phase filtered ei -',
+                 'Resultant ev out of bottom {extremes} ei +', 'Resultant ev out of top {extremes} ei +',
+                 'Resultant ev out of phase filtered ei +',
+                 'Average Resultant', 'Unit Circle'],
+                x_tag = 'x',
+                y_tag = 'y',
+                lw=1,
+                export=True,
+                path=path,
+                filename=f'Diachronic Resultant filtered pm inout eps={round(eps,3):.3f} dim={dim} extremes={extremes} roll{pha_dev:.3f}.png',
+                close=True)
+        
+        qtplot(f'Diachronic resultant +/- for eps={round(eps,3):.3f} with dim={dim} with 4 layers',
+                [bm_in_res[0], tm_in_res[0], dm_in_res[0], bp_in_res[0], tp_in_res[0], dp_in_res[0],
+                 t_res[0], np.cos(space)],
+                [bm_in_res[1], tm_in_res[1], dm_in_res[1], bp_in_res[1], tp_in_res[1], dp_in_res[1],
+                 t_res[1], np.sin(space)],
+                ['Resultant ev in of bottom {extremes} ei -', 'Resultant ev in of top {extremes} ei -',
+                 'Resultant ev in of phase filtered ei -',
+                 'Resultant ev in of bottom {extremes} ei +', 'Resultant ev in of top {extremes} ei +',
+                 'Resultant ev in of phase filtered ei +',
+                 'Average Resultant', 'Unit Circle'],
+                x_tag = 'x',
+                y_tag = 'y',
+                lw=1,
+                export=True,
+                path=path,
+                filename=f'Evolution of Resultant filtered pm in eps={round(eps,3):.3f} dim={dim} extremes={extremes} roll{pha_dev:.3f}.png',
+                close=True)
+        
+        qtplot(f'Diachronic resultant +/- for eps={round(eps,3):.3f} with dim={dim} with 4 layers',
+                [bm_out_res[0], tm_out_res[0], dm_out_res[0], bp_out_res[0], tp_out_res[0], dp_out_res[0],
+                 t_res[0], np.cos(space)],
+                [bm_out_res[1], tm_out_res[1], dm_out_res[1], bp_out_res[1], tp_out_res[1], dp_out_res[1],
+                 t_res[1], np.sin(space)],
+                ['Resultant ev out of bottom {extremes} ei -', 'Resultant ev out of top {extremes} ei -',
+                 'Resultant ev out of phase filtered ei -',
+                 'Resultant ev out of bottom {extremes} ei +', 'Resultant ev out of top {extremes} ei +',
+                 'Resultant ev out of phase filtered ei +',
+                 'Average Resultant', 'Unit Circle'],
+                x_tag = 'x',
+                y_tag = 'y',
+                lw=1,
+                export=True,
+                path=path,
+                filename=f'Evolution of Resultant filtered pm out eps={round(eps,3):.3f} dim={dim} extremes={extremes} roll{pha_dev:.3f}.png',
+                close=True)
+        # qtplot(f'Diachronic ei +/- for eps={round(eps,3):.3f} with 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,
+        #         filename=f'Diachronic EI filtered 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)
+        
+        

evaluation/4Layer Resultant Phase Evolution.py

@@ -0,0 +1,181 @@
+#!/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 = lambda x,y: (m.atan2(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 [7]:
+    for eps in eps_space:
+        path=f'/cloud/Public/_data/neuropercolation/4lay/cons=7-knight_steps=1000100/dim=07/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 = 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(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.pi<abs(phase_diff[i])<0.58*m.pi else from_sync(i-1)
+        
+        dev_ind = sorted([i for i,val in enumerate(ei[:-maxdt]) if (pha_center-pha_dev)<=abs(phase_diff[i])<=(pha_center+pha_dev)], key = lambda i: ei[i])
+        dev_in_ind = [i for i in dev_ind if not from_sync(i)]
+        dev_out_ind = [i for i in dev_ind if from_sync(i)]
+        
+        devm_in_ind = [i for i in dev_in_ind if phase_diff[i]<0]
+        devp_in_ind = [i for i in dev_in_ind if phase_diff[i]>0]
+        devm_out_ind = [i for i in dev_out_ind if phase_diff[i]<0]
+        devp_out_ind = [i for i in dev_out_ind if phase_diff[i]>0]
+        
+        if extremes is None:
+            extremes = min(len(dev_in_ind),len(dev_out_ind))//1000*100
+
+        print(len(dev_out_ind), len(dev_in_ind))
+        #print(all_res, av_diff)
+        
+        
+        bot_in_ind = dev_in_ind[         :extremes]
+        top_in_ind = dev_in_ind[-extremes:        ]
+        bot_out_ind = dev_out_ind[         :extremes]
+        top_out_ind = dev_out_ind[-extremes:        ]
+        
+        bot_ind = bot_in_ind + bot_out_ind
+        top_ind = top_in_ind + top_out_ind
+        
+        bot_res = []
+        top_res = []
+        dev_res = []
+        tot_res = []
+        
+        bot_ph = []
+        top_ph = []
+        dev_ph = []
+        tot_ph = []
+        
+        bot_ei = []
+        top_ei = []
+        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]
+            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)) )
+            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)) )
+            tot_ph.append( phase(*resultant(tot_pha)) )
+            
+            bot_ei.append( np.average([ei[i+dt] for i in bot_ind]) )
+            top_ei.append( np.average([ei[i+dt] for i in top_ind]) )
+            dev_ei.append( np.average([ei[i+dt] for i in dev_ind]) )
+            tot_ei.append( np.average(ei[dt:dt-maxdt]) )
+            
+            if dt%100==99:
+                print(f'Done dt={dt}')
+            
+        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'],
+                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',
+                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'],
+                x_tag = 'dt',
+                y_tag = 'concentration',
+                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',
+                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)
+        
+        

evaluation/plot.py

@@ -78,11 +78,12 @@ def plot_execute():
 def qtplot(titletext, spaces, vals, names, x_tag=f'noise level {chr(949)}', y_tag=None, colors=None, export=False, path=None, filename=None, lw=4, close=False):
     linewidth = lw
     
-    #app = qt.mkQApp() 
+    if not close:
+        app = qt.mkQApp()
     
     ph = qt.plot()
     ph.showGrid(x=True,y=True)
-    ph.setXRange(np.min(spaces), np.max(spaces))
+    ph.setXRange(min([min(sp) for sp in spaces]), max([max(sp) for sp in spaces]))
     # ph.setYRange(0.0, 6)
     
     #ph.setTitle(title='Susceptibility density evolution for different automaton sizes', offset=(1000,500))#.format(dim))
@@ -100,7 +101,10 @@ def qtplot(titletext, spaces, vals, names, x_tag=f'noise level {chr(949)}', y_ta
     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)]
+        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)]
     
     
     for i in range(len(vals)):
@@ -134,7 +138,7 @@ def qtplot(titletext, spaces, vals, names, x_tag=f'noise level {chr(949)}', y_ta
     
     if close:
         ph.close()
-        
+    else:
-    # app.aboutToQuit.connect(handleAppExit)
+        app.aboutToQuit.connect(handleAppExit)
-    # app.exec_()
+        app.exec_()
     

examples/Simulate1Layer.py

@@ -15,7 +15,7 @@ eps_space = np.linspace(0.135,0.15,4)
 for dim in [49,100]:
     for eps in eps_space:
         eps = round(eps,3)
-        Simulate1Layer(dim,
+        Simulate2Layers(dim,
                        eps,
                        res=2,
                        path=f'/cloud/Public/_data/neuropercolation/1lay/steps=100000/dim={dim:02}/',

neuropercolation/automaton.py

@@ -152,8 +152,9 @@ class Neurolattice(CellularAutomatonCreator, abc.ABC):
 class Neuropercolation(CellularAutomatonCreator, abc.ABC):
     """ Cellular automaton with the evolution rules of conways game of life """
     
-    def __init__(self, dim, eps, *args, **kwargs):
+    def __init__(self, dim, eps, init, *args, **kwargs):
         super().__init__(dimension=[dim, dim, 2],
+                         init=init,
                          neighborhood=VonNeumannNeighborhood(EdgeRule.FIRST_AND_LAST_CELL_OF_DIMENSION_ARE_NEIGHBORS))
         self._evolution_step = 0
         self.epsilon = eps