Started working on the orientation maps and added uniformirty correction to the correlation length measure

scripts/interneuron_placement.py

@@ -7,7 +7,7 @@ from brian2.units import *
 from matplotlib.patches import Ellipse
 from tqdm import tqdm
 
-from scripts.spatial_maps.orientation_map import OrientationMap
+from scripts.spatial_maps.orientation_maps.orientation_map import OrientationMap
 
 
 class Pickle:
@@ -82,8 +82,11 @@ def get_correct_position_mesh(ex_positions):
     ymin = np.min(y)
     ymax = np.max(y)
     dy = (ymax - ymin) / (y.shape[0] - 1)
+    # print(np.arange(xmin, xmax + 2 * dx, dx) - dx / 2.)
+    # X, Y = np.meshgrid(np.arange(xmin, xmax + 2 * dx, dx) - dx / 2., np.arange(ymin, ymax + 2 * dy, dy) - dy / 2.)
+    # print(np.linspace(xmin, xmax + dx, n_ex + 1, endpoint=True) - dx / 2)
+    X, Y = np.meshgrid(np.linspace(xmin, xmax + dx, n_ex + 1, endpoint=True) - dx / 2, np.linspace(ymin, ymax + dy, n_ex + 1, endpoint=True) - dy / 2)
 
-    X, Y = np.meshgrid(np.arange(xmin, xmax + 2 * dx, dx) - dx / 2., np.arange(ymin, ymax + 2 * dy, dy) - dy / 2.)
 
     return X, Y
 

scripts/spatial_maps/orientation_map_generator_pypet.py

@@ -1,114 +0,0 @@
-import matplotlib.pyplot as plt
-import numpy as np
-import pandas as pd
-
-from scripts.interneuron_placement import create_grid_of_excitatory_neurons, plot_neural_sheet
-from scripts.spatial_maps.orientation_map import OrientationMap
-from pypet import Environment, cartesian_product, Trajectory, Result
-
-DATA_FOLDER = "../../data/"
-LOG_FOLDER = "../../logs/"
-
-TRAJ_NAME_ORIENTATION_MAPS = "precalculated_orientation_maps"
-
-def generate_and_save_seeded_map(traj):
-    corr_len = traj.corr_len
-    seed = traj.seed
-    print(corr_len,seed)
-    N_E = 900
-
-    sheet_x = 450
-    sheet_y = 450
-
-    number_of_excitatory_neurons_per_row = int(np.sqrt(N_E))
-
-    map = OrientationMap(number_of_excitatory_neurons_per_row + 1, number_of_excitatory_neurons_per_row + 1,
-                         corr_len, sheet_x, sheet_y, seed)
-
-    loaded = False
-    try:
-        map.load_orientation_map()
-        loaded = True
-    except:
-        print(
-            'No map yet with {}x{} pixels and {} pixel correllation length and {} seed'.format(map.x_dim, map.y_dim, map.corr_len, map.rnd_seed))
-    if not loaded:
-        map.improve(10)
-        traj.f_add_result('map', map.angle_grid, comment='map with {}x{} pixels and {} pixel correllation length and seed {}'.format(map.x_dim, map.y_dim, map.corr_len, map.rnd_seed))
-
-        print('map_{}x{}_pixels_{}_corr_pix_{}_seed is saved'.format(map.x_dim, map.y_dim,map.corr_len,map.rnd_seed))
-
-    return map.angle_grid
-
-
-def create_map_frame(traj, result_list):
-    print(len(result_list))
-    corr_len_range = traj.par.f_get('corr_len').f_get_range()
-    seed_range = traj.par.f_get('seed').f_get_range()
-
-    corr_len_index = sorted(set(corr_len_range))
-    seed_index = sorted(set(seed_range))
-    map_frame = pd.DataFrame(columns=corr_len_index, index=seed_index)
-    print(map_frame)
-    # This frame is basically a two dimensional table that we can index with our
-    # parameters
-
-    # Now iterate over the results. The result list is a list of tuples, with the
-    # run index at first position and our result at the second
-    for result_tuple in result_list:
-        run_idx = result_tuple[0]
-        print(run_idx)
-        result_map = result_tuple[1]
-        corr_len_val = corr_len_range[run_idx]
-        print(corr_len_range[run_idx])
-        seed_val = seed_range[run_idx]
-        # print(corr_len_val, seed_val, result_map)
-        map_frame.loc[corr_len_val, seed_val] = corr_len_val *seed_val
-
-    print(map_frame)
-
-    # Finally we going to store our new firing rate table into the trajectory
-    traj.f_add_result('summary.map_frames', map_frame=map_frame,
-                      comment='Contains precalculated orientation maps')
-
-def create_seeded_maps_in_scale_range():
-    import multiprocessing
-    import itertools
-
-    env = Environment(trajectory=TRAJ_NAME_ORIENTATION_MAPS, filename=DATA_FOLDER,
-                      file_title='orientation_map_set',
-                      comment='A number of precalculated orientation maps!',
-                      large_overview_tables=True,
-                      overwrite_file=True,
-                      multiproc=True,
-                      log_folder=LOG_FOLDER,
-                      ncores=0)
-    traj = env.trajectory
-
-    traj.f_add_parameter('corr_len', 200.0, comment='Correlation length')
-    traj.f_add_parameter('seed', 1, comment='Seed for random')
-
-    correlation_length_range = np.linspace(0.0, 400.0, 30).tolist()
-    # correlation_length_range = [150.0,200.0]
-    seed_range = range(15)
-    # seed_range = [1,2]
-    
-    traj.f_explore(cartesian_product({'corr_len': correlation_length_range, 'seed': seed_range}))
-
-    # env.add_postprocessing(create_map_frame)
-
-    env.run(generate_and_save_seeded_map)
-
-    env.disable_logging()
-
-    # corr_len_range = range(0, 451, 15)
-    # seed_range = range(10)
-    # pool_arguments = itertools.product(corr_len_range,seed_range)
-    # # print([*pool_arguments])
-    #
-    # pool = multiprocessing.Pool()
-    # pool.starmap(generate_and_save_seeded_map,[*pool_arguments])
-
-
-if __name__ == "__main__":
-    create_seeded_maps_in_scale_range()

scripts/spatial_maps/orientation_map.py

@@ -66,20 +66,23 @@ class OrientationMap:
         delta_z_list = []
         for vec_i in vec_list:
             (i_x, i_y), z_i = vec_i
-            # print(i_x,i_y,z_i)
-            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:
-                    r = np.sqrt((i_x - j_x) ** 2 + (i_y - j_y) ** 2)
-                    d_z = z_j * self.w(r) * self.bound(z_i)
-                    # print(d_z)
-                    d_z_i += d_z
-            delta_z_list.append(d_z_i)
-            # print(d_z_i)
+            if self.bound(z_i) == 1.0:
+                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:
+                        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)
+                        d_z_i += d_z
+                delta_z_list.append(d_z_i)
+            else:
+                delta_z_list.append(np.array([0., 0.]))
+
         for i, vec_i in enumerate(vec_list):
             (i_x, i_y), z_i = vec_i
-            vec_grid[i_x, i_y] += delta_z_list[i]
+            if self.bound(z_i) == 1.0:
+                vec_grid[i_x, i_y] += delta_z_list[i]
 
     def improve(self, iterations):
         for i in range(iterations):

scripts/spatial_maps/orientation_map_generator.py

@@ -1,9 +1,8 @@
 import matplotlib.pyplot as plt
 import numpy as np
-from numpy.linalg import norm
 from scripts.interneuron_placement import create_grid_of_excitatory_neurons, \
         plot_neural_sheet
-from scripts.spatial_maps.orientation_map import OrientationMap
+from scripts.spatial_maps.orientation_maps.orientation_map import OrientationMap
 
 def plotting(X, Y, Z, o_map, ax=None, field=True):
     for y_idx in range(Z.shape[1]):
@@ -107,17 +106,14 @@ def create_seeded_maps_in_scale_range():
 
 def test():
 
-    N_E = 900
+    dim = 30
 
-    sheet_x = 450
-    sheet_y = 450
+    size = 450
 
     corr_len = 150
 
-    number_of_excitatory_neurons_per_row = int(np.sqrt(N_E))
-
-    map = OrientationMap(number_of_excitatory_neurons_per_row + 1, number_of_excitatory_neurons_per_row + 1,
-                         corr_len, sheet_x, sheet_y)
+    map = OrientationMap(dim, dim,
+                         corr_len, size, size)
 
     map.improve(5)
 
@@ -135,5 +131,5 @@ def test():
 
 if __name__ == "__main__":
     # create_maps_in_scale_range()
-    # test()
-    create_seeded_maps_in_scale_range()
+    test()
+    # create_seeded_maps_in_scale_range()

scripts/spatial_maps/orientation_maps/orientation_map_generator_pypet.py

@@ -0,0 +1,65 @@
+import numpy as np
+import pandas as pd
+
+from scripts.spatial_maps.orientation_maps.orientation_map import OrientationMap
+from pypet import Environment, cartesian_product
+
+DATA_FOLDER = "../../../data/"
+LOG_FOLDER = "../../../logs/"
+
+TRAJ_NAME_ORIENTATION_MAPS = "precalculated_orientation_maps_test"
+
+def generate_and_save_seeded_map(traj):
+    corr_len = traj.corr_len
+    seed = traj.seed
+    print(corr_len,seed)
+
+    dim = 60
+    size = 900
+
+    map = OrientationMap(dim, dim, corr_len, size, size, seed)
+
+    loaded = False
+    try:
+        map.load_orientation_map()
+        loaded = True
+    except:
+        print(
+            'No map yet with {}x{} pixels and {} pixel correllation length and {} seed'.format(map.x_dim, map.y_dim, map.corr_len, map.rnd_seed))
+    if not loaded:
+        map.improve(10)
+        traj.f_add_result('map', map.angle_grid, comment='map with {}x{} pixels and {} pixel correllation length and seed {}'.format(map.x_dim, map.y_dim, map.corr_len, map.rnd_seed))
+
+        print('map_{}x{}_pixels_{}_corr_pix_{}_seed is saved'.format(map.x_dim, map.y_dim,map.corr_len,map.rnd_seed))
+
+    return map.angle_grid
+
+
+def create_seeded_maps_in_scale_range():
+    env = Environment(trajectory=TRAJ_NAME_ORIENTATION_MAPS, filename=DATA_FOLDER,
+                      file_title='orientation_map_set',
+                      comment='A number of precalculated orientation maps!',
+                      large_overview_tables=True,
+                      overwrite_file=True,
+                      multiproc=True,
+                      log_folder=LOG_FOLDER,
+                      ncores=0)
+    traj = env.trajectory
+
+    traj.f_add_parameter('corr_len', 200.0, comment='Correlation length')
+    traj.f_add_parameter('seed', 1, comment='Seed for random')
+
+    correlation_length_range = np.linspace(1.0, 800.0, 12, endpoint=True).tolist()
+    # correlation_length_range = [200.0, 400.0, 800.0]
+    seed_range = range(10)
+    # seed_range = [1]
+    
+    traj.f_explore(cartesian_product({'corr_len': correlation_length_range, 'seed': seed_range}))
+
+    env.run(generate_and_save_seeded_map)
+
+    env.disable_logging()
+
+
+if __name__ == "__main__":
+    create_seeded_maps_in_scale_range()

scripts/spatial_maps/orientation_map_plotter.py

@@ -1,7 +1,7 @@
 import numpy as np
 
 import matplotlib.pyplot as plt
-from scripts.spatial_maps.orientation_map import OrientationMap
+from scripts.spatial_maps.orientation_maps.orientation_map import OrientationMap
 from scripts.interneuron_placement import create_grid_of_excitatory_neurons, \
     plot_neural_sheet
 

scripts/spatial_maps/orientation_map_test_plot_pypet.py

@@ -1,19 +1,18 @@
 import matplotlib.pyplot as plt
 import numpy as np
-import pandas as pd
 
 from scripts.interneuron_placement import create_grid_of_excitatory_neurons, plot_neural_sheet
-from scripts.spatial_maps.orientation_map import OrientationMap
-from pypet import Environment, cartesian_product, Trajectory, Result
-
+from scripts.spatial_maps.orientation_maps.orientation_map import OrientationMap
+from pypet import Trajectory
+from scripts.spatial_maps.orientation_maps.orientation_map_generator_pypet import DATA_FOLDER, TRAJ_NAME_ORIENTATION_MAPS
 
 if __name__ == "__main__":
 
-    traj = Trajectory(filename='../../noise_maps/pypet_test/HDF/orientation_map_set.hdf5')
+    traj = Trajectory(filename=DATA_FOLDER + TRAJ_NAME_ORIENTATION_MAPS + ".hdf5")
 
     traj.f_load(index=-1, load_parameters=2, load_results=2)
 
-    corr_len = 200.0
+    corr_len = 800.0
     seed = 1
 
     map_by_params = lambda x, y: x == corr_len and y == seed
@@ -25,14 +24,14 @@ if __name__ == "__main__":
         traj.v_idx = idx
         map_angle_grid = traj.crun.map
 
-    N_E = 900
+    N_E = 3600
 
-    sheet_x = 450
-    sheet_y = 450
+    sheet_x = 900
+    sheet_y = 900
 
     number_of_excitatory_neurons_per_row = int(np.sqrt(N_E))
 
-    map = OrientationMap(number_of_excitatory_neurons_per_row + 1, number_of_excitatory_neurons_per_row + 1,
+    map = OrientationMap(number_of_excitatory_neurons_per_row, number_of_excitatory_neurons_per_row,
                          corr_len, sheet_x, sheet_y, seed)
     map.angle_grid = map_angle_grid
 

scripts/spatial_maps/orientation_perlin_comparison.py

@@ -1,13 +1,6 @@
-import numpy as np
 import matplotlib.pyplot as plt
-from opensimplex import OpenSimplex
-from pypet import Trajectory
-from scipy.signal import correlate2d
-from scipy.optimize import leastsq
-
-from scripts.spatial_maps.orientation_map import OrientationMap
-from scripts.spatial_maps.orientation_map_generator_pypet import TRAJ_NAME_ORIENTATION_MAPS
-from scripts.spatial_maps.uniform_perlin_map import UniformPerlinMap, plot_example, get_orientation_map
+
+from scripts.spatial_maps.uniform_perlin_map import UniformPerlinMap, get_orientation_map
 
 dim = 30
 size = 450

scripts/spatial_maps/placement_jitter_test.py

@@ -1,16 +1,10 @@
 import numpy as np
 import matplotlib.pyplot as plt
-from opensimplex import OpenSimplex
-from pypet import Trajectory
-from scipy.signal import correlate2d
-from scipy.optimize import leastsq
 
 from scripts.interneuron_placement import create_grid_of_excitatory_neurons, \
     create_interneuron_sheet_entropy_max_orientation, get_correct_position_mesh, \
     get_excitatory_neurons_in_inhibitory_axonal_clouds
-from scripts.spatial_maps.orientation_map import OrientationMap
-from scripts.spatial_maps.orientation_map_generator_pypet import TRAJ_NAME_ORIENTATION_MAPS
-from scripts.spatial_maps.uniform_perlin_map import UniformPerlinMap, plot_example, get_orientation_map
+from scripts.spatial_maps.uniform_perlin_map import get_orientation_map
 
 dim = 30
 size = 450

scripts/spatial_maps/simplex_input_tuning_correlation/correlation_length_distance_perlin.py

@@ -1,12 +1,7 @@
-import os
-import matplotlib.pyplot as plt
 import numpy as np
-from pypet import Trajectory
-from pypet.brian2 import Brian2MonitorResult
 from scripts.spatial_network.perlin.run_entropy_maximisation_perlin_map import DATA_FOLDER
-from scripts.interneuron_placement import create_grid_of_excitatory_neurons
-from scripts.spatial_maps.uniform_perlin_map import UniformPerlinMap
 from scipy.optimize import curve_fit
+import pingouin as pg
 
 
 def filter_run_names_by_par_dict(traj, par_dict):
@@ -34,53 +29,82 @@ def correlation_length_fit_dict(traj, load=True):
     corr_len_expl = traj.f_get('correlation_length').f_get_range()
     corr_len_list = np.unique(corr_len_expl)
 
+    seed_expl = traj.f_get('seed').f_get_range()
+    seed_list = np.unique(seed_expl)
+
     dim = 60
     size = 900
-    n_exc = dim * dim
 
     if not load:
         fit_correlation_lengths_dict = {}
 
-        ex_positions = traj.results.runs[traj.f_get_run_names()[0]].ex_positions
         for corr_len in corr_len_list:
+            fit_corr_lens_different_seeds = []
+            for map_seed in seed_list:
+                par_dict = {'seed': map_seed, 'correlation_length': corr_len, 'long_axis': 100.0}
+
+                run_name = filter_run_names_by_par_dict(traj, par_dict)
+
+                ex_tunings = traj.results.runs[run_name].ex_tunings
+
+                ex_tun_array = np.array(ex_tunings).reshape((dim, dim))
+
+                x_displacements = range(-dim + 1, dim)
+                y_displacements = range(-dim + 1, dim)
+                r_c_array = np.ndarray((len(x_displacements), len(y_displacements)))
+                dist_array = np.ndarray((len(x_displacements), len(y_displacements)))
+
+                for x_displacement in x_displacements:
+                    for y_displacement in y_displacements:
+                        a_1 = ex_tun_array[max(0, x_displacement): min(dim, dim + x_displacement),
+                              max(0, y_displacement): min(dim, dim + y_displacement)]
+
+                        a_2 = ex_tun_array[max(0, -x_displacement): min(dim, dim - x_displacement),
+                              max(0, -y_displacement): min(dim, dim - y_displacement)]
+
+                        if a_1.shape == (1, 1):
+                            # TODO: Still not sure which value to put here.
+                            r_c = np.nan
+                        else:
+                            # TODO: For some reason, some values come out as NaN or +/-inf
+                            r_c, p_c = pg.circ_corrcc(a_1.flatten(), a_2.flatten(), correction_uniform=True)
+
+                        r_c_array[x_displacement + dim - 1, y_displacement + dim - 1] = r_c
+                        dist_array[x_displacement + dim - 1, y_displacement + dim - 1] = np.sqrt(
+                            x_displacement ** 2 + y_displacement ** 2)
+
+                # correlations over distances
+
+                dist_array *= size / dim
+
+                distances = dist_array.flatten()
+                correlations = r_c_array.flatten()
 
-            par_dict = {'seed': 1, 'correlation_length': corr_len, 'long_axis': 100.0}
+                distances = distances[~np.isnan(correlations)]
+                correlations = correlations[~np.isnan(correlations)]
 
-            run_name = filter_run_names_by_par_dict(traj, par_dict)
+                # Binned mean and exponential fit
 
-            ex_tunings = traj.results.runs[run_name].ex_tunings
+                bins = np.linspace(0.0, size, 301)
 
-            bins = np.linspace(0.0, 800, 301)
+                binned_ids = np.digitize(distances, bins)
+                binned_correlations = [[] for _ in range(len(bins))]
+                for id, bin_id in enumerate(binned_ids):
+                    bin_corr = correlations[id]
+                    binned_correlations[bin_id - 1].append(bin_corr)
 
-            dist_list = []
-            phi_list = []
-            d_phi_list = []
-            for i in range(n_exc):
-                phi_list.append(ex_tunings[i])
-                for j in range(n_exc):
-                    xi, yi = ex_positions[i]
-                    xj, yj = ex_positions[j]
-                    dist = np.sqrt((xi - xj) ** 2 + (yi - yj) ** 2)
-                    dist_list.append(dist)
-                    d_phi = np.abs(phi_diff(ex_tunings[i] - ex_tunings[j]))
-                    d_phi_list.append(d_phi)
+                binned_mean_correlations = np.array([np.mean(bin) for bin in binned_correlations])
 
-            binned_ids = np.digitize(dist_list, bins)
-            binned_d_phi = [[] for _ in range(len(bins))]
-            for id, bin_id in enumerate(binned_ids):
-                current_d_phi = d_phi_list[id]
-                binned_d_phi[bin_id - 1].append(current_d_phi)
+                bins = bins[np.isfinite(binned_mean_correlations)]
+                binned_mean_correlations = binned_mean_correlations[np.isfinite(binned_mean_correlations)]
 
-            binned_mean = np.array([np.mean(bin) for bin in binned_d_phi])
+                # TODO: bins is wrong here, as it is the bin borders
 
-            nan_ids = np.isnan(binned_mean)
-            binned_mean = binned_mean[~nan_ids]
-            bins = bins[~nan_ids]
+                params, _ = curve_fit(exp_fit_func, bins, binned_mean_correlations, [0.5])
+                fit_corr_len = 1 / params[0]
+                fit_corr_lens_different_seeds.append(fit_corr_len)
 
-            params, _ = curve_fit(exp_fit_func, bins, np.pi / 2 - np.array(binned_mean), [1])
-            fit_corr_len = 1. / params[0]
-            print(corr_len, fit_corr_len)
-            fit_correlation_lengths_dict[corr_len] = fit_corr_len
+            fit_correlation_lengths_dict[corr_len] = np.mean(fit_corr_lens_different_seeds)
 
         np.save(DATA_FOLDER + 'fit_correlation_lengths_dict.npy', fit_correlation_lengths_dict)
         return fit_correlation_lengths_dict

scripts/spatial_maps/simplex_input_tuning_correlation/correlation_length_with_circular_correlation_formula.py

@@ -3,6 +3,8 @@ import numpy as np
 
 from scripts.interneuron_placement import create_grid_of_excitatory_neurons, get_correct_position_mesh
 from scripts.spatial_maps.uniform_perlin_map import UniformPerlinMap
+import pingouin as pg
+from scipy.optimize import curve_fit
 
 
 def phi_diff(d_phi):
@@ -10,21 +12,20 @@ def phi_diff(d_phi):
 
 
 def exp_fit_func(x, corr_len):
-    return np.pi / 2. * np.exp(-corr_len * x)
+    return np.exp(-corr_len * x)
 
 
 use_saved_array = False
-dim = 60
-size = 900
+dim = 120
+size = 1800
 n_exc = dim * dim
 
 corr_len_list = [200, 400, 800]
 
-fig, axes = plt.subplots(3, 2)
+fig, axes = plt.subplots(3, 4)
 
 for idx, (corr_len, axes_line) in enumerate(zip(corr_len_list, axes)):
 
-    bins = np.linspace(0.0, 800, 301)
 
     perlin_map = UniformPerlinMap(dim, dim, corr_len, size, size, 3)
 
@@ -35,60 +36,47 @@ for idx, (corr_len, axes_line) in enumerate(zip(corr_len_list, axes)):
     ex_tun_array = np.array(ex_tunings).reshape((dim, dim))
     # TODO: Why was this transposed for plotting? (now changed)
 
-    x_index = [pos[0] for pos in ex_positions]
-    y_index = [pos[1] for pos in ex_positions]
-
     x_displacements = range(-dim + 1, dim)
     y_displacements = range(-dim + 1, dim)
+    print(len(x_displacements))
     r_c_array = np.ndarray((len(x_displacements), len(y_displacements)))
+    dist_array = np.ndarray((len(x_displacements), len(y_displacements)))
 
     # For test purposes:
     # for i in range(dim):
     #     for j in range(dim):
-    #         ex_tun_array[i,j] = np.sin(2*np.pi*i/dim)
+    #         ex_tun_array[i,j] = np.sin(2*np.pi*np.sqrt(i**2+j**2)/(0.5*dim))
+    #         # ex_tun_array[i,j] = -np.pi+2*np.pi/dim*np.sqrt(i**2+j**2)
+
 
     for x_displacement in x_displacements:
         for y_displacement in y_displacements:
             a_1 = ex_tun_array[max(0, x_displacement): min(dim, dim + x_displacement), max(0, y_displacement): min(dim, dim + y_displacement)]
-            C_11 = np.sum(np.cos(a_1)) / len(a_1)
-            S_11 = np.sum(np.sin(a_1)) / len(a_1)
-            T_11 = np.arctan2(S_11, C_11)
 
             a_2 = ex_tun_array[max(0, -x_displacement): min(dim, dim - x_displacement), max(0, -y_displacement): min(dim, dim - y_displacement)]
-            C_21 = np.sum(np.cos(a_2)) / len(a_2)
-            S_21 = np.sum(np.sin(a_2)) / len(a_2)
-            T_21 = np.arctan2(S_21, C_21)
-
-            # TODO: hack to avoid arbitrary mean directions in effectively uniform data
-            T_11 = 0
-            T_21 = 0
-            r_c_denominator = np.sin(a_1 - T_11) * np.sin(a_2 - T_21)
-            r_c_norm = np.sqrt(np.sum(np.sin(a_1 - T_11) ** 2) * np.sum(np.sin(a_2 - T_21) ** 2))
-            if r_c_norm != 0:
-                r_c = np.sum(r_c_denominator) / r_c_norm
+
+            if a_1.shape == (1,1):
+                # TODO: Still not sure which value to put here.
+                r_c = np.nan
             else:
-                # TODO: Not sure which value to put here.
-                r_c = 0
-            if False and (x_displacement == 0 and y_displacement == 1):
-                plt.figure()
-                plt.suptitle("{:.2f}".format(r_c))
-                plt.hist2d(a_1.flatten() - T_11, a_2.flatten() - T_21)
-                print(T_11)
-                print(T_21)
-                print()
-                plt.show()
-                plt.close()
+                # TODO: For some reason, some values come out as NaN or +/-inf
+                r_c, p_c = pg.circ_corrcc(a_1.flatten(), a_2.flatten(), correction_uniform=True)
 
             r_c_array[x_displacement + dim - 1, y_displacement + dim - 1] = r_c
-    print(r_c_array)
+            dist_array[x_displacement + dim - 1, y_displacement + dim - 1] = np.sqrt(x_displacement**2 + y_displacement**2)
+
+    # Noise map
 
     X, Y = get_correct_position_mesh(ex_positions)
 
     head_dir_preference = np.array(ex_tunings).reshape((dim, dim))
-    c = axes_line[0].pcolor(X, Y, head_dir_preference, vmin=-np.pi, vmax=np.pi, cmap='hsv')
+    print(X.shape, Y.shape, ex_tun_array.shape)
+    c = axes_line[0].pcolor(X, Y, ex_tun_array, vmin=-np.pi, vmax=np.pi, cmap='hsv')
     axes_line[0].set_title('Scale factor: {}'.format(corr_len))
     axes_line[0].set_aspect('equal')
 
+    # Correlation map
+
     xmin = np.min(x_displacements)
     xmax = np.max(x_displacements)
     dx = (xmax - xmin) / (len(x_displacements) - 1)
@@ -98,9 +86,45 @@ for idx, (corr_len, axes_line) in enumerate(zip(corr_len_list, axes)):
     dy = (ymax - ymin) / (len(x_displacements) - 1)
 
     DX, DY = np.meshgrid(np.arange(xmin, xmax + 2 * dx, dx) - dx / 2., np.arange(ymin, ymax + 2 * dy, dy) - dy / 2.)
-    d = axes_line[1].pcolor(DX, DY, r_c_array)
+    d = axes_line[1].pcolor(DX, DY, r_c_array, vmin=-1., vmax=1.)
     axes_line[1].set_aspect('equal')
 
-    fig.colorbar(d, ax=axes_line[1], label="correlation")
+    # Scatter plot correlations over distances
+
+    dist_array *= size/dim
+
+    distances = dist_array.flatten()
+    correlations = r_c_array.flatten()
+
+    distances = distances[~np.isnan(correlations)]
+    correlations = correlations[~np.isnan(correlations)]
+
+    axes_line[2].scatter(distances, correlations, color='blue', alpha=0.05)
+
+    # Binned mean and exponential fit
+
+    bins = np.linspace(0.0, size, 301)
+
+    binned_ids = np.digitize(distances, bins)
+    binned_correlations = [[] for _ in range(len(bins))]
+    for id, bin_id in enumerate(binned_ids):
+        bin_corr = correlations[id]
+        binned_correlations[bin_id - 1].append(bin_corr)
+
+    binned_mean_correlations = np.array([np.mean(bin) for bin in binned_correlations])
+
+    bins = bins[np.isfinite(binned_mean_correlations)]
+    binned_mean_correlations = binned_mean_correlations[np.isfinite(binned_mean_correlations)]
+
+    # TODO: bins is wrong here, as it is the bin borders
+
+    axes_line[3].scatter(bins, binned_mean_correlations, color='blue')
+
+    params, _ = curve_fit(exp_fit_func, bins, binned_mean_correlations, [0.5])
+    corr_len = 1 / params[0]
+    axes_line[3].plot(bins, exp_fit_func(bins, params[0]), "r--")
+    axes_line[3].set_title('Circ. corr. and exp. fit. Corr. Len.: {:.1f}'.format(corr_len))
+
+
 
 plt.show()

scripts/spatial_maps/uniform_perlin_map.py

@@ -7,8 +7,8 @@ from scipy.optimize import leastsq
 
 from scripts.interneuron_placement import create_grid_of_excitatory_neurons, \
     create_interneuron_sheet_entropy_max_orientation, get_correct_position_mesh
-from scripts.spatial_maps.orientation_map import OrientationMap
-from scripts.spatial_maps.orientation_map_generator_pypet import TRAJ_NAME_ORIENTATION_MAPS
+from scripts.spatial_maps.orientation_maps.orientation_map import OrientationMap
+from scripts.spatial_maps.orientation_maps.orientation_map_generator_pypet import TRAJ_NAME_ORIENTATION_MAPS
 
 DATA_FOLDER = "../../data/"
 

scripts/spatial_network/analyse_hdi_optimization_orientation_map.py

@@ -1,16 +1,12 @@
 import numpy as np
 from brian2.units import *
-from pypet import Environment, cartesian_product, Trajectory
-from pypet.brian2.parameter import Brian2MonitorResult
-from scipy.optimize import differential_evolution
+from pypet import Environment
 import matplotlib.pyplot as plt
 
 from scripts.interneuron_placement import create_grid_of_excitatory_neurons, \
     create_interneuron_sheet_entropy_max_orientation, get_excitatory_neurons_in_inhibitory_axonal_clouds, Pickle, \
-    plot_neural_sheet, get_position_mesh
+    get_position_mesh
 from scripts.ring_network.head_direction import ex_in_network
-from scripts.spatial_maps.orientation_map import OrientationMap
-from scripts.spatial_maps.orientation_map_generator_pypet import TRAJ_NAME_ORIENTATION_MAPS
 from scripts.spatial_network.head_direction_index_over_noise_scale 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, \

scripts/spatial_network/hdi_optimization_via_syn_strength.py

@@ -2,7 +2,7 @@ import matplotlib.pyplot as plt
 import noise
 import numpy as np
 from brian2.units import *
-from scripts.spatial_maps.orientation_map import OrientationMap
+from scripts.spatial_maps.orientation_maps.orientation_map import OrientationMap
 from scripts.interneuron_placement import create_interneuron_sheet_entropy_max_orientation
 import scripts.models as modellib
 from scripts.ring_network.head_direction import get_head_direction_input, \

scripts/spatial_network/head_direction_index_over_noise_scale.py

@@ -13,7 +13,7 @@ from scripts.interneuron_placement import create_grid_of_excitatory_neurons, \
 from scripts.interneuron_placement import create_interneuron_sheet_entropy_max_orientation
 from scripts.ring_network.head_direction import get_head_direction_input, \
     ex_in_network
-from scripts.spatial_maps.orientation_map import OrientationMap
+from scripts.spatial_maps.orientation_maps.orientation_map import OrientationMap
 
 trials_per_scale = 1
 

scripts/spatial_network/inhibitory_input/run_entropy_maximisation_orientation_map_interneuron_input.py

@@ -6,8 +6,8 @@ from pypet.brian2.parameter import Brian2MonitorResult
 from scripts.interneuron_placement 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_maps.orientation_map import OrientationMap
-from scripts.spatial_maps.orientation_map_generator_pypet import TRAJ_NAME_ORIENTATION_MAPS
+from scripts.spatial_maps.orientation_maps.orientation_map import OrientationMap
+from scripts.spatial_maps.orientation_maps.orientation_map_generator_pypet import TRAJ_NAME_ORIENTATION_MAPS
 from scripts.spatial_network.head_direction_index_over_noise_scale 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, \

scripts/spatial_network/orientation_map/analyse_orientation_map.py

@@ -0,0 +1,195 @@
+import multiprocessing
+
+import numpy as np
+from pypet import Trajectory
+from pypet.brian2 import Brian2MonitorResult
+
+from scripts.spatial_network.head_direction_index_over_noise_scale import calculate_rates
+from scripts.spatial_network.orientation_map.run_orientation_map import DATA_FOLDER, TRAJ_NAME
+
+traj = None
+directions = None
+def get_spike_train_dictionary(number_of_neurons, spike_times, neuron_indices):
+    spike_train_dict = {}
+    for neuron_idx in range(number_of_neurons):
+        spike_train_dict[neuron_idx] = []
+
+    for neuron_idx, t in zip(neuron_indices, spike_times):
+        spike_train_dict[neuron_idx].append(t)
+    return spike_train_dict
+
+
+def get_firing_rate_dict_per_cell_and_direction(traj, run_name):
+    traj.f_set_crun(run_name)
+    firing_rate_dict = {}
+    direction_names = ["dir{:d}".format(idx) for idx in range(traj.input.number_of_directions)]
+    for idx in range(traj.N_E):
+        firing_rate_dict[idx] = []
+    for direction in direction_names:
+        number_of_neurons = traj.N_E
+        all_spike_times = traj.results.runs[run_name][direction].spikes.e.t
+        neuron_indices = traj.results.runs[run_name][direction].spikes.e.i
+        ex_spike_trains = get_spike_train_dictionary(number_of_neurons, all_spike_times,
+                                                     neuron_indices)
+        ex_spike_rates = calculate_rates(ex_spike_trains.values())
+        for idx, spike_rate in enumerate(ex_spike_rates):
+            firing_rate_dict[idx].append(spike_rate)
+    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))
+    direction_names = ["dir{:d}".format(idx) for idx in range(traj.input.number_of_directions)]
+    for dir_idx, direction in enumerate(direction_names):
+        number_of_neurons = traj.N_E
+        all_spike_times = traj.results.runs[run_name][direction].spikes.e.t
+        neuron_indices = traj.results.runs[run_name][direction].spikes.e.i
+        ex_spike_trains = get_spike_train_dictionary(number_of_neurons, all_spike_times,
+                                                     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?
+            firing_rate_array[n_idx, dir_idx] = spike_rate
+    traj.f_restore_default()
+    return firing_rate_array
+
+
+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))
+    rate_sums = np.zeros((n_exc_neurons,))
+    for ex_id, ex_rates in enumerate(firing_rate_array):
+        rate_sum = 0.
+        for dir_id, dir in enumerate(directions):
+            tuning_vectors[ex_id, dir_id] = np.array([np.cos(dir), np.sin(dir)]) * ex_rates[dir_id]
+            rate_sum += ex_rates[dir_id]
+        rate_sums[ex_id] = rate_sum
+    tuning_vectors_return = tuning_vectors.copy()
+    for ex_id in range(n_exc_neurons):
+        if rate_sums[ex_id] != 0.:
+            tuning_vectors[ex_id, :, :] /= rate_sums[ex_id]
+    tuning_vectors_summed = np.sum(tuning_vectors, axis=1)
+
+    head_direction_indices = np.array([np.linalg.norm(v) for v in tuning_vectors_summed])
+
+    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)
+    firing_rate_array = np.ndarray((number_of_neurons, traj.input.number_of_directions))
+    direction_names = ["dir{:d}".format(idx) for idx in range(traj.input.number_of_directions)]
+
+    try:
+        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))
+        return np.zeros((number_of_neurons, traj.input.number_of_directions))
+
+    for dir_idx, direction in enumerate(direction_names):
+        all_spike_times = traj.results.runs[run_name][direction].spikes.i.t
+        neuron_indices = traj.results.runs[run_name][direction].spikes.i.i
+        inh_spike_trains = get_spike_train_dictionary(number_of_neurons, all_spike_times,
+                                                      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?
+            firing_rate_array[n_idx, dir_idx] = spike_rate
+    traj.f_restore_default()
+    return firing_rate_array
+
+
+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))
+    rate_sums = np.zeros((n_inh_neurons,))
+    for inh_id, inh_rates in enumerate(firing_rate_array):
+        rate_sum = 0.
+        for dir_id, dir in enumerate(directions):
+            tuning_vectors[inh_id, dir_id] = np.array([np.cos(dir), np.sin(dir)]) * inh_rates[dir_id]
+            rate_sum += inh_rates[dir_id]
+        rate_sums[inh_id] = rate_sum
+    tuning_vectors_return = tuning_vectors.copy()
+    for inh_id in range(n_inh_neurons):
+        if rate_sums[inh_id] != 0.:
+            tuning_vectors[inh_id, :, :] /= rate_sums[inh_id]
+    tuning_vectors_summed = np.sum(tuning_vectors, axis=1)
+
+    head_direction_indices = np.array([np.linalg.norm(v) for v in tuning_vectors_summed])
+
+    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)
+
+    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)
+    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))
+
+    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)
+    return exc_firing_rate_array, exc_head_direction_indices, exc_tuning_vectors, inh_firing_rate_array, inh_head_direction_indices, inh_tuning_vectors
+
+
+def main():
+    global traj, directions
+    traj = Trajectory(TRAJ_NAME, add_time=False, dynamic_imports=Brian2MonitorResult)
+    NO_LOADING = 0
+    FULL_LOAD = 2
+    traj.f_load(filename=DATA_FOLDER + TRAJ_NAME + ".hdf5", load_parameters=FULL_LOAD, load_results=FULL_LOAD)
+
+    # traj.v_auto_load = True
+
+    directions = np.linspace(-np.pi, np.pi, traj.input.number_of_directions, endpoint=False)
+
+    run_names = traj.f_get_run_names()[::-1]
+
+    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)
+        traj.f_add_result('runs.$.firing_rate_array', multi_proc_result[idx][0],
+                          comment='The firing rates of the excitatory population')
+        traj.f_add_result('runs.$.head_direction_indices', multi_proc_result[idx][1],
+                          comment='The HDIs of the excitatory population')
+        traj.f_add_result('runs.$.tuning_vectors', multi_proc_result[idx][2],
+                          comment='The tuning vectors of the excitatory population')
+
+        traj.f_add_result('runs.$.inh_firing_rate_array', multi_proc_result[idx][3],
+                          comment='The firing rates of the inhibitory population')
+        traj.f_add_result('runs.$.inh_head_direction_indices', multi_proc_result[idx][4],
+                          comment='The HDIs of the inhibitory population')
+        traj.f_add_result('runs.$.inh_tuning_vectors', multi_proc_result[idx][5],
+                          comment='The tuning vectors of the inhibitory population')
+
+        traj.f_restore_default()
+
+    traj.f_store()
+    print('Analysis done')
+
+if __name__ == "__main__":
+    main()

scripts/spatial_network/orientation_map/orientation_map_full_figure_runner.py

@@ -0,0 +1,7 @@
+import scripts.spatial_network.perlin.run_entropy_maximisation_perlin_map as runner
+import scripts.spatial_network.perlin.analyse_entropy_maximisation_perlin_map as analyser
+import scripts.spatial_network.perlin.paper_figures_spatial_head_direction_network_perlin_map as plotter
+
+runner.main()
+analyser.main()
+# plotter.main()

scripts/spatial_network/orientation_map/run_orientation_map.py

@@ -0,0 +1,269 @@
+import json
+import os
+
+import numpy as np
+from brian2.units import *
+from pypet import Environment, cartesian_product, Trajectory
+from pypet.brian2.parameter import Brian2MonitorResult
+
+from scripts.interneuron_placement 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_maps.orientation_maps.orientation_map import OrientationMap
+from scripts.spatial_maps.orientation_maps.orientation_map_generator_pypet import TRAJ_NAME_ORIENTATION_MAPS
+from scripts.spatial_maps.uniform_perlin_map import UniformPerlinMap
+from scripts.spatial_network.head_direction_index_over_noise_scale 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/"
+
+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.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
+
+    number_of_excitatory_neurons_per_row = int(np.sqrt(N_E))
+
+    map = OrientationMap(number_of_excitatory_neurons_per_row + 1, number_of_excitatory_neurons_per_row + 1,
+                         corr_len, sheet_size, sheet_size, seed)
+    map.angle_grid = map_angle_grid
+
+    return map.tuning
+
+
+def spatial_network_with_entropy_maximisation(traj):
+    sheet_size = traj.map.sheet_size
+
+    N_E = traj.network.N_E
+    N_I = traj.network.N_I
+
+    orientation_map = get_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)
+
+    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
+
+    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)
+
+    ie_connections = get_excitatory_neurons_in_inhibitory_axonal_clouds(ex_positions, inhibitory_axonal_clouds)
+
+    inhibitory_synapse_strength = traj.synapse.inhibitory * nS
+    excitatory_synapse_strength = traj.synapse.excitatory * mV
+
+    if inhibitory_synapse_strength != 0.0 * nS and excitatory_synapse_strength != 0.0 * mV \
+            and inhibitory_axon_long_axis == inhibitory_axon_short_axis:
+        traj.f_add_derived_parameter("morphology.morph_label", CIRCULAR,
+                                     comment="Interneuron morphology of this run is circular")
+    elif inhibitory_synapse_strength != 0.0 * nS and excitatory_synapse_strength != 0.0 * mV:
+        traj.f_add_derived_parameter("morphology.morph_label", POLARIZED,
+                                     comment="Interneuron morphology of this run is ellipsoid")
+    else:
+        traj.f_add_derived_parameter("morphology.morph_label", NO_SYNAPSES,
+                                     comment="There are no interneurons")
+
+    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):
+        # 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,
+                            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)
+        input_to_excitatory_population = create_head_direction_input(traj.input.baseline * nA, ex_tunings,
+                                                                     sharpness,
+                                                                     traj.input.amplitude * nA, dir)
+
+        excitatory_neurons = net["excitatory_neurons"]
+        excitatory_neurons.I = input_to_excitatory_population
+
+        inhibitory_neurons = net["interneurons"]
+        inhibitory_neurons.u_ext = traj.inh_input.baseline * mV
+        inhibitory_neurons.tau = traj.interneuron.tau * ms
+
+        net.run(traj.simulation.duration * ms)
+
+        direction_id = 'dir{:d}'.format(idx)
+        traj.f_add_result(Brian2MonitorResult, '{:s}.spikes.e'.format(direction_id), net["excitatory_spike_monitor"],
+                          comment='The spiketimes of the excitatory population')
+        traj.f_add_result(Brian2MonitorResult, '{:s}.spikes.i'.format(direction_id), net["inhibitory_spike_monitor"],
+                          comment='The spiketimes of the inhibitory population')
+    traj.f_add_result('ex_positions', np.array(ex_positions),
+                      comment='The positions of the excitatory neurons on the sheet')
+    traj.f_add_result('ex_tunings', np.array(ex_tunings),
+                      comment='The input tunings of the excitatory neurons')
+
+    ie_connections_save_array = np.zeros((N_I, N_E))
+    for i_idx, ie_conn in enumerate(ie_connections):
+        for e_idx in ie_conn:
+            ie_connections_save_array[i_idx, e_idx] = 1
+    traj.f_add_result('ie_adjacency', ie_connections_save_array,
+                      comment='Recurrent connection adjacency matrix')
+
+    axon_cloud_save_list = [[p.x, p.y, p.phi] for p in inhibitory_axonal_clouds]
+    axon_cloud_save_array = np.array(axon_cloud_save_list)
+    traj.f_add_result('inhibitory_axonal_cloud_array', axon_cloud_save_array,
+                      comment='The inhibitory axonal clouds')
+
+    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():
+    env = Environment(trajectory=TRAJ_NAME,
+                      comment="Compare the head direction tuning for circular and ellipsoid interneuron morphology, "
+                              "when tuning orientations to maximise entropy of connected excitatory tunings.",
+                      multiproc=True, filename=DATA_FOLDER, ncores=30, 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("network")
+    traj.f_add_parameter("network.N_E", 3600, comment="Number of excitatory neurons")
+    traj.f_add_parameter("network.N_I", 400, comment="Number of inhibitory neurons")
+
+    traj.f_add_parameter_group("interneuron")
+    traj.f_add_parameter("interneuron.tau", 7., comment="Interneuron timescale in ms")
+
+    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.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")
+
+    traj.f_add_parameter_group("inh_input")
+    traj.f_add_parameter("inh_input.baseline", -50., comment="Head direction input baseline")
+    traj.f_add_parameter("inh_input.amplitude", 0., comment="Head direction input amplitude")
+
+    traj.f_add_parameter_group("morphology")
+    traj.f_add_parameter("morphology.long_axis", 100.0, comment="Long axis of axon ellipsoid")
+    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 = np.linspace(1.0, 800.0, 12, endpoint=True).tolist()
+    # correlation_length_range = [200.0]
+    seed_range = range(10)
+    # 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,
+        "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[
+                                                      "morphology.long_axis"][0] * ellipsoid_parameter_exploration[
+                                                      "morphology.short_axis"][0]))
+
+    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"],
+        "synapse.inhibitory": ellipsoid_parameter_exploration["synapse.inhibitory"],
+        "synapse.excitatory": ellipsoid_parameter_exploration["synapse.excitatory"]
+    }
+
+    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"],
+        "synapse.inhibitory": [0.],
+        "synapse.excitatory": [0.]
+    }
+
+    expanded_dicts = [cartesian_product(dict) for dict in [ellipsoid_parameter_exploration,
+                                                           circle_parameter_exploration,
+                                                           no_conn_parameter_exploration]]
+    final_dict = {}
+
+    for key in expanded_dicts[0].keys():
+        list_of_parameter_lists = [dict[key] for dict in expanded_dicts]
+        final_dict[key] = sum(list_of_parameter_lists, [])
+
+    traj.f_explore(final_dict)
+    env.run(spatial_network_with_entropy_maximisation)
+
+    env.disable_logging()
+
+if __name__ == "__main__":
+    main()

scripts/spatial_network/perlin/paper_figures_spatial_head_direction_network_perlin_map.py

@@ -1006,7 +1006,7 @@ def plot_hdi_histogram_combined_and_overlayed(traj, plot_run_names, ex_polar_plo
         ax.axvline(np.mean(exc_hdis[i]), color='r')
 
         ax.axvline(np.mean(inh_hdis[i]), color='b')
-        ax.axvline(no_conn_hdi, color='dimgrey', linestyle='--', linewidth=2.5)
+        ax.axvline(no_conn_hdi, color='dimgrey', linestyle='--', linewidth=1.5)
 
         ax.set_ylabel(labels[i], rotation='vertical')
         # plt.axis('on')
@@ -1032,8 +1032,10 @@ def plot_hdi_histogram_combined_and_overlayed(traj, plot_run_names, ex_polar_plo
                     textcoords="offset points",
                     va="top", ha="left", color="r")
 
-        i_ha = "left" if i == 1 else "right"
-        i_offset = 2 if i == 1 else -2
+        # i_ha = "left" if i == 1 else "right"
+        # i_offset = 2 if i == 1 else -2
+        i_ha = "right"
+        i_offset = -1
         ax.annotate("{:.2f}".format(np.mean(inh_hdis[i])), xy=(np.mean(inh_hdis[i]), max_density),
                     xytext=(i_offset, 0), xycoords="data",
                     textcoords="offset points",
@@ -1221,7 +1223,7 @@ def plot_exc_and_inh_hdi_over_simplex_grid_scale(traj, plot_run_names):
     # plt.legend()
 
     if save_figs:
-        plt.savefig(FIGURE_SAVE_PATH + 'F_hdi_over_corr_len_scaled.png')
+        plt.savefig(FIGURE_SAVE_PATH + 'F_hdi_over_grid_scale.png')
         plt.close(fig)
 
 def plot_exc_and_inh_hdi_over_fit_corr_len(traj, plot_run_names):
@@ -1327,7 +1329,7 @@ def plot_exc_and_inh_hdi_over_fit_corr_len(traj, plot_run_names):
     # plt.legend()
 
     if save_figs:
-        plt.savefig(FIGURE_SAVE_PATH + 'hdi_over_fit_corr_len.png')
+        plt.savefig(FIGURE_SAVE_PATH + 'F_hdi_over_corr_len_scaled.png')
         plt.close(fig)
 
 def plot_in_degree_map(traj, plot_run_names):
@@ -1547,7 +1549,7 @@ if __name__ == "__main__":
     map_seed = 1
     exemplary_head_direction = 0
 
-    # corr_len_fit_dict = correlation_length_fit_dict(traj, load=True)
+    # corr_len_fit_dict = correlation_length_fit_dict(traj, load=False)
     # plt.plot(corr_len_fit_dict.keys(),corr_len_fit_dict.values())
     # plt.show()
     # abbrechen
@@ -1577,8 +1579,8 @@ if __name__ == "__main__":
     directions = get_input_head_directions(traj)
     direction_idx = np.argmin(np.abs(np.array(directions) - np.deg2rad(exemplary_head_direction)))
 
-    selected_neuron_excitatory = 1716
-    selected_inhibitory_neuron = 144
+    selected_neuron_excitatory = 1052
+    selected_inhibitory_neuron = 28
 
     print("## Figure specification")
     print("\tpanel size: {:.2f} cm".format(panel_size * cm_per_inch))

scripts/spatial_network/placement_jitter/run_entropy_maximisation_orientation_map_placement_jitter.py

@@ -6,8 +6,8 @@ from pypet.brian2.parameter import Brian2MonitorResult
 from scripts.interneuron_placement 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_maps.orientation_map import OrientationMap
-from scripts.spatial_maps.orientation_map_generator_pypet import TRAJ_NAME_ORIENTATION_MAPS
+from scripts.spatial_maps.orientation_maps.orientation_map import OrientationMap
+from scripts.spatial_maps.orientation_maps.orientation_map_generator_pypet import TRAJ_NAME_ORIENTATION_MAPS
 from scripts.spatial_network.head_direction_index_over_noise_scale 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, \

scripts/spatial_network/run_entropy_maximisation_orientation_map.py

@@ -6,8 +6,8 @@ from pypet.brian2.parameter import Brian2MonitorResult
 from scripts.interneuron_placement 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_maps.orientation_map import OrientationMap
-from scripts.spatial_maps.orientation_map_generator_pypet import TRAJ_NAME_ORIENTATION_MAPS
+from scripts.spatial_maps.orientation_maps.orientation_map import OrientationMap
+from scripts.spatial_maps.orientation_maps.orientation_map_generator_pypet import TRAJ_NAME_ORIENTATION_MAPS
 from scripts.spatial_network.head_direction_index_over_noise_scale 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, \

scripts/spatial_network/run_synaptic_strength_scan_orientation_map.py

@@ -6,8 +6,8 @@ from pypet.brian2.parameter import Brian2MonitorResult
 from scripts.interneuron_placement 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_maps.orientation_map import OrientationMap
-from scripts.spatial_maps.orientation_map_generator_pypet import TRAJ_NAME_ORIENTATION_MAPS
+from scripts.spatial_maps.orientation_maps.orientation_map import OrientationMap
+from scripts.spatial_maps.orientation_maps.orientation_map_generator_pypet import TRAJ_NAME_ORIENTATION_MAPS
 from scripts.spatial_network.head_direction_index_over_noise_scale 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, \

scripts/spatial_network/run_synaptic_strength_scan_orientation_map_small_scale.py

@@ -6,8 +6,8 @@ from pypet.brian2.parameter import Brian2MonitorResult
 from scripts.interneuron_placement 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_maps.orientation_map import OrientationMap
-from scripts.spatial_maps.orientation_map_generator_pypet import TRAJ_NAME_ORIENTATION_MAPS
+from scripts.spatial_maps.orientation_maps.orientation_map import OrientationMap
+from scripts.spatial_maps.orientation_maps.orientation_map_generator_pypet import TRAJ_NAME_ORIENTATION_MAPS
 from scripts.spatial_network.head_direction_index_over_noise_scale 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, \

scripts/spatial_network/sharpening_over_noise_scale.py

@@ -2,7 +2,7 @@ import matplotlib.pyplot as plt
 import noise
 import numpy as np
 from brian2.units import *
-from scripts.spatial_maps.orientation_map import OrientationMap
+from scripts.spatial_maps.orientation_maps.orientation_map import OrientationMap
 from scripts.interneuron_placement import create_interneuron_sheet_entropy_max_orientation
 from tqdm import tqdm
 import scripts.models as modellib

scripts/spatial_network/sharpening_over_noise_scale_multiprocessing.py

@@ -2,7 +2,7 @@ import matplotlib.pyplot as plt
 import noise
 import numpy as np
 from brian2.units import *
-from scripts.spatial_maps.orientation_map import OrientationMap
+from scripts.spatial_maps.orientation_maps.orientation_map import OrientationMap
 from scripts.interneuron_placement import create_interneuron_sheet_entropy_max_orientation
 import scripts.models as modellib
 from scripts.ring_network.head_direction import get_head_direction_input, \

scripts/spatial_network/spatial_head_dir_network.py

@@ -2,7 +2,7 @@ import matplotlib.pyplot as plt
 import noise
 import numpy as np
 from brian2.units import *
-from scripts.spatial_maps.orientation_map import OrientationMap
+from scripts.spatial_maps.orientation_maps.orientation_map import OrientationMap
 from scripts.interneuron_placement import create_interneuron_sheet_entropy_max_orientation, \
     create_grid_of_excitatory_neurons, plot_neural_sheet, create_interneuron_sheet_by_repulsive_force, \
     get_excitatory_neurons_in_inhibitory_axonal_clouds, get_position_mesh

scripts/spatial_packing_entropy/entropy_over_noise_scale.py

@@ -2,7 +2,7 @@ import matplotlib.pyplot as plt
 import noise
 import numpy as np
 from brian2.units import *
-from scripts.spatial_maps.orientation_map import OrientationMap
+from scripts.spatial_maps.orientation_maps.orientation_map import OrientationMap
 from scripts.interneuron_placement import create_interneuron_sheet_entropy_max_orientation
 from tqdm import tqdm
 

scripts/spatial_packing_entropy/entropy_over_noise_scale_multiprocessing.py

@@ -2,7 +2,7 @@ import matplotlib.pyplot as plt
 import noise
 import numpy as np
 from brian2.units import *
-from scripts.spatial_maps.orientation_map import OrientationMap
+from scripts.spatial_maps.orientation_maps.orientation_map import OrientationMap
 from scripts.interneuron_placement import create_interneuron_sheet_entropy_max_orientation
 
 from scripts.interneuron_placement import create_grid_of_excitatory_neurons, \