Merge remote-tracking branch 'origin/publication'

.gitignore

@@ -49,4 +49,5 @@ MANIFEST
 .venv*/
 
 # Output data
+data
 logs

scripts/spatial_maps/correlation_length_fit/correlation_length_fit.py

@@ -40,7 +40,7 @@ def correlation_length_fit_dict(traj, map_type='perlin_map', load=True):
     if load:
         try:
             fit_correlation_lengths_dict = np.load(DATA_FOLDER + map_type + traj_name + '_fit_correlation_lengths_dict.npy')
-            all_values_there = [scale in fit_correlation_lengths_dict for scale in scale_list]
+            all_values_there = [scale in fit_correlation_lengths_dict.item().keys() for scale in scale_list]
             if all(all_values_there):
                 return fit_correlation_lengths_dict.item()
             print('Some scale values missing in corr. len. fit dictionary, new one will be generated')

scripts/spatial_maps/supplement_max_entropy_rule/figure_max_entropy_rule.py

@@ -1,12 +1,31 @@
 import matplotlib.pyplot as plt
 import numpy as np
 # from scripts.spatial_network.perlin_map.paper_figures_spatial_head_direction_network_perlin_map import FIGURE_SAVE_PATH
-FIGURE_SAVE_PATH = '../../../figures/figure_4_paper_perlin/'
+FIGURE_SAVE_PATH = '../../../figures/supplement_max_entropy/'
 
 from scripts.spatial_network.perlin_map.run_simulation_perlin_map import get_perlin_map
 from scripts.spatial_maps.supplement_max_entropy_rule.spatial_layout import SpatialLayout
-from scripts.spatial_maps.spatial_network_layout import Interneuron, get_excitatory_phases_in_inhibitory_axon, get_position_mesh, plot_neural_sheet
+from scripts.spatial_maps.spatial_network_layout import Interneuron, get_excitatory_phases_in_inhibitory_axon, get_position_mesh
 from scripts.spatial_maps.supplement_max_entropy_rule.spatial_layout import get_entropy
+
+def plot_neural_sheet(ex_positions, ex_tunings, axonal_clouds=None):
+    X, Y = get_position_mesh(ex_positions)
+
+    n_ex = int(np.sqrt(len(ex_positions)))
+    head_dir_preference = np.array(ex_tunings).reshape((n_ex, n_ex))
+
+    fig = plt.figure(figsize=(4, 4))
+    ax = fig.add_subplot(111)
+
+    c = ax.pcolor(X, Y, head_dir_preference, vmin=-np.pi, vmax=np.pi, cmap="twilight")
+    fig.colorbar(c, ax=ax, label="Orientation")
+    if axonal_clouds is not None:
+        for i, p in enumerate(axonal_clouds):
+            ell = p.get_ellipse()
+            # print(ell)
+            ax.add_artist(ell)
+
+
 # Get input map
 
 seed = 1
@@ -70,7 +89,7 @@ import matplotlib
 # In[20]:
 
 
-plt.style.use('../../spatial_network/perlin_map/figures.mplstyle')
+plt.style.use('../../model_figure/figures.mplstyle')
 
 # In[21]:
 

scripts/spatial_maps/supplement_max_entropy_rule/spatial_layout.py

@@ -72,7 +72,7 @@ class SpatialLayout(object):
         self.short_axis = short_axis
         self.sheet_size = sheet_size
 
-        ex_positions, ex_tunings = create_grid_of_excitatory_neurons(sheet_size, sheet_size, int(np.sqrt(NE)),
+        ex_positions, ex_tunings = create_grid_of_excitatory_neurons(sheet_size, int(np.sqrt(NE)),
                                                                      orientation_map)
         self.ex_positions = ex_positions
         self.ex_tunings = ex_tunings

scripts/spatial_maps/supplement_pinwheel_map/pinwheel_map.py

@@ -14,23 +14,23 @@ class PinwheelMap:
     # def w(self,r):
     #     return self.A*np.exp(-self.lambda_1*r**2)-self.B*np.exp(-self.lambda_2*r**2)
     def w(self, r):
-        corr_len = self.corr_len
-        return self.w_scale * np.piecewise(r, [r < int(corr_len / 2), r >= int(corr_len / 2) and r < corr_len,
-                                               r >= corr_len], [1.0, -1.0, 0.0])
+        scale = self.scale
+        return self.w_scale * np.piecewise(r, [r < int(scale / 2), r >= int(scale / 2) and r < scale,
+                                               r >= scale], [1.0, -1.0, 0.0])
 
     def bound(self, z_i):
         return 1.0 if norm(z_i) < 1.0 else 0.0
 
-    def __init__(self,x_dim,y_dim, corr_len, sheet_x=0, sheet_y=0, rnd_seed = -1):
+    def __init__(self,x_dim,y_dim, scale, sheet_x=0, sheet_y=0, seed = -1):
         self.x_dim = x_dim
         self.y_dim = y_dim
-        self.corr_len = int((corr_len/sheet_x)*x_dim) # in number of neurons?!?
+        self.scale = int((scale/sheet_x)*x_dim) # in number of neurons?!?
         self.sheet_x = sheet_x
         self.sheet_y = sheet_y
 
-        if rnd_seed != -1:
-            np.random.seed(rnd_seed)
-            self.rnd_seed = rnd_seed
+        if seed != -1:
+            np.random.seed(seed)
+            self.seed = seed
 
         self.angle_grid = np.random.rand(x_dim,y_dim)*2.0*np.pi-np.pi
         self.vec_grid = np.ndarray((x_dim,y_dim,2))
@@ -42,10 +42,10 @@ class PinwheelMap:
                 self.vec_grid[idx, idy] = vec
 
     def pinwheel_map_by_id(self, id_x, id_y):
-        return self.angle_grid[id_x,id_y]
+        return self.angle_grid[id_x, id_y]
 
-    def get_tuning_by_id(self,id_x,id_y):
-        return self.angle_grid[id_x,id_y]
+    def get_tuning_by_id(self, id_x, id_y):
+        return self.angle_grid[id_x, id_y]
 
     def tuning(self, x, y):
         id_x = int(x * (self.x_dim - 1) / self.sheet_x)
@@ -62,7 +62,7 @@ class PinwheelMap:
                 d_z_i = np.array([0., 0.])
                 for vec_j in vec_list:
                     (j_x, j_y), z_j = vec_j
-                    if np.abs(i_x - j_x) < self.corr_len and np.abs(i_y - j_y) < self.corr_len and vec_i != vec_j:
+                    if np.abs(i_x - j_x) < self.scale and np.abs(i_y - j_y) < self.scale and vec_i != vec_j:
                         r = np.sqrt((i_x - j_x) * (i_x - j_x) + (i_y - j_y) * (i_y - j_y))
                         d_z = z_j * self.w(r)
                         # print(d_z)

scripts/spatial_maps/supplement_pinwheel_map/pinwheel_map_generator.py

@@ -1,13 +1,13 @@
-import numpy as np
-
 from scripts.spatial_maps.supplement_pinwheel_map.pinwheel_map import PinwheelMap
+import warnings
+warnings.simplefilter(action='ignore', category=FutureWarning) # Otherwise pypet's FutureWarning spams the console output
 from pypet import Environment, cartesian_product
-from scripts.spatial_network.supplement_pinwheel_map.run_simulation_pinwheel_map import SCALE_RANGE, SEED_RANGE
+from scripts.spatial_network.supplement_pinwheel_map.run_simulation_pinwheel_map import SCALE_RANGE, SEED_RANGE, \
+    TRAJ_NAME_PINWHEEL_MAPS
 
 DATA_FOLDER = "../../../data/"
 LOG_FOLDER = "../../../logs/"
 
-TRAJ_NAME_PINWHEEL_MAPS = "precalculated_pinwheel_maps"
 
 def generate_and_save_seeded_map(traj):
     scale = traj.scale

scripts/spatial_network/perlin_map/final_analysis_and_plotting_perlin_map.py

@@ -356,7 +356,7 @@ def plot_example_input_maps(traj, figsize=(panel_size, panel_size), figname='per
 
     for idx, (ax, (fit_corr_len, seed)) in enumerate(zip(axes, corr_and_seed)):
         corresp_scale = get_closest_fit_correlation_length(traj, fit_corr_len)
-        par_dict = {'seed': seed, 'correlation_length': corresp_scale}
+        par_dict = {'seed': seed, 'scale': corresp_scale}
 
         run_name = filter_run_names_by_par_dict(traj, par_dict)[0]
         traj.f_set_crun(run_name)
@@ -410,7 +410,6 @@ def plot_example_input_maps(traj, figsize=(panel_size, panel_size), figname='per
             ell.set_linewidth(linewidth)
             ax.add_artist(ell)
 
-
         if idx == 8:
             start_scale_x = 100
             end_scale_x = start_scale_x + scale_length
@@ -543,70 +542,6 @@ def plot_orientation_maps_diff_scales(traj):
         plt.close(fig)
 
 
-def plot_orientation_maps_diff_scales_with_ellipse(traj):
-    n_ex = int(np.sqrt(traj.N_E))
-
-    scale_run_names = []
-    plot_scales = [0.0, 400.0, 800.0]
-    real_scales = [get_closest_scale(traj, scale) for scale in plot_scales]
-    for scale in plot_scales:
-        par_dict = {'seed': 1, 'correlation_length': get_closest_scale(traj, scale), 'long_axis': 100.}
-        scale_run_names.append(*filter_run_names_by_par_dict(traj, par_dict))
-
-    inhibitory_positions = traj.results.runs[scale_run_names[1]].inhibitory_axonal_cloud_array[:, :2]
-    selected_polar_neuron = get_neurons_close_to_given_position((300, 300), 1, inhibitory_positions)[0]
-    selected_circular_neuron = get_neurons_close_to_given_position((600, 600), 1, inhibitory_positions)[0]
-
-    width = panel_size
-    height = 1.5 * panel_size
-    fig = plt.figure(figsize=(width, height))
-
-    axes = ImageGrid(fig, 111, axes_pad=0.15,
-                     nrows_ncols=(len(plot_scales), 1))
-    for ax, run_name, scale in zip(axes, scale_run_names, real_scales):
-        traj.f_set_crun(run_name)
-
-        X, Y = get_position_mesh(traj.results.runs[run_name].ex_positions)
-
-        inhibitory_axonal_cloud_array = traj.results.runs[run_name].inhibitory_axonal_cloud_array
-        axonal_clouds = [Interneuron(p[0], p[1], traj.morphology.long_axis, traj.morphology.short_axis, p[2]) for p in
-                         inhibitory_axonal_cloud_array]
-
-        head_dir_preference = np.array(traj.results.runs[run_name].ex_tunings).reshape((n_ex, n_ex))
-        # TODO: Why was this transposed for plotting? (now changed)
-        c = ax.pcolor(X, Y, head_dir_preference, vmin=-np.pi, vmax=np.pi, cmap=head_direction_input_colormap)
-        # ax.set_title('Correlation length: {}'.format(scale))
-        # fig.colorbar(c, ax=ax, label="Tuning")
-        ax.set_xticks([])
-        ax.set_yticks([])
-
-        p1 = axonal_clouds[selected_polar_neuron]
-        ell = p1.get_ellipse()
-        ell._linewidth = 2.
-        ax.add_artist(ell)
-
-        p2 = axonal_clouds[selected_circular_neuron]
-        circ_r = 2 * np.sqrt(p2.a * p2.b)
-        circ = Ellipse((p2.x, p2.y), circ_r, circ_r, fill=False, zorder=2, edgecolor='k')
-        circ._linewidth = 2.
-
-        ax.add_artist(circ)
-        ax.annotate("{:.0f}".format(scale), xy=(1.0, 0.5), xytext=(2, 0), xycoords="axes fraction", textcoords="offset "
-                                                                                                               "points",
-                    va="center", ha="left")
-        remove_frame(ax)
-
-    # fig.suptitle('axonal cloud', fontsize=16)
-    traj.f_restore_default()
-    add_length_scale(axes[1], 200, -300, -100, 50, 50)
-    axes[0].annotate("input maps", xy=(-0.5, 1.0), xycoords="axes fraction", rotation=90, ha="left", va="top")
-    fig.tight_layout()
-
-    if save_figs:
-        plt.savefig(FIGURE_SAVE_PATH + 'F_orientation_maps_diff_scales_with_ellipse' + FIGURE_SAVE_FORMAT)
-        plt.close(fig)
-
-
 def plot_polar_plot_excitatory(traj, plot_run_names, selected_neuron_idx,
                                figname=FIGURE_SAVE_PATH + 'D_polar_plot_excitatory' + FIGURE_SAVE_FORMAT):
     directions = np.linspace(-np.pi, np.pi, traj.input.number_of_directions, endpoint=False)
@@ -630,7 +565,6 @@ def plot_polar_plot_excitatory(traj, plot_run_names, selected_neuron_idx,
     ax.plot([], [], color='white',label='          ')
     hdis = []
     for run_idx, run_name in enumerate(plot_run_names):
-        # label = traj.derived_parameters.runs[run_name].morphology.morph_label
         label = labels[run_idx]
         hdi = traj.results.runs[run_name].head_direction_indices[selected_neuron_idx]
         hdis.append(hdi)
@@ -640,13 +574,9 @@ def plot_polar_plot_excitatory(traj, plot_run_names, selected_neuron_idx,
         if run_max_rate > max_rate:
             max_rate = run_max_rate
         rate_plot.append(rate_plot[0])
-        # plt_label = '{:s} {:.2f}'.format(short_labels(label), hdi)
         plt_label = '{:s}'.format(short_labels(label))
         ax.plot(directions_plt, rate_plot, linewidth=line_widths[run_idx],
                 label=plt_label, color=colors[run_idx], linestyle=line_styles[run_idx], zorder=zorders[run_idx])
-    # ax.set_title('Firing Rate')
-    # ax.plot([0.0, 0.0], [0.0, 1.05 * max_rate], color='red', alpha=0.25, linewidth=4.)
-    # TODO: Set ticks for polar
     ticks = [40., 80.]
     ax.set_rgrids(ticks, labels=["{:.0f} Hz".format(ticklabel) if idx == len(ticks) - 1 else "" for idx, ticklabel in
                                  enumerate(ticks)], angle=60)
@@ -692,8 +622,6 @@ def plot_polar_plot_inhibitory(traj, plot_run_names, selected_neuron_idx, fignam
     ax.plot([], [], color='white',label='          ')
     hdis = []
     for run_idx, run_name in enumerate(plot_run_names[:2]):
-        # ax = axes[max_hdi_idx, run_idx]
-        # label = traj.derived_parameters.runs[run_name].morphology.morph_label
         label = labels[run_idx]
         hdi = traj.results.runs[run_name].inh_head_direction_indices[selected_neuron_idx]
         hdis.append(hdi)
@@ -701,14 +629,10 @@ def plot_polar_plot_inhibitory(traj, plot_run_names, selected_neuron_idx, fignam
         tuning_vectors = traj.results.runs[run_name].inh_tuning_vectors
         rate_plot = [np.linalg.norm(v) for v in tuning_vectors[plot_n_idx]]
         rate_plot.append(rate_plot[0])
-        # plt_label = '{:s} {:.2f}'.format(short_labels(label), hdi)
         plt_label = '{:s}'.format(short_labels(label))
         ax.plot(directions_plt, rate_plot, linewidth=line_widths[run_idx],
                 label=plt_label,
                 color=colors[run_idx], linestyle=line_styles[run_idx], zorder=zorders[run_idx])
-    # ax.set_title('Inh. Firing Rate')
-    # TODO: Set ticks for polar
-    # ticks = [np.round(max_rate / 3.), np.round(max_rate * 2. / 3.), np.round(max_rate)]
     ticks = [40., 80., 120.]
     ax.set_rgrids(ticks, labels=["{:.0f} Hz".format(ticklabel) if idx == len(ticks) - 1 else "" for idx, ticklabel in
                                  enumerate(ticks)], angle=60)
@@ -1095,7 +1019,6 @@ def plot_exc_and_inh_hdi_over_fit_corr_len(traj, plot_run_names, cut_off_dist):
     ax.set_ylabel('HDI')
     ax.axvline(corr_len_fit_dict[get_closest_scale(traj, 200.0)], color='k', linewidth=0.5, zorder=0)
     ax.set_ylim(0.0, 1.0)
-    # ax.set_xlim(0.0, np.max(corr_len_range))
     remove_frame(ax, ["right", "top"])
 
     tablelegend(ax, ncol=2, bbox_to_anchor=(1.1, 1.1), loc="upper right",
@@ -1109,6 +1032,7 @@ def plot_exc_and_inh_hdi_over_fit_corr_len(traj, plot_run_names, cut_off_dist):
         plt.savefig(FIGURE_SAVE_PATH + 'F_hdi_over_corr_len_scaled' + FIGURE_SAVE_FORMAT)
         plt.close(fig)
 
+
 def get_phase_difference(total_difference):
     """
     Map accumulated phase difference to shortest possible difference
@@ -1152,14 +1076,10 @@ def plot_firing_rate_similar_vs_diff_tuning(traj, plot_run_names, figsize=(9, 9)
 
     plot_fr_array = np.array(plot_fr_array)
 
-    # rects1 = ax.bar(x - width / 2, plot_fr_array[:, 0], width,
-    #                 label=r'$\theta_{pref} \pm 30°$')
-    # rects2 = ax.bar(x + width / 2, plot_fr_array[:, 1], width,
-    #                 label=r'$\theta_{opp} \pm 30°$')
     rects1 = ax.bar(x - width / 2, plot_fr_array[:, 0], width,
-                    label=r'sim.')
+                    label=r'preferred')
     rects2 = ax.bar(x + width / 2, plot_fr_array[:, 1], width,
-                    label=r'diff.')
+                    label=r'opposite')
 
     ax.set_xticks(x)
     ax.set_xticklabels(labels)
@@ -1290,16 +1210,15 @@ if __name__ == "__main__":
     print("\tpanel size: {:.2f} cm".format(panel_size * cm_per_inch))
     print()
 
-
     plot_colorbar(figsize=(0.8 * panel_size, 0.8 * panel_size), figname=FIGURE_SAVE_PATH + "A_i_colormap.svg")
-    plot_input_map(traj, run_name_dict[POLARIZED], figname="A_i_exemplary_input_map"+FIGURE_SAVE_FORMAT,
+
+    plot_input_map(traj, run_name_dict[POLARIZED], figname="A_i_exemplary_input_map" + FIGURE_SAVE_FORMAT,
                    figsize=(panel_size, panel_size))
 
-    # plot_example_input_maps(traj, figsize=(2 * panel_size, 2 * panel_size))
+    plot_example_input_maps(traj, figsize=(2 * panel_size, 2 * panel_size))
 
     plot_axonal_clouds(traj, plot_run_names)
 
-    #
     plot_firing_rate_map_excitatory(traj, direction_idx, plot_run_names, selected_neuron_excitatory)
 
     in_max_rate = plot_firing_rate_map_inhibitory(traj, direction_idx, plot_run_names, selected_inhibitory_neuron)

scripts/spatial_network/perlin_map/run_simulation_perlin_map.py

@@ -5,10 +5,10 @@ from pypet.brian2.parameter import Brian2MonitorResult
 
 from scripts.spatial_maps.spatial_network_layout import create_grid_of_excitatory_neurons, \
     create_interneuron_sheet_entropy_max_orientation, get_excitatory_neurons_in_inhibitory_axonal_clouds
-from scripts.spatial_maps.perlin_map.uniform_perlin_map import UniformPerlinMap
 from scripts.spatial_network.spatial_network_setup import get_synaptic_weights, \
     create_head_direction_input, ex_in_network
 from scripts.spatial_network.models import *
+from scripts.spatial_maps.perlin_map.uniform_perlin_map import UniformPerlinMap
 
 POLARIZED = 'ellipsoid'
 CIRCULAR = 'circular'
@@ -39,8 +39,7 @@ def spatial_network_with_entropy_maximisation(traj):
     N_E = traj.network.N_E
     N_I = traj.network.N_I
 
-    perlin_map = get_perlin_map(traj.input_map.scale, traj.input_map.seed,
-                                sheet_size, N_E)
+    perlin_map = get_perlin_map(traj.input_map.scale, traj.input_map.seed, sheet_size, N_E)
 
     ex_positions, ex_tunings = create_grid_of_excitatory_neurons(sheet_size, int(np.sqrt(N_E)), perlin_map)
 
@@ -137,7 +136,7 @@ def main():
     env = Environment(trajectory=TRAJ_NAME,
                       comment="Compare the head direction tuning for circular and polarized interneuron morphology, "
                               "when tuning orientations to maximise entropy of connected excitatory tunings.",
-                      multiproc=True, filename=DATA_FOLDER, ncores=3, overwrite_file=True, log_folder=LOG_FOLDER)
+                      multiproc=True, filename=DATA_FOLDER, ncores=24, overwrite_file=True, log_folder=LOG_FOLDER)
 
     traj = env.trajectory
 
@@ -179,13 +178,13 @@ def main():
 
     # To test if the simulation, analysis and plotting work as intended, use this setup to have a lower overall runtime
 
-    scale_range = [200.0]
-    seed_range = [1]
+    # scale_range = [200.0]
+    # seed_range = [1]
 
     # TODO: This is the parameter range used for the results shown in Fig 4. Uncomment for the full run.
 
-    # scale_range = np.linspace(1.0, 800.0, 24, endpoint=True).tolist()
-    # seed_range = range(10)
+    scale_range = np.linspace(1.0, 800.0, 24, endpoint=True).tolist()
+    seed_range = range(10)
 
     ellipsoid_parameter_exploration = {
         "morphology.long_axis": [100.0],

scripts/spatial_network/supplement_pinwheel_map/orientation_map_full_figure_runner.py

@@ -1,5 +0,0 @@
-import scripts.spatial_network.supplement_pinwheel_map.run_simulation_pinwheel_map as runner
-import scripts.spatial_network.supplement_pinwheel_map.analyse_orientation_map as analyser
-
-runner.main()
-analyser.main()

scripts/spatial_network/supplement_pinwheel_map/analyse_orientation_map.py

@@ -3,12 +3,21 @@ import multiprocessing
 import numpy as np
 from pypet import Trajectory
 from pypet.brian2 import Brian2MonitorResult
+from brian2.units import ms, khertz
 
-from scripts.spatial_network.spatial_network_setup import calculate_rates
 from scripts.spatial_network.supplement_pinwheel_map.run_simulation_pinwheel_map import DATA_FOLDER, TRAJ_NAME
 
 traj = None
 directions = None
+
+
+def calculate_rates(list_of_spike_times, min_time=0*ms):
+    isis = [np.ediff1d(np.extract(spike_times / ms > min_time / ms, spike_times / ms)) * ms for spike_times in
+            list_of_spike_times]
+    rates = np.array([1.0 / np.mean(isi / ms) if isi.shape[0] != 0 else 0 for isi in isis]) * khertz
+    return rates
+
+
 def get_spike_train_dictionary(number_of_neurons, spike_times, neuron_indices):
     spike_train_dict = {}
     for neuron_idx in range(number_of_neurons):
@@ -37,6 +46,7 @@ def get_firing_rate_dict_per_cell_and_direction(traj, run_name):
     traj.f_restore_default()
     return firing_rate_dict
 
+
 def get_firing_rate_array_per_cell_and_direction(traj, run_name):
     traj.f_set_crun(run_name)
     firing_rate_array = np.ndarray((traj.N_E, traj.input.number_of_directions))
@@ -49,7 +59,7 @@ def get_firing_rate_array_per_cell_and_direction(traj, run_name):
                                                      neuron_indices)
         ex_spike_rates = calculate_rates(ex_spike_trains.values())
         for n_idx, spike_rate in enumerate(ex_spike_rates):
-            #TODO: Why on earth does the unit vanish?
+            # TODO: Why on earth does the unit vanish?
             firing_rate_array[n_idx, dir_idx] = spike_rate
     traj.f_restore_default()
     return firing_rate_array
@@ -59,7 +69,7 @@ def get_head_direction_indices(directions, firing_rate_array):
     n_exc_neurons = firing_rate_array.shape[0]
     n_directions = len(directions)
 
-    tuning_vectors = np.zeros((n_exc_neurons,n_directions,2))
+    tuning_vectors = np.zeros((n_exc_neurons, n_directions, 2))
     rate_sums = np.zeros((n_exc_neurons,))
     for ex_id, ex_rates in enumerate(firing_rate_array):
         rate_sum = 0.
@@ -78,7 +88,6 @@ def get_head_direction_indices(directions, firing_rate_array):
     return head_direction_indices, tuning_vectors_return
 
 
-
 def get_inhibitory_firing_rate_array_per_cell_and_direction(traj, run_name):
     number_of_neurons = traj.N_I
     traj.f_set_crun(run_name)
@@ -89,7 +98,7 @@ def get_inhibitory_firing_rate_array_per_cell_and_direction(traj, run_name):
         traj.results.runs[run_name]['dir0'].spikes.i.t
     except:
         label = traj.derived_parameters.runs[run_name].morphology.morph_label
-        print('Cant find t for run {} with label {}'.format(run_name,label))
+        print('Cant find inhibitory spike times for run {}, probably because there where no spikes'.format(run_name, label))
         return np.zeros((number_of_neurons, traj.input.number_of_directions))
 
     for dir_idx, direction in enumerate(direction_names):
@@ -99,7 +108,7 @@ def get_inhibitory_firing_rate_array_per_cell_and_direction(traj, run_name):
                                                       neuron_indices)
         inh_spike_rates = calculate_rates(inh_spike_trains.values())
         for n_idx, spike_rate in enumerate(inh_spike_rates):
-            #TODO: Why on earth does the unit vanish?
+            # TODO: Why on earth does the unit vanish?
             firing_rate_array[n_idx, dir_idx] = spike_rate
     traj.f_restore_default()
     return firing_rate_array
@@ -109,7 +118,7 @@ def get_inhibitory_head_direction_indices(directions, firing_rate_array):
     n_inh_neurons = firing_rate_array.shape[0]
     n_directions = len(directions)
 
-    tuning_vectors = np.zeros((n_inh_neurons,n_directions,2))
+    tuning_vectors = np.zeros((n_inh_neurons, n_directions, 2))
     rate_sums = np.zeros((n_inh_neurons,))
     for inh_id, inh_rates in enumerate(firing_rate_array):
         rate_sum = 0.
@@ -128,28 +137,26 @@ def get_inhibitory_head_direction_indices(directions, firing_rate_array):
     return head_direction_indices, tuning_vectors_return
 
 
-def get_runs_with_circular_morphology(traj):
-    filtered_indices = traj.f_find_idx(('parameters.long_axis', 'parameters.short_axis'), lambda r1, r2: r1 == r2)
-    return filtered_indices
-
 def analyse_single_run(run_name):
     traj.f_set_crun(run_name)
     label = traj.derived_parameters.runs[run_name].morphology.morph_label
-    print(run_name, label)
+    print('Starting analysis of run {}'.format(run_name))
 
     if label != 'no conn':
-        print(run_name,'inh')
         inh_firing_rate_array = get_inhibitory_firing_rate_array_per_cell_and_direction(traj, run_name)
-        inh_head_direction_indices, inh_tuning_vectors = get_inhibitory_head_direction_indices(directions, inh_firing_rate_array)
+        inh_head_direction_indices, inh_tuning_vectors = get_inhibitory_head_direction_indices(directions,
+                                                                                               inh_firing_rate_array)
     else:
         n_inh = traj.N_I
         n_dir = traj.input.number_of_directions
         inh_firing_rate_array = np.zeros((n_inh, n_dir))
         inh_head_direction_indices = np.zeros(n_inh)
-        inh_tuning_vectors = np.zeros((n_inh,n_dir,2))
+        inh_tuning_vectors = np.zeros((n_inh, n_dir, 2))
 
     exc_firing_rate_array = get_firing_rate_array_per_cell_and_direction(traj, run_name)
     exc_head_direction_indices, exc_tuning_vectors = get_head_direction_indices(directions, exc_firing_rate_array)
+    print('Finishing analysis of run {}'.format(run_name))
+
     return exc_firing_rate_array, exc_head_direction_indices, exc_tuning_vectors, inh_firing_rate_array, inh_head_direction_indices, inh_tuning_vectors
 
 
@@ -160,6 +167,7 @@ def main():
     FULL_LOAD = 2
     traj.f_load(filename=DATA_FOLDER + TRAJ_NAME + ".hdf5", load_parameters=FULL_LOAD, load_results=FULL_LOAD)
 
+    # Use in conjunction with NO_LOADING to save on memory. Beware, that it might not always work as intended.
     # traj.v_auto_load = True
 
     directions = np.linspace(-np.pi, np.pi, traj.input.number_of_directions, endpoint=False)
@@ -168,7 +176,6 @@ def main():
 
     pool = multiprocessing.Pool()
     multi_proc_result = pool.map(analyse_single_run, run_names)
-    print(type(multi_proc_result), len(multi_proc_result))
 
     for idx, run_name in enumerate(run_names):
         traj.f_set_crun(run_name)
@@ -189,7 +196,7 @@ def main():
         traj.f_restore_default()
 
     traj.f_store()
-    print('Analysis done')
+
 
 if __name__ == "__main__":
     main()

scripts/spatial_network/supplement_pinwheel_map/run_and_analyze_pinwheel.py

@@ -0,0 +1,7 @@
+import scripts.spatial_network.supplement_pinwheel_map.run_simulation_pinwheel_map as runner
+import scripts.spatial_network.supplement_pinwheel_map.preliminary_analysis_pinwheel_map as analyser
+import scripts.spatial_maps.supplement_pinwheel_map.pinwheel_map_generator as map_generator
+
+map_generator.create_seeded_maps_in_scale_range()
+runner.main()
+analyser.main()

scripts/spatial_network/supplement_pinwheel_map/run_simulation_pinwheel_map.py

@@ -1,113 +1,65 @@
-import json
-import os
-
 import numpy as np
+from brian2 import BrianLogger
 from brian2.units import *
+import warnings
+warnings.simplefilter(action='ignore', category=FutureWarning) # Otherwise pypet's FutureWarning spams the console output
 from pypet import Environment, cartesian_product, Trajectory
 from pypet.brian2.parameter import Brian2MonitorResult
 
-from scripts.spatial_network import create_grid_of_excitatory_neurons, \
+from scripts.spatial_maps.spatial_network_layout import create_grid_of_excitatory_neurons, \
     create_interneuron_sheet_entropy_max_orientation, get_excitatory_neurons_in_inhibitory_axonal_clouds
-from scripts.ring_network.head_direction import ex_in_network
+from scripts.spatial_network.spatial_network_setup import get_synaptic_weights, \
+    create_head_direction_input, ex_in_network
+from scripts.spatial_network.models import *
 from scripts.spatial_maps.supplement_pinwheel_map.pinwheel_map import PinwheelMap
-from scripts.spatial_maps.supplement_pinwheel_map.pinwheel_map_generator import TRAJ_NAME_ORIENTATION_MAPS
-from scripts.spatial_network.spatial_network_setup import excitatory_eqs, excitatory_params, \
-    lif_interneuron_eqs, lif_interneuron_params, lif_interneuron_options, ei_synapse_model, ei_synapse_on_pre, \
-    ei_synapse_param, ie_synapse_model, ie_synapse_on_pre, ie_synapse_param, get_synaptic_weights, \
-    create_head_direction_input
+
 
 POLARIZED = 'ellipsoid'
 CIRCULAR = 'circular'
 NO_SYNAPSES = 'no conn'
-
-def get_local_data_folder():
-    data_folder = "../../../data/"
-    config_file_name = ".config.json"
-    if os.path.isfile(config_file_name):
-        with open(config_file_name)  as config_file:
-            config_dict = json.load(config_file)
-            data_folder = os.path.abspath(config_dict["data"])
-            print(data_folder)
-
-    return data_folder
-
-
 DATA_FOLDER = "../../../data/"
 LOG_FOLDER = "../../../logs/"
-SEED_RANGE = range(10)
-SCALE_RANGE = np.linspace(1.0, 650.0, 24, endpoint=True).tolist()
-
-TRAJ_NAME = "full_figure_orientation_map"
-
-
-def get_orientation_map(correlation_length, seed, sheet_size, N_E, data_folder=None):
-    if data_folder is None:
-        data_folder = DATA_FOLDER
-
-    traj = Trajectory(filename=data_folder + TRAJ_NAME_ORIENTATION_MAPS + ".hdf5")
+TRAJ_NAME = "spatial_network_pinwheel_map"
+TRAJ_NAME_PINWHEEL_MAPS = "precalculated_pinwheel_maps"
 
-    traj.f_load(index=-1, load_parameters=2, load_results=2)
-
-    available_lengths = sorted(list(set(traj.f_get("corr_len").f_get_range())))
-    closest_length = available_lengths[np.argmin(np.abs(np.array(available_lengths)-correlation_length))]
-    if closest_length!=correlation_length:
-        print("Warning: desired correlation length {:.1f} not available. Taking {:.1f} instead".format(
-            correlation_length, closest_length))
-    corr_len = closest_length
-    seed = seed
-
-    map_by_params = lambda x, y: x == corr_len and y == seed
-
-    idx_iterator = traj.f_find_idx(['corr_len', 'seed'], map_by_params)
-
-    # TODO: Since it has only one entry, maybe iterator can be replaced
-    for idx in idx_iterator:
-        traj.v_idx = idx
-        map_angle_grid = traj.crun.map
+# SCALE_RANGE = [200.0]
+# SEED_RANGE = [1]
 
-    number_of_excitatory_neurons_per_row = int(np.sqrt(N_E))
+SCALE_RANGE = np.linspace(1.0, 650.0, 24, endpoint=True).tolist()
+SEED_RANGE = range(10)
 
-    map = PinwheelMap(number_of_excitatory_neurons_per_row, number_of_excitatory_neurons_per_row,
-                      corr_len, sheet_size, sheet_size, seed)
-    map.angle_grid = map_angle_grid
 
-    return map.tuning
 
-def get_uniform_orientation_map(correlation_length, seed, sheet_size, N_E, data_folder=None):
+def get_uniform_pinwheel_map(scale, seed, sheet_size, N_E, data_folder=None):
     if data_folder is None:
         data_folder = DATA_FOLDER
 
-    traj = Trajectory(filename=data_folder + TRAJ_NAME_ORIENTATION_MAPS + ".hdf5")
+    traj = Trajectory(filename=data_folder + TRAJ_NAME_PINWHEEL_MAPS + ".hdf5")
 
     traj.f_load(index=-1, load_parameters=2, load_results=2)
 
-    available_lengths = sorted(list(set(traj.f_get("corr_len").f_get_range())))
-    closest_length = available_lengths[np.argmin(np.abs(np.array(available_lengths)-correlation_length))]
-    if closest_length!=correlation_length:
+    available_scales = sorted(list(set(traj.f_get("scale").f_get_range())))
+    closest_scale = available_scales[np.argmin(np.abs(np.array(available_scales)-scale))]
+    if closest_scale != scale:
         print("Warning: desired correlation length {:.1f} not available. Taking {:.1f} instead".format(
-            correlation_length, closest_length))
-    corr_len = closest_length
-    seed = seed
+            scale, closest_scale))
 
-    map_by_params = lambda x, y: x == corr_len and y == seed
+    map_by_params = lambda x, y: x == scale and y == seed
 
-    idx_iterator = traj.f_find_idx(['corr_len', 'seed'], map_by_params)
+    idx_iterator = traj.f_find_idx(['scale', 'seed'], map_by_params)
 
-    # TODO: Since it has only one entry, maybe iterator can be replaced
     for idx in idx_iterator:
         traj.v_idx = idx
-        map_angle_grid = traj.crun.map
+        map_angle_grid = traj.crun.pinwheel_map
 
     number_of_excitatory_neurons_per_row = int(np.sqrt(N_E))
 
     map = PinwheelMap(number_of_excitatory_neurons_per_row, number_of_excitatory_neurons_per_row,
-                      corr_len, sheet_size, sheet_size, seed)
+                      scale, sheet_size, sheet_size, seed)
 
-    # Uniformize orientation map
+    # Uniformize pinwheel map
 
     nrow = number_of_excitatory_neurons_per_row
-    size = sheet_size
-    scale = corr_len  # TODO: Probably this needs to be a linear interpolation
 
     n = map_angle_grid / np.pi
     m = np.concatenate(n)
@@ -118,34 +70,35 @@ def get_uniform_orientation_map(correlation_length, seed, sheet_size, N_E, data_
         m[sorted_idx[ii * idx:(ii + 1) * idx]] = val
     p_map = (m - nrow) / nrow
     map.angle_grid = p_map.reshape(nrow, -1) * np.pi
-    # self.map *= np.pi
-
-    # map.angle_grid = map_angle_grid
 
     return map.tuning
 
 
+def get_input_head_directions(traj):
+    directions = np.linspace(-np.pi, np.pi, traj.input.number_of_directions, endpoint=False)
+    return directions
+
+
 def spatial_network_with_entropy_maximisation(traj):
-    sheet_size = traj.map.sheet_size
+    sheet_size = traj.input_map.sheet_size
 
     N_E = traj.network.N_E
     N_I = traj.network.N_I
 
-    orientation_map = get_uniform_orientation_map(traj.map.correlation_length, traj.map.seed, sheet_size, N_E)
-    ex_positions, ex_tunings = create_grid_of_excitatory_neurons(sheet_size,
-                                                                 sheet_size,
-                                                                 int(np.sqrt(N_E)), orientation_map)
+    pinwheel_map = get_uniform_pinwheel_map(traj.input_map.scale, traj.input_map.seed, sheet_size, N_E)
+
+    ex_positions, ex_tunings = create_grid_of_excitatory_neurons(sheet_size, int(np.sqrt(N_E)), pinwheel_map)
 
     inhibitory_axon_long_axis = traj.morphology.long_axis
     inhibitory_axon_short_axis = traj.morphology.short_axis
 
-    entropy_maximisation_steps = traj.simulation.entropy_maximisation.steps if inhibitory_axon_long_axis != \
-                                                                               inhibitory_axon_short_axis else 1
+    entropy_maximisation_trial_orientations = traj.interneuron.entropy_maximisation.trial_orientations if \
+        inhibitory_axon_long_axis != inhibitory_axon_short_axis else 0
 
     inhibitory_axonal_clouds, ellipse_single_trial_entropy = create_interneuron_sheet_entropy_max_orientation(
         ex_positions, ex_tunings, N_I, inhibitory_axon_long_axis,
         inhibitory_axon_short_axis, sheet_size,
-        sheet_size, trial_orientations=entropy_maximisation_steps)
+        sheet_size, trial_orientations=entropy_maximisation_trial_orientations)
 
     ie_connections = get_excitatory_neurons_in_inhibitory_axonal_clouds(ex_positions, inhibitory_axonal_clouds)
 
@@ -166,21 +119,29 @@ def spatial_network_with_entropy_maximisation(traj):
     ex_in_weights, in_ex_weights = get_synaptic_weights(N_E, N_I, ie_connections, excitatory_synapse_strength,
                                                         inhibitory_synapse_strength)
 
-    sharpness = 1.0 / (traj.input.width) ** 2
-
-    directions = get_input_head_directions(traj)
-    for idx, dir in enumerate(directions):
+    input_directions = get_input_head_directions(traj)
+    for idx, direction in enumerate(input_directions):
         # We recreate the network here for every dir, which slows down the simulation quite considerably. Otherwise,
         # we get a problem with saving and restoring the spike times (0s spike for neuron 0)
-        net = ex_in_network(N_E, N_I, excitatory_eqs, excitatory_params, lif_interneuron_eqs,
+
+        net = ex_in_network(N_E, N_I,
+                            hodgkin_huxley_eqs_with_synaptic_conductance,
+                            hodgkin_huxley_params,
+                            lif_interneuron_eqs,
                             lif_interneuron_params,
-                            lif_interneuron_options, ei_synapse_model, ei_synapse_on_pre,
-                            ei_synapse_param,
-                            ex_in_weights, ie_synapse_model, ie_synapse_on_pre,
-                            ie_synapse_param, in_ex_weights, random_seed=2)
+                            lif_interneuron_options,
+                            delta_synapse_model,
+                            delta_synapse_on_pre,
+                            delta_synapse_param,
+                            ex_in_weights,
+                            exponential_synapse,
+                            exponential_synapse_on_pre,
+                            exponential_synapse_params,
+                            in_ex_weights,
+                            random_seed=traj.input_map.seed)
         input_to_excitatory_population = create_head_direction_input(traj.input.baseline * nA, ex_tunings,
-                                                                     sharpness,
-                                                                     traj.input.amplitude * nA, dir)
+                                                                     traj.input.phase_dispersion,
+                                                                     traj.input.amplitude * nA, direction)
 
         excitatory_neurons = net["excitatory_neurons"]
         excitatory_neurons.I = input_to_excitatory_population
@@ -216,25 +177,21 @@ def spatial_network_with_entropy_maximisation(traj):
     return 1
 
 
-def get_input_head_directions(traj):
-    directions = np.linspace(-np.pi, np.pi, traj.input.number_of_directions, endpoint=False)
-    return directions
-
-
 def main():
+    BrianLogger.suppress_name('method_choice')
+
+    # TODO: Set ncores to the desired number of processes to use by pypet
     env = Environment(trajectory=TRAJ_NAME,
-                      comment="Compare the head direction tuning for circular and ellipsoid interneuron morphology, "
+                      comment="Compare the head direction tuning for circular and polarized interneuron morphology, "
                               "when tuning orientations to maximise entropy of connected excitatory tunings.",
-                      multiproc=True, filename=DATA_FOLDER, ncores=0, overwrite_file=True, log_folder=LOG_FOLDER)
+                      multiproc=True, filename=DATA_FOLDER, ncores=24, overwrite_file=True, log_folder=LOG_FOLDER)
 
     traj = env.trajectory
 
-    traj.f_add_parameter_group("map")
-
-    traj.f_add_parameter("map.correlation_length", 200.0,
-                         comment="Correlation length of orientations in um")
-    traj.f_add_parameter("map.seed", 1, comment="Random seed for map generation.")
-    traj.f_add_parameter("map.sheet_size", 900, comment="Sheet size in um")
+    traj.f_add_parameter_group("input_map")
+    traj.f_add_parameter("input_map.scale", 200.0, comment="Scaling factor of the input map")
+    traj.f_add_parameter("input_map.seed", 1, comment="Random seed for input map generation.")
+    traj.f_add_parameter("input_map.sheet_size", 900, comment="Sheet size in um")
 
     traj.f_add_parameter_group("network")
     traj.f_add_parameter("network.N_E", 3600, comment="Number of excitatory neurons")
@@ -242,13 +199,15 @@ def main():
 
     traj.f_add_parameter_group("interneuron")
     traj.f_add_parameter("interneuron.tau", 7., comment="Interneuron timescale in ms")
+    traj.f_add_parameter("interneuron.entropy_maximisation.trial_orientations", 30,
+                         comment="Steps for entropy maximisation")
 
     traj.f_add_parameter_group("synapse")
     traj.f_add_parameter("synapse.inhibitory", 30.0, "Strength of conductance-based inhibitory synapse in nS.")
     traj.f_add_parameter("synapse.excitatory", 2.5, "Strength of conductance-based inhibitory synapse in mV.")
 
     traj.f_add_parameter_group("input")
-    traj.f_add_parameter("input.width", 1. / np.sqrt(2.5), comment="Standard deviation of incoming head direction input.")
+    traj.f_add_parameter("input.phase_dispersion", 2.5, comment="Standard deviation of incoming head direction input.")
     traj.f_add_parameter("input.baseline", 0.05, comment="Head direction input baseline")
     traj.f_add_parameter("input.amplitude", 0.6, comment="Head direction input amplitude")
     traj.f_add_parameter("input.number_of_directions", 12, comment="Number of probed directions")
@@ -262,23 +221,19 @@ def main():
     traj.f_add_parameter("morphology.short_axis", 25.0, comment="Short axis of axon ellipsoid")
 
     traj.f_add_parameter_group("simulation")
-    traj.f_add_parameter("simulation.entropy_maximisation.steps", 30, comment="Steps for entropy maximisation")
     traj.f_add_parameter("simulation.dt", 0.1, comment="Network simulation time step in ms")
     traj.f_add_parameter("simulation.duration", 1000, comment="Network simulation duration in ms")
 
-    correlation_length_range = SCALE_RANGE
-    # correlation_length_range = [200.0]
+    scale_range = SCALE_RANGE
     seed_range = SEED_RANGE
-    # seed_range = [1]
 
     ellipsoid_parameter_exploration = {
         "morphology.long_axis": [100.0],
         "morphology.short_axis": [25.0],
-        "map.correlation_length": correlation_length_range,
-        "map.seed": seed_range,
+        "input_map.scale": scale_range,
+        "input_map.seed": seed_range,
         "synapse.inhibitory": [30.],
         "synapse.excitatory": [2.5]
-        # "map.correlation_length": np.arange(0.0, 200.0, 50).tolist()
     }
 
     corresponding_circular_radius = float(np.sqrt(ellipsoid_parameter_exploration[
@@ -288,8 +243,8 @@ def main():
     circle_parameter_exploration = {
         "morphology.long_axis": [corresponding_circular_radius],
         "morphology.short_axis": [corresponding_circular_radius],
-        "map.correlation_length": ellipsoid_parameter_exploration["map.correlation_length"],
-        "map.seed": ellipsoid_parameter_exploration["map.seed"],
+        "input_map.scale": ellipsoid_parameter_exploration["input_map.scale"],
+        "input_map.seed": ellipsoid_parameter_exploration["input_map.seed"],
         "synapse.inhibitory": ellipsoid_parameter_exploration["synapse.inhibitory"],
         "synapse.excitatory": ellipsoid_parameter_exploration["synapse.excitatory"]
     }
@@ -297,8 +252,8 @@ def main():
     no_conn_parameter_exploration = {
         "morphology.long_axis": [corresponding_circular_radius],
         "morphology.short_axis": [corresponding_circular_radius],
-        "map.correlation_length": ellipsoid_parameter_exploration["map.correlation_length"],
-        "map.seed": ellipsoid_parameter_exploration["map.seed"],
+        "input_map.scale": ellipsoid_parameter_exploration["input_map.scale"],
+        "input_map.seed": ellipsoid_parameter_exploration["input_map.seed"],
         "synapse.inhibitory": [0.],
         "synapse.excitatory": [0.]
     }
@@ -317,5 +272,6 @@ def main():
 
     env.disable_logging()
 
+
 if __name__ == "__main__":
-    main()
+    main()