More work

diff --git a/README b/README
index 35378eb881c0a4316f6db8de65947d44b3a54fba..312e3c3411de67902e784c830c55a16a5e2b2960 100644 (file)
--- a/README
+++ b/README
@@ -28,7 +28,7 @@ Modify conf.py as needed:
 
        vi conf.py
 
-Run simulation (use the -h flag to see more options):
+Run simulation:
 
        python3 model.py
 

diff --git a/conf.py b/conf.py
index 02c981c9e2ba9f967d5c63319917601232adec08..b9664989ef3d3b51eecf46e6e93b318393967e41 100644 (file)
--- a/conf.py
+++ b/conf.py
@@ -1,9 +1,12 @@
+import numpy as np
+
 dt = 0.01              # Time step
 t_cue = 1.0            # Duration of cue presentation
+t_delay = 8.0          # Duration of delay period between cue and decision
 cue_scale = 1.0        # How strong the cuelus is from the visual system
-perceived = 0          # ???
+misperceive = 0.1      # ???
 time_scale = 0.4       # ???
-steps = 100            # How fine to run the outer loop of the simulation
+steps = np.arange(750) # Steps to use
 noise_wm = 0.005       # Standard deviation of white noise added to WM
 noise_decision = 0.005 # Standard deviation of white noise added to decision
 neurons_decide = 100   # Number of neurons for decision
@@ -11,4 +14,3 @@ neurons_inputs = 100   # Number of neurons for inputs ensemble
 neurons_wm = 100       # Number of neurons for working memory ensemble
 tau_wm = 0.1           # Synapse on recurrent connection in wm
 tau = 0.01             # Synaptic time constant between ensembles
-t_delay = 8.0          # Duration of delay period between cue and decision

diff --git a/model.py b/model.py
index 0ea66c629cfde8c0818d0cc5db8bd4f97fe5862d..db3c69156877bbc37b53828b891ed987b764cfdc 100644 (file)
--- a/model.py
+++ b/model.py
@@ -1,15 +1,17 @@
 from datetime import datetime
 from os import mkdir
 import logging
+import pickle
 
 import matplotlib.pyplot as plt
 import matplotlib.ticker as mtick
 import nengo
 import numpy as np
+import pandas as pd
+from tqdm import tqdm
 
 exec(open("conf.py").read())
 
-
 def fmt_num(num, width=18):
     """
     Format number to string.
@@ -35,7 +37,14 @@ def time_function(t):
 
 
 def decision_function(x):
-    return 1.0 if x[0] + x[1] > 0.0 else -1.0
+    output = 0.0
+    value = x[0] + x[1]
+    if value > 0.0:
+        output = 1.0
+    elif value < 0.0:
+        output = -1.0
+    return output
+    #return 1.0 if x[0] + x[1] > 0.0 else -1.0
 
 
 class Alpha():
@@ -44,16 +53,15 @@ class Alpha():
     """
 
     def __init__(self):
-        self.x = np.logspace(0, 3, steps)
+        self.x = steps
         self.y = 1 / (1 + (999 * np.exp(-0.1233 * (self.x / self.offset))))
 
         self.gains = []
         self.biass = []
 
-        for i in range(steps):
-            y = self.y[i]
-            self.gains.append(1 + self.gaind * y)
-            self.biass.append(1 + self.biasd * y)
+        for i in range(len(steps)):
+            self.gains.append(1 + self.gaind * self.y[i])
+            self.biass.append(1 + self.biasd * self.y[i])
 
     def plot(self):
         out = f"./out/{self.__class__.__name__}"
@@ -145,49 +153,89 @@ class Simulation():
     def __init__(self):
         self.a1 = Alpha1()
         self.a2 = Alpha2()
+        self.num_spikes = np.ones(len(steps))
+        self.biass = np.ones(len(steps))
+        self.gains = np.ones(len(steps))
+        self.out = np.ones(3)
+        self.trial = 0
+
+        self.perceived = np.ones(3)  # correctly perceived (not necessarily remembered) cues
+        rng=np.random.RandomState(seed=111)
+        self.cues=2 * rng.randint(2, size=3)-1  # whether the cues is on the left or right
+        for n in range(len(self.perceived)):
+            if rng.rand() < misperceive:
+                self.perceived[n] = 0
+
+    def plot(self):
+        title = "Norepinephrine Concentration vs Spiking Rate"
+        logging.info("Plotting " + title)
+        plt.figure()
+        plt.plot(steps, self.num_spikes)
+
+        plt.xscale("log")
+
+        plt.xlabel("Norepinephrine concentration (nM)")
+        plt.ylabel("Spiking rate (spikes/time step)")
+        plt.title(title)
+
+        plt.draw()
+        plt.savefig("./out/concentration-spiking.png")
+
+    def cue_function(self, t):
+        if t < t_cue and self.perceived[self.trial] != 0:
+            return cue_scale * self.cues[self.trial]
+        else:
+            return 0
 
     def run(self):
-        for i in range(steps):
-            gain = self.a1.gains[i] + self.a2.gains[i] - 1
-            bias = self.a1.biass[i] + self.a2.biass[i] - 1
-            logging.info(f"gain: {fmt_num(gain)}, bias: {fmt_num(bias)}")
-
-            with nengo.Network() as net:
-                # Nodes
-                time_node = nengo.Node(output=time_function)
-                noise_wm_node = nengo.Node(output=noise_bias_function)
-                noise_decision_node = nengo.Node(
-                    output=noise_decision_function)
-
-                # Ensembles
-                wm = nengo.Ensemble(neurons_wm, 2)
-                wm.gain = np.full(wm.n_neurons, gain)
-                wm.bias = np.full(wm.n_neurons, bias)
-                decision = nengo.Ensemble(neurons_decide, 2)
-                inputs = nengo.Ensemble(neurons_inputs, 2)
-                output = nengo.Ensemble(neurons_decide, 1)
-
-                # Connections
-                nengo.Connection(time_node, inputs[1], synapse=None)
-                nengo.Connection(inputs, wm, synapse=tau_wm,
-                                 function=inputs_function)
-                wm_recurrent = nengo.Connection(wm, wm, synapse=tau_wm)
-                nengo.Connection(noise_wm_node, wm.neurons, synapse=tau_wm,
-                                 transform=np.ones((neurons_wm, 1)) * tau_wm)
-                wm_to_decision = nengo.Connection(
-                    wm[0], decision[0], synapse=tau)
-                nengo.Connection(noise_decision_node,
-                                 decision[1], synapse=None)
-                nengo.Connection(decision, output, function=decision_function)
-
-                # Probes
-                probes_wm = nengo.Probe(wm[0], synapse=0.01)
-                probe_output = nengo.Probe(output, synapse=None)
+        with nengo.Network() as net:
+            # Nodes
+            cue_node = nengo.Node(output=self.cue_function)
+            time_node = nengo.Node(output=time_function)
+            noise_wm_node = nengo.Node(output=noise_bias_function)
+            noise_decision_node = nengo.Node(
+                output=noise_decision_function)
+
+            # Ensembles
+            wm = nengo.Ensemble(neurons_wm, 2)
+            decision = nengo.Ensemble(neurons_decide, 2)
+            inputs = nengo.Ensemble(neurons_inputs, 2)
+            output = nengo.Ensemble(neurons_decide, 1)
+
+            # Connections
+            nengo.Connection(cue_node, inputs[0], synapse=None)
+            nengo.Connection(time_node, inputs[1], synapse=None)
+            nengo.Connection(inputs, wm, synapse=tau_wm,
+                             function=inputs_function)
+            wm_recurrent = nengo.Connection(wm, wm, synapse=tau_wm)
+            nengo.Connection(noise_wm_node, wm.neurons, synapse=tau_wm,
+                             transform=np.ones((neurons_wm, 1)) * tau_wm)
+            wm_to_decision = nengo.Connection(
+                wm[0], decision[0], synapse=tau)
+            nengo.Connection(noise_decision_node,
+                             decision[1], synapse=None)
+            nengo.Connection(decision, output, function=decision_function)
+
+            # Probes
+            probes_wm = nengo.Probe(wm[0], synapse=0.01)
+            probe_spikes = nengo.Probe(wm.neurons)
+            probe_output = nengo.Probe(output, synapse=None)
 
             # Run simulation
             with nengo.Simulator(net, dt=dt, progress_bar=False) as sim:
-                sim.run(t_cue + t_delay)
+                for i, _ in tqdm(enumerate(steps), total=len(steps)):
+                    wm.gain = (self.a1.gains[i] + self.a2.gains[i]) * sim.data[wm].gain
+                    wm.bias = (self.a1.biass[i] + self.a2.biass[i]) * sim.data[wm].gain
+                    for self.trial in range(3):
+                        logging.debug(f"Simulating: trial: {self.trial}, gain: {fmt_num(self.gains[i])}, bias: {fmt_num(self.biass[i])}")
+                        sim.run(t_cue + t_delay)
+                        self.out[self.trial] = np.count_nonzero(sim.data[probe_spikes])
+                    self.num_spikes[i] = np.average(self.out)
+
+        with open(f"out/{datetime.now().isoformat()}-spikes.pkl", "wb") as pout:
+            pickle.dump(self.num_spikes, pout)
 
+        self.plot()
 
 def main():
     logging.info("Initializing simulation")
@@ -203,8 +251,5 @@ if __name__ == "__main__":
 
     logging.basicConfig(filename=f"out/{datetime.now().isoformat()}.log",
                         level=logging.DEBUG)
-    console = logging.StreamHandler()
-    console.setLevel(logging.INFO)
-    logging.getLogger("").addHandler(console)
 
     main()

diff --git a/requirements.txt b/requirements.txt
index e2c606f57f0e6e184ab684fba7e2a8bdb84b4a74..0f714c8c805bfcdb82419a587716542570486a23 100644 (file)
--- a/requirements.txt
+++ b/requirements.txt
@@ -13,3 +13,4 @@ pytz==2021.1
 scipy==1.5.3
 six==1.16.0
 toml==0.10.2
+tqdm==4.62.0