Firing rate still incorrect and thus HDI turns out to be NAN.

environment.yml

@@ -186,5 +186,5 @@ dependencies:
     - numexpr==2.7.1
     - pypet==0.4.3
     - tables==3.6.1
-prefix: /home/pfeiffer/Applications/miniconda2/envs/interneuron_polarity
+
 

scripts/spatial_maps/orientation_map_generator_pypet.py

@@ -0,0 +1,109 @@
+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=3)
+    traj = env.trajectory
+
+    traj.f_add_parameter('corr_len', 200.0, comment='Correlation length')
+    traj.f_add_parameter('seed', 1, comment='Seed for random')
+
+    traj.f_explore(cartesian_product({'corr_len': [200.0], 'seed': [1]}))
+
+    # 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_network/analyse_entropy_maximisation_orientation_map.py

@@ -1,6 +1,7 @@
 import numpy as np
 from pypet import Trajectory
 from pypet.brian2 import Brian2MonitorResult
+from tqdm import tqdm
 
 from scripts.spatial_network.head_direction_index_over_noise_scale import calculate_rates
 from scripts.spatial_network.run_entropy_maximisation_orientation_map import DATA_FOLDER, TRAJ_NAME
@@ -13,6 +14,7 @@ def get_spike_train_dictionary(number_of_neurons, spike_times, neuron_indices):
 
     for neuron_idx, t in zip(neuron_indices, spike_times):
         spike_train_dict[neuron_idx].append(t)
+    print(spike_train_dict)
     return spike_train_dict
 
 
@@ -34,7 +36,23 @@ def get_firing_rate_per_cell_and_direction(traj, run_name):
 
 
 def get_head_direction_indices(directions, firing_rate_dict):
-    return np.zeros(len(firing_rate_dict.keys()))
+    n_exc_neurons = len(firing_rate_dict)
+    n_directions = len(directions)
+
+    tuning_vectors = np.ndarray((n_exc_neurons,n_directions,2))
+
+    for ex_id, ex_rates in firing_rate_dict.items():
+        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]
+        tuning_vectors[ex_id, :] /= rate_sum
+
+    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
 
 
 def get_runs_with_circular_morphology(traj):
@@ -46,7 +64,7 @@ if __name__ == "__main__":
     traj = Trajectory(TRAJ_NAME, add_time=False, dynamic_imports=Brian2MonitorResult)
     NO_LOADING = 0
     FULL_LOAD = 2
-
+    # TODO: Why no loading of results? Because of dynamic loading?
     traj.f_load(filename=DATA_FOLDER + TRAJ_NAME + ".hdf5", load_parameters=FULL_LOAD, load_results=NO_LOADING)
 
     traj.v_auto_load = True
@@ -54,13 +72,14 @@ if __name__ == "__main__":
     correlation_lengths = traj.f_get('correlation_length').f_get_range()
     long_axis = traj.f_get('long_axis').f_get_range()
     short_axis = traj.f_get('short_axis').f_get_range()
-
+    # TODO: Again, maybe directions as parameters
     directions = np.linspace(-np.pi, np.pi, traj.input.number_of_directions)
 
     circular_indices = list(get_runs_with_circular_morphology(traj))
 
     for idx, run_name in enumerate(traj.f_get_run_names()):
         firing_rate_dict = get_firing_rate_per_cell_and_direction(traj, run_name)
+        print(firing_rate_dict)
         head_direction_indices = get_head_direction_indices(directions, firing_rate_dict)
         traj.f_set_crun(run_name)
         print("Circle" if idx in circular_indices else "Ellipsoid" )

scripts/spatial_network/head_direction_index_over_noise_scale.py

@@ -240,7 +240,7 @@ def get_synaptic_weights(N_E, N_I, ie_connections, excitatory_synapse_strength,
     return ex_in_weights, in_ex_weights
 
 
-if __name__ == "main":
+if __name__ == "__main__":
     corr_len_range = range(0, 451, 15)
     corr_len_range_len = len(corr_len_range)
     tuning_range_len = 12
@@ -261,7 +261,7 @@ if __name__ == "main":
         np.save('../../simulations/2020_02_27_head_direction_index_over_noise_scale/data.npy', np.array(data))
     else:
         data = np.load('../../simulations/2020_02_27_head_direction_index_over_noise_scale/data.npy', allow_pickle=True)
-
+    print('Calculating HDI')
     no_conn_trial_hdi_array = np.array((corr_len_range_len, seed_range_len, tuning_range_len))
     ellipse_trial_hdi_array = np.array((corr_len_range_len, seed_range_len, tuning_range_len))
     circle_trial_hdi_array = np.array((corr_len_range_len, seed_range_len, tuning_range_len))

scripts/spatial_network/run_entropy_maximisation_orientation_map.py

@@ -1,11 +1,13 @@
 import numpy as np
 from brian2.units import *
-from pypet import Environment, cartesian_product
+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_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, \
@@ -17,10 +19,33 @@ LOG_FOLDER = "../../logs/"
 TRAJ_NAME = "entropy_maximisation_over_different_input_correlation_lengths"
 
 
-def get_orientation_map(correlation_length, seed, sheet_size):
-    return lambda x, y: 0.0
+def get_orientation_map(correlation_length, seed, sheet_size, N_E):
+    traj = Trajectory(filename=DATA_FOLDER + TRAJ_NAME_ORIENTATION_MAPS + ".hdf5")
 
+    traj.f_load(index=-1, load_parameters=2, load_results=2)
 
+    corr_len = correlation_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
+
+    tuning_map = lambda x, y: map.orientation_map(x, y)
+    return tuning_map
+
+# TODO: Rename to avoid namespace conflicts
 def experiment(traj):
     sheet_size = traj.orientation_map.sheet_size
 
@@ -28,7 +53,7 @@ def experiment(traj):
     N_I = traj.network.N_I
 
     orientation_map = get_orientation_map(traj.orientation_map.correlation_length, traj.orientation_map.seed,
-                                          sheet_size)
+                                          sheet_size, N_E)
     ex_positions, ex_tunings = create_grid_of_excitatory_neurons(sheet_size,
                                                                  sheet_size,
                                                                  int(np.sqrt(N_E)), orientation_map)
@@ -36,6 +61,7 @@ def experiment(traj):
     inhibitory_axon_long_axis = traj.morphology.long_axis
     inhibitory_axon_short_axis = traj.morphology.short_axis
 
+    # TODO: Remove one step for circles as it's unnecessary
     entropy_maximisation_steps = traj.simulation.entropy_maximisation.steps if inhibitory_axon_long_axis != \
                                                                                inhibitory_axon_short_axis else 1
 
@@ -62,8 +88,10 @@ def experiment(traj):
 
     sharpness = 1.0 / (traj.input.width) ** 2
 
+    # TODO: Parameterize input direction?
     directions = np.linspace(-np.pi, np.pi, traj.input.number_of_directions)
 
+    # TODO: Maybe smarter to use pypet explore instead of loop
     for idx, dir in enumerate(directions):
         net.restore()
         input_to_excitatory_population = create_head_direction_input(traj.input.baseline * nA, ex_tunings,
@@ -94,7 +122,7 @@ if __name__ == "__main__":
 
     traj.f_add_parameter_group("orientation_map")
 
-    traj.f_add_parameter("orientation_map.correlation_length", 100, comment="Correlation length of orientations in um")
+    traj.f_add_parameter("orientation_map.correlation_length", 200.0, comment="Correlation length of orientations in um")
     traj.f_add_parameter("orientation_map.seed", 1, comment="Random seed for map generation.")
     traj.f_add_parameter("orientation_map.sheet_size", 450, comment="Sheet size in um")
 
@@ -124,10 +152,13 @@ if __name__ == "__main__":
     traj.f_add_parameter("simulation.dt", 0.1, comment="Network simulation time step in ms")
     traj.f_add_parameter("simulation.duration", 10, comment="Network simulation duration in ms")
 
+    # TODO: The Orientation Maps need to be precalculated for these values
+    # TODO: Random seeds need to be added to exploration
     ellipsoid_parameter_exploration = {
         "morphology.long_axis": [200.0],
         "morphology.short_axis": [50.0],
-        "orientation_map.correlation_length": np.arange(1, 40, 50).tolist()
+        "orientation_map.correlation_length": [200.0]
+        # "orientation_map.correlation_length": np.arange(0.0, 200.0, 50).tolist()
     }
 
     corresponding_circular_radius = float(np.sqrt(ellipsoid_parameter_exploration[