diff --git a/examples/notebooks/README.md b/examples/notebooks/README.md index 4074e606b..4e6f82ed6 100644 --- a/examples/notebooks/README.md +++ b/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. diff --git a/examples/notebooks/bend_morphologies.ipynb b/examples/notebooks/bend_morphologies.ipynb index 53805f18d..2c37bacc6 100644 --- a/examples/notebooks/bend_morphologies.ipynb +++ b/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": [] diff --git a/snudda/place/bend_morphologies.py b/snudda/place/bend_morphologies.py index 66aa6f324..17df4dc24 100644 --- a/snudda/place/bend_morphologies.py +++ b/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,15 +96,14 @@ 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): @@ -105,14 +111,9 @@ def bend_morphology(self, morphology: NeuronMorphologyExtended, k=30e-6, n_rando 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) diff --git a/snudda/place/place.py b/snudda/place/place.py index 08791b2d5..ef32dab14 100644 --- a/snudda/place/place.py +++ b/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