spike-based homeostasis: averaging over trials introduced

STDP/compare.py

@@ -103,7 +103,7 @@ def main(variant):
                       label='Arbor', marker='None')
     axes.flat[5].plot(brian_classical_data[:, 0], brian_classical_data[:, 1],
                       label='Brian', linestyle='dashed', marker='None')
-    axes.flat[5].set_xlabel(r'Time difference $\Delta t$ (ms)')
+    axes.flat[5].set_xlabel(r'Time delay $\Delta t$ (ms)')
     axes.flat[5].set_ylabel(r'Weight change $\Delta w$ (µS)')
 
     fig.legend(*axes[0, 0].get_legend_handles_labels(), loc="lower right")

spike_based_homeostasis/Makefile

@@ -2,18 +2,19 @@
 
 BUILD_CATALOGUE_SCRIPT := arbor-build-catalogue
 CONFIG := config_arbor.json
+NUM_TRIALS := 100
 
 homeostasis-catalogue.so: $(wildcard mechanisms/*.mod)
 	$(BUILD_CATALOGUE_SCRIPT) homeostasis mechanisms
 
-arbor_spikes.dat arbor_traces.dat: homeostasis-catalogue.so $(CONFIG) arbor_homeostasis.py
-	./arbor_homeostasis.py $(CONFIG) --catalogue ./$<
+arbor_input.dat arbor_traces.dat arbor_spikes.dat: homeostasis-catalogue.so $(CONFIG) arbor_homeostasis.py
+	./run_arbor.sh $(CONFIG) $(NUM_TRIALS) --catalogue ./$<
 
-brian2_spikes.dat brian2_traces.dat: config_brian2.json brian2_homeostasis.py
-	./brian2_homeostasis.py config_brian2.json
+brian2_input.dat brian2_traces.dat brian2_spikes.dat: config_brian2.json brian2_homeostasis.py
+	./run_brian2.sh config_brian2.json $(NUM_TRIALS)
 
-homeostasis.svg: arbor_input.dat arbor_traces.dat brian2_input.dat brian2_traces.dat compare.py
-	./compare.py $(CONFIG) --save homeostasis.svg
+homeostasis.svg: arbor_input.dat arbor_traces.dat arbor_spikes.dat brian2_input.dat brian2_traces.dat brian2_spikes.dat compare.py
+	./compare.py $(CONFIG) $(NUM_TRIALS) --save homeostasis.svg
 
 .PHONY: clean
 clean:

spike_based_homeostasis/arbor_homeostasis.py

@@ -29,13 +29,16 @@ def generate_poisson_spike_train(rate, duration):
 class SingleRecipe(arbor.recipe):
     """Implementation of Arbor simulation recipe."""
 
-    def __init__(self, config, catalogue):
+    def __init__(self, config, cat_file):
         """Initialize the recipe from config."""
 
         # The base C++ class constructor must be called first, to ensure that
         # all memory in the C++ class is initialized correctly.
         arbor.recipe.__init__(self)
 
+        # load custom catalogue
+        catalogue = arbor.load_catalogue(cat_file)
+
         self.the_props = arbor.neuron_cable_properties()
         self.the_cat = catalogue
         self.the_cat.extend(arbor.default_catalogue(), "")
@@ -192,19 +195,28 @@ def main(config_file, catalogue):
 
     spike_times = sorted([s[1] for s in sim.spikes()])
 
-    numpy.savetxt(f'arbor_traces.dat', data_stacked)
-    numpy.savetxt(f'arbor_spikes.dat', spike_times)
+    return data_stacked, spike_times
 
 
 if __name__ == '__main__':
 
     import argparse
     parser = argparse.ArgumentParser()
     parser.add_argument('config', help="name of config file")
+    parser.add_argument('num_trials', type=int, help="number of trials to consider")
     parser.add_argument('--catalogue', help="name of catalogue file library", default="homeostasis-catalogue.so")
     args = parser.parse_args()
 
-    # load custom catalogue
-    catalogue = arbor.load_catalogue(args.catalogue)
-
-    main(args.config, catalogue)
+    num_trials = args.num_trials
+    data_stacked_sum = numpy.array([])
+    spike_times_all = []
+    for i in range(num_trials):
+        data_stacked, spike_times = main(args.config, args.catalogue)
+        if data_stacked_sum.size == 0:
+            data_stacked_sum = numpy.array(data_stacked)
+        else:
+            data_stacked_sum += data_stacked
+        spike_times_all.extend(spike_times)
+
+    numpy.savetxt(f'arbor_traces.dat', data_stacked_sum / num_trials)
+    numpy.savetxt(f'arbor_spikes.dat', spike_times_all)

spike_based_homeostasis/brian2_homeostasis.py

@@ -88,7 +88,7 @@ def main(config_file):
     defaultclock.dt = config["simulation"]["dt"]*ms
     run(duration, report='text')
 
-    # collect data and store
+    # collect and store data
     input_dts = np.arange(0., len(varying_stimulus.values)*varying_stimulus.dt/ms, varying_stimulus.dt/ms)
     input_x = list(itertools.chain(*zip(input_dts, input_dts)))
     input_y = [0] + list(itertools.chain(*zip(varying_stimulus.values, varying_stimulus.values)))[:-1]
@@ -97,16 +97,31 @@ def main(config_file):
     spike_times_with = SM_with.t/ms
     spike_times_without = SM_without.t/ms
     np.savetxt(f'brian2_input.dat', input_stacked)
-    np.savetxt(f'brian2_traces.dat', results_stacked)
-    np.savetxt(f'brian2_spikes.dat', spike_times_with)
-    np.savetxt(f'brian2_spikes_without_homestasis.dat', spike_times_without)
+    #np.savetxt(f'brian2_traces.dat', results_stacked)
+    #np.savetxt(f'brian2_spikes.dat', spike_times_with)
+    #np.savetxt(f'brian2_spikes_without_homestasis.dat', spike_times_without)
+
+    return results_stacked, spike_times_with, spike_times_without
 
 
 if __name__ == '__main__':
 
     import argparse
     parser = argparse.ArgumentParser()
     parser.add_argument('config', help="name of config file")
+    parser.add_argument('num_trials', type=int, help="number of trials to consider")
     args = parser.parse_args()
 
-    main(args.config)
+    num_trials = args.num_trials
+    data_stacked_sum = np.array([])
+    spike_times_all = []
+    for i in range(num_trials):
+        data_stacked, spike_times_with, _ = main(args.config)
+        if data_stacked_sum.size == 0:
+            data_stacked_sum = np.array(data_stacked)
+        else:
+            data_stacked_sum += data_stacked
+        spike_times_all.extend(spike_times_with)
+
+    np.savetxt(f'brian2_traces.dat', data_stacked_sum / num_trials)
+    np.savetxt(f'brian2_spikes.dat', spike_times_all)

spike_based_homeostasis/compare.py

@@ -20,7 +20,7 @@ def compute_and_print_goodness(data_1, data_2):
     print(f"  RMSE = {rmse}")
 
 
-def plot(config,
+def plot(config, num_trials,
          arbor_input, arbor_traces, arbor_spikes,
          brian2_input, brian2_traces, brian2_spikes):
     """
@@ -47,10 +47,10 @@ def plot(config,
     times = np.arange(0, config["simulation"]["runtime"], dt)
     idxs = np.searchsorted(arbor_spikes, times)
     counts = np.array([len(arbor_spikes[start:stop]) for start, stop in zip(idxs, idxs[1:])])
-    arbor_firing_rate = counts / (dt/1000) # convert to Hz
+    arbor_firing_rate = counts / num_trials / (dt/1000) # convert to Hz
     idxs = np.searchsorted(brian2_spikes, times)
     counts = np.array([len(brian2_spikes[start:stop]) for start, stop in zip(idxs, idxs[1:])])
-    brian2_firing_rate = counts / (dt/1000) # convert to Hz
+    brian2_firing_rate = counts / num_trials / (dt/1000) # convert to Hz
     axes[2].plot(times[:-1], arbor_firing_rate, label='Arbor', color='#1f77b4')
     axes[2].plot(times[:-1], brian2_firing_rate, label='Brian', linestyle='dashed', color='#ff7f0e')
     axes[2].set_ylabel("Target neuron rate (Hz)")
@@ -70,6 +70,7 @@ def plot(config,
     import argparse
     parser = argparse.ArgumentParser()
     parser.add_argument('config', help="name of config file")
+    parser.add_argument('num_trials', type=int, help="number of trials in the data")
     parser.add_argument('--show', help="show plot",
                         action="store_true", default=False)
     parser.add_argument('--save', help="save to given file name", default=None)
@@ -89,7 +90,7 @@ def plot(config,
     brian2_spikes = np.loadtxt("brian2_spikes.dat")
     config = json.load(open(args.config, 'r'))
 
-    fig = plot(config, arbor_input, arbor_traces, arbor_spikes, brian2_input, brian2_traces, brian2_spikes)
+    fig = plot(config, args.num_trials, arbor_input, arbor_traces, arbor_spikes, brian2_input, brian2_traces, brian2_spikes)
 
     if args.save:
         fig.savefig(args.save)