From: Armaan Bhojwani <me@armaanb.net>
Date: Sat, 7 Aug 2021 02:52:52 +0000 (-0400)
Subject: Restructure repository
X-Git-Url: https://git.armaanb.net/?p=norepinephrine_wm.git;a=commitdiff_plain;h=34c295265845573ded0eb3d96d81d9a1adea1609

Restructure repository
---

diff --git a/.gitignore b/.gitignore
index 0ef23ea..8609bed 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,5 +1,3 @@
-*.png
-*.pkl
 *.pyc
-params.json
 out/
+venv
diff --git a/README b/README
index 9b1367a..35378eb 100644
--- a/README
+++ b/README
@@ -8,5 +8,28 @@ Marianne Bezaire.
 
 The work in this repository was in-part derived from "Effects of Guanfacine and
 Phenylephrine on a Spiking Neuron Model of Working Memory" by Peter Duggins,
-Terry Stewart, Xuan Choo, Chris Eliasmith.
+Terry Stewart, Xuan Choo, Chris Eliasmith. Some of those original files are
+preserved in old/.
 <https://github.com/psipeter/drugs_and_working_memory>
+
+Usage
+=====
+
+Create and activate a new virtual environment:
+
+    python3 -m venv venv
+    . ./venv/bin/activate
+
+Install necesary dependencies:
+
+	pip install -r requirements.txt
+
+Modify conf.py as needed:
+
+	vi conf.py
+
+Run simulation (use the -h flag to see more options):
+
+	python3 model.py
+
+Results are stored in ./out/.
diff --git a/model.py b/model.py
new file mode 100644
index 0000000..4a4b54b
--- /dev/null
+++ b/model.py
@@ -0,0 +1,112 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from os import mkdir
+
+# Refer to parameters.txt for documentation
+dt = 0.001
+t_cue = 1.0
+cue_scale = 1.0
+perceived = 0 # ???
+time_scale = 0.4
+steps = 100
+
+class Alpha(object):
+    def __init__(self):
+        self.x = np.logspace(0, 3, steps)
+        self.y = 1 / (1 + (999 * np.exp(-0.1233 * (self.x / self.offset))))
+        
+        self.gain = []
+        self.bias = []
+        
+    def calcgb(self, gaind, biasd):
+        for i in range(steps):
+            y = self.y[i]
+            self.gain.append(1 + gaind * y)
+            self.bias.append(1 + biasd * y)
+
+    def plot(self):
+        try:
+            mkdir("./out")
+        except FileExistsError:
+            pass
+
+        out = f"./out/{self.__class__.__name__}"
+        plt.plot(self.x, self.y)
+
+        plt.xlabel("Norepinephrine concentration (nM)")
+        plt.ylabel("Activity (%)")
+        plt.title("Norepinepherine Concentration vs Neuron Activity in " +
+            self.pretty)
+
+        plt.vlines(self.ki, 0, 1, linestyles="dashed")
+        plt.text(1.1 * self.ki, 0.1, "Affinity")
+
+        plt.hlines(0.5, 0, 1000, linestyles="dashed")
+        plt.text(1, 0.51, "50%")
+
+        plt.xscale("log")
+        gc = plt.gca()
+        gc.set_yticklabels(['{:.0f}%'.format(x * 100) for x in gc.get_yticks()])
+
+        plt.draw()
+        plt.savefig(f"{out}-norep-activity.png", dpi=1000)
+        
+        #######################################################################
+        
+        plt.plot(self.x, self.gain)
+        
+        plt.xlabel("Norepinephrine concentration (nM)")
+        plt.ylabel("Gain")
+        plt.title(f"Concentration vs Gain in {self.pretty}")
+
+        plt.draw()
+        plt.savefig(f"{out}-concentration-gain.png", dpi=1000)
+        
+        #######################################################################
+        
+        plt.plot(self.x, self.bias)
+        
+        plt.xlabel("Norepinephrine concentration (nM)")
+        plt.ylabel("Bias")
+        plt.title("Concentration vs Bias in " + self.pretty)
+        
+        plt.draw()
+        plt.savefig(f"{out}-concentration-bias.png", dpi=1000)
+
+class Alpha1(Alpha):
+    def __init__(self):
+        self.ki = 330
+        self.offset = 5.895
+        self.pretty = "α1 Receptor"
+        self.gaind = 0.1
+        self.biasd = 0.1
+        super().__init__()
+        
+    def calcgb(self):
+        super().calcgb(self.gaind, self.biasd)
+
+class Alpha2(Alpha):
+    def __init__(self):
+        self.ki = 56
+        self.offset = 1
+        self.pretty = "α2 Receptor"
+        self.gaind = -0.04
+        self.biasd = -0.02
+        super().__init__()
+        
+    def calcgb(self):
+        super().calcgb(self.gaind, self.biasd)
+
+def main():
+    plt.style.use("ggplot")
+    
+    a1 = Alpha1()
+    a1.calcgb()
+    a1.plot()
+
+    a2 = Alpha2()
+    a2.calcgb()
+    a2.plot()
+
+if __name__ == "__main__":
+    main()
diff --git a/norepinephrine_wm/helper.py b/norepinephrine_wm/helper.py
deleted file mode 100644
index 078b2f1..0000000
--- a/norepinephrine_wm/helper.py
+++ /dev/null
@@ -1,131 +0,0 @@
-def make_addon(N):
-    import string
-    import random
-    addon=str(''.join(random.choice(string.ascii_uppercase+string.digits) for _ in range(N)))
-    return addon
-
-def ch_dir():
-    #change directory for data and plot outputs
-    import os
-    import sys
-    root=os.getcwd()
-    addon=make_addon(9)
-    datadir=''
-    if sys.platform == "linux" or sys.platform == "linux2" or sys.platform == "darwin":
-        datadir=root+'/data/'+addon+'/' #linux or mac
-    elif sys.platform == "win32":
-        datadir=root+'\\data\\'+addon+'\\' #windows
-    os.makedirs(datadir)
-    os.chdir(datadir)
-    return datadir
-
-def empirical_dataframe():
-    import numpy as np
-    import pandas as pd
-    columns=('time','drug','accuracy','trial')
-    emp_times = [3.0,5.0,7.0,9.0]
-    emp_dataframe = pd.DataFrame(columns=columns,index=np.arange(0, 12))
-    pre_PHE=[0.972, 0.947, 0.913, 0.798]
-    pre_GFC=[0.970, 0.942, 0.882, 0.766]
-    post_GFC=[0.966, 0.928, 0.906, 0.838]
-    post_PHE=[0.972, 0.938, 0.847, 0.666]
-    q=0
-    for t in range(len(emp_times)):
-        emp_dataframe.loc[q]=[emp_times[t],'control (empirical)',np.average([pre_GFC[t],pre_PHE[t]]),0]
-        emp_dataframe.loc[q+1]=[emp_times[t],'PHE (empirical)',post_PHE[t],0]
-        emp_dataframe.loc[q+2]=[emp_times[t],'GFC (empirical)',post_GFC[t],0]
-        q+=3
-    return emp_dataframe
-
-def make_cues(P):
-    import numpy as np
-    trials=np.arange(P['n_trials'])
-    perceived=np.ones(P['n_trials']) #list of correctly perceived (not necessarily remembered) cues
-    rng=np.random.RandomState(seed=P['seed'])
-    cues=2*rng.randint(2,size=P['n_trials'])-1 #whether the cues is on the left or right
-    for n in range(len(perceived)):
-        if rng.rand()<P['misperceive']: perceived[n]=0
-    return trials, perceived, cues
-
-
-'''drug approximations'''
-import nengo
-class MySolver(nengo.solvers.Solver):
-    #When the simulator builds the network, it looks for a solver to calculate the decoders
-    #instead of the normal least-squares solver, we define our own, so that we can return
-    #the old decoders
-    def __init__(self,weights): #feed in old decoders
-        self.weights=False #they are not weights but decoders
-        self.my_weights=weights
-    def __call__(self,A,Y,rng=None,E=None): #the function that gets called by the builder
-        return self.my_weights.T, dict()
-
-def reset_gain_bias(P,model,sim,wm,wm_recurrent,wm_to_decision,drug):
-    #set gains and biases as a constant multiple of the old values
-    wm.gain = sim.data[wm].gain * P['drug_effect_biophysical'][drug][0]
-    wm.bias = sim.data[wm].bias * P['drug_effect_biophysical'][drug][1]
-    #set the solver of each of the connections coming out of wm using the custom MySolver class
-    #with input equal to the old decoders. We use the old decoders because we don't want the builder
-    #to optimize the decoders to the new gainbias/bias values, otherwise it would "adapt" to the drug
-    wm_recurrent.solver = MySolver(sim.model.params[wm_recurrent].weights)
-    wm_to_decision.solver=MySolver(sim.model.params[wm_to_decision].weights)
-    #rebuild the network to affect the gain/bias change
-    sim=nengo.Simulator(model,dt=P['dt'])
-    return sim
-
-'''dataframe initialization'''
-def primary_dataframe(P,sim,drug,trial,probe_wm,probe_output):
-    import numpy as np
-    import pandas as pd
-    columns=('time','drug','wm','output','correct','trial')
-    df_primary = pd.DataFrame(columns=columns, index=np.arange(0,len(P['timesteps'])))
-    i=0
-    for t in P['timesteps']:
-        wm_val=np.abs(sim.data[probe_wm][t][0])
-        output_val=sim.data[probe_output][t][0]
-        correct=get_correct(P['cues'][trial],output_val)
-        rt=t*P['dt_sample']
-        df_primary.loc[i]=[rt,drug,wm_val,output_val,correct,trial]
-        i+=1
-    return df_primary
-
-def firing_dataframe(P,sim,drug,trial,sim_wm,probe_spikes):
-    import numpy as np
-    import pandas as pd
-    columns=('time','drug','neuron-trial','tuning','firing_rate')
-    df_firing = pd.DataFrame(columns=columns, index=np.arange(0,len(P['timesteps'])*\
-        int(P['neurons_wm']*P['frac'])))
-    t_width = 0.2
-    t_h = np.arange(t_width / P['dt']) * P['dt'] - t_width / 2.0
-    h = np.exp(-t_h ** 2 / (2 * P['sigma_smoothing'] ** 2))
-    h = h / np.linalg.norm(h, 1)
-    j=0
-    for nrn in range(int(P['neurons_wm']*P['frac'])):
-        enc = sim_wm.encoders[nrn]
-        tuning = get_tuning(P,trial,enc)
-        spikes = sim.data[probe_spikes][:,nrn]
-        firing_rate = np.convolve(spikes,h,mode='same')
-        for t in P['timesteps']:
-            rt=t*P['dt_sample']
-            df_firing.loc[j]=[rt,drug,nrn+trial*P['neurons_wm'],tuning,firing_rate[t]]
-            j+=1
-        # print 'appending dataframe for neuron %s' %f
-    return df_firing
-
-def get_correct(cue,output_val):
-    if (cue > 0.0 and output_val > 0.0) or (cue < 0.0 and output_val < 0.0): correct=1
-    else: correct=0
-    return correct
-
-def get_tuning(P,trial,enc):
-    cue=P['cues'][trial]
-    enc_min_cutoff=P['enc_min_cutoff']
-    enc_max_cutoff=P['enc_max_cutoff']
-    if (cue > 0.0 and 0.0 < enc[0] < enc_min_cutoff) or \
-        (cue < 0.0 and 0.0 > enc[0] > -1.0*enc_min_cutoff): tuning='superweak'
-    if (cue > 0.0 and enc_min_cutoff < enc[0] < enc_max_cutoff) or \
-        (cue < 0.0 and -1.0*enc_max_cutoff < enc[0] < -1.0*enc_min_cutoff): tuning='weak'
-    elif (cue > 0.0 and enc[0] > enc_max_cutoff) or \
-        (cue < 0.0 and enc[0] < -1.0*enc_max_cutoff): tuning='strong'
-    else: tuning='nonpreferred'
-    return tuning
\ No newline at end of file
diff --git a/norepinephrine_wm/model.py b/norepinephrine_wm/model.py
deleted file mode 100644
index a81eb76..0000000
--- a/norepinephrine_wm/model.py
+++ /dev/null
@@ -1,208 +0,0 @@
-'''
-Peter Duggins, Terry Stewart, Xuan Choo, Chris Eliasmith
-Effects of Guanfacine and Phenylephrine on a Spiking Neuron Model of Working Memory
-June-August 2016
-Main Model File
-'''
-
-def run(params):
-    import nengo
-    from nengo.rc import rc
-    import numpy as np
-    import pandas as pd
-    from helper import reset_gain_bias, primary_dataframe, firing_dataframe
-
-    decision_type=params[0]
-    drug_type=params[1]
-    drug = params[2]
-    trial = params[3]
-    seed = params[4]
-    P = params[5]
-    dt=P['dt']
-    dt_sample=P['dt_sample']
-    t_cue=P['t_cue']
-    t_delay=P['t_delay']
-    drug_effect_neural=P['drug_effect_neural']
-    drug_effect_functional=P['drug_effect_functional']
-    drug_effect_biophysical=P['drug_effect_biophysical']
-    enc_min_cutoff=P['enc_min_cutoff']
-    enc_max_cutoff=P['enc_max_cutoff']
-    sigma_smoothing=P['sigma_smoothing']
-    frac=P['frac']
-    neurons_inputs=P['neurons_inputs']
-    neurons_wm=P['neurons_wm']
-    neurons_decide=P['neurons_decide']
-    time_scale=P['time_scale']
-    cue_scale=P['cue_scale']
-    tau=P['tau']
-    tau_wm=P['tau_wm']
-    noise_wm=P['noise_wm']
-    noise_decision=P['noise_decision']
-    perceived=P['perceived']
-    cues=P['cues']
-
-    if drug_type == 'biophysical': rc.set("decoder_cache", "enabled", "False") #don't try to remember old decoders
-    else: rc.set("decoder_cache", "enabled", "True")
-
-    def cue_function(t):
-        if t < t_cue and perceived[trial]!=0:
-            return cue_scale * cues[trial]
-        else: return 0
-
-    def time_function(t):
-        if t > t_cue:
-            return time_scale
-        else: return 0
-
-    def noise_bias_function(t):
-        import numpy as np
-        if drug_type=='neural':
-            return np.random.normal(drug_effect_neural[drug],noise_wm)
-        else:
-            return np.random.normal(0.0,noise_wm)
-
-    def noise_decision_function(t):
-        import numpy as np
-        if decision_type == 'default':
-            return np.random.normal(0.0,noise_decision)
-        elif decision_type == 'basal_ganglia':
-            return np.random.normal(0.0,noise_decision,size=2)
-
-    def inputs_function(x):
-        return x * tau_wm
-
-    def wm_recurrent_function(x):
-        if drug_type == 'functional':
-            return x * drug_effect_functional[drug]
-        else:
-            return x
-
-    def decision_function(x):
-        output=0.0
-        if decision_type=='default':
-            value=x[0]+x[1]
-            if value > 0.0: output = 1.0
-            elif value < 0.0: output = -1.0
-        elif decision_type=='basal_ganglia':
-            if x[0] > x[1]: output = 1.0
-            elif x[0] < x[1]: output = -1.0
-        return output
-
-    def BG_rescale(x): #rescales -1 to 1 into 0.3 to 1, makes 2-dimensional
-        pos_x = 0.5 * (x + 1)
-        rescaled = 0.4 + 0.6 * pos_x, 0.4 + 0.6 * (1 - pos_x)
-        return rescaled
-
-    '''model definition'''
-    with nengo.Network(seed=seed+trial) as model:
-
-        #Ensembles
-        cue = nengo.Node(output=cue_function)
-        time = nengo.Node(output=time_function)
-        inputs = nengo.Ensemble(neurons_inputs,2)
-        noise_wm_node = nengo.Node(output=noise_bias_function)
-        noise_decision_node = nengo.Node(output=noise_decision_function)
-        wm = nengo.Ensemble(neurons_wm,2)
-        decision = nengo.Ensemble(neurons_decide,2)
-        output = nengo.Ensemble(neurons_decide,1)
-
-        #Connections
-        nengo.Connection(cue,inputs[0],synapse=None)
-        nengo.Connection(time,inputs[1],synapse=None)
-        nengo.Connection(inputs,wm,synapse=tau_wm,function=inputs_function)
-        wm_recurrent=nengo.Connection(wm,wm,synapse=tau_wm,function=wm_recurrent_function)
-        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
-        probe_wm=nengo.Probe(wm[0],synapse=0.01,sample_every=dt_sample)
-        probe_spikes=nengo.Probe(wm.neurons, 'spikes', sample_every=dt_sample)
-        probe_output=nengo.Probe(output,synapse=None,sample_every=dt_sample)
-
-
-
-
-    '''SIMULATION'''
-    print 'Running drug \"%s\", trial %s...' %(drug,trial+1)
-    with nengo.Simulator(model,dt=dt) as sim:
-        if drug_type == 'biophysical': sim=reset_gain_bias(
-                P,model,sim,wm,wm_recurrent,wm_to_decision,drug)
-        sim.run(t_cue+t_delay)
-        df_primary=primary_dataframe(P,sim,drug,trial,probe_wm,probe_output)
-        df_firing=firing_dataframe(P,sim,drug,trial,sim.data[wm],probe_spikes)
-    return [df_primary, df_firing]
-
-
-
-'''MAIN'''
-def main():
-    import matplotlib.pyplot as plt
-    import seaborn as sns
-    import pandas as pd
-    import numpy as np
-    from helper import make_cues, empirical_dataframe
-    from pathos.helpers import freeze_support #for Windows
-    # import ipdb
-
-    '''Import Parameters from File'''
-    P=eval(open('parameters.txt').read()) #parameter dictionary
-    seed=P['seed'] #sets tuning curves equal to control before drug application
-    n_trials=P['n_trials']
-    drug_type=str(P['drug_type'])
-    decision_type=str(P['decision_type'])
-    drugs=P['drugs']
-    trials, perceived, cues = make_cues(P)
-    P['timesteps']=np.arange(0,int((P['t_cue']+P['t_delay'])/P['dt_sample']))
-    P['cues']=cues
-    P['perceived']=perceived
-
-    '''Multiprocessing'''
-    print "drug_type=%s, decision_type=%s, trials=%s..." %(drug_type,decision_type,n_trials)
-    freeze_support()
-    exp_params=[]
-    for drug in drugs:
-        for trial in trials:
-            exp_params.append([decision_type, drug_type, drug, trial, seed, P])
-    df_list=[run(exp_params[0]),run(exp_params[-1])]
-    primary_dataframe = pd.concat([df_list[i][0] for i in range(len(df_list))], ignore_index=True)
-    firing_dataframe = pd.concat([df_list[i][1] for i in range(len(df_list))], ignore_index=True)
-
-    '''Plot and Export'''
-    print 'Exporting Data...'
-    primary_dataframe.to_pickle('primary_data.pkl')
-    firing_dataframe.to_pickle('firing_data.pkl')
-    param_df=pd.DataFrame([P])
-    param_df.reset_index().to_json('params.json',orient='records')
-
-    print 'Plotting...'
-    emp_dataframe=empirical_dataframe()
-    sns.set(context='poster')
-    figure, (ax1, ax2) = plt.subplots(2, 1)
-    sns.tsplot(time="time",value="wm",data=primary_dataframe,unit="trial",condition='drug',ax=ax1,ci=95)
-    sns.tsplot(time="time",value="correct",data=primary_dataframe,unit="trial",condition='drug',ax=ax2,ci=95)
-    sns.tsplot(time="time",value="accuracy",data=emp_dataframe,unit='trial',condition='drug',
-                interpolate=False,ax=ax2)
-    sns.tsplot(time="time",value="accuracy",data=emp_dataframe, unit='trial',condition='drug',
-                interpolate=True,ax=ax2, legend=False)
-    ax1.set(xlabel='',ylabel='decoded $\hat{cue}$',xlim=(0,9.5),ylim=(0,1),
-                title="drug_type=%s, decision_type=%s, trials=%s" %(drug_type,decision_type,n_trials))
-    ax2.set(xlabel='time (s)',xlim=(0,9.5),ylim=(0.5,1),ylabel='DRT accuracy')
-    figure.savefig('primary_plots.png')
-
-    figure2, (ax3, ax4) = plt.subplots(1, 2)
-    if len(firing_dataframe.query("tuning=='weak'"))>0:
-        sns.tsplot(time="time",value="firing_rate",unit="neuron-trial",condition='drug',ax=ax3,ci=95,
-                data=firing_dataframe.query("tuning=='weak'").reset_index())
-    if len(firing_dataframe.query("tuning=='nonpreferred'"))>0:
-        sns.tsplot(time="time",value="firing_rate",unit="neuron-trial",condition='drug',ax=ax4,ci=95,
-                data=firing_dataframe.query("tuning=='nonpreferred'").reset_index())
-    ax3.set(xlabel='time (s)',xlim=(0.0,9.5),ylim=(0,250),ylabel='Normalized Firing Rate',title='Preferred Direction')
-    ax4.set(xlabel='time (s)',xlim=(0.0,9.5),ylim=(0,250),ylabel='',title='Nonpreferred Direction')
-    figure2.savefig('firing_plots.png')
-
-    plt.show()
-
-if __name__=='__main__':
-    main()
\ No newline at end of file
diff --git a/norepinephrine_wm/newmodel.py b/norepinephrine_wm/newmodel.py
deleted file mode 100644
index 4a4b54b..0000000
--- a/norepinephrine_wm/newmodel.py
+++ /dev/null
@@ -1,112 +0,0 @@
-import numpy as np
-import matplotlib.pyplot as plt
-from os import mkdir
-
-# Refer to parameters.txt for documentation
-dt = 0.001
-t_cue = 1.0
-cue_scale = 1.0
-perceived = 0 # ???
-time_scale = 0.4
-steps = 100
-
-class Alpha(object):
-    def __init__(self):
-        self.x = np.logspace(0, 3, steps)
-        self.y = 1 / (1 + (999 * np.exp(-0.1233 * (self.x / self.offset))))
-        
-        self.gain = []
-        self.bias = []
-        
-    def calcgb(self, gaind, biasd):
-        for i in range(steps):
-            y = self.y[i]
-            self.gain.append(1 + gaind * y)
-            self.bias.append(1 + biasd * y)
-
-    def plot(self):
-        try:
-            mkdir("./out")
-        except FileExistsError:
-            pass
-
-        out = f"./out/{self.__class__.__name__}"
-        plt.plot(self.x, self.y)
-
-        plt.xlabel("Norepinephrine concentration (nM)")
-        plt.ylabel("Activity (%)")
-        plt.title("Norepinepherine Concentration vs Neuron Activity in " +
-            self.pretty)
-
-        plt.vlines(self.ki, 0, 1, linestyles="dashed")
-        plt.text(1.1 * self.ki, 0.1, "Affinity")
-
-        plt.hlines(0.5, 0, 1000, linestyles="dashed")
-        plt.text(1, 0.51, "50%")
-
-        plt.xscale("log")
-        gc = plt.gca()
-        gc.set_yticklabels(['{:.0f}%'.format(x * 100) for x in gc.get_yticks()])
-
-        plt.draw()
-        plt.savefig(f"{out}-norep-activity.png", dpi=1000)
-        
-        #######################################################################
-        
-        plt.plot(self.x, self.gain)
-        
-        plt.xlabel("Norepinephrine concentration (nM)")
-        plt.ylabel("Gain")
-        plt.title(f"Concentration vs Gain in {self.pretty}")
-
-        plt.draw()
-        plt.savefig(f"{out}-concentration-gain.png", dpi=1000)
-        
-        #######################################################################
-        
-        plt.plot(self.x, self.bias)
-        
-        plt.xlabel("Norepinephrine concentration (nM)")
-        plt.ylabel("Bias")
-        plt.title("Concentration vs Bias in " + self.pretty)
-        
-        plt.draw()
-        plt.savefig(f"{out}-concentration-bias.png", dpi=1000)
-
-class Alpha1(Alpha):
-    def __init__(self):
-        self.ki = 330
-        self.offset = 5.895
-        self.pretty = "α1 Receptor"
-        self.gaind = 0.1
-        self.biasd = 0.1
-        super().__init__()
-        
-    def calcgb(self):
-        super().calcgb(self.gaind, self.biasd)
-
-class Alpha2(Alpha):
-    def __init__(self):
-        self.ki = 56
-        self.offset = 1
-        self.pretty = "α2 Receptor"
-        self.gaind = -0.04
-        self.biasd = -0.02
-        super().__init__()
-        
-    def calcgb(self):
-        super().calcgb(self.gaind, self.biasd)
-
-def main():
-    plt.style.use("ggplot")
-    
-    a1 = Alpha1()
-    a1.calcgb()
-    a1.plot()
-
-    a2 = Alpha2()
-    a2.calcgb()
-    a2.plot()
-
-if __name__ == "__main__":
-    main()
diff --git a/norepinephrine_wm/parameters.txt b/norepinephrine_wm/parameters.txt
deleted file mode 100644
index d36549b..0000000
--- a/norepinephrine_wm/parameters.txt
+++ /dev/null
@@ -1,32 +0,0 @@
-{
-'seed':3, #for the simulator build process
-'n_trials':3,
-'n_processes':9,
-'dt':0.001, #timestep
-'dt_sample':0.01, #timestep for data recording through probes
-'t_cue':1.0, #duration of cue presentation
-'t_delay':8.0, #duration of delay period between cue and decision
-'misperceive':0.1, #chance of failing to perceive the cue, causing no info to go into WM
-
-'neurons_inputs':100, #neurons for the inputs ensemble
-'neurons_wm':1000, #neurons for workimg memory ensemble
-'neurons_decide':100, #neurons for decision or basal ganglia
-'time_scale':0.4, #how fast does the 'time' dimension accumulate in WM neurons
-'cue_scale':1.0, #how strong is the cueulus from the visual system
-'tau':0.01, #synaptic time constant between ensembles
-'tau_wm':0.1, #synapse on recurrent connection in wm
-'noise_wm':0.005, #std of full-spectrum white noise inject into wm.neurons
-'noise_decision':0.3, #std of added gaussian noise in the default decision procedure; 
-
-'decision_type':'default', #decision procedure: 'default', 'basal_ganglia'
-'drug_type':'biophysical', #drug simulation: 'neural','functional','biophysical'				
-'drugs':['control','PHE','GFC'], #names of the of drugs to simulate, see dicts below
-'drug_effect_neural':{'control':0.0,'PHE':-0.2,'GFC':0.5}, #mean of cue onto wm.neurons
-'drug_effect_biophysical':{'control':[1.0,1,0],'PHE':[0.99,1.02],'GFC':[1.05,0.95]}, #k_gain, k_bias
-'drug_effect_functional':{'control':1.0,'PHE':0.985,'GFC':1.03}, #k_multiply
-
-'enc_min_cutoff':0.3, #minimum cutoff for "weak" encoders in preferred directions
-'enc_max_cutoff':0.6, #maximum cutoff for "weak" encoders in preferred directions
-'sigma_smoothing':0.005, #gaussian smoothing applied to spike data to calculate firing rate
-'frac':0.01, #fraction of neurons in WM to add to firing rate dataframe and plot
-}
\ No newline at end of file
diff --git a/old/.gitignore b/old/.gitignore
new file mode 100644
index 0000000..12579e0
--- /dev/null
+++ b/old/.gitignore
@@ -0,0 +1,3 @@
+*
+!.py
+!parameters.txt
diff --git a/old/helper.py b/old/helper.py
new file mode 100644
index 0000000..078b2f1
--- /dev/null
+++ b/old/helper.py
@@ -0,0 +1,131 @@
+def make_addon(N):
+    import string
+    import random
+    addon=str(''.join(random.choice(string.ascii_uppercase+string.digits) for _ in range(N)))
+    return addon
+
+def ch_dir():
+    #change directory for data and plot outputs
+    import os
+    import sys
+    root=os.getcwd()
+    addon=make_addon(9)
+    datadir=''
+    if sys.platform == "linux" or sys.platform == "linux2" or sys.platform == "darwin":
+        datadir=root+'/data/'+addon+'/' #linux or mac
+    elif sys.platform == "win32":
+        datadir=root+'\\data\\'+addon+'\\' #windows
+    os.makedirs(datadir)
+    os.chdir(datadir)
+    return datadir
+
+def empirical_dataframe():
+    import numpy as np
+    import pandas as pd
+    columns=('time','drug','accuracy','trial')
+    emp_times = [3.0,5.0,7.0,9.0]
+    emp_dataframe = pd.DataFrame(columns=columns,index=np.arange(0, 12))
+    pre_PHE=[0.972, 0.947, 0.913, 0.798]
+    pre_GFC=[0.970, 0.942, 0.882, 0.766]
+    post_GFC=[0.966, 0.928, 0.906, 0.838]
+    post_PHE=[0.972, 0.938, 0.847, 0.666]
+    q=0
+    for t in range(len(emp_times)):
+        emp_dataframe.loc[q]=[emp_times[t],'control (empirical)',np.average([pre_GFC[t],pre_PHE[t]]),0]
+        emp_dataframe.loc[q+1]=[emp_times[t],'PHE (empirical)',post_PHE[t],0]
+        emp_dataframe.loc[q+2]=[emp_times[t],'GFC (empirical)',post_GFC[t],0]
+        q+=3
+    return emp_dataframe
+
+def make_cues(P):
+    import numpy as np
+    trials=np.arange(P['n_trials'])
+    perceived=np.ones(P['n_trials']) #list of correctly perceived (not necessarily remembered) cues
+    rng=np.random.RandomState(seed=P['seed'])
+    cues=2*rng.randint(2,size=P['n_trials'])-1 #whether the cues is on the left or right
+    for n in range(len(perceived)):
+        if rng.rand()<P['misperceive']: perceived[n]=0
+    return trials, perceived, cues
+
+
+'''drug approximations'''
+import nengo
+class MySolver(nengo.solvers.Solver):
+    #When the simulator builds the network, it looks for a solver to calculate the decoders
+    #instead of the normal least-squares solver, we define our own, so that we can return
+    #the old decoders
+    def __init__(self,weights): #feed in old decoders
+        self.weights=False #they are not weights but decoders
+        self.my_weights=weights
+    def __call__(self,A,Y,rng=None,E=None): #the function that gets called by the builder
+        return self.my_weights.T, dict()
+
+def reset_gain_bias(P,model,sim,wm,wm_recurrent,wm_to_decision,drug):
+    #set gains and biases as a constant multiple of the old values
+    wm.gain = sim.data[wm].gain * P['drug_effect_biophysical'][drug][0]
+    wm.bias = sim.data[wm].bias * P['drug_effect_biophysical'][drug][1]
+    #set the solver of each of the connections coming out of wm using the custom MySolver class
+    #with input equal to the old decoders. We use the old decoders because we don't want the builder
+    #to optimize the decoders to the new gainbias/bias values, otherwise it would "adapt" to the drug
+    wm_recurrent.solver = MySolver(sim.model.params[wm_recurrent].weights)
+    wm_to_decision.solver=MySolver(sim.model.params[wm_to_decision].weights)
+    #rebuild the network to affect the gain/bias change
+    sim=nengo.Simulator(model,dt=P['dt'])
+    return sim
+
+'''dataframe initialization'''
+def primary_dataframe(P,sim,drug,trial,probe_wm,probe_output):
+    import numpy as np
+    import pandas as pd
+    columns=('time','drug','wm','output','correct','trial')
+    df_primary = pd.DataFrame(columns=columns, index=np.arange(0,len(P['timesteps'])))
+    i=0
+    for t in P['timesteps']:
+        wm_val=np.abs(sim.data[probe_wm][t][0])
+        output_val=sim.data[probe_output][t][0]
+        correct=get_correct(P['cues'][trial],output_val)
+        rt=t*P['dt_sample']
+        df_primary.loc[i]=[rt,drug,wm_val,output_val,correct,trial]
+        i+=1
+    return df_primary
+
+def firing_dataframe(P,sim,drug,trial,sim_wm,probe_spikes):
+    import numpy as np
+    import pandas as pd
+    columns=('time','drug','neuron-trial','tuning','firing_rate')
+    df_firing = pd.DataFrame(columns=columns, index=np.arange(0,len(P['timesteps'])*\
+        int(P['neurons_wm']*P['frac'])))
+    t_width = 0.2
+    t_h = np.arange(t_width / P['dt']) * P['dt'] - t_width / 2.0
+    h = np.exp(-t_h ** 2 / (2 * P['sigma_smoothing'] ** 2))
+    h = h / np.linalg.norm(h, 1)
+    j=0
+    for nrn in range(int(P['neurons_wm']*P['frac'])):
+        enc = sim_wm.encoders[nrn]
+        tuning = get_tuning(P,trial,enc)
+        spikes = sim.data[probe_spikes][:,nrn]
+        firing_rate = np.convolve(spikes,h,mode='same')
+        for t in P['timesteps']:
+            rt=t*P['dt_sample']
+            df_firing.loc[j]=[rt,drug,nrn+trial*P['neurons_wm'],tuning,firing_rate[t]]
+            j+=1
+        # print 'appending dataframe for neuron %s' %f
+    return df_firing
+
+def get_correct(cue,output_val):
+    if (cue > 0.0 and output_val > 0.0) or (cue < 0.0 and output_val < 0.0): correct=1
+    else: correct=0
+    return correct
+
+def get_tuning(P,trial,enc):
+    cue=P['cues'][trial]
+    enc_min_cutoff=P['enc_min_cutoff']
+    enc_max_cutoff=P['enc_max_cutoff']
+    if (cue > 0.0 and 0.0 < enc[0] < enc_min_cutoff) or \
+        (cue < 0.0 and 0.0 > enc[0] > -1.0*enc_min_cutoff): tuning='superweak'
+    if (cue > 0.0 and enc_min_cutoff < enc[0] < enc_max_cutoff) or \
+        (cue < 0.0 and -1.0*enc_max_cutoff < enc[0] < -1.0*enc_min_cutoff): tuning='weak'
+    elif (cue > 0.0 and enc[0] > enc_max_cutoff) or \
+        (cue < 0.0 and enc[0] < -1.0*enc_max_cutoff): tuning='strong'
+    else: tuning='nonpreferred'
+    return tuning
\ No newline at end of file
diff --git a/old/model.py b/old/model.py
new file mode 100644
index 0000000..a81eb76
--- /dev/null
+++ b/old/model.py
@@ -0,0 +1,208 @@
+'''
+Peter Duggins, Terry Stewart, Xuan Choo, Chris Eliasmith
+Effects of Guanfacine and Phenylephrine on a Spiking Neuron Model of Working Memory
+June-August 2016
+Main Model File
+'''
+
+def run(params):
+    import nengo
+    from nengo.rc import rc
+    import numpy as np
+    import pandas as pd
+    from helper import reset_gain_bias, primary_dataframe, firing_dataframe
+
+    decision_type=params[0]
+    drug_type=params[1]
+    drug = params[2]
+    trial = params[3]
+    seed = params[4]
+    P = params[5]
+    dt=P['dt']
+    dt_sample=P['dt_sample']
+    t_cue=P['t_cue']
+    t_delay=P['t_delay']
+    drug_effect_neural=P['drug_effect_neural']
+    drug_effect_functional=P['drug_effect_functional']
+    drug_effect_biophysical=P['drug_effect_biophysical']
+    enc_min_cutoff=P['enc_min_cutoff']
+    enc_max_cutoff=P['enc_max_cutoff']
+    sigma_smoothing=P['sigma_smoothing']
+    frac=P['frac']
+    neurons_inputs=P['neurons_inputs']
+    neurons_wm=P['neurons_wm']
+    neurons_decide=P['neurons_decide']
+    time_scale=P['time_scale']
+    cue_scale=P['cue_scale']
+    tau=P['tau']
+    tau_wm=P['tau_wm']
+    noise_wm=P['noise_wm']
+    noise_decision=P['noise_decision']
+    perceived=P['perceived']
+    cues=P['cues']
+
+    if drug_type == 'biophysical': rc.set("decoder_cache", "enabled", "False") #don't try to remember old decoders
+    else: rc.set("decoder_cache", "enabled", "True")
+
+    def cue_function(t):
+        if t < t_cue and perceived[trial]!=0:
+            return cue_scale * cues[trial]
+        else: return 0
+
+    def time_function(t):
+        if t > t_cue:
+            return time_scale
+        else: return 0
+
+    def noise_bias_function(t):
+        import numpy as np
+        if drug_type=='neural':
+            return np.random.normal(drug_effect_neural[drug],noise_wm)
+        else:
+            return np.random.normal(0.0,noise_wm)
+
+    def noise_decision_function(t):
+        import numpy as np
+        if decision_type == 'default':
+            return np.random.normal(0.0,noise_decision)
+        elif decision_type == 'basal_ganglia':
+            return np.random.normal(0.0,noise_decision,size=2)
+
+    def inputs_function(x):
+        return x * tau_wm
+
+    def wm_recurrent_function(x):
+        if drug_type == 'functional':
+            return x * drug_effect_functional[drug]
+        else:
+            return x
+
+    def decision_function(x):
+        output=0.0
+        if decision_type=='default':
+            value=x[0]+x[1]
+            if value > 0.0: output = 1.0
+            elif value < 0.0: output = -1.0
+        elif decision_type=='basal_ganglia':
+            if x[0] > x[1]: output = 1.0
+            elif x[0] < x[1]: output = -1.0
+        return output
+
+    def BG_rescale(x): #rescales -1 to 1 into 0.3 to 1, makes 2-dimensional
+        pos_x = 0.5 * (x + 1)
+        rescaled = 0.4 + 0.6 * pos_x, 0.4 + 0.6 * (1 - pos_x)
+        return rescaled
+
+    '''model definition'''
+    with nengo.Network(seed=seed+trial) as model:
+
+        #Ensembles
+        cue = nengo.Node(output=cue_function)
+        time = nengo.Node(output=time_function)
+        inputs = nengo.Ensemble(neurons_inputs,2)
+        noise_wm_node = nengo.Node(output=noise_bias_function)
+        noise_decision_node = nengo.Node(output=noise_decision_function)
+        wm = nengo.Ensemble(neurons_wm,2)
+        decision = nengo.Ensemble(neurons_decide,2)
+        output = nengo.Ensemble(neurons_decide,1)
+
+        #Connections
+        nengo.Connection(cue,inputs[0],synapse=None)
+        nengo.Connection(time,inputs[1],synapse=None)
+        nengo.Connection(inputs,wm,synapse=tau_wm,function=inputs_function)
+        wm_recurrent=nengo.Connection(wm,wm,synapse=tau_wm,function=wm_recurrent_function)
+        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
+        probe_wm=nengo.Probe(wm[0],synapse=0.01,sample_every=dt_sample)
+        probe_spikes=nengo.Probe(wm.neurons, 'spikes', sample_every=dt_sample)
+        probe_output=nengo.Probe(output,synapse=None,sample_every=dt_sample)
+
+
+
+
+    '''SIMULATION'''
+    print 'Running drug \"%s\", trial %s...' %(drug,trial+1)
+    with nengo.Simulator(model,dt=dt) as sim:
+        if drug_type == 'biophysical': sim=reset_gain_bias(
+                P,model,sim,wm,wm_recurrent,wm_to_decision,drug)
+        sim.run(t_cue+t_delay)
+        df_primary=primary_dataframe(P,sim,drug,trial,probe_wm,probe_output)
+        df_firing=firing_dataframe(P,sim,drug,trial,sim.data[wm],probe_spikes)
+    return [df_primary, df_firing]
+
+
+
+'''MAIN'''
+def main():
+    import matplotlib.pyplot as plt
+    import seaborn as sns
+    import pandas as pd
+    import numpy as np
+    from helper import make_cues, empirical_dataframe
+    from pathos.helpers import freeze_support #for Windows
+    # import ipdb
+
+    '''Import Parameters from File'''
+    P=eval(open('parameters.txt').read()) #parameter dictionary
+    seed=P['seed'] #sets tuning curves equal to control before drug application
+    n_trials=P['n_trials']
+    drug_type=str(P['drug_type'])
+    decision_type=str(P['decision_type'])
+    drugs=P['drugs']
+    trials, perceived, cues = make_cues(P)
+    P['timesteps']=np.arange(0,int((P['t_cue']+P['t_delay'])/P['dt_sample']))
+    P['cues']=cues
+    P['perceived']=perceived
+
+    '''Multiprocessing'''
+    print "drug_type=%s, decision_type=%s, trials=%s..." %(drug_type,decision_type,n_trials)
+    freeze_support()
+    exp_params=[]
+    for drug in drugs:
+        for trial in trials:
+            exp_params.append([decision_type, drug_type, drug, trial, seed, P])
+    df_list=[run(exp_params[0]),run(exp_params[-1])]
+    primary_dataframe = pd.concat([df_list[i][0] for i in range(len(df_list))], ignore_index=True)
+    firing_dataframe = pd.concat([df_list[i][1] for i in range(len(df_list))], ignore_index=True)
+
+    '''Plot and Export'''
+    print 'Exporting Data...'
+    primary_dataframe.to_pickle('primary_data.pkl')
+    firing_dataframe.to_pickle('firing_data.pkl')
+    param_df=pd.DataFrame([P])
+    param_df.reset_index().to_json('params.json',orient='records')
+
+    print 'Plotting...'
+    emp_dataframe=empirical_dataframe()
+    sns.set(context='poster')
+    figure, (ax1, ax2) = plt.subplots(2, 1)
+    sns.tsplot(time="time",value="wm",data=primary_dataframe,unit="trial",condition='drug',ax=ax1,ci=95)
+    sns.tsplot(time="time",value="correct",data=primary_dataframe,unit="trial",condition='drug',ax=ax2,ci=95)
+    sns.tsplot(time="time",value="accuracy",data=emp_dataframe,unit='trial',condition='drug',
+                interpolate=False,ax=ax2)
+    sns.tsplot(time="time",value="accuracy",data=emp_dataframe, unit='trial',condition='drug',
+                interpolate=True,ax=ax2, legend=False)
+    ax1.set(xlabel='',ylabel='decoded $\hat{cue}$',xlim=(0,9.5),ylim=(0,1),
+                title="drug_type=%s, decision_type=%s, trials=%s" %(drug_type,decision_type,n_trials))
+    ax2.set(xlabel='time (s)',xlim=(0,9.5),ylim=(0.5,1),ylabel='DRT accuracy')
+    figure.savefig('primary_plots.png')
+
+    figure2, (ax3, ax4) = plt.subplots(1, 2)
+    if len(firing_dataframe.query("tuning=='weak'"))>0:
+        sns.tsplot(time="time",value="firing_rate",unit="neuron-trial",condition='drug',ax=ax3,ci=95,
+                data=firing_dataframe.query("tuning=='weak'").reset_index())
+    if len(firing_dataframe.query("tuning=='nonpreferred'"))>0:
+        sns.tsplot(time="time",value="firing_rate",unit="neuron-trial",condition='drug',ax=ax4,ci=95,
+                data=firing_dataframe.query("tuning=='nonpreferred'").reset_index())
+    ax3.set(xlabel='time (s)',xlim=(0.0,9.5),ylim=(0,250),ylabel='Normalized Firing Rate',title='Preferred Direction')
+    ax4.set(xlabel='time (s)',xlim=(0.0,9.5),ylim=(0,250),ylabel='',title='Nonpreferred Direction')
+    figure2.savefig('firing_plots.png')
+
+    plt.show()
+
+if __name__=='__main__':
+    main()
\ No newline at end of file
diff --git a/old/parameters.txt b/old/parameters.txt
new file mode 100644
index 0000000..d36549b
--- /dev/null
+++ b/old/parameters.txt
@@ -0,0 +1,32 @@
+{
+'seed':3, #for the simulator build process
+'n_trials':3,
+'n_processes':9,
+'dt':0.001, #timestep
+'dt_sample':0.01, #timestep for data recording through probes
+'t_cue':1.0, #duration of cue presentation
+'t_delay':8.0, #duration of delay period between cue and decision
+'misperceive':0.1, #chance of failing to perceive the cue, causing no info to go into WM
+
+'neurons_inputs':100, #neurons for the inputs ensemble
+'neurons_wm':1000, #neurons for workimg memory ensemble
+'neurons_decide':100, #neurons for decision or basal ganglia
+'time_scale':0.4, #how fast does the 'time' dimension accumulate in WM neurons
+'cue_scale':1.0, #how strong is the cueulus from the visual system
+'tau':0.01, #synaptic time constant between ensembles
+'tau_wm':0.1, #synapse on recurrent connection in wm
+'noise_wm':0.005, #std of full-spectrum white noise inject into wm.neurons
+'noise_decision':0.3, #std of added gaussian noise in the default decision procedure; 
+
+'decision_type':'default', #decision procedure: 'default', 'basal_ganglia'
+'drug_type':'biophysical', #drug simulation: 'neural','functional','biophysical'				
+'drugs':['control','PHE','GFC'], #names of the of drugs to simulate, see dicts below
+'drug_effect_neural':{'control':0.0,'PHE':-0.2,'GFC':0.5}, #mean of cue onto wm.neurons
+'drug_effect_biophysical':{'control':[1.0,1,0],'PHE':[0.99,1.02],'GFC':[1.05,0.95]}, #k_gain, k_bias
+'drug_effect_functional':{'control':1.0,'PHE':0.985,'GFC':1.03}, #k_multiply
+
+'enc_min_cutoff':0.3, #minimum cutoff for "weak" encoders in preferred directions
+'enc_max_cutoff':0.6, #maximum cutoff for "weak" encoders in preferred directions
+'sigma_smoothing':0.005, #gaussian smoothing applied to spike data to calculate firing rate
+'frac':0.01, #fraction of neurons in WM to add to firing rate dataframe and plot
+}
\ No newline at end of file
diff --git a/requirements.txt b/requirements.txt
new file mode 100644
index 0000000..e2c606f
--- /dev/null
+++ b/requirements.txt
@@ -0,0 +1,15 @@
+autopep8==1.5.7
+cycler==0.10.0
+kiwisolver==1.1.0
+matplotlib==3.4.2
+nengo==3.1.0
+numpy==1.21.1
+pandas==1.3.1
+Pillow==8.3.1
+pycodestyle==2.7.0
+pyparsing==2.4.7
+python-dateutil==2.8.2
+pytz==2021.1
+scipy==1.5.3
+six==1.16.0
+toml==0.10.2