diff --git a/evaluation/4Layer Bootstrap.py b/evaluation/4Layer Bootstrap.py
new file mode 100644
index 0000000..8562d81
--- /dev/null
+++ b/evaluation/4Layer Bootstrap.py	
@@ -0,0 +1,111 @@
+#!/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 = 50000
+maxdt = 500
+
+for dim in [7]:
+    for eps in eps_space:
+        path=f'/cloud/Public/_data/neuropercolation/4lay/steps=1000100/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)
+        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))
+        
+        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)
+        
+        bot_ei = sorted(enumerate(ei[:-maxdt]), key = lambda x: x[1])[ extremes]
+        top_ei = sorted(enumerate(ei[:-maxdt]), key = lambda x: x[1])[-extremes]
+        
+        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_eis]
+        top_ind = [enum[0] for enum in top_eis]
+        
+        bot_res = []
+        top_res = []
+        for dt in range(maxdt):
+            bot_pha = [phase_diff[i+dt] for i in bot_ind]
+            top_pha = [phase_diff[i+dt] for i in top_ind]
+            
+            bot_res.append( norm(resultant(bot_pha)) )
+            top_res.append( norm(resultant(top_pha)) )
+            if dt%100==99:
+                print(f'Done dt={dt}')
+            
+            
+        qtplot(f'Diachronic resultant for dim={dim} with 4 layers',
+                [np.array(range(maxdt))]*3,
+                [bot_res, top_res, [all_res]*maxdt],
+                ['Resultant ev of bottom 100 ei', 'Resultant ev of top 100 ei', 'Average Resultant'],
+                x_tag = 'dt',
+                y_tag = 'concentration',
+                export=True,
+                path=path,
+                filename=f'Diachronic Resultant for eps={round(eps,3)} dim={dim} extremes={extremes}.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
new file mode 100644
index 0000000..1f143f9
--- /dev/null
+++ b/evaluation/4Layer Res+EI.py	
@@ -0,0 +1,60 @@
+#!/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 datetime import datetime
+
+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)
+
+for dim in [7]:
+    diff_res = [0]
+    av_ei = [0]
+    for eps in eps_space:
+        path=f'/cloud/Public/_data/neuropercolation/4lay/steps=100000/dim={dim:02}/'
+        
+        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
+        
+        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))
+        
+        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)
+        
+        
\ No newline at end of file
diff --git a/evaluation/4Layer Survey.py b/evaluation/4Layer Survey.py
new file mode 100644
index 0000000..5b1c8a8
--- /dev/null
+++ b/evaluation/4Layer Survey.py	
@@ -0,0 +1,94 @@
+#!/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 = 50000
+maxdt = 500
+
+for dim in [7]:
+    for eps in eps_space:
+        path=f'/cloud/Public/_data/neuropercolation/4lay/steps=1000100/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)
+        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))
+        
+        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_divs = 8
+        pha_sort = [0]*pha_divs
+        
+        pha_sort = [[pha for pha in list(enumerate(phase_diff))[:-maxdt] if i*m.pi/pha_divs <= abs(pha[1]) <= (i+1)*m.pi/pha_divs] for i in range(pha_divs)]
+        dia_ei = [[np.average([ei[pha[0]+dt] for pha in pha_div]) for dt in range(maxdt)] for pha_div in pha_sort]
+            
+        qtplot(f'Diachronic EI for dim={dim} with 4 layers',
+                [list(range(maxdt))]*pha_divs,
+                dia_ei,
+                [f'EI evolution for initial phase in [{div}pi,{div+1}pi]' for div in range(pha_divs)],
+                x_tag = 'dt',
+                y_tag = 'average EI',
+                colors = 'rgb',
+                export=True,
+                path=path,
+                filename=f'Diachronic EI for eps={round(eps,3)} dim={dim} pha_divs={pha_divs}.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/plot.py b/evaluation/plot.py
new file mode 100644
index 0000000..ccbb090
--- /dev/null
+++ b/evaluation/plot.py
@@ -0,0 +1,138 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Mon Aug 21 16:12:53 2023
+
+@author: astral
+"""
+
+import sys
+import os
+
+import pyqtgraph as qt
+import pyqtgraph.exporters
+from PyQt5.QtGui import QFont
+
+
+import math as m
+import numpy as np
+
+def __get_color(factor, gamma=0.8):
+    frequency=380+factor*(750-380)
+    
+    lightfrequency = 0.4*(3*np.log10(frequency/2)-2)/4
+    
+    wavelength = 300/lightfrequency
+    
+    '''steps of 10Hz: 22 means 220Hz'''
+
+    '''This converts a given wavelength of light to an 
+    approximate RGB color value. The wavelength must be given
+    in values between 0 and 1 for 0=380nm through 1=750nm
+    (789 THz through 400 THz).
+
+    Based on code by Dan Bruton
+    http://www.physics.sfasu.edu/astro/color/spectra.html
+    '''
+
+    wavelength = float(wavelength)
+    if wavelength >= 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, 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() 
+    
+    ph = qt.plot()
+    ph.showGrid(x=True,y=True)
+    ph.setXRange(np.min(spaces), np.max(spaces))
+    # ph.setYRange(0.0, 6)
+    
+    #ph.setTitle(title='Susceptibility density evolution for different automaton sizes', offset=(1000,500))#.format(dim))
+    ph.setLabel('left', y_tag)
+    ph.setLabel('bottom', x_tag)
+    #ph.setXRange(0, 0.15)
+    
+    ph.addLegend(offset=(400, 30))
+    
+    
+    #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)]
+    
+    
+    for i in range(len(vals)):
+        s = ph.plot(spaces[i], vals[i], 
+                              name=names[i], pen=qt.mkPen(colors[i], width=linewidth))
+    
+    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:
+        if not os.path.exists(path):
+            os.makedirs(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)
+    
+    def handleAppExit():
+        # Add any cleanup tasks here
+        print("closing")
+    
+    if close:
+        ph.close()
+        
+    # app.aboutToQuit.connect(handleAppExit)
+    # app.exec_()
+    
\ No newline at end of file
diff --git a/examples/Simulate1Layer.py b/examples/Simulate1Layer.py
index 67adc50..535dc43 100644
--- a/examples/Simulate1Layer.py
+++ b/examples/Simulate1Layer.py
@@ -5,7 +5,7 @@ Created on Fri Aug 18 19:05:04 2023
 
 @author: astral
 """
-
+from datetime import datetime
 import numpy as np
 from neuropercolation import Simulate1Layer
 
@@ -23,4 +23,4 @@ for dim in [49,100]:
                        ).run(evolutions_per_second=30,
                              last_evolution_step=100000,
                              save_states=False)
-        print(f'Done eps={eps:.3f} at dim={dim}')
\ No newline at end of file
+        print(f'Done eps={eps:.3f} with dim={dim} at {datetime.now()}')
\ No newline at end of file
diff --git a/examples/Simulate2Layers.py b/examples/Simulate2Layers.py
new file mode 100644
index 0000000..7f71aee
--- /dev/null
+++ b/examples/Simulate2Layers.py
@@ -0,0 +1,37 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Fri Aug 18 19:05:04 2023
+
+@author: astral
+"""
+
+from datetime import datetime
+import numpy as np
+from neuropercolation import Simulate2Layers, Simulate4Layers
+
+eps_space = np.linspace(0.005,0.5,100)
+#eps_space = np.linspace(0.135,0.15,4)
+
+for dim in [8]:
+    for eps in eps_space:
+        eps = round(eps,3)
+        sim = Simulate4Layers(dim,
+                       eps,
+                       steps=100100,
+                       draw=None,
+                       save='all',
+                       path=f'/cloud/Public/_data/neuropercolation/4lay/steps=100100/dim={dim:02}/',
+                       )
+        print(f'Done eps={eps:.3f} with dim={dim} at {datetime.now()}')
+    
+    for eps in eps_space:
+        eps = round(eps,3)
+        sim = Simulate2Layers(dim,
+                       eps,
+                       steps=100100,
+                       draw=None,
+                       save='all',
+                       path=f'/cloud/Public/_data/neuropercolation/2lay/steps=100100/dim={dim:02}/',
+                       )
+        print(f'Done eps={eps:.3f} with dim={dim} at {datetime.now()}')
\ No newline at end of file
diff --git a/examples/Simulate4Layers.py b/examples/Simulate4Layers.py
new file mode 100644
index 0000000..6d7db01
--- /dev/null
+++ b/examples/Simulate4Layers.py
@@ -0,0 +1,32 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Fri Aug 18 19:05:04 2023
+
+@author: astral
+"""
+
+from datetime import datetime
+import numpy as np
+from neuropercolation import Simulate4Layers
+
+eps_space = np.linspace(0.005,0.5,100)
+#eps_space = np.linspace(0.135,0.15,4)
+
+stp = 1000100
+
+for batch in range(1,5):
+    for dim in [7]:
+        for eps in eps_space[1:41:2]:
+            eps = round(eps,3)
+            cons = [(n,n) for n in range(dim)]+[(n,(n+3)%7) for n in range(dim)]
+            sim = Simulate4Layers(dim,
+                           eps,
+                           coupling=cons,
+                           steps=stp,
+                           draw=None,
+                           res=2,
+                           save='simple',
+                           path=f'/cloud/Public/_data/neuropercolation/4lay/cons={len(cons)}_steps={stp}/dim={dim:02}/batch={batch}/',
+                           )
+            print(f'Done eps={eps:.3f} with dim={dim} at {datetime.now()}')
\ No newline at end of file