Fixed seeds

examples/notebooks/README.md

@@ -10,6 +10,7 @@ Here is a collection of Jupyter Notebooks, some of the workflows are split over
 * [population_unit_network](population_unit_network.ipynb) how to define population units.
 * [example_of_density_function](example_of_density_function.ipynb) how to specify density variations using a function of (x,y,z) in a volume.
 * [example_of_neuron_rotations](example_of_neuron_rotations.ipynb) shows how to rotate neurons based on position.
+* [bend_morphologies](bend_morphologies.ipynb) shows how to make the neurons bend the axons and dendrites at the edge of the mesh, to keep them constrained to the volume.
 * [connect_structures_example](connect_structures_example.ipynb) shows how to create neuron projections between volumes when no-axon data is available ([parallel version](connect_structures_example_parallel.ipynb)). There is also an [alternative version](connect_structures_example_projection_detection.ipynb) that places axon voxels randomly within the projection zone and then applies touch detection.
 * [virtual_neurons](VirtualNeurons/VirtualNeurons.ipynb) shows how to only simulate the core of a volume of neurons, by turning the outer neurons to virtual neurons that spike at predetermined times. This is useful to avoid edge effects in your simulations.
 

examples/notebooks/bend_morphologies.ipynb

@@ -2,7 +2,7 @@
  "cells": [
   {
    "cell_type": "markdown",
-   "id": "86b06c63",
+   "id": "858c7f9d",
    "metadata": {},
    "source": [
     "# Bend Morphologies to stay inside Mesh\n",
@@ -13,7 +13,7 @@
   {
    "cell_type": "code",
    "execution_count": 1,
-   "id": "240433fe",
+   "id": "c18817cc",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -29,7 +29,7 @@
   {
    "cell_type": "code",
    "execution_count": 2,
-   "id": "389368d3",
+   "id": "8f1a7c56",
    "metadata": {},
    "outputs": [
     {
@@ -60,7 +60,7 @@
   {
    "cell_type": "code",
    "execution_count": 3,
-   "id": "6c2be32a",
+   "id": "89a08a73",
    "metadata": {},
    "outputs": [
     {
@@ -119,7 +119,7 @@
   {
    "cell_type": "code",
    "execution_count": 4,
-   "id": "d4f90195",
+   "id": "2e073869",
    "metadata": {},
    "outputs": [
     {
@@ -153,7 +153,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "b2c629fc",
+   "id": "f854ad33",
    "metadata": {},
    "outputs": [],
    "source": []

snudda/place/bend_morphologies.py

@@ -24,10 +24,17 @@ def check_if_inside(self, morphology: NeuronMorphologyExtended):
 
         return inside_flag
 
-    def bend_morphology(self, morphology: NeuronMorphologyExtended, k=30e-6, n_random=5):
+    def bend_morphology(self, morphology: NeuronMorphologyExtended, k=30e-6, n_random=5, random_seed=None):
 
         # k -- how early will the neuron start bending when it approaches the border
 
+        # print(f"random_seed = {random_seed}")
+
+        if random_seed is not None:
+            rng = np.random.default_rng(random_seed)
+        else:
+            rng = self.rng
+
         # Check distance to border for all points, negative distance means inside
         all_original_dist = self.region_mesh.distance_to_border(morphology.geometry[:, :3])
         if (all_original_dist < 0).all():
@@ -89,30 +96,24 @@ def bend_morphology(self, morphology: NeuronMorphologyExtended, k=30e-6, n_rando
                 P_move = 1 / (1 + np.exp(-dist/k))
 
                 # Cache the random numbers for segments in the section...
-                if dist > parent_dist and self.rng.uniform() < P_move:
-                #  if self.rng.uniform() < P_move:
+                if dist > parent_dist and rng.uniform() < P_move:
 
                     morphology_changed = True
                     parent_moved = True
 
                     # We need to randomize new rotation matrix
                     # https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.transform.Rotation.html
-                    angles = self.rng.uniform(size=(n_random, 3), low=-0.1, high=0.1)  # Angles in radians
+                    angles = rng.uniform(size=(n_random, 3), low=-0.1, high=0.1)  # Angles in radians
                     avoidance_rotations = Rotation.from_euler(seq="XYZ", angles=angles)
 
                     for idx2, av_rot in enumerate(avoidance_rotations):
                         candidate_pos[idx2, :] = parent_point + length * (av_rot*rotation).apply(vectors=parent_dir)
 
                     candidate_dist = self.region_mesh.distance_to_border(points=candidate_pos)
 
-                    if False:
-                        P_candidate = np.divide(1, 1 + np.exp(-candidate_dist/k))
-                        P_candidate = P_candidate / np.sum(P_candidate)
-                        picked_idx = self.rng.choice(n_random, p=P_candidate)
-                    else:
-                        # We want the smallest (or most negative) distance
-                        picked_idx = np.argsort(candidate_dist)[0]
-                        dist = candidate_dist[picked_idx]
+                    # We want the smallest (or most negative) distance
+                    picked_idx = np.argsort(candidate_dist)[0]
+                    dist = candidate_dist[picked_idx]
 
                     new_rot = avoidance_rotations[picked_idx]*rotation
                     new_rot_rep.append((new_rot, length))
@@ -341,11 +342,11 @@ def write_swc(self, morphology: MorphologyData, output_file, comment=None):
 
         print(f"Wrote {output_file}")
 
-    def edge_avoiding_morphology(self, swc_file, new_file, original_position, original_rotation):
+    def edge_avoiding_morphology(self, swc_file, new_file, original_position, original_rotation, random_seed=None):
 
         md = MorphologyData(swc_file=swc_file)
         md.place(rotation=original_rotation, position=original_position)
-        rot_rep, morphology_changed = self.bend_morphology(md)
+        rot_rep, morphology_changed = self.bend_morphology(md, random_seed=random_seed)
 
         if morphology_changed:
             new_coord = self.apply_rotation(md, rot_rep)

snudda/place/place.py

@@ -567,25 +567,26 @@ def parse_config(self, config_file=None, resort_neurons=True):
 
     def avoid_edges_parallel(self):
 
+        ss = np.random.SeedSequence(self.random_seed + 100)
+        neuron_random_seed = ss.generate_state(len(self.neurons))
+
         if self.d_view is None:
-            return self.avoid_edges()
+            # Make sure we use the same random seeds if we run in serial, as would have been used in parallel
+            return self.avoid_edges(neuron_random_seeds=neuron_random_seed)
 
         # Make random permutation of neurons, to spread out the edge neurons
         unsorted_neuron_id = self.random_generator.permutation(len(self.neurons))
-        ss = self.random_generator.SeedSequence()
-
-        worker_random_seed = ss.generate_state(len(self.d_view))
 
         self.d_view.scatter("worker_neuron_id", unsorted_neuron_id, block=True)
-        self.d_view.scatter("worker_random_seed", worker_random_seed, block=True)
         self.dview.push({"config": self.config,
-                         "neurons": self.neurons})
+                         "neurons": self.neurons,
+                         "neuron_random_seeds": neuron_random_seed})
 
         cmd_str = f"sp = SnuddaPlace(config_file={self.config_file},network_path={self.network_path},snudda_data={self.snudda_data}, random_seed=worker_random_seed[0])"
         self.d_view.execute(cmd_str)
         cmd_str2 = f"sp.config = config; sp.neurons = neurons"
         self.d_view.execute(cmd_str2)
-        cmd_str3 = f"modified_neurons = sp.avoid_edges(neuron_id=unsorted_neuron_id)"
+        cmd_str3 = f"modified_neurons = sp.avoid_edges(neuron_id=unsorted_neuron_id, random_seed=neuron_random_seeds[unsorted_neuron_id])"
         self.d_view.execute(cmd_str3)
         modified_neurons = self.d_view.gather("modified_neurons", block=True)
 
@@ -594,7 +595,7 @@ def avoid_edges_parallel(self):
             self.neurons[neuron_id].swc_filename = new_morphology
             self.neurons[neuron_id].rotation = np.eye(3)
 
-    def avoid_edges(self, neuron_id=None):
+    def avoid_edges(self, neuron_id=None, neuron_random_seeds=None):
 
         from snudda.place.bend_morphologies import BendMorphologies
 
@@ -609,9 +610,16 @@ def avoid_edges(self, neuron_id=None):
         else:
             neurons = self.neurons[neuron_id]
 
+        if neuron_random_seeds is None:
+            neuron_random_seeds = [None for n in neurons]
+        else:
+            neuron_random_seeds = neuron_random_seeds[neuron_id]
+
         modified_morphologies = []
 
-        for neuron in neurons:
+        for neuron, random_seed in zip(neurons, neuron_random_seeds.flatten()):
+
+            # print(f"neuron.name = {neuron.name}, random_seed = {random_seed}")
             config = self.config["Neurons"][neuron.name]
 
             if "stayInsideMesh" in config and config["stayInsideMesh"]:
@@ -626,7 +634,8 @@ def avoid_edges(self, neuron_id=None):
                 new_morphology = bend_morph[volume_id].edge_avoiding_morphology(swc_file=neuron.swc_filename,
                                                                                 new_file=new_morph_name,
                                                                                 original_position=neuron.position,
-                                                                                original_rotation=neuron.rotation)
+                                                                                original_rotation=neuron.rotation,
+                                                                                random_seed=random_seed)
 
                 if new_morphology:
                     # Replace the original morphology with the warped morphology, morphology includes rotation