generation of the connectivity matrix

src/interneuron_polarity/actions/create_connectivity.py

@@ -1,12 +1,30 @@
+import numpy
+from interneuron_polarity.io.io import IO
+from interneuron_polarity.model.morphology.shapes import Ellipsoid
+from tqdm import tqdm
+
+
 def create_connectivity(simulation_directory):
-    # morphology = load_morphology(simulation_directory)
-    #
-    # %actual additional step
-    # excitatory_neuron_positions
-    # inhibitory_neuron_positions
-    # inhibitory_axonal_shapes
-    #
-    #
-    # in_ex_connectivity
-    #
-    pass
+    my_io = IO(simulation_directory)
+    morphology = my_io.load_morphology_dict()
+    ex_positions, in_positions, in_orientations = my_io.load_positions_and_orientations()
+
+    n_ex = ex_positions.shape[1]
+    n_in = in_positions.shape[1]
+
+    long_axis = morphology["inhibitory_neurons"]["long_axis"]
+    short_axis = morphology["inhibitory_neurons"]["long_axis"]
+    inhibitory_axon_shape = Ellipsoid(long_axis, short_axis)
+
+    in_ex_connectivity = numpy.zeros((n_in, n_ex))
+    t = tqdm(total=n_in, desc="Identify morphological connection of interneurons to excitatory neurons")
+    for inhibitory_index, inhibitory_specification in enumerate(zip(in_positions.T, in_orientations.T)):
+        in_pos, in_orientation = inhibitory_specification
+        contained_excitatory_neurons = inhibitory_axon_shape.contains(ex_positions, in_pos,
+                                                                               in_orientation)
+
+        in_ex_connectivity[inhibitory_index, contained_excitatory_neurons] = 1
+        t.update(1)
+    t.close()
+
+    my_io.save_in_ex_connectivity(in_ex_connectivity)

src/interneuron_polarity/application.py

@@ -105,7 +105,7 @@ def main(args):
 
     if args.action == "space":
         create_spatial_arrangement(path_to_simulation_dir)
-    if args.action == "space":
+    if args.action == "conn":
         create_connectivity(path_to_simulation_dir)
 
 

src/interneuron_polarity/io/io.py

@@ -26,6 +26,23 @@ class IO:
         save_array(self.output_folder, "in_positions", inhibitory_somata_positions)
         save_array(self.output_folder, "in_orientations", inhibitory_axonal_orientations)
 
+    def load_positions_and_orientations(self):
+        excitatory_somata_positions = load_array(self.output_folder, "ex_positions")
+        inhibitory_somata_positions = load_array(self.output_folder, "in_positions")
+        inhibitory_axonal_orientations = load_array(self.output_folder, "in_orientations")
+        return excitatory_somata_positions, inhibitory_somata_positions, inhibitory_axonal_orientations
+
+    def save_in_ex_connectivity(self, in_ex_connectivity):
+        save_array(self.output_folder, "in_ex_connectivity", in_ex_connectivity)
+
 
 def save_array(dir, filename, array):
-    np.save(os.path.join(dir, "{:s}.npy".format(filename)), array)
+    full_path = os.path.join(dir, "{:s}.npy".format(filename))
+    print("Save {:s}".format(full_path))
+    np.save(full_path, array)
+
+
+def load_array(dir, filename):
+    full_path = os.path.join(dir, "{:s}.npy".format(filename))
+    print("Load {:s}".format(full_path))
+    return np.load(full_path)

conmorph/shapes.py

@@ -1,19 +1,22 @@
 import numpy as np
 
-from conmorph.volume import NeuralVolume
+from interneuron_polarity.model.morphology.volume import NeuralVolume
 
 
 def uniform(distance_to_foci):
     return np.ones(distance_to_foci.shape)
 
+
 def get_gaussian(length_scale):
     def gaussian(distance_to_foci):
-        return np.exp(-distance_to_foci**2/length_scale)
+        return np.exp(-distance_to_foci ** 2 / length_scale)
+
     return gaussian
 
+
 class Ellipsoid:
     def __init__(self, r_long, r_small):
-        self.center = np.array([0,0,0])
+        self.center = np.array([0, 0, 0])
         self.r_long = r_long
         self.r_small = r_small
 
@@ -21,18 +24,20 @@ class Ellipsoid:
         self.distance_to_foci = 2 * r_long
 
     def get_containing_volume(self, center, long_axis, voxel_size):
-        first_direction = tuple(center + long_axis*self.r_long) #is this big enough?, not quite add orthogonal direction with r_small
-        second_direction = tuple(center - long_axis*self.r_long)
-        x_min, y_min, z_min = tuple([min(f,s) for f,s in zip(first_direction, second_direction)])
-        x_max, y_max, z_max = tuple([max(f,s) for f,s in zip(first_direction, second_direction)])
+        first_direction = tuple(
+            center + long_axis * self.r_long)  # is this big enough?, not quite add orthogonal direction with r_small
+        second_direction = tuple(center - long_axis * self.r_long)
+        x_min, y_min, z_min = tuple([min(f, s) for f, s in zip(first_direction, second_direction)])
+        x_max, y_max, z_max = tuple([max(f, s) for f, s in zip(first_direction, second_direction)])
         return NeuralVolume(x_max, y_max, z_max, voxel_size, x_min, y_min, z_min)
 
     def get_density(self, center, long_axis, volume, density=uniform):
         coordinates = volume.get_coordinates()
         distance_to_foci = self.get_distance_to_foci(coordinates, center, long_axis)
 
-        density = np.where(distance_to_foci<self.distance_to_foci, density(distance_to_foci), np.zeros(distance_to_foci.shape))
-        return density.reshape(volume.mesh_shape)/np.linalg.norm(density)
+        density = np.where(distance_to_foci < self.distance_to_foci, density(distance_to_foci),
+                           np.zeros(distance_to_foci.shape))
+        return density.reshape(volume.mesh_shape) / np.linalg.norm(density)
 
     def get_distance_to_foci(self, coordinates, center, long_axis):
         focus_1 = center + long_axis * self.distance_between_foci / 2.0
@@ -48,14 +53,16 @@ class Ellipsoid:
     def contains(self, coordinates, center, long_axis):
         return np.where(self.get_distance_to_foci(coordinates, center, long_axis) < self.distance_to_foci)[0]
 
+
 def get_overlap(density_1, density_2):
     return np.sum(density_1 * density_2)
 
-def distance(x,y):
-    return np.linalg.norm(x-y)
+
+def distance(x, y):
+    return np.linalg.norm(x - y)
 
 
 class NeuralProcess:
     def __init__(self, form, orientation):
         self.form = form
-        self.orientation = orientation
+        self.orientation = orientation