From c6fdfefb6b06ffe006ee73e59d96ddd84d6bac5e Mon Sep 17 00:00:00 2001 From: arnaudon Date: Wed, 19 Apr 2023 14:30:00 +0200 Subject: [PATCH 01/63] squashed --- .pre-commit-config.yaml | 2 +- region_grower/atlas_helper.py | 2 +- region_grower/context.py | 100 ++++++++++++++++++++--- region_grower/synthesize_morphologies.py | 21 ++++- setup.py | 1 + tests/data_factories.py | 2 +- tests/test_context.py | 37 ++++----- tox.ini | 4 +- 8 files changed, 132 insertions(+), 37 deletions(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 6fa9e94..92d0ae0 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -1,5 +1,5 @@ default_language_version: - python: python3.8 + python: python3 repos: - repo: https://github.com/pre-commit/pre-commit-hooks rev: v4.4.0 diff --git a/region_grower/atlas_helper.py b/region_grower/atlas_helper.py index f41b8da..98c37fb 100644 --- a/region_grower/atlas_helper.py +++ b/region_grower/atlas_helper.py @@ -32,7 +32,7 @@ def __init__(self, atlas: Atlas, region_structure_path: str): with open(region_structure_path, "r", encoding="utf-8") as region_file: self.region_structure = yaml.safe_load(region_file) else: - raise ValueError("Please provide an existing region_structure file.") + raise ValueError(f"region_structure file at {region_structure_path} not found.") self.regions = list(self.region_structure.keys()) self.layers = {} diff --git a/region_grower/context.py b/region_grower/context.py index 5a46284..62e802c 100644 --- a/region_grower/context.py +++ b/region_grower/context.py @@ -7,12 +7,14 @@ the placement_algorithm package to synthesize circuit morphologies. """ import os +import json from collections import defaultdict from copy import deepcopy from typing import Dict from typing import List from typing import Optional from typing import Union +import trimesh os.environ.setdefault("NEURON_MODULE_OPTIONS", "-nogui") # suppress NEURON warning # This environment variable must be set before 'NEURON' is loaded in 'morph_tool.nrnhines'. @@ -83,6 +85,86 @@ class SpaceContext: layer_depths: List cortical_depths: List + boundaries: Dict = None + atlas_info: Dict = {} # voxel_dimnesions, offset and shape from atlas for indices conversion + + def indices_to_positions(self, indices): + """Local version of voxcel's function to prevent atlas loading.""" + return indices * np.array(self.atlas_info["voxel_dimensions"]) + np.array( + self.atlas_info["offset"] + ) + + def positions_to_indices(self, positions): + """Local and reduced version of voxcel's function to prevent atlas loading.""" + result = (positions - np.array(self.atlas_info["offset"])) / np.array( + self.atlas_info["voxel_dimensions"] + ) + result[np.abs(result) < 1e-7] = 0.0 # suppress rounding errors around 0 + result = np.floor(result).astype(int) + result[result < 0] = -1 + result[result >= self.atlas_info["shape"]] = -1 + return result + + def get_boundaries(self): + """Returns a dict with boundaries data for NeuroTS.""" + # add here necessary logic to convert raw config data to NeuroTS context data + self.boundaries = json.loads(self.boundaries) + mesh = trimesh.load_mesh(self.boundaries["path"]) + + # to save points for plotting while growing to debug the prob function + debug = True + + def get_distance_to_mesh(mesh, ray_origin, ray_direction): + """Compute distances from point/directions to a mesh.""" + vox_ray_origin = self.positions_to_indices(ray_origin) + locations = mesh.ray.intersects_location( + ray_origins=[vox_ray_origin], ray_directions=[ray_direction] + )[0] + if len(locations): + intersect = self.indices_to_positions(locations[0]) + dist = np.linalg.norm(intersect - ray_origin) + + if debug: + x, y, z = ray_origin + vx, vy, vz = ray_direction + with open("data.csv", "a", encoding="utf8") as file: + print(f"{x}, {y}, {z}, {vx}, {vy}, {vz}, {dist}", file=file) + return dist + if debug: + x, y, z = ray_origin + vx, vy, vz = ray_direction + with open("data.csv", "a", encoding="utf8") as file: + print(f"{x}, {y}, {z}, {vx}, {vy}, {vz}, {-100}", file=file) + return np.inf + + def _prob(dist, params): + """Probability function from distance to boundary.""" + p = (dist - params.get("d_min", 0)) / ( + params.get("d_max", 100) - params.get("d_min", 0) + ) ** params.get("power", 1.5) + + # TODO: other variants for this function coud be included here + return np.clip(p, 0, 1) + + if "params_section" in self.boundaries: + + def section_prob(direction, current_point): + """Probability function for the sections.""" + dist = get_distance_to_mesh(mesh, current_point, direction) + return _prob(dist, self.boundaries["params_section"]) + + self.boundaries["section_prob"] = section_prob + + if "params_trunk" in self.boundaries: + + def trunk_prob(direction, current_point): + """Probability function for the trunks.""" + dist = get_distance_to_mesh(mesh, current_point, direction) + return _prob(dist, self.boundaries["params_trunk"]) + + self.boundaries["trunk_prob"] = trunk_prob + + return self.boundaries def layer_fraction_to_position(self, layer: int, layer_fraction: float) -> float: """Returns an absolute position from a layer and a fraction of the layer. @@ -141,7 +223,6 @@ class SynthesisParameters: axon_morph_scale: Optional[float] = None rotational_jitter_std: float = None scaling_jitter_std: float = None - recenter: bool = True seed: int = 0 min_hard_scale: float = 0 @@ -313,14 +394,7 @@ def completion(self, synthesized): def _synthesize_once(self, rng) -> SynthesisResult: """One try to synthesize the cell.""" params = self._correct_position_orientation_scaling(self.params.tmd_parameters) - - # Today we don't use the atlas during the synthesis (we just use it to - # generate the parameters) so we can - # grow the cell as if it was in [0, 0, 0] - # But the day we use it during the actual growth, we will need to grow the cell at its - # absolute position and translate to [0, 0, 0] after the growth - if self.params.recenter: - params["origin"] = [0, 0, 0] + params["origin"] = self.cell.position.tolist() if self.params.tmd_parameters["diameter_params"]["method"] == "external": external_diametrizer = build_diameters.build @@ -332,16 +406,20 @@ def _synthesize_once(self, rng) -> SynthesisResult: else: axon_morph = None + context = {"debug_data": self.debug_infos["input_scaling"]} + if self.context.boundaries is not None: + context.update(self.context.get_boundaries()) + grower = NeuronGrower( input_parameters=params, input_distributions=self.params.tmd_distributions, external_diametrizer=external_diametrizer, skip_preprocessing=True, - context={"debug_data": self.debug_infos["input_scaling"]}, + context=context, rng_or_seed=rng, ) grower.grow() - + mt.translate(grower.neuron, -self.cell.position) self._post_growth_rescaling(grower, params) if axon_morph is not None: diff --git a/region_grower/synthesize_morphologies.py b/region_grower/synthesize_morphologies.py index 2d7b772..6cda8d3 100644 --- a/region_grower/synthesize_morphologies.py +++ b/region_grower/synthesize_morphologies.py @@ -126,7 +126,9 @@ def _parallel_wrapper( ) current_space_context = SpaceContext( layer_depths=row["layer_depths"], - cortical_depths=cortical_depths[row["synthesis_region"]], + cortical_depths=cortical_depths, + boundaries=row["boundaries"] if "boundaries" in row else None, + atlas_info=json.loads(row["atlas_info"]) if "atlas_info" in row else None, ) axon_scale = row.get("axon_scale", None) @@ -144,7 +146,6 @@ def _parallel_wrapper( axon_morph_scale=axon_scale, rotational_jitter_std=rotational_jitter_std, scaling_jitter_std=scaling_jitter_std, - recenter=True, seed=row["seed"], min_hard_scale=min_hard_scale, ) @@ -326,12 +327,28 @@ def __init__( self.cells_data["synthesis_region"] = self.cells_data["region"].apply( lambda region: self.region_mapper[region] ) + self.cells_data["current_depth"] = current_depths + self.cells_data["layer_depths"] = layer_depths.T.tolist() + + if "boundaries" in self.atlas.region_structure: + self.cells_data["boundaries"] = json.dumps(self.atlas.region_structure["boundaries"]) + self.cells_data["atlas_info"] = json.dumps( + { + "voxel_dimensions": self.atlas.depths.voxel_dimensions.tolist(), + "offset": self.atlas.depths.offset.tolist(), + "shape": self.atlas.depths.shape, + } + ) + + if not self.cells.orientations: # pragma: no cover + self.cells_data["orientation"] = orientations.tolist() LOGGER.info("Checking TMD parameters and distributions according to cells mtypes") self.verify() LOGGER.info("Fetching atlas data from %s", atlas) self.assign_atlas_data(min_depth, max_depth) + if morph_axon is not None: LOGGER.info("Loading axon morphologies from %s", morph_axon) self.axon_morph_list = load_morphology_list(morph_axon, self.task_ids) diff --git a/setup.py b/setup.py index 383b085..a20d4e8 100644 --- a/setup.py +++ b/setup.py @@ -28,6 +28,7 @@ "tqdm>=4.28.1", "urllib3>=1.26,<2; python_version < '3.9'", "voxcell>=3.1.1,<4", + "trimesh", ] mpi_extras = [ diff --git a/tests/data_factories.py b/tests/data_factories.py index a496ba5..565ce2e 100644 --- a/tests/data_factories.py +++ b/tests/data_factories.py @@ -132,7 +132,7 @@ def get_tmd_distributions(filename): def get_cell_position(): """The cell position.""" - return [0, 500, 0] + return np.array([0, 500, 0]) def get_cell_mtype(): diff --git a/tests/test_context.py b/tests/test_context.py index bc1736e..934e585 100644 --- a/tests/test_context.py +++ b/tests/test_context.py @@ -116,19 +116,18 @@ def test_context( assert_array_almost_equal( synthesis_parameters.tmd_parameters["apical_dendrite"]["orientation"], [[0.0, 1.0, 0.0]] ) - assert_array_almost_equal( result.neuron.soma.points, np.array( [ - [-5.785500526428223, 4.9841227531433105, 0.0], - [-7.5740227699279785, -0.9734848141670227, 0.0], - [-1.966903805732727, -7.378671169281006, 0.0], - [3.565324068069458, -6.752922534942627, 0.0], - [7.266839027404785, -2.346604108810425, 0.0], - [7.384983062744141, 1.059786081314087, 0.0], - [6.818241119384766, 3.4387624263763428, 0.0], - [4.675901919111924e-16, 7.636327266693115, 0.0], + [-5.785500526428223, 4.984130859375, 0.0], + [-7.5740227699279785, -0.973480224609375, 0.0], + [-1.966903805732727, -7.378662109375, 0.0], + [3.565324068069458, -6.7529296875, 0.0], + [7.266839027404785, -2.34661865234375, 0.0], + [7.384983062744141, 1.059783935546875, 0.0], + [6.818241119384766, 3.438751220703125, 0.0], + [4.675901919111924e-16, 7.636322021484375, 0.0], ], dtype=np.float32, ), @@ -189,25 +188,25 @@ def test_context_external_diametrizer( assert_array_almost_equal( synthesis_parameters.tmd_parameters["apical_dendrite"]["orientation"], [[0.0, 1.0, 0.0]] ) - assert_array_almost_equal( result.neuron.soma.points, np.array( [ - [-5.785500526428223, 4.9841227531433105, 0.0], - [-7.5740227699279785, -0.9734848141670227, 0.0], - [-1.966903805732727, -7.378671169281006, 0.0], - [3.565324068069458, -6.752922534942627, 0.0], - [7.266839027404785, -2.346604108810425, 0.0], - [7.384983062744141, 1.059786081314087, 0.0], - [6.818241119384766, 3.4387624263763428, 0.0], - [4.675901919111924e-16, 7.636327266693115, 0.0], + [-5.785500526428223, 4.984130859375, 0.0], + [-7.5740227699279785, -0.973480224609375, 0.0], + [-1.966903805732727, -7.378662109375, 0.0], + [3.565324068069458, -6.7529296875, 0.0], + [7.266839027404785, -2.34661865234375, 0.0], + [7.384983062744141, 1.059783935546875, 0.0], + [6.818241119384766, 3.438751220703125, 0.0], + [4.675901919111924e-16, 7.636322021484375, 0.0], ], dtype=np.float32, ), ) assert len(result.neuron.root_sections) == 4 + assert_array_almost_equal( next(result.neuron.iter()).points, np.array( @@ -258,7 +257,7 @@ def test_scale(self, small_context_worker, tmd_parameters): assert_array_almost_equal( [ # Check only first and last points of neurites - np.around(np.array([neu.points[0], neu.points[-1]]), 6) + np.around(np.array([neu.points[0], neu.points[-1]]), 5) for neu in result.neuron.root_sections ], [ diff --git a/tox.ini b/tox.ini index 5227512..649008e 100644 --- a/tox.ini +++ b/tox.ini @@ -81,7 +81,7 @@ commands = basepython = python3.8 [testenv:lint] -basepython = python3.8 +basepython = python3 deps = pre-commit pylint @@ -90,7 +90,7 @@ commands = pylint -j {env:PYLINT_NPROCS:1} {[base]files} [testenv:format] -basepython = python3.8 +basepython = python3 skip_install = true deps = codespell From aee7a0aad6f00f8a5a988cfb3b76ca0296a51a30 Mon Sep 17 00:00:00 2001 From: arnaudon Date: Mon, 15 May 2023 16:13:50 +0200 Subject: [PATCH 02/63] improve --- region_grower/synthesize_morphologies.py | 31 ++++++++++-------------- 1 file changed, 13 insertions(+), 18 deletions(-) diff --git a/region_grower/synthesize_morphologies.py b/region_grower/synthesize_morphologies.py index 6cda8d3..4dce458 100644 --- a/region_grower/synthesize_morphologies.py +++ b/region_grower/synthesize_morphologies.py @@ -109,7 +109,6 @@ def inverse_mapper(self): def _parallel_wrapper( row, computation_parameters, - cortical_depths, rotational_jitter_std, scaling_jitter_std, min_hard_scale, @@ -126,7 +125,7 @@ def _parallel_wrapper( ) current_space_context = SpaceContext( layer_depths=row["layer_depths"], - cortical_depths=cortical_depths, + cortical_depths=row["cortical_depths"], boundaries=row["boundaries"] if "boundaries" in row else None, atlas_info=json.loads(row["atlas_info"]) if "atlas_info" in row else None, ) @@ -327,21 +326,6 @@ def __init__( self.cells_data["synthesis_region"] = self.cells_data["region"].apply( lambda region: self.region_mapper[region] ) - self.cells_data["current_depth"] = current_depths - self.cells_data["layer_depths"] = layer_depths.T.tolist() - - if "boundaries" in self.atlas.region_structure: - self.cells_data["boundaries"] = json.dumps(self.atlas.region_structure["boundaries"]) - self.cells_data["atlas_info"] = json.dumps( - { - "voxel_dimensions": self.atlas.depths.voxel_dimensions.tolist(), - "offset": self.atlas.depths.offset.tolist(), - "shape": self.atlas.depths.shape, - } - ) - - if not self.cells.orientations: # pragma: no cover - self.cells_data["orientation"] = orientations.tolist() LOGGER.info("Checking TMD parameters and distributions according to cells mtypes") self.verify() @@ -538,6 +522,18 @@ def assign_atlas_data(self, min_depth=25, max_depth=5000): dtype=object, ) + if "boundaries" in self.atlas.region_structure[_region]: + self.cells_data.loc[region_mask, "boundaries"] = json.dumps( + self.atlas.region_structure[_region]["boundaries"] + ) + self.cells_data.loc[region_mask, "atlas_info"] = json.dumps( + { + "voxel_dimensions": self.atlas.depths[_region].voxel_dimensions.tolist(), + "offset": self.atlas.depths[_region].offset.tolist(), + "shape": self.atlas.depths[_region].shape, + } + ) + @property def task_ids(self): """Task IDs (= CellCollection IDs).""" @@ -601,7 +597,6 @@ def compute(self): func_kwargs = { "computation_parameters": computation_parameters, - "cortical_depths": self.cortical_depths, "rotational_jitter_std": self.rotational_jitter_std, "scaling_jitter_std": self.scaling_jitter_std, "min_hard_scale": self.min_hard_scale, From 828cd00bd9076cbd96dbfc947ec47ab8b51a36af Mon Sep 17 00:00:00 2001 From: arnaudon Date: Mon, 15 May 2023 18:50:47 +0200 Subject: [PATCH 03/63] fix CI --- region_grower/context.py | 2 +- setup.py | 3 +- tests/conftest.py | 31 ++- tests/data/apical_boundary.yaml | 16 ++ tests/data_factories.py | 22 ++ tests/test_context.py | 14 +- .../test_synthesize_morphologies_boundary.py | 212 ++++++++++++++++++ tox.ini | 4 +- 8 files changed, 297 insertions(+), 7 deletions(-) create mode 100644 tests/data/apical_boundary.yaml create mode 100644 tests/test_synthesize_morphologies_boundary.py diff --git a/region_grower/context.py b/region_grower/context.py index 62e802c..cb25055 100644 --- a/region_grower/context.py +++ b/region_grower/context.py @@ -143,7 +143,7 @@ def _prob(dist, params): params.get("d_max", 100) - params.get("d_min", 0) ) ** params.get("power", 1.5) - # TODO: other variants for this function coud be included here + # TODO: other variants for this function could be included here return np.clip(p, 0, 1) if "params_section" in self.boundaries: diff --git a/setup.py b/setup.py index a20d4e8..40e1358 100644 --- a/setup.py +++ b/setup.py @@ -28,7 +28,7 @@ "tqdm>=4.28.1", "urllib3>=1.26,<2; python_version < '3.9'", "voxcell>=3.1.1,<4", - "trimesh", + "trimesh>=3.21", ] mpi_extras = [ @@ -54,6 +54,7 @@ "pytest-html>=2", "pytest-mock>=3.5", "pytest-xdist>=3.0.2", + "neurocollage>=0.3.0", ] setup( diff --git a/tests/conftest.py b/tests/conftest.py index ad40e30..75a12e9 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -19,6 +19,8 @@ from .data_factories import generate_cell_collection from .data_factories import generate_cells_df from .data_factories import generate_input_cells +from .data_factories import generate_mesh +from .data_factories import generate_region_structure_boundary from .data_factories import generate_small_O1 from .data_factories import get_cell_mtype from .data_factories import get_cell_orientation @@ -37,6 +39,25 @@ def small_O1_path(tmpdir_factory): return atlas_directory +@pytest.fixture(scope="session") +def mesh(small_O1_path, tmpdir_factory): + """Generate mesh from atlas.""" + mesh_path = str(tmpdir_factory.mktemp("mesh") / "mesh.obj") + + atlas = {"atlas": small_O1_path, "structure": DATA / "region_structure.yaml"} + generate_mesh(atlas, mesh_path) + return mesh_path + + +@pytest.fixture(scope="session") +def region_structure_boundary(tmpdir_factory, mesh): + """Generate region_structure file with boundary entries.""" + out_path = str(tmpdir_factory.mktemp("structure") / "region_structure.yaml") + region_structure_path = DATA / "region_structure.yaml" + generate_region_structure_boundary(region_structure_path, out_path, mesh) + return out_path + + @pytest.fixture(scope="session") def small_O1(small_O1_path): """Open the atlas.""" @@ -113,7 +134,15 @@ def space_context(): """A space context object.""" layer_depth = [0.0, 200.0, 300.0, 400.0, 500.0, 600.0, 800.0] thicknesses = [165, 149, 353, 190, 525, 700] - return SpaceContext(layer_depths=layer_depth, cortical_depths=np.cumsum(thicknesses)) + atlas_info = { + "voxel_dimensions": [100.0, 100.0, 100.0], + "offset": [-1100.0, -100.0, -1000.0], + "shape": (22, 10, 20), + } + + return SpaceContext( + layer_depths=layer_depth, cortical_depths=np.cumsum(thicknesses), atlas_info=atlas_info + ) @pytest.fixture(scope="function") diff --git a/tests/data/apical_boundary.yaml b/tests/data/apical_boundary.yaml new file mode 100644 index 0000000..6663c73 --- /dev/null +++ b/tests/data/apical_boundary.yaml @@ -0,0 +1,16 @@ +14a03569d26b949692e5dfe8cb1855fe: +- 2.271697998046875 +- 114.72030639648438 +- -5.1029510498046875 +216363698b529b4a97b750923ceb3ffd: +- -3.509429931640625 +- 115.28182983398438 +- 5.8795013427734375 +4462ebfc5f915ef09cfbac6e7687a66e: +- 12.257843017578125 +- 85.43954467773438 +- -19.165298461914062 +e3e70682c2094cac629f6fbed82c07cd: +- -38.110015869140625 +- 152.74227905273438 +- -20.128509521484375 diff --git a/tests/data_factories.py b/tests/data_factories.py index 565ce2e..b2566e6 100644 --- a/tests/data_factories.py +++ b/tests/data_factories.py @@ -8,6 +8,8 @@ import numpy as np import pandas as pd +import yaml +from neurocollage.mesh_helper import MeshHelper from voxcell import CellCollection DF_SIZE = 12 @@ -33,6 +35,26 @@ def generate_small_O1(directory): return str(directory) +def generate_mesh(atlas, mesh_path): + """Generate a mesh from atlas to test boundary code.""" + mesh_helper = MeshHelper(atlas, "O0") + mesh = mesh_helper.get_boundary_mesh() + mesh.export(mesh_path) # pylint: disable=no-member + + +def generate_region_structure_boundary(region_structure_path, out_path, mesh): + """Generate region_structure file with boundary entries.""" + with open(region_structure_path, encoding="utf-8") as f: + structure = yaml.safe_load(f) + structure["O0"]["boundaries"] = { + "params_section": {"d_min": 5, "d_max": 50}, + "params_trunk": {"d_min": 5, "d_max": 500}, + "path": mesh, + } + with open(out_path, "w", encoding="utf-8") as f: + yaml.dump(structure, f) + + def generate_cells_df(): """Raw data for the cell collection.""" x = [200] * DF_SIZE diff --git a/tests/test_context.py b/tests/test_context.py index 934e585..f3f16b8 100644 --- a/tests/test_context.py +++ b/tests/test_context.py @@ -78,9 +78,17 @@ def test_distance_to_constraint(self, space_context): assert space_context.distance_to_constraint(-100, {"layer": 6, "fraction": 0}) == -900 assert space_context.distance_to_constraint(1000, {"layer": 6, "fraction": 0}) == 200 - for i in [np.nan, None, []]: - space_context.layer_depths = i - assert space_context.distance_to_constraint(1000, {"layer": 6, "fraction": 0}) is None + def test_indices_to_positions(self, space_context): + pos = space_context.indices_to_positions([0, 0, 0]) + assert_array_equal(pos, [-1100.0, -100.0, -1000.0]) + pos = space_context.indices_to_positions([100, 100, 100]) + assert_array_equal(pos, [8900.0, 9900.0, 9000.0]) + + def test_positions_to_indices(self, space_context): + indices = space_context.positions_to_indices([-1100.0, -100.0, -1000.0]) + assert_array_equal(indices, [0, 0, 0]) + indices = space_context.positions_to_indices([8900.0, 9900.0, -9000.0]) + assert_array_equal(indices, [-1, -1, -1]) class TestSpaceWorker: diff --git a/tests/test_synthesize_morphologies_boundary.py b/tests/test_synthesize_morphologies_boundary.py new file mode 100644 index 0000000..ef5b74b --- /dev/null +++ b/tests/test_synthesize_morphologies_boundary.py @@ -0,0 +1,212 @@ +"""Test the region_grower.synthesize_morphologies module.""" +# pylint: disable=missing-function-docstring +import logging +import shutil +from pathlib import Path +from uuid import uuid4 + +from morph_tool.utils import iter_morphology_files +from numpy.testing import assert_allclose + +from region_grower.synthesize_morphologies import SynthesizeMorphologies + +from .test_synthesize_morphologies import check_yaml + +DATA = Path(__file__).parent / "data" + + +def create_args( + with_mpi, + tmp_folder, + input_cells, + atlas_path, + axon_morph_tsv, + out_apical_NRN_sections, + min_depth, + region_structure_boundary, +): + """Create the arguments used for tests.""" + args = {} + + # Circuit + args["input_cells"] = input_cells + + # Atlas + args["atlas"] = atlas_path + + # Parameters + args["tmd_distributions"] = DATA / "distributions.json" + args["tmd_parameters"] = DATA / "parameters.json" + args["seed"] = 0 + args["min_depth"] = min_depth + + # Internals + args["overwrite"] = True + args["out_morph_ext"] = ["h5", "swc", "asc"] + args["out_morph_dir"] = tmp_folder + args["out_apical"] = tmp_folder / "apical.yaml" + args["out_cells"] = str(tmp_folder / "test_cells.mvd3") + if out_apical_NRN_sections: + args["out_apical_nrn_sections"] = tmp_folder / out_apical_NRN_sections + else: + args["out_apical_nrn_sections"] = None + if with_mpi: + args["with_mpi"] = with_mpi + else: + args["nb_processes"] = 2 + + # Axons + args["base_morph_dir"] = str(DATA / "input-cells") + args["morph_axon"] = axon_morph_tsv + args["max_drop_ratio"] = 0.5 + args["rotational_jitter_std"] = 10 + args["scaling_jitter_std"] = 0.5 + args["region_structure"] = region_structure_boundary + + return args + + +def test_synthesize( + tmpdir, + small_O1_path, + input_cells, + axon_morph_tsv, + mesh, + region_structure_boundary, +): # pylint: disable=unused-argument + """Test morphology synthesis but skip write step.""" + with_axon = True + with_NRN = True + min_depth = 25 + tmp_folder = Path(tmpdir) + + args = create_args( + False, + tmp_folder, + input_cells, + small_O1_path, + axon_morph_tsv if with_axon else None, + "apical_NRN_sections.yaml" if with_NRN else None, + min_depth, + region_structure_boundary, + ) + args["skip_write"] = True + + synthesizer = SynthesizeMorphologies(**args) + res = synthesizer.synthesize() + + assert (res["x"] == 200).all() + assert (res["y"] == 200).all() + assert (res["z"] == 200).all() + assert res["name"].tolist() == [ + "e3e70682c2094cac629f6fbed82c07cd", + None, + "216363698b529b4a97b750923ceb3ffd", + None, + "14a03569d26b949692e5dfe8cb1855fe", + None, + "4462ebfc5f915ef09cfbac6e7687a66e", + None, + ] + print(res["apical_points"]) + assert_allclose(res["apical_points"][0], [[-38.110016, 152.74228, -20.12851]]) + assert res["apical_points"][1] is None + + assert_allclose(res["apical_points"][3], [[-3.50943, 115.28183, 5.8795013]]) + assert res["apical_points"][4] is None + + assert_allclose(res["apical_points"][6], [[2.271698, 114.72031, -5.102951]]) + assert res["apical_points"][7] is None + + assert_allclose(res["apical_points"][9], [[12.257843, 85.439545, -19.165298]]) + assert res["apical_points"][10] is None + + # Check that the morphologies were not written + res_files = tmpdir.listdir() + assert len(res_files) == 5 + assert sorted(i.basename for i in res_files) == [ + "apical.yaml", + "apical_NRN_sections.yaml", + "axon_morphs.tsv", + "input_cells.mvd3", + "test_cells.mvd3", + ] + + +def run_with_mpi(): + """Test morphology synthesis with MPI.""" + # pylint: disable=import-outside-toplevel, too-many-locals, import-error + from data_factories import generate_axon_morph_tsv + from data_factories import generate_cell_collection + from data_factories import generate_cells_df + from data_factories import generate_input_cells + from data_factories import generate_mesh + from data_factories import generate_region_structure_boundary + from data_factories import generate_small_O1 + from data_factories import input_cells_path + from mpi4py import MPI + + from region_grower.utils import setup_logger + + COMM = MPI.COMM_WORLD # pylint: disable=c-extension-no-member + rank = COMM.Get_rank() + MASTER_RANK = 0 + is_master = rank == MASTER_RANK + + tmp_folder = Path("/tmp/test-run-synthesis_" + str(uuid4())) + tmp_folder = COMM.bcast(tmp_folder, root=MASTER_RANK) + input_cells = input_cells_path(tmp_folder) + small_O1_path = str(tmp_folder / "atlas") + region_structure_path = str(tmp_folder / "region_structure_boundary.yaml") + args = create_args( + True, + tmp_folder, + input_cells, + small_O1_path, + tmp_folder / "axon_morphs.tsv", + "apical_NRN_sections.yaml", + min_depth=25, + region_structure_boundary=region_structure_path, + ) + + setup_logger("info") + logging.getLogger("distributed").setLevel(logging.ERROR) + + if is_master: + tmp_folder.mkdir(exist_ok=True) + print(f"============= #{rank}: Create data") + cells_raw_data = generate_cells_df() + cell_collection = generate_cell_collection(cells_raw_data) + generate_input_cells(cell_collection, tmp_folder) + generate_small_O1(small_O1_path) + atlas = {"atlas": small_O1_path, "structure": DATA / "region_structure.yaml"} + generate_mesh(atlas, tmp_folder / "mesh.obj") + region_structure_path = generate_region_structure_boundary( + DATA / "region_structure.yaml", region_structure_path, str(tmp_folder / "mesh.obj") + ) + + generate_axon_morph_tsv(tmp_folder) + for dest in range(1, COMM.Get_size()): + req = COMM.isend("done", dest=dest) + else: + req = COMM.irecv(source=0) + req.wait() + + synthesizer = SynthesizeMorphologies(**args) + try: + print(f"============= #{rank}: Start synthesize") + synthesizer.synthesize() + + # Check results + print(f"============= #{rank}: Checking results") + expected_size = 12 + assert len(list(iter_morphology_files(tmp_folder))) == expected_size + check_yaml(DATA / "apical_boundary.yaml", args["out_apical"]) + finally: + # Clean the directory + print(f"============= #{rank}: Cleaning") + shutil.rmtree(tmp_folder) + + +if __name__ == "__main__": # pragma: no cover + run_with_mpi() diff --git a/tox.ini b/tox.ini index 649008e..5669993 100644 --- a/tox.ini +++ b/tox.ini @@ -49,7 +49,9 @@ deps = pytest allowlist_externals = mpiexec -commands = mpiexec -n 4 {envbindir}/python tests/test_synthesize_morphologies.py +commands = + mpiexec -n 4 {envbindir}/python tests/test_synthesize_morphologies.py + mpiexec -n 4 {envbindir}/python tests/test_synthesize_morphologies_boundary.py [testenv:coverage] skip_install = true From 09d037f10d68f7c72fe416d6fce740da29f4b1a6 Mon Sep 17 00:00:00 2001 From: arnaudon Date: Mon, 15 May 2023 21:09:10 +0200 Subject: [PATCH 04/63] set debug off --- region_grower/context.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/region_grower/context.py b/region_grower/context.py index cb25055..c970e93 100644 --- a/region_grower/context.py +++ b/region_grower/context.py @@ -112,7 +112,7 @@ def get_boundaries(self): mesh = trimesh.load_mesh(self.boundaries["path"]) # to save points for plotting while growing to debug the prob function - debug = True + debug = False def get_distance_to_mesh(mesh, ray_origin, ray_direction): """Compute distances from point/directions to a mesh.""" @@ -163,7 +163,6 @@ def trunk_prob(direction, current_point): return _prob(dist, self.boundaries["params_trunk"]) self.boundaries["trunk_prob"] = trunk_prob - return self.boundaries def layer_fraction_to_position(self, layer: int, layer_fraction: float) -> float: From de73ace463cbe1e014272614dd5c7b83c2ae4412 Mon Sep 17 00:00:00 2001 From: arnaudon Date: Tue, 16 May 2023 11:57:38 +0200 Subject: [PATCH 05/63] get cov to 100 --- region_grower/context.py | 7 +- setup.py | 1 + tests/conftest.py | 10 - tests/data/apical_boundary.yaml | 24 +- tests/data/apical_no_axon.yaml | 30 +-- tests/data/parameters.json | 58 +++-- tests/data_factories.py | 17 +- tests/test_context.py | 54 ++--- tests/test_synthesize_morphologies.py | 195 +++++++++++++++- .../test_synthesize_morphologies_boundary.py | 212 ------------------ tox.ini | 2 +- 11 files changed, 302 insertions(+), 308 deletions(-) delete mode 100644 tests/test_synthesize_morphologies_boundary.py diff --git a/region_grower/context.py b/region_grower/context.py index c970e93..eed8bf9 100644 --- a/region_grower/context.py +++ b/region_grower/context.py @@ -124,18 +124,18 @@ def get_distance_to_mesh(mesh, ray_origin, ray_direction): intersect = self.indices_to_positions(locations[0]) dist = np.linalg.norm(intersect - ray_origin) - if debug: + if debug: # pragma: no cover x, y, z = ray_origin vx, vy, vz = ray_direction with open("data.csv", "a", encoding="utf8") as file: print(f"{x}, {y}, {z}, {vx}, {vy}, {vz}, {dist}", file=file) return dist - if debug: + if debug: # pragma: no cover x, y, z = ray_origin vx, vy, vz = ray_direction with open("data.csv", "a", encoding="utf8") as file: print(f"{x}, {y}, {z}, {vx}, {vy}, {vz}, {-100}", file=file) - return np.inf + return np.inf # pragma: no cover def _prob(dist, params): """Probability function from distance to boundary.""" @@ -322,7 +322,6 @@ def _post_growth_rescaling(self, grower: NeuronGrower, params: Dict) -> None: target_min_length=target_min_length, target_max_length=target_max_length, ) - if scale > self.params.min_hard_scale: scale_section(root_section, ScaleParameters(mean=scale), recursive=True) is_deleted = False diff --git a/setup.py b/setup.py index 40e1358..e663045 100644 --- a/setup.py +++ b/setup.py @@ -29,6 +29,7 @@ "urllib3>=1.26,<2; python_version < '3.9'", "voxcell>=3.1.1,<4", "trimesh>=3.21", + "rtree>=1.0.1", ] mpi_extras = [ diff --git a/tests/conftest.py b/tests/conftest.py index 75a12e9..7eaf84e 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -20,7 +20,6 @@ from .data_factories import generate_cells_df from .data_factories import generate_input_cells from .data_factories import generate_mesh -from .data_factories import generate_region_structure_boundary from .data_factories import generate_small_O1 from .data_factories import get_cell_mtype from .data_factories import get_cell_orientation @@ -49,15 +48,6 @@ def mesh(small_O1_path, tmpdir_factory): return mesh_path -@pytest.fixture(scope="session") -def region_structure_boundary(tmpdir_factory, mesh): - """Generate region_structure file with boundary entries.""" - out_path = str(tmpdir_factory.mktemp("structure") / "region_structure.yaml") - region_structure_path = DATA / "region_structure.yaml" - generate_region_structure_boundary(region_structure_path, out_path, mesh) - return out_path - - @pytest.fixture(scope="session") def small_O1(small_O1_path): """Open the atlas.""" diff --git a/tests/data/apical_boundary.yaml b/tests/data/apical_boundary.yaml index 6663c73..fe76a95 100644 --- a/tests/data/apical_boundary.yaml +++ b/tests/data/apical_boundary.yaml @@ -1,16 +1,16 @@ 14a03569d26b949692e5dfe8cb1855fe: -- 2.271697998046875 -- 114.72030639648438 -- -5.1029510498046875 +- 4.369110107421875 +- 55.535736083984375 +- -6.2556915283203125 216363698b529b4a97b750923ceb3ffd: -- -3.509429931640625 -- 115.28182983398438 -- 5.8795013427734375 +- 5.216094970703125 +- 114.55657958984375 +- -12.4947509765625 4462ebfc5f915ef09cfbac6e7687a66e: -- 12.257843017578125 -- 85.43954467773438 -- -19.165298461914062 +- -11.084274291992188 +- 111.49124145507812 +- 0.980438232421875 e3e70682c2094cac629f6fbed82c07cd: -- -38.110015869140625 -- 152.74227905273438 -- -20.128509521484375 +- 16.126876831054688 +- 160.52789306640625 +- 22.863800048828125 diff --git a/tests/data/apical_no_axon.yaml b/tests/data/apical_no_axon.yaml index bfd6aa6..424b442 100644 --- a/tests/data/apical_no_axon.yaml +++ b/tests/data/apical_no_axon.yaml @@ -19,9 +19,9 @@ - 127.51496124267578 - 16.63180923461914 6513270e269e0d37f2a74de452e6b438: -- 1.8200105428695679 -- 113.30069732666016 -- 0.27308279275894165 +- -8.12664794921875 +- 114.60183715820312 +- 1.4770965576171875 7b89296c6dcbac5008577eb1924770d3: - -2.5671095848083496 - 58.11355972290039 @@ -31,21 +31,21 @@ - 371.3540344238281 - 6.843407154083252 b8a1abcd1a6916c74da4f9fc3c6da5d7: -- 2.9995338916778564 -- 114.8524169921875 -- 12.829068183898926 +- -13.87677001953125 +- 162.44619750976562 +- 1.412567138671875 cd613e30d8f16adf91b7584a2265b1f5: -- 14.439825057983398 -- 114.19844818115234 -- -7.658063888549805 +- 1.1701507568359375 +- 97.20172119140625 +- 1.061187744140625 d95bafc8f2a4d27bdcf4bb99f4bea973: -- 3.320021867752075 -- 59.75845718383789 -- -2.885958433151245 +- -2.2343902587890625 +- 59.40399169921875 +- -0.331756591796875 db5b5fab8f4d3e27dda1494c73cf256d: -- 10.5548677444458 -- 93.02537536621094 -- 0.9866926074028015 +- -1.606964111328125 +- 58.9842529296875 +- 1.571624755859375 e3e70682c2094cac629f6fbed82c07cd: - -4.8529953956604 - 158.9239959716797 diff --git a/tests/data/parameters.json b/tests/data/parameters.json index bb571ac..3a024fb 100644 --- a/tests/data/parameters.json +++ b/tests/data/parameters.json @@ -7,13 +7,16 @@ "has_apical_tuft": true, "metric": "path_distances", "modify": null, - "orientation": [ - [ - 0.0, - 1.0, - 0.0 - ] - ], + "orientation": { + "mode": "normal_pia_constraint", + "values": { + "direction": { + "mean": 0.0, + "std": 0.0 + } + } + }, + "radius": 0.3, "randomness": 0.15, "step_size": { "norm": { @@ -30,7 +33,17 @@ "growth_method": "tmd", "metric": "path_distances", "modify": null, - "orientation": null, + "orientation": { + "mode": "pia_constraint", + "values": { + "form": "step", + "params": [ + 0.44484279339093336, + 0.5461595884090537 + ] + } + }, + "radius": 0.3, "randomness": 0.15, "step_size": { "norm": { @@ -66,13 +79,16 @@ "has_apical_tuft": true, "metric": "path_distances", "modify": null, - "orientation": [ - [ - 0.0, - 1.0, - 0.0 - ] - ], + "orientation": { + "mode": "normal_pia_constraint", + "values": { + "direction": { + "mean": 0.0, + "std": 0.0 + } + } + }, + "radius": 0.3, "randomness": 0.15, "step_size": { "norm": { @@ -89,7 +105,17 @@ "growth_method": "tmd", "metric": "path_distances", "modify": null, - "orientation": null, + "orientation": { + "mode": "pia_constraint", + "values": { + "form": "step", + "params": [ + 0.44484279339093336, + 0.5461595884090537 + ] + } + }, + "radius": 0.3, "randomness": 0.15, "step_size": { "norm": { diff --git a/tests/data_factories.py b/tests/data_factories.py index b2566e6..457f65d 100644 --- a/tests/data_factories.py +++ b/tests/data_factories.py @@ -42,15 +42,20 @@ def generate_mesh(atlas, mesh_path): mesh.export(mesh_path) # pylint: disable=no-member -def generate_region_structure_boundary(region_structure_path, out_path, mesh): +def generate_region_structure_boundary( + region_structure_path, out_path, mesh, with_sections=True, with_trunks=True +): """Generate region_structure file with boundary entries.""" with open(region_structure_path, encoding="utf-8") as f: structure = yaml.safe_load(f) - structure["O0"]["boundaries"] = { - "params_section": {"d_min": 5, "d_max": 50}, - "params_trunk": {"d_min": 5, "d_max": 500}, - "path": mesh, - } + structure["O0"]["boundaries"] = {"path": mesh} + + if with_sections: + structure["O0"]["boundaries"]["params_section"] = {"d_min": 5, "d_max": 50} + + if with_trunks: + structure["O0"]["boundaries"]["params_trunk"] = {"d_min": 5, "d_max": 5000} + with open(out_path, "w", encoding="utf-8") as f: yaml.dump(structure, f) diff --git a/tests/test_context.py b/tests/test_context.py index f3f16b8..bc8237a 100644 --- a/tests/test_context.py +++ b/tests/test_context.py @@ -121,21 +121,30 @@ def test_context( ) # This tests that input orientations are not mutated by the synthesize() call - assert_array_almost_equal( - synthesis_parameters.tmd_parameters["apical_dendrite"]["orientation"], [[0.0, 1.0, 0.0]] + assert ( + list( + dictdiffer.diff( + synthesis_parameters.tmd_parameters["apical_dendrite"]["orientation"], + { + "mode": "normal_pia_constraint", + "values": {"direction": {"mean": 0.0, "std": 0.0}}, + }, + ) + ) + == [] ) assert_array_almost_equal( result.neuron.soma.points, np.array( [ - [-5.785500526428223, 4.984130859375, 0.0], - [-7.5740227699279785, -0.973480224609375, 0.0], - [-1.966903805732727, -7.378662109375, 0.0], - [3.565324068069458, -6.7529296875, 0.0], - [7.266839027404785, -2.34661865234375, 0.0], - [7.384983062744141, 1.059783935546875, 0.0], - [6.818241119384766, 3.438751220703125, 0.0], - [4.675901919111924e-16, 7.636322021484375, 0.0], + [7.0957971, -2.8218994, 0.0], + [0.0, 7.6363220, 0.0], + [-3.0748143, 6.9899292, 0.0], + [-4.8864875, 5.8681946, 0.0], + [-7.5141201, -1.3606873, 0.0], + [-6.1610136, -3.6210632, 0.0], + [-1.2113731, -7.5396423, 0.0], + [0.67872918, -7.6061096, 0.0], ], dtype=np.float32, ), @@ -447,6 +456,11 @@ def test_scale(self, small_context_worker, tmd_parameters): } } result = small_context_worker.synthesize() + assert [i.type for i in result.neuron.root_sections] == expected_types[1:] + + # test removing basals if scale is too small + small_context_worker.params.min_hard_scale = 2.2 + result = small_context_worker.synthesize() assert [i.type for i in result.neuron.root_sections] == [expected_types[-1]] params["basal_dendrite"]["orientation"] = {} @@ -491,18 +505,9 @@ def test_debug_scales(self, small_context_worker, tmd_parameters): "min_target_thickness": 1.0, }, "scaling": [ - { - "max_ph": 126.96580088057513, - "scaling_ratio": 0.9554140127388535, - }, - { - "max_ph": 139.93140451880743, - "scaling_ratio": 0.9554140127388535, - }, - { - "max_ph": 143.28779114733433, - "scaling_ratio": 0.9554140127388535, - }, + {"max_ph": 78.53590605172101, "scaling_ratio": 0.9554140127388535}, + {"max_ph": 126.96580088057513, "scaling_ratio": 0.9554140127388535}, + {"max_ph": 126.96580088057513, "scaling_ratio": 0.9554140127388535}, ], }, "target_func": { @@ -514,10 +519,7 @@ def test_debug_scales(self, small_context_worker, tmd_parameters): "min_target_path_distance": 1.0, }, "scaling": [ - { - "max_ph": 420.70512274498446, - "scaling_ratio": 0.18064909574697127, - } + {"max_ph": 255.45843100194077, "scaling_ratio": 0.2975043716581138} ], }, }, diff --git a/tests/test_synthesize_morphologies.py b/tests/test_synthesize_morphologies.py index d1d72e2..6e12098 100644 --- a/tests/test_synthesize_morphologies.py +++ b/tests/test_synthesize_morphologies.py @@ -1,6 +1,7 @@ """Test the region_grower.synthesize_morphologies module.""" # pylint: disable=missing-function-docstring import json +import sys import logging import os import shutil @@ -26,6 +27,7 @@ from region_grower.synthesize_morphologies import RegionMapper from region_grower.synthesize_morphologies import SynthesizeMorphologies + DATA = Path(__file__).parent / "data" @@ -53,6 +55,7 @@ def create_args( axon_morph_tsv, out_apical_NRN_sections, min_depth, + region_structure, ): """Create the arguments used for tests.""" args = {} @@ -90,7 +93,7 @@ def create_args( args["max_drop_ratio"] = 0.5 args["rotational_jitter_std"] = 10 args["scaling_jitter_std"] = 0.5 - args["region_structure"] = DATA / "region_structure.yaml" + args["region_structure"] = region_structure return args @@ -118,6 +121,7 @@ def test_synthesize( axon_morph_tsv if with_axon else None, "apical_NRN_sections.yaml" if with_NRN else None, min_depth, + DATA / "region_structure.yaml", ) synthesizer = SynthesizeMorphologies(**args) @@ -233,6 +237,7 @@ def test_synthesize_skip_write( axon_morph_tsv if with_axon else None, "apical_NRN_sections.yaml" if with_NRN else None, min_depth, + DATA / "region_structure.yaml", ) args["skip_write"] = True args["nb_processes"] = nb_processes @@ -285,6 +290,109 @@ def test_synthesize_skip_write( ] +@pytest.mark.parametrize("with_sections", [True, False]) +@pytest.mark.parametrize("with_trunks", [True, False]) +def test_synthesize_boundary( + tmpdir, + small_O1_path, + input_cells, + axon_morph_tsv, + mesh, + with_sections, + with_trunks, +): # pylint: disable=unused-argument + """Test morphology synthesis but skip write step.""" + with_axon = True + with_NRN = True + min_depth = 25 + tmp_folder = Path(tmpdir) + + # pylint: disable=import-outside-toplevel + from .data_factories import generate_region_structure_boundary + + region_structure = "region_structure.yaml" + generate_region_structure_boundary( + DATA / "region_structure.yaml", region_structure, mesh, with_sections, with_trunks + ) + args = create_args( + False, + tmp_folder, + input_cells, + small_O1_path, + axon_morph_tsv if with_axon else None, + "apical_NRN_sections.yaml" if with_NRN else None, + min_depth, + region_structure, + ) + args["skip_write"] = True + + synthesizer = SynthesizeMorphologies(**args) + res = synthesizer.synthesize() + + assert (res["x"] == 200).all() + assert (res["y"] == 200).all() + assert (res["z"] == 200).all() + assert res["name"].tolist() == [ + "e3e70682c2094cac629f6fbed82c07cd", + None, + "216363698b529b4a97b750923ceb3ffd", + None, + "14a03569d26b949692e5dfe8cb1855fe", + None, + "4462ebfc5f915ef09cfbac6e7687a66e", + None, + ] + if with_sections and with_trunks: + assert_allclose(res["apical_points"][0], [[16.126877, 160.5279, 22.8638]]) + assert res["apical_points"][1] is None + + assert_allclose(res["apical_points"][3], [[5.216095, 114.55658, -12.494751]]) + assert res["apical_points"][4] is None + + assert_allclose(res["apical_points"][6], [[4.36911, 55.535736, -6.2556915]]) + assert res["apical_points"][7] is None + + assert_allclose(res["apical_points"][9], [[-11.084274, 111.49124, 0.98043823]]) + assert res["apical_points"][10] is None + + if with_sections and not with_trunks: + assert_allclose(res["apical_points"][0], [[-0.29234314, 58.81488, 0.5384369]]) + assert res["apical_points"][1] is None + + assert_allclose(res["apical_points"][3], [[8.230438, 116.570435, -7.345169]]) + assert res["apical_points"][4] is None + + assert_allclose(res["apical_points"][6], [[11.181992, 96.11371, -12.706863]]) + assert res["apical_points"][7] is None + + assert_allclose(res["apical_points"][9], [[3.5267944, 164.42618, 1.2018433]]) + assert res["apical_points"][10] is None + + if not with_sections and with_trunks: + assert_allclose(res["apical_points"][0], [[8.792313, 131.91104, -4.2198334]]) + assert res["apical_points"][1] is None + + assert_allclose(res["apical_points"][3], [[2.0048676, 116.672, -5.334656]]) + assert res["apical_points"][4] is None + + assert_allclose(res["apical_points"][6], [[-2.347641, 58.229614, -0.56985474]]) + assert res["apical_points"][7] is None + + assert_allclose(res["apical_points"][9], [[6.861679, 112.31445, -13.528854]]) + assert res["apical_points"][10] is None + + # Check that the morphologies were not written + res_files = tmpdir.listdir() + assert len(res_files) == 5 + assert sorted(i.basename for i in res_files) == [ + "apical.yaml", + "apical_NRN_sections.yaml", + "axon_morphs.tsv", + "input_cells.mvd3", + "test_cells.mvd3", + ] + + def run_with_mpi(): """Test morphology synthesis with MPI.""" # pylint: disable=import-outside-toplevel, too-many-locals, import-error @@ -316,6 +424,7 @@ def run_with_mpi(): tmp_folder / "axon_morphs.tsv", "apical_NRN_sections.yaml", min_depth=25, + region_structure=DATA / "region_structure.yaml", ) setup_logger("debug", prefix=f"Rank = {rank} - ") @@ -346,12 +455,83 @@ def run_with_mpi(): print(f"============= #{rank}: Checking results") expected_size = 18 assert len(list(iter_morphology_files(tmp_folder))) == expected_size + check_yaml(DATA / "apical.yaml", args["out_apical"]) + check_yaml(DATA / "apical_NRN_sections.yaml", args["out_apical_nrn_sections"]) + finally: + # Clean the directory + print(f"============= #{rank}: Cleaning") + shutil.rmtree(tmp_folder) + + +def run_with_mpi_boundary(): + """Test morphology synthesis with MPI.""" + # pylint: disable=import-outside-toplevel, too-many-locals, import-error + from data_factories import generate_axon_morph_tsv + from data_factories import generate_cell_collection + from data_factories import generate_cells_df + from data_factories import generate_input_cells + from data_factories import generate_mesh + from data_factories import generate_region_structure_boundary + from data_factories import generate_small_O1 + from data_factories import input_cells_path + from mpi4py import MPI + + from region_grower.utils import setup_logger + + COMM = MPI.COMM_WORLD # pylint: disable=c-extension-no-member + rank = COMM.Get_rank() + MASTER_RANK = 0 + is_master = rank == MASTER_RANK + + tmp_folder = Path("/tmp/test-run-synthesis_" + str(uuid4())) + tmp_folder = COMM.bcast(tmp_folder, root=MASTER_RANK) + input_cells = input_cells_path(tmp_folder) + small_O1_path = str(tmp_folder / "atlas") + region_structure_path = str(tmp_folder / "region_structure_boundary.yaml") + args = create_args( + True, + tmp_folder, + input_cells, + small_O1_path, + tmp_folder / "axon_morphs.tsv", + "apical_NRN_sections.yaml", + min_depth=25, + region_structure=region_structure_path, + ) + + setup_logger("info") + logging.getLogger("distributed").setLevel(logging.ERROR) - check_yaml(DATA / ("apical.yaml"), args["out_apical"]) - check_yaml( - DATA / ("apical_NRN_sections.yaml"), - args["out_apical_nrn_sections"], + if is_master: + tmp_folder.mkdir(exist_ok=True) + print(f"============= #{rank}: Create data") + cells_raw_data = generate_cells_df() + cell_collection = generate_cell_collection(cells_raw_data) + generate_input_cells(cell_collection, tmp_folder) + generate_small_O1(small_O1_path) + atlas = {"atlas": small_O1_path, "structure": DATA / "region_structure.yaml"} + generate_mesh(atlas, tmp_folder / "mesh.obj") + region_structure_path = generate_region_structure_boundary( + DATA / "region_structure.yaml", region_structure_path, str(tmp_folder / "mesh.obj") ) + + generate_axon_morph_tsv(tmp_folder) + for dest in range(1, COMM.Get_size()): + req = COMM.isend("done", dest=dest) + else: + req = COMM.irecv(source=0) + req.wait() + + synthesizer = SynthesizeMorphologies(**args) + try: + print(f"============= #{rank}: Start synthesize") + synthesizer.synthesize() + + # Check results + print(f"============= #{rank}: Checking results") + expected_size = 12 + assert len(list(iter_morphology_files(tmp_folder))) == expected_size + check_yaml(DATA / "apical_boundary.yaml", args["out_apical"]) finally: # Clean the directory print(f"============= #{rank}: Cleaning") @@ -603,4 +783,7 @@ def test_inconsistent_context( if __name__ == "__main__": # pragma: no cover - run_with_mpi() + if sys.argv[-1] == "boundary": + run_with_mpi_boundary() + else: + run_with_mpi() diff --git a/tests/test_synthesize_morphologies_boundary.py b/tests/test_synthesize_morphologies_boundary.py deleted file mode 100644 index ef5b74b..0000000 --- a/tests/test_synthesize_morphologies_boundary.py +++ /dev/null @@ -1,212 +0,0 @@ -"""Test the region_grower.synthesize_morphologies module.""" -# pylint: disable=missing-function-docstring -import logging -import shutil -from pathlib import Path -from uuid import uuid4 - -from morph_tool.utils import iter_morphology_files -from numpy.testing import assert_allclose - -from region_grower.synthesize_morphologies import SynthesizeMorphologies - -from .test_synthesize_morphologies import check_yaml - -DATA = Path(__file__).parent / "data" - - -def create_args( - with_mpi, - tmp_folder, - input_cells, - atlas_path, - axon_morph_tsv, - out_apical_NRN_sections, - min_depth, - region_structure_boundary, -): - """Create the arguments used for tests.""" - args = {} - - # Circuit - args["input_cells"] = input_cells - - # Atlas - args["atlas"] = atlas_path - - # Parameters - args["tmd_distributions"] = DATA / "distributions.json" - args["tmd_parameters"] = DATA / "parameters.json" - args["seed"] = 0 - args["min_depth"] = min_depth - - # Internals - args["overwrite"] = True - args["out_morph_ext"] = ["h5", "swc", "asc"] - args["out_morph_dir"] = tmp_folder - args["out_apical"] = tmp_folder / "apical.yaml" - args["out_cells"] = str(tmp_folder / "test_cells.mvd3") - if out_apical_NRN_sections: - args["out_apical_nrn_sections"] = tmp_folder / out_apical_NRN_sections - else: - args["out_apical_nrn_sections"] = None - if with_mpi: - args["with_mpi"] = with_mpi - else: - args["nb_processes"] = 2 - - # Axons - args["base_morph_dir"] = str(DATA / "input-cells") - args["morph_axon"] = axon_morph_tsv - args["max_drop_ratio"] = 0.5 - args["rotational_jitter_std"] = 10 - args["scaling_jitter_std"] = 0.5 - args["region_structure"] = region_structure_boundary - - return args - - -def test_synthesize( - tmpdir, - small_O1_path, - input_cells, - axon_morph_tsv, - mesh, - region_structure_boundary, -): # pylint: disable=unused-argument - """Test morphology synthesis but skip write step.""" - with_axon = True - with_NRN = True - min_depth = 25 - tmp_folder = Path(tmpdir) - - args = create_args( - False, - tmp_folder, - input_cells, - small_O1_path, - axon_morph_tsv if with_axon else None, - "apical_NRN_sections.yaml" if with_NRN else None, - min_depth, - region_structure_boundary, - ) - args["skip_write"] = True - - synthesizer = SynthesizeMorphologies(**args) - res = synthesizer.synthesize() - - assert (res["x"] == 200).all() - assert (res["y"] == 200).all() - assert (res["z"] == 200).all() - assert res["name"].tolist() == [ - "e3e70682c2094cac629f6fbed82c07cd", - None, - "216363698b529b4a97b750923ceb3ffd", - None, - "14a03569d26b949692e5dfe8cb1855fe", - None, - "4462ebfc5f915ef09cfbac6e7687a66e", - None, - ] - print(res["apical_points"]) - assert_allclose(res["apical_points"][0], [[-38.110016, 152.74228, -20.12851]]) - assert res["apical_points"][1] is None - - assert_allclose(res["apical_points"][3], [[-3.50943, 115.28183, 5.8795013]]) - assert res["apical_points"][4] is None - - assert_allclose(res["apical_points"][6], [[2.271698, 114.72031, -5.102951]]) - assert res["apical_points"][7] is None - - assert_allclose(res["apical_points"][9], [[12.257843, 85.439545, -19.165298]]) - assert res["apical_points"][10] is None - - # Check that the morphologies were not written - res_files = tmpdir.listdir() - assert len(res_files) == 5 - assert sorted(i.basename for i in res_files) == [ - "apical.yaml", - "apical_NRN_sections.yaml", - "axon_morphs.tsv", - "input_cells.mvd3", - "test_cells.mvd3", - ] - - -def run_with_mpi(): - """Test morphology synthesis with MPI.""" - # pylint: disable=import-outside-toplevel, too-many-locals, import-error - from data_factories import generate_axon_morph_tsv - from data_factories import generate_cell_collection - from data_factories import generate_cells_df - from data_factories import generate_input_cells - from data_factories import generate_mesh - from data_factories import generate_region_structure_boundary - from data_factories import generate_small_O1 - from data_factories import input_cells_path - from mpi4py import MPI - - from region_grower.utils import setup_logger - - COMM = MPI.COMM_WORLD # pylint: disable=c-extension-no-member - rank = COMM.Get_rank() - MASTER_RANK = 0 - is_master = rank == MASTER_RANK - - tmp_folder = Path("/tmp/test-run-synthesis_" + str(uuid4())) - tmp_folder = COMM.bcast(tmp_folder, root=MASTER_RANK) - input_cells = input_cells_path(tmp_folder) - small_O1_path = str(tmp_folder / "atlas") - region_structure_path = str(tmp_folder / "region_structure_boundary.yaml") - args = create_args( - True, - tmp_folder, - input_cells, - small_O1_path, - tmp_folder / "axon_morphs.tsv", - "apical_NRN_sections.yaml", - min_depth=25, - region_structure_boundary=region_structure_path, - ) - - setup_logger("info") - logging.getLogger("distributed").setLevel(logging.ERROR) - - if is_master: - tmp_folder.mkdir(exist_ok=True) - print(f"============= #{rank}: Create data") - cells_raw_data = generate_cells_df() - cell_collection = generate_cell_collection(cells_raw_data) - generate_input_cells(cell_collection, tmp_folder) - generate_small_O1(small_O1_path) - atlas = {"atlas": small_O1_path, "structure": DATA / "region_structure.yaml"} - generate_mesh(atlas, tmp_folder / "mesh.obj") - region_structure_path = generate_region_structure_boundary( - DATA / "region_structure.yaml", region_structure_path, str(tmp_folder / "mesh.obj") - ) - - generate_axon_morph_tsv(tmp_folder) - for dest in range(1, COMM.Get_size()): - req = COMM.isend("done", dest=dest) - else: - req = COMM.irecv(source=0) - req.wait() - - synthesizer = SynthesizeMorphologies(**args) - try: - print(f"============= #{rank}: Start synthesize") - synthesizer.synthesize() - - # Check results - print(f"============= #{rank}: Checking results") - expected_size = 12 - assert len(list(iter_morphology_files(tmp_folder))) == expected_size - check_yaml(DATA / "apical_boundary.yaml", args["out_apical"]) - finally: - # Clean the directory - print(f"============= #{rank}: Cleaning") - shutil.rmtree(tmp_folder) - - -if __name__ == "__main__": # pragma: no cover - run_with_mpi() diff --git a/tox.ini b/tox.ini index 5669993..ced097a 100644 --- a/tox.ini +++ b/tox.ini @@ -51,7 +51,7 @@ allowlist_externals = mpiexec commands = mpiexec -n 4 {envbindir}/python tests/test_synthesize_morphologies.py - mpiexec -n 4 {envbindir}/python tests/test_synthesize_morphologies_boundary.py + mpiexec -n 4 {envbindir}/python tests/test_synthesize_morphologies.py boundary [testenv:coverage] skip_install = true From ddc8eb1f7584298deb493f7df2834619aba0895d Mon Sep 17 00:00:00 2001 From: Adrien Berchet Date: Mon, 15 May 2023 19:26:02 +0000 Subject: [PATCH 06/63] Apply 1 suggestion(s) to 1 file(s) --- region_grower/context.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/region_grower/context.py b/region_grower/context.py index eed8bf9..ae75502 100644 --- a/region_grower/context.py +++ b/region_grower/context.py @@ -86,7 +86,7 @@ class SpaceContext: layer_depths: List cortical_depths: List boundaries: Dict = None - atlas_info: Dict = {} # voxel_dimnesions, offset and shape from atlas for indices conversion + atlas_info: Dict = {} # voxel_dimensions, offset and shape from atlas for indices conversion def indices_to_positions(self, indices): """Local version of voxcel's function to prevent atlas loading.""" From 21e984c0aa743cac9b5310d7f61d226c8c34e67b Mon Sep 17 00:00:00 2001 From: Adrien Berchet Date: Mon, 15 May 2023 19:28:15 +0000 Subject: [PATCH 07/63] Apply 1 suggestion(s) to 1 file(s) --- region_grower/atlas_helper.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/region_grower/atlas_helper.py b/region_grower/atlas_helper.py index 98c37fb..cdac2af 100644 --- a/region_grower/atlas_helper.py +++ b/region_grower/atlas_helper.py @@ -32,7 +32,7 @@ def __init__(self, atlas: Atlas, region_structure_path: str): with open(region_structure_path, "r", encoding="utf-8") as region_file: self.region_structure = yaml.safe_load(region_file) else: - raise ValueError(f"region_structure file at {region_structure_path} not found.") + raise ValueError(f"region_structure file not found at {region_structure_path}.") self.regions = list(self.region_structure.keys()) self.layers = {} From 38ac8af7cb1cb555357b935971d57d0caf074f0d Mon Sep 17 00:00:00 2001 From: arnaudon Date: Tue, 16 May 2023 13:34:39 +0200 Subject: [PATCH 08/63] lint --- tests/test_synthesize_morphologies.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/tests/test_synthesize_morphologies.py b/tests/test_synthesize_morphologies.py index 6e12098..e41fdd6 100644 --- a/tests/test_synthesize_morphologies.py +++ b/tests/test_synthesize_morphologies.py @@ -1,10 +1,10 @@ """Test the region_grower.synthesize_morphologies module.""" # pylint: disable=missing-function-docstring import json -import sys import logging import os import shutil +import sys from copy import deepcopy from itertools import combinations from pathlib import Path @@ -27,7 +27,6 @@ from region_grower.synthesize_morphologies import RegionMapper from region_grower.synthesize_morphologies import SynthesizeMorphologies - DATA = Path(__file__).parent / "data" From a00b4ba2bab536e7c2566b7bbfd1b621b7e04b29 Mon Sep 17 00:00:00 2001 From: arnaudon Date: Thu, 22 Jun 2023 16:32:13 +0200 Subject: [PATCH 09/63] more --- region_grower/atlas_helper.py | 1 + region_grower/context.py | 52 ++++++++++++++---------- region_grower/synthesize_morphologies.py | 10 +++-- 3 files changed, 38 insertions(+), 25 deletions(-) diff --git a/region_grower/atlas_helper.py b/region_grower/atlas_helper.py index cdac2af..63ddbb1 100644 --- a/region_grower/atlas_helper.py +++ b/region_grower/atlas_helper.py @@ -31,6 +31,7 @@ def __init__(self, atlas: Atlas, region_structure_path: str): if region_structure_path is not None and Path(region_structure_path).exists(): with open(region_structure_path, "r", encoding="utf-8") as region_file: self.region_structure = yaml.safe_load(region_file) + self.region_structure_base_path = Path(region_structure_path).parent else: raise ValueError(f"region_structure file not found at {region_structure_path}.") diff --git a/region_grower/context.py b/region_grower/context.py index ae75502..825a928 100644 --- a/region_grower/context.py +++ b/region_grower/context.py @@ -85,7 +85,7 @@ class SpaceContext: layer_depths: List cortical_depths: List - boundaries: Dict = None + boundaries: List = None atlas_info: Dict = {} # voxel_dimensions, offset and shape from atlas for indices conversion def indices_to_positions(self, indices): @@ -105,16 +105,10 @@ def positions_to_indices(self, positions): result[result >= self.atlas_info["shape"]] = -1 return result - def get_boundaries(self): + def get_boundaries(self, mtype, neurite_types): """Returns a dict with boundaries data for NeuroTS.""" - # add here necessary logic to convert raw config data to NeuroTS context data - self.boundaries = json.loads(self.boundaries) - mesh = trimesh.load_mesh(self.boundaries["path"]) - - # to save points for plotting while growing to debug the prob function - debug = False - def get_distance_to_mesh(mesh, ray_origin, ray_direction): + def get_distance_to_mesh(mesh, ray_origin, ray_direction, debug=False): """Compute distances from point/directions to a mesh.""" vox_ray_origin = self.positions_to_indices(ray_origin) locations = mesh.ray.intersects_location( @@ -146,23 +140,35 @@ def _prob(dist, params): # TODO: other variants for this function could be included here return np.clip(p, 0, 1) - if "params_section" in self.boundaries: + # add here necessary logic to convert raw config data to NeuroTS context data + self.boundaries = json.loads(self.boundaries) + for i, boundary in enumerate(self.boundaries): + + mtypes = boundary.pop('mtypes', None) + if mtypes is not None and mtype not in mtypes: + continue + # TODO: add check a similar check on neurite_types + + mesh = trimesh.load_mesh(boundary["path"]) + + if "params_section" in boundary: - def section_prob(direction, current_point): - """Probability function for the sections.""" - dist = get_distance_to_mesh(mesh, current_point, direction) - return _prob(dist, self.boundaries["params_section"]) + def section_prob(direction, current_point): + """Probability function for the sections.""" + dist = get_distance_to_mesh(mesh, current_point, direction) + return _prob(dist, boundary["params_section"]) - self.boundaries["section_prob"] = section_prob + self.boundaries[i]["section_prob"] = section_prob - if "params_trunk" in self.boundaries: + if "params_trunk" in boundary: - def trunk_prob(direction, current_point): - """Probability function for the trunks.""" - dist = get_distance_to_mesh(mesh, current_point, direction) - return _prob(dist, self.boundaries["params_trunk"]) + def trunk_prob(direction, current_point): + """Probability function for the trunks.""" + dist = get_distance_to_mesh(mesh, current_point, direction) + return _prob(dist, boundary["params_trunk"]) + + self.boundaries[i]["trunk_prob"] = trunk_prob - self.boundaries["trunk_prob"] = trunk_prob return self.boundaries def layer_fraction_to_position(self, layer: int, layer_fraction: float) -> float: @@ -406,7 +412,9 @@ def _synthesize_once(self, rng) -> SynthesisResult: context = {"debug_data": self.debug_infos["input_scaling"]} if self.context.boundaries is not None: - context.update(self.context.get_boundaries()) + context["boundaries"] = self.context.get_boundaries( + mtype=self.cell.mtype, neurite_types=params["grow_types"] + ) grower = NeuronGrower( input_parameters=params, diff --git a/region_grower/synthesize_morphologies.py b/region_grower/synthesize_morphologies.py index 4dce458..9461382 100644 --- a/region_grower/synthesize_morphologies.py +++ b/region_grower/synthesize_morphologies.py @@ -523,9 +523,13 @@ def assign_atlas_data(self, min_depth=25, max_depth=5000): ) if "boundaries" in self.atlas.region_structure[_region]: - self.cells_data.loc[region_mask, "boundaries"] = json.dumps( - self.atlas.region_structure[_region]["boundaries"] - ) + boundaries = self.atlas.region_structure[_region]["boundaries"] + for boundary in boundaries: + if not Path(boundary["path"]).is_absolute(): + boundary["path"] = str( + self.atlas.region_structure_base_path / boundary["path"] + ) + self.cells_data.loc[region_mask, "boundaries"] = json.dumps(boundaries) self.cells_data.loc[region_mask, "atlas_info"] = json.dumps( { "voxel_dimensions": self.atlas.depths[_region].voxel_dimensions.tolist(), From dad096bf99c9808b61fe5735a8dcea74e9b6f536 Mon Sep 17 00:00:00 2001 From: arnaudon Date: Thu, 22 Jun 2023 21:27:11 +0200 Subject: [PATCH 10/63] small fix --- region_grower/synthesize_morphologies.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/region_grower/synthesize_morphologies.py b/region_grower/synthesize_morphologies.py index 9461382..9330e33 100644 --- a/region_grower/synthesize_morphologies.py +++ b/region_grower/synthesize_morphologies.py @@ -527,7 +527,7 @@ def assign_atlas_data(self, min_depth=25, max_depth=5000): for boundary in boundaries: if not Path(boundary["path"]).is_absolute(): boundary["path"] = str( - self.atlas.region_structure_base_path / boundary["path"] + (self.atlas.region_structure_base_path / boundary["path"]).absolute() ) self.cells_data.loc[region_mask, "boundaries"] = json.dumps(boundaries) self.cells_data.loc[region_mask, "atlas_info"] = json.dumps( From 321b5513d41dd01351e49517397904c7c5c29942 Mon Sep 17 00:00:00 2001 From: arnaudon Date: Fri, 30 Jun 2023 11:51:33 +0200 Subject: [PATCH 11/63] more --- region_grower/context.py | 29 +++++++++++++++++++++-------- 1 file changed, 21 insertions(+), 8 deletions(-) diff --git a/region_grower/context.py b/region_grower/context.py index 825a928..a436685 100644 --- a/region_grower/context.py +++ b/region_grower/context.py @@ -108,14 +108,22 @@ def positions_to_indices(self, positions): def get_boundaries(self, mtype, neurite_types): """Returns a dict with boundaries data for NeuroTS.""" - def get_distance_to_mesh(mesh, ray_origin, ray_direction, debug=False): + def get_distance_to_mesh(mesh, ray_origin, ray_direction, mesh_type): """Compute distances from point/directions to a mesh.""" - vox_ray_origin = self.positions_to_indices(ray_origin) + debug = os.environ.get("REGION_GROWER_BOUNDARY_DEBUG", 0) + + if mesh_type == "voxel": + ray_origin = self.positions_to_indices(ray_origin) + locations = mesh.ray.intersects_location( - ray_origins=[vox_ray_origin], ray_directions=[ray_direction] + ray_origins=[ray_origin], ray_directions=[ray_direction] )[0] + if len(locations): - intersect = self.indices_to_positions(locations[0]) + intersect = locations[0] + if mesh_type == "voxel": + intersect = self.indices_to_positions(intersect) + dist = np.linalg.norm(intersect - ray_origin) if debug: # pragma: no cover @@ -124,11 +132,13 @@ def get_distance_to_mesh(mesh, ray_origin, ray_direction, debug=False): with open("data.csv", "a", encoding="utf8") as file: print(f"{x}, {y}, {z}, {vx}, {vy}, {vz}, {dist}", file=file) return dist + if debug: # pragma: no cover x, y, z = ray_origin vx, vy, vz = ray_direction with open("data.csv", "a", encoding="utf8") as file: print(f"{x}, {y}, {z}, {vx}, {vy}, {vz}, {-100}", file=file) + return np.inf # pragma: no cover def _prob(dist, params): @@ -143,8 +153,7 @@ def _prob(dist, params): # add here necessary logic to convert raw config data to NeuroTS context data self.boundaries = json.loads(self.boundaries) for i, boundary in enumerate(self.boundaries): - - mtypes = boundary.pop('mtypes', None) + mtypes = boundary.pop("mtypes", None) if mtypes is not None and mtype not in mtypes: continue # TODO: add check a similar check on neurite_types @@ -155,7 +164,9 @@ def _prob(dist, params): def section_prob(direction, current_point): """Probability function for the sections.""" - dist = get_distance_to_mesh(mesh, current_point, direction) + dist = get_distance_to_mesh( + mesh, current_point, direction, mesh_type=boundary.get("mesh_type", "voxel") + ) return _prob(dist, boundary["params_section"]) self.boundaries[i]["section_prob"] = section_prob @@ -164,7 +175,9 @@ def section_prob(direction, current_point): def trunk_prob(direction, current_point): """Probability function for the trunks.""" - dist = get_distance_to_mesh(mesh, current_point, direction) + dist = get_distance_to_mesh( + mesh, current_point, direction, mesh_type=boundary.get("mesh_type", "voxel") + ) return _prob(dist, boundary["params_trunk"]) self.boundaries[i]["trunk_prob"] = trunk_prob From ff089d99a157387d4554dbe01e1cb719ec4ee48e Mon Sep 17 00:00:00 2001 From: arnaudon Date: Tue, 5 Sep 2023 16:14:52 +0200 Subject: [PATCH 12/63] more --- region_grower/context.py | 48 +++++++++++++++++++----- region_grower/synthesize_morphologies.py | 5 +++ 2 files changed, 43 insertions(+), 10 deletions(-) diff --git a/region_grower/context.py b/region_grower/context.py index a436685..42fafe3 100644 --- a/region_grower/context.py +++ b/region_grower/context.py @@ -86,6 +86,7 @@ class SpaceContext: layer_depths: List cortical_depths: List boundaries: List = None + directions: Dict = None atlas_info: Dict = {} # voxel_dimensions, offset and shape from atlas for indices conversion def indices_to_positions(self, indices): @@ -105,7 +106,25 @@ def positions_to_indices(self, positions): result[result >= self.atlas_info["shape"]] = -1 return result - def get_boundaries(self, mtype, neurite_types): + def get_direction(self, mtype): + for i, direction in enumerate(self.directions): + target_direction = direction["params"]["direction"] + strength = direction["params"]["strength"] + + def section_prob(seg_direction, current_point): + """Probability function for the sections.""" + p = np.clip( + 1.0 - np.arccos(seg_direction.dot(target_direction)) * strength / np.pi, + 0, + 1, + ) + return p + + self.directions[i]["section_prob"] = section_prob + + return self.directions + + def get_boundaries(self, mtype): """Returns a dict with boundaries data for NeuroTS.""" def get_distance_to_mesh(mesh, ray_origin, ray_direction, mesh_type): @@ -151,12 +170,11 @@ def _prob(dist, params): return np.clip(p, 0, 1) # add here necessary logic to convert raw config data to NeuroTS context data - self.boundaries = json.loads(self.boundaries) - for i, boundary in enumerate(self.boundaries): + boundaries = json.loads(self.boundaries) + for i, boundary in enumerate(boundaries): mtypes = boundary.pop("mtypes", None) if mtypes is not None and mtype not in mtypes: continue - # TODO: add check a similar check on neurite_types mesh = trimesh.load_mesh(boundary["path"]) @@ -169,7 +187,7 @@ def section_prob(direction, current_point): ) return _prob(dist, boundary["params_section"]) - self.boundaries[i]["section_prob"] = section_prob + boundaries[i]["section_prob"] = section_prob if "params_trunk" in boundary: @@ -180,9 +198,9 @@ def trunk_prob(direction, current_point): ) return _prob(dist, boundary["params_trunk"]) - self.boundaries[i]["trunk_prob"] = trunk_prob + boundaries[i]["trunk_prob"] = trunk_prob - return self.boundaries + return boundaries def layer_fraction_to_position(self, layer: int, layer_fraction: float) -> float: """Returns an absolute position from a layer and a fraction of the layer. @@ -282,6 +300,12 @@ def __init__( def _correct_position_orientation_scaling(self, params_orig: Dict) -> Dict: """Return a copy of the parameter with the correct orientation and scaling.""" params = deepcopy(params_orig) + self.context.directions = json.loads(self.context.directions) + if self.context.directions is not None: + for i, direction in enumerate(self.context.directions): + self.context.directions[i]["params"]["direction"] = self.cell.lookup_orientation( + direction["params"]["direction"] + ) for neurite_type in params["grow_types"]: if isinstance(params[neurite_type]["orientation"], dict): @@ -424,10 +448,14 @@ def _synthesize_once(self, rng) -> SynthesisResult: axon_morph = None context = {"debug_data": self.debug_infos["input_scaling"]} + if self.context.directions is not None: + if "constraints" not in context: + context["constraints"] = [] + context["constraints"] += self.context.get_direction(mtype=self.cell.mtype) if self.context.boundaries is not None: - context["boundaries"] = self.context.get_boundaries( - mtype=self.cell.mtype, neurite_types=params["grow_types"] - ) + if "constraints" not in context: + context["constraints"] = [] + context["constraints"] += self.context.get_boundaries(mtype=self.cell.mtype) grower = NeuronGrower( input_parameters=params, diff --git a/region_grower/synthesize_morphologies.py b/region_grower/synthesize_morphologies.py index 9330e33..28e01b7 100644 --- a/region_grower/synthesize_morphologies.py +++ b/region_grower/synthesize_morphologies.py @@ -126,6 +126,7 @@ def _parallel_wrapper( current_space_context = SpaceContext( layer_depths=row["layer_depths"], cortical_depths=row["cortical_depths"], + directions=row["directions"] if "directions" in row else None, boundaries=row["boundaries"] if "boundaries" in row else None, atlas_info=json.loads(row["atlas_info"]) if "atlas_info" in row else None, ) @@ -522,6 +523,10 @@ def assign_atlas_data(self, min_depth=25, max_depth=5000): dtype=object, ) + if "directions" in self.atlas.region_structure[_region]: + self.cells_data.loc[region_mask, "directions"] = json.dumps( + self.atlas.region_structure[_region]["directions"] + ) if "boundaries" in self.atlas.region_structure[_region]: boundaries = self.atlas.region_structure[_region]["boundaries"] for boundary in boundaries: From 38ac87161134294748e2e8c676ac010536e6a983 Mon Sep 17 00:00:00 2001 From: arnaudon Date: Thu, 21 Sep 2023 15:23:17 +0200 Subject: [PATCH 13/63] more --- region_grower/context.py | 3 ++- region_grower/synthesize_morphologies.py | 6 +++--- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/region_grower/context.py b/region_grower/context.py index 42fafe3..14ffb97 100644 --- a/region_grower/context.py +++ b/region_grower/context.py @@ -142,6 +142,7 @@ def get_distance_to_mesh(mesh, ray_origin, ray_direction, mesh_type): intersect = locations[0] if mesh_type == "voxel": intersect = self.indices_to_positions(intersect) + ray_origin = self.indices_to_positions(ray_origin) dist = np.linalg.norm(intersect - ray_origin) @@ -300,8 +301,8 @@ def __init__( def _correct_position_orientation_scaling(self, params_orig: Dict) -> Dict: """Return a copy of the parameter with the correct orientation and scaling.""" params = deepcopy(params_orig) - self.context.directions = json.loads(self.context.directions) if self.context.directions is not None: + self.context.directions = json.loads(self.context.directions) for i, direction in enumerate(self.context.directions): self.context.directions[i]["params"]["direction"] = self.cell.lookup_orientation( direction["params"]["direction"] diff --git a/region_grower/synthesize_morphologies.py b/region_grower/synthesize_morphologies.py index 28e01b7..0eaa44f 100644 --- a/region_grower/synthesize_morphologies.py +++ b/region_grower/synthesize_morphologies.py @@ -537,9 +537,9 @@ def assign_atlas_data(self, min_depth=25, max_depth=5000): self.cells_data.loc[region_mask, "boundaries"] = json.dumps(boundaries) self.cells_data.loc[region_mask, "atlas_info"] = json.dumps( { - "voxel_dimensions": self.atlas.depths[_region].voxel_dimensions.tolist(), - "offset": self.atlas.depths[_region].offset.tolist(), - "shape": self.atlas.depths[_region].shape, + "voxel_dimensions": self.atlas.brain_regions.voxel_dimensions.tolist(), + "offset": self.atlas.brain_regions.offset.tolist(), + "shape": self.atlas.brain_regions.shape, } ) From 8b01201929a0beb35b9c83ba5fa9524b17b21fb7 Mon Sep 17 00:00:00 2001 From: arnaudon Date: Mon, 25 Sep 2023 15:32:20 +0200 Subject: [PATCH 14/63] more --- region_grower/context.py | 15 +++++++++------ region_grower/modify.py | 2 +- 2 files changed, 10 insertions(+), 7 deletions(-) diff --git a/region_grower/context.py b/region_grower/context.py index 14ffb97..6f3418e 100644 --- a/region_grower/context.py +++ b/region_grower/context.py @@ -110,15 +110,17 @@ def get_direction(self, mtype): for i, direction in enumerate(self.directions): target_direction = direction["params"]["direction"] strength = direction["params"]["strength"] + mode = direction["params"]["mode"] def section_prob(seg_direction, current_point): """Probability function for the sections.""" - p = np.clip( - 1.0 - np.arccos(seg_direction.dot(target_direction)) * strength / np.pi, - 0, - 1, - ) - return p + if mode == "parallel": + p = 1.0 - np.arccos(seg_direction.dot(target_direction)) * strength / np.pi + if mode == "perpendicular": + p = 1.0 - strength * abs( + np.arccos(seg_direction.dot(target_direction)) / np.pi * 2.0 - 1.0 + ) + return np.clip(p, 0, 1) self.directions[i]["section_prob"] = section_prob @@ -402,6 +404,7 @@ def synthesize(self) -> SynthesisResult: """ rng = np.random.default_rng(self.params.seed) + return self._synthesize_once(rng) for _ in range(self.internals.retries): try: return self._synthesize_once(rng) diff --git a/region_grower/modify.py b/region_grower/modify.py index 4bf42a0..618c786 100644 --- a/region_grower/modify.py +++ b/region_grower/modify.py @@ -134,7 +134,7 @@ def input_scaling( "with_debug_info": debug_info is not None, }, } - else: + elif neurite_type != "axon": if debug_info is not None: debug_info.update( { From 7eafd856905994cd6c56f3de5cf6be7acbcdf24b Mon Sep 17 00:00:00 2001 From: arnaudon Date: Wed, 27 Sep 2023 15:51:00 +0200 Subject: [PATCH 15/63] more --- region_grower/context.py | 60 ++++++++++++++++-------- region_grower/synthesize_morphologies.py | 2 + 2 files changed, 43 insertions(+), 19 deletions(-) diff --git a/region_grower/context.py b/region_grower/context.py index 6f3418e..563ddc6 100644 --- a/region_grower/context.py +++ b/region_grower/context.py @@ -10,6 +10,7 @@ import json from collections import defaultdict from copy import deepcopy +from functools import partial from typing import Dict from typing import List from typing import Optional @@ -85,6 +86,8 @@ class SpaceContext: layer_depths: List cortical_depths: List + soma_position: float + soma_depth: float boundaries: List = None directions: Dict = None atlas_info: Dict = {} # voxel_dimensions, offset and shape from atlas for indices conversion @@ -107,13 +110,30 @@ def positions_to_indices(self, positions): return result def get_direction(self, mtype): + directions = [] for i, direction in enumerate(self.directions): - target_direction = direction["params"]["direction"] - strength = direction["params"]["strength"] - mode = direction["params"]["mode"] + mtypes = direction.pop("mtypes", None) + if mtypes is not None and mtype not in mtypes: + continue - def section_prob(seg_direction, current_point): + def section_prob(seg_direction, current_point, direction=None): """Probability function for the sections.""" + + target_direction = direction["params"]["direction"] + strength = direction["params"]["strength"] + mode = direction["params"]["mode"] + layers = direction["params"].get("layers", []) + if layers: + depth = self.soma_depth + (self.soma_position - current_point).dot( + direction["pia_direction"] + ) + in_layer = False + for layer in layers: + if self.layer_depths[layer - 1] < depth < self.layer_depths[layer]: + in_layer = True + if not in_layer: + return 1.0 + if mode == "parallel": p = 1.0 - np.arccos(seg_direction.dot(target_direction)) * strength / np.pi if mode == "perpendicular": @@ -122,9 +142,9 @@ def section_prob(seg_direction, current_point): ) return np.clip(p, 0, 1) - self.directions[i]["section_prob"] = section_prob - - return self.directions + direction["section_prob"] = partial(section_prob, direction=direction) + directions.append(direction) + return directions def get_boundaries(self, mtype): """Returns a dict with boundaries data for NeuroTS.""" @@ -173,35 +193,34 @@ def _prob(dist, params): return np.clip(p, 0, 1) # add here necessary logic to convert raw config data to NeuroTS context data - boundaries = json.loads(self.boundaries) - for i, boundary in enumerate(boundaries): + all_boundaries = json.loads(self.boundaries) + boundaries = [] + for i, boundary in enumerate(all_boundaries): mtypes = boundary.pop("mtypes", None) if mtypes is not None and mtype not in mtypes: continue mesh = trimesh.load_mesh(boundary["path"]) + mesh_type = boundary.get("mesh_type", "voxel") if "params_section" in boundary: - def section_prob(direction, current_point): + def section_prob(direction, current_point, mesh=None, mesh_type=None): """Probability function for the sections.""" - dist = get_distance_to_mesh( - mesh, current_point, direction, mesh_type=boundary.get("mesh_type", "voxel") - ) + dist = get_distance_to_mesh(mesh, current_point, direction, mesh_type=mesh_type) return _prob(dist, boundary["params_section"]) - boundaries[i]["section_prob"] = section_prob + boundary["section_prob"] = partial(section_prob, mesh=mesh, mesh_type=mesh_type) if "params_trunk" in boundary: - def trunk_prob(direction, current_point): + def trunk_prob(direction, current_point, mesh=None, mesh_type=None): """Probability function for the trunks.""" - dist = get_distance_to_mesh( - mesh, current_point, direction, mesh_type=boundary.get("mesh_type", "voxel") - ) + dist = get_distance_to_mesh(mesh, current_point, direction, mesh_type) return _prob(dist, boundary["params_trunk"]) - boundaries[i]["trunk_prob"] = trunk_prob + boundary["trunk_prob"] = partial(trunk_prob, mesh=mesh, mesh_type=mesh_type) + boundaries.append(boundary) return boundaries @@ -309,6 +328,9 @@ def _correct_position_orientation_scaling(self, params_orig: Dict) -> Dict: self.context.directions[i]["params"]["direction"] = self.cell.lookup_orientation( direction["params"]["direction"] ) + self.context.directions[i]["pia_direction"] = self.cell.lookup_orientation( + PIA_DIRECTION + ).tolist() for neurite_type in params["grow_types"]: if isinstance(params[neurite_type]["orientation"], dict): diff --git a/region_grower/synthesize_morphologies.py b/region_grower/synthesize_morphologies.py index 0eaa44f..63b9713 100644 --- a/region_grower/synthesize_morphologies.py +++ b/region_grower/synthesize_morphologies.py @@ -129,6 +129,8 @@ def _parallel_wrapper( directions=row["directions"] if "directions" in row else None, boundaries=row["boundaries"] if "boundaries" in row else None, atlas_info=json.loads(row["atlas_info"]) if "atlas_info" in row else None, + soma_position=current_cell.position, + soma_depth=row["current_depth"], ) axon_scale = row.get("axon_scale", None) From 639991789a74ebb4d6dd0fdc739be6077e94184b Mon Sep 17 00:00:00 2001 From: arnaudon Date: Tue, 3 Oct 2023 14:49:00 +0200 Subject: [PATCH 16/63] more --- region_grower/context.py | 82 ++++++++++++++++++++++++++++------------ 1 file changed, 57 insertions(+), 25 deletions(-) diff --git a/region_grower/context.py b/region_grower/context.py index 563ddc6..2f5c04b 100644 --- a/region_grower/context.py +++ b/region_grower/context.py @@ -10,6 +10,7 @@ import json from collections import defaultdict from copy import deepcopy +from pathlib import Path from functools import partial from typing import Dict from typing import List @@ -120,7 +121,7 @@ def section_prob(seg_direction, current_point, direction=None): """Probability function for the sections.""" target_direction = direction["params"]["direction"] - strength = direction["params"]["strength"] + power = direction["params"].get("power", 1.0) mode = direction["params"]["mode"] layers = direction["params"].get("layers", []) if layers: @@ -135,11 +136,12 @@ def section_prob(seg_direction, current_point, direction=None): return 1.0 if mode == "parallel": - p = 1.0 - np.arccos(seg_direction.dot(target_direction)) * strength / np.pi + p = (1.0 - np.arccos(seg_direction.dot(target_direction)) / np.pi) ** power if mode == "perpendicular": - p = 1.0 - strength * abs( - np.arccos(seg_direction.dot(target_direction)) / np.pi * 2.0 - 1.0 - ) + p = ( + 1.0 + - abs(np.arccos(seg_direction.dot(target_direction)) / np.pi * 2.0 - 1.0) + ) ** power return np.clip(p, 0, 1) direction["section_prob"] = partial(section_prob, direction=direction) @@ -183,15 +185,6 @@ def get_distance_to_mesh(mesh, ray_origin, ray_direction, mesh_type): return np.inf # pragma: no cover - def _prob(dist, params): - """Probability function from distance to boundary.""" - p = (dist - params.get("d_min", 0)) / ( - params.get("d_max", 100) - params.get("d_min", 0) - ) ** params.get("power", 1.5) - - # TODO: other variants for this function could be included here - return np.clip(p, 0, 1) - # add here necessary logic to convert raw config data to NeuroTS context data all_boundaries = json.loads(self.boundaries) boundaries = [] @@ -200,26 +193,65 @@ def _prob(dist, params): if mtypes is not None and mtype not in mtypes: continue - mesh = trimesh.load_mesh(boundary["path"]) + meshes = [] + if Path(boundary["path"]).is_dir(): + for mesh_path in Path(boundary["path"]).iterdir(): + meshes.append(trimesh.load_mesh(mesh_path)) + distances = [ + trimesh.proximity.closest_point(mesh, [self.soma_position])[1][0] + for mesh in meshes + ] + mesh = meshes[np.argmin(distances)] + + else: + mesh = trimesh.load_mesh(boundary["path"]) + mesh_type = boundary.get("mesh_type", "voxel") + mode = boundary.get("mode", "repulsive") - if "params_section" in boundary: + if mode == "repulsive": - def section_prob(direction, current_point, mesh=None, mesh_type=None): + def prob(direction, current_point, mesh=None, mesh_type=None, params=None): """Probability function for the sections.""" dist = get_distance_to_mesh(mesh, current_point, direction, mesh_type=mesh_type) - return _prob(dist, boundary["params_section"]) + p = (dist - params.get("d_min", 0)) / ( + params.get("d_max", 100) - params.get("d_min", 0) + ) ** params.get("power", 1.5) + return np.clip(p, 0, 1) - boundary["section_prob"] = partial(section_prob, mesh=mesh, mesh_type=mesh_type) + elif mode == "attractive": - if "params_trunk" in boundary: + def prob(direction, current_point, mesh=None, mesh_type=None, params=None): + """Probability function for the sections.""" + closest_point = trimesh.proximity.closest_point(mesh, [current_point])[0][0] + + if mesh_type == "voxel": + closest_point = self.indices_to_positions(closest_point) + + mesh_direction = closest_point - current_point + distance = np.linalg.norm(mesh_direction) + if params.get("d_min", 0) < distance < params.get("d_max", 100): + if mesh_type == "voxel": + current_point = self.positions_to_indices(current_point) + p = ( + 1 - np.arccos(direction.dot(mesh_direction) / distance) / np.pi + ) ** params.get("power", 1.0) + return np.clip(p, 0, 1) + return 1.0 + + else: + raise ValueError(f"boundary mode {mode} not understood!") + + if "params_section" in boundary: + boundary["section_prob"] = partial( + prob, mesh=mesh, mesh_type=mesh_type, params=boundary["params_section"] + ) - def trunk_prob(direction, current_point, mesh=None, mesh_type=None): - """Probability function for the trunks.""" - dist = get_distance_to_mesh(mesh, current_point, direction, mesh_type) - return _prob(dist, boundary["params_trunk"]) + if "params_trunk" in boundary: + boundary["trunk_prob"] = partial( + prob, mesh=mesh, mesh_type=mesh_type, params=boundary["params_ trunk"] + ) - boundary["trunk_prob"] = partial(trunk_prob, mesh=mesh, mesh_type=mesh_type) boundaries.append(boundary) return boundaries From 030c961b7834308cc60cdee837197df5d519d796 Mon Sep 17 00:00:00 2001 From: arnaudon Date: Tue, 10 Oct 2023 11:38:00 +0200 Subject: [PATCH 17/63] option to synthesize axons --- region_grower/cli.py | 6 ++++++ region_grower/synthesize_morphologies.py | 11 ++++++++++- 2 files changed, 16 insertions(+), 1 deletion(-) diff --git a/region_grower/cli.py b/region_grower/cli.py index 3c0f93f..b7842fe 100644 --- a/region_grower/cli.py +++ b/region_grower/cli.py @@ -284,6 +284,12 @@ def generate_distributions( is_flag=True, default=False, ) +@click.option( + "--synthesize-axons", + help="Display the versions of all the accessible modules in a logger entry", + is_flag=True, + default=False, +) def synthesize_morphologies(**kwargs): # pylint: disable=too-many-arguments, too-many-locals """Synthesize morphologies.""" if mpi_enabled and kwargs.get("with_mpi", False): # pragma: no cover diff --git a/region_grower/synthesize_morphologies.py b/region_grower/synthesize_morphologies.py index 63b9713..873c088 100644 --- a/region_grower/synthesize_morphologies.py +++ b/region_grower/synthesize_morphologies.py @@ -201,6 +201,7 @@ class SynthesizeMorphologies: should be used to graft the axon on each synthesized morphology. base_morph_dir: the path containing the morphologies listed in the TSV file given in ``morph_axon``. + synthesize_axons: set to True to synthesize axons instead of grafting atlas_cache: the path to the directory used for the atlas cache. seed: the starting seed to use (note that the GID of each cell is added to this seed to ensure all cells has different seeds). @@ -254,6 +255,7 @@ def __init__( hide_progress_bar=False, dask_config=None, chunksize=None, + synthesize_axons=False, ): # pylint: disable=too-many-arguments, too-many-locals, too-many-statements self.seed = seed self.scaling_jitter_std = scaling_jitter_std @@ -303,6 +305,13 @@ def __init__( with open(tmd_distributions, "r", encoding="utf-8") as f: self.tmd_distributions = convert_from_legacy_neurite_type(json.load(f)) + for region, params in self.tmd_parameters.items(): + for param in params.values(): + if synthesize_axons: + assert "axon" in param["grow_types"] + elif "axon" in param["grow_types"]: + param["grow_types"].remove("axon") + # Set default values to tmd_parameters and tmd_distributions self.set_default_params_and_distrs() @@ -336,7 +345,7 @@ def __init__( LOGGER.info("Fetching atlas data from %s", atlas) self.assign_atlas_data(min_depth, max_depth) - if morph_axon is not None: + if morph_axon is not None and not synthesize_axons: LOGGER.info("Loading axon morphologies from %s", morph_axon) self.axon_morph_list = load_morphology_list(morph_axon, self.task_ids) check_na_morphologies( From 0f5e8d9154478cf1bf9f3f0a188ddac9bae4d926 Mon Sep 17 00:00:00 2001 From: arnaudon Date: Wed, 18 Oct 2023 15:44:54 +0200 Subject: [PATCH 18/63] add option for barrel --- region_grower/context.py | 54 +++++++++++++++++++++++++--------------- 1 file changed, 34 insertions(+), 20 deletions(-) diff --git a/region_grower/context.py b/region_grower/context.py index 2f5c04b..2291b07 100644 --- a/region_grower/context.py +++ b/region_grower/context.py @@ -193,26 +193,39 @@ def get_distance_to_mesh(mesh, ray_origin, ray_direction, mesh_type): if mtypes is not None and mtype not in mtypes: continue + mesh_type = boundary.get("mesh_type", "voxel") + mode = boundary.get("mode", "repulsive") + meshes = [] if Path(boundary["path"]).is_dir(): - for mesh_path in Path(boundary["path"]).iterdir(): - meshes.append(trimesh.load_mesh(mesh_path)) - distances = [ - trimesh.proximity.closest_point(mesh, [self.soma_position])[1][0] - for mesh in meshes - ] - mesh = meshes[np.argmin(distances)] + soma_position = self.soma_position + if mesh_type == "voxel": + soma_position = self.positions_to_indices(soma_position) + + if boundary.get("multimesh_mode", "closest") == "closest": + for mesh_path in Path(boundary["path"]).iterdir(): + meshes.append(trimesh.load_mesh(mesh_path)) + distances = [ + trimesh.proximity.closest_point(mesh, [soma_position])[1][0] + for mesh in meshes + ] + mesh = meshes[np.argmin(distances)] + + if boundary.get("multimesh_mode", "closest") == "inside": + mesh = None + for mesh_path in Path(boundary["path"]).iterdir(): + _mesh = trimesh.load_mesh(mesh_path) + if _mesh.contains([soma_position])[0]: + mesh = _mesh + break else: mesh = trimesh.load_mesh(boundary["path"]) - mesh_type = boundary.get("mesh_type", "voxel") - mode = boundary.get("mode", "repulsive") - if mode == "repulsive": def prob(direction, current_point, mesh=None, mesh_type=None, params=None): - """Probability function for the sections.""" + """Probability function for repulsive mode.""" dist = get_distance_to_mesh(mesh, current_point, direction, mesh_type=mesh_type) p = (dist - params.get("d_min", 0)) / ( params.get("d_max", 100) - params.get("d_min", 0) @@ -222,7 +235,7 @@ def prob(direction, current_point, mesh=None, mesh_type=None, params=None): elif mode == "attractive": def prob(direction, current_point, mesh=None, mesh_type=None, params=None): - """Probability function for the sections.""" + """Probability function for attractive mode.""" closest_point = trimesh.proximity.closest_point(mesh, [current_point])[0][0] if mesh_type == "voxel": @@ -242,15 +255,16 @@ def prob(direction, current_point, mesh=None, mesh_type=None, params=None): else: raise ValueError(f"boundary mode {mode} not understood!") - if "params_section" in boundary: - boundary["section_prob"] = partial( - prob, mesh=mesh, mesh_type=mesh_type, params=boundary["params_section"] - ) + if mesh is not None: + if "params_section" in boundary: + boundary["section_prob"] = partial( + prob, mesh=mesh, mesh_type=mesh_type, params=boundary["params_section"] + ) - if "params_trunk" in boundary: - boundary["trunk_prob"] = partial( - prob, mesh=mesh, mesh_type=mesh_type, params=boundary["params_ trunk"] - ) + if "params_trunk" in boundary: + boundary["trunk_prob"] = partial( + prob, mesh=mesh, mesh_type=mesh_type, params=boundary["params_trunk"] + ) boundaries.append(boundary) From 777ac401f75e3d473011c4b8af6a390fd60c2628 Mon Sep 17 00:00:00 2001 From: arnaudon Date: Thu, 9 Nov 2023 15:24:41 +0100 Subject: [PATCH 19/63] more --- region_grower/context.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/region_grower/context.py b/region_grower/context.py index 2291b07..1f8a171 100644 --- a/region_grower/context.py +++ b/region_grower/context.py @@ -137,11 +137,13 @@ def section_prob(seg_direction, current_point, direction=None): if mode == "parallel": p = (1.0 - np.arccos(seg_direction.dot(target_direction)) / np.pi) ** power - if mode == "perpendicular": + elif mode == "perpendicular": p = ( 1.0 - abs(np.arccos(seg_direction.dot(target_direction)) / np.pi * 2.0 - 1.0) ) ** power + else: + raise ValueError(f"mode {mode} not understood for direction probability") return np.clip(p, 0, 1) direction["section_prob"] = partial(section_prob, direction=direction) @@ -228,7 +230,7 @@ def prob(direction, current_point, mesh=None, mesh_type=None, params=None): """Probability function for repulsive mode.""" dist = get_distance_to_mesh(mesh, current_point, direction, mesh_type=mesh_type) p = (dist - params.get("d_min", 0)) / ( - params.get("d_max", 100) - params.get("d_min", 0) + params.get("d_max", 100) - max(0, params.get("d_min", 0)) ) ** params.get("power", 1.5) return np.clip(p, 0, 1) @@ -246,6 +248,7 @@ def prob(direction, current_point, mesh=None, mesh_type=None, params=None): if params.get("d_min", 0) < distance < params.get("d_max", 100): if mesh_type == "voxel": current_point = self.positions_to_indices(current_point) + p = ( 1 - np.arccos(direction.dot(mesh_direction) / distance) / np.pi ) ** params.get("power", 1.0) From 3154bc79e9bbb668d15c5b5f8b0146dd1aedb516 Mon Sep 17 00:00:00 2001 From: arnaudon Date: Tue, 21 Nov 2023 13:25:41 +0100 Subject: [PATCH 20/63] more --- region_grower/context.py | 1 - region_grower/synthesize_morphologies.py | 14 ++++++++------ 2 files changed, 8 insertions(+), 7 deletions(-) diff --git a/region_grower/context.py b/region_grower/context.py index 1f8a171..028e177 100644 --- a/region_grower/context.py +++ b/region_grower/context.py @@ -475,7 +475,6 @@ def synthesize(self) -> SynthesisResult: """ rng = np.random.default_rng(self.params.seed) - return self._synthesize_once(rng) for _ in range(self.internals.retries): try: return self._synthesize_once(rng) diff --git a/region_grower/synthesize_morphologies.py b/region_grower/synthesize_morphologies.py index 873c088..0699965 100644 --- a/region_grower/synthesize_morphologies.py +++ b/region_grower/synthesize_morphologies.py @@ -308,7 +308,10 @@ def __init__( for region, params in self.tmd_parameters.items(): for param in params.values(): if synthesize_axons: - assert "axon" in param["grow_types"] + if "axon" not in param["grow_types"]: + LOGGER.warning( + "No axon data, but axon synthesis requested, you will not have axons" + ) elif "axon" in param["grow_types"]: param["grow_types"].remove("axon") @@ -368,7 +371,7 @@ def set_cortical_depths(self): for region in self.regions: if ( region not in self.atlas.region_structure - or self.atlas.region_structure[region]["thicknesses"] is None + or self.atlas.region_structure[region].get("thicknesses") is None ): # pragma: no cover self.cortical_depths[region] = self.cortical_depths["default"] else: @@ -503,8 +506,8 @@ def assign_atlas_data(self, min_depth=25, max_depth=5000): LOGGER.debug("Extract atlas data for %s region", _region) if ( _region in self.atlas.regions - and self.atlas.region_structure[_region].get("thicknesses", None) is not None - and self.atlas.region_structure[_region].get("layers", None) is not None + and self.atlas.region_structure[_region].get("thicknesses") is not None + and self.atlas.region_structure[_region].get("layers") is not None ): layers = self.atlas.layers[_region] thicknesses = [self.atlas.layer_thickness(layer) for layer in layers] @@ -515,8 +518,7 @@ def assign_atlas_data(self, min_depth=25, max_depth=5000): current_depths = np.clip(depths.lookup(positions), min_depth, max_depth) else: LOGGER.warning( - "We are not able to synthesize the region %s, we fallback to 'default' region", - _region, + "We are not able to synthesize the region %s with thicknesses", _region ) layer_depths = None current_depths = None From ccc33884cf741510d3e4e8c94c6a9e94c2957a57 Mon Sep 17 00:00:00 2001 From: arnaudon Date: Thu, 21 Dec 2023 14:55:45 +0100 Subject: [PATCH 21/63] more --- region_grower/context.py | 2 +- region_grower/synthesize_morphologies.py | 8 ++++---- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/region_grower/context.py b/region_grower/context.py index 028e177..b132d00 100644 --- a/region_grower/context.py +++ b/region_grower/context.py @@ -122,7 +122,7 @@ def section_prob(seg_direction, current_point, direction=None): target_direction = direction["params"]["direction"] power = direction["params"].get("power", 1.0) - mode = direction["params"]["mode"] + mode = direction["params"].get("mode", "parallel") layers = direction["params"].get("layers", []) if layers: depth = self.soma_depth + (self.soma_position - current_point).dot( diff --git a/region_grower/synthesize_morphologies.py b/region_grower/synthesize_morphologies.py index 0699965..beff644 100644 --- a/region_grower/synthesize_morphologies.py +++ b/region_grower/synthesize_morphologies.py @@ -486,11 +486,11 @@ def _init_parallel(self, mpi_only=False): def _close_parallel(self): if self._parallel_client is not None: LOGGER.debug("Closing the Dask client") - self._parallel_client.retire_workers() - time.sleep(1) - self._parallel_client.shutdown() + #self._parallel_client.retire_workers() + #time.sleep(1) + #self._parallel_client.shutdown() self._parallel_client.close() - self._parallel_client = None + #self._parallel_client = None def assign_atlas_data(self, min_depth=25, max_depth=5000): """Open an Atlas and compute depths and orientations according to the given positions.""" From 347559f4ade2b7a53a1bc01dce936f0e19a84d85 Mon Sep 17 00:00:00 2001 From: arnaudon Date: Thu, 11 Jan 2024 14:40:35 +0100 Subject: [PATCH 22/63] save mesh_name --- region_grower/context.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/region_grower/context.py b/region_grower/context.py index b132d00..0e84427 100644 --- a/region_grower/context.py +++ b/region_grower/context.py @@ -205,13 +205,16 @@ def get_distance_to_mesh(mesh, ray_origin, ray_direction, mesh_type): soma_position = self.positions_to_indices(soma_position) if boundary.get("multimesh_mode", "closest") == "closest": - for mesh_path in Path(boundary["path"]).iterdir(): + mesh_paths = list(Path(boundary["path"]).iterdir()) + for mesh_path in mesh_paths: meshes.append(trimesh.load_mesh(mesh_path)) distances = [ trimesh.proximity.closest_point(mesh, [soma_position])[1][0] for mesh in meshes ] - mesh = meshes[np.argmin(distances)] + mesh_id = np.argmin(distances) + boundary["mesh_name"] = Path(mesh_paths[mesh_id]).stem + mesh = meshes[mesh_id] if boundary.get("multimesh_mode", "closest") == "inside": mesh = None @@ -530,6 +533,8 @@ def _synthesize_once(self, rng) -> SynthesisResult: if "constraints" not in context: context["constraints"] = [] context["constraints"] += self.context.get_boundaries(mtype=self.cell.mtype) + if context["constraints"] and "mesh_name" in context["constraints"][-1]: + self.debug_infos["mesh_name"] = context["constraints"][-1]["mesh_name"] grower = NeuronGrower( input_parameters=params, From 3d5fa1eae80e024e40cf18c514f71e791cdabb92 Mon Sep 17 00:00:00 2001 From: arnaudon Date: Fri, 12 Jan 2024 11:47:18 +0100 Subject: [PATCH 23/63] shuffle rows --- region_grower/synthesize_morphologies.py | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/region_grower/synthesize_morphologies.py b/region_grower/synthesize_morphologies.py index beff644..ae5993b 100644 --- a/region_grower/synthesize_morphologies.py +++ b/region_grower/synthesize_morphologies.py @@ -6,11 +6,11 @@ - assign identity cell rotations to MVD3/sonata - optional axon grafting "on-the-fly" """ + import json import logging import os import subprocess -import time from pathlib import Path from shutil import which from typing import Optional @@ -486,11 +486,11 @@ def _init_parallel(self, mpi_only=False): def _close_parallel(self): if self._parallel_client is not None: LOGGER.debug("Closing the Dask client") - #self._parallel_client.retire_workers() - #time.sleep(1) - #self._parallel_client.shutdown() + # self._parallel_client.retire_workers() + # time.sleep(1) + # self._parallel_client.shutdown() self._parallel_client.close() - #self._parallel_client = None + # self._parallel_client = None def assign_atlas_data(self, min_depth=25, max_depth=5000): """Open an Atlas and compute depths and orientations according to the given positions.""" @@ -639,6 +639,9 @@ def compute(self): } ) + # shuffle rows to get more even loads on tasks + self.cells_data = self.cells_data.sample(frac=1.0).reset_index() + if self.nb_processes == 0: LOGGER.info("Start computation") computed = self.cells_data.apply( @@ -658,8 +661,7 @@ def compute(self): computed = future.compute() LOGGER.info("Format results") - res = self.cells_data.join(computed) - return res + return self.cells_data.join(computed).sort_values(by="index").set_index("index") def finalize(self, result: pd.DataFrame): """Finalize master work. From e327baae54e1165597bee730508f8f7c16dc9b20 Mon Sep 17 00:00:00 2001 From: arnaudon Date: Thu, 1 Feb 2024 10:21:24 +0100 Subject: [PATCH 24/63] closest_y option --- region_grower/context.py | 25 ++++++++++++++++++++++++- 1 file changed, 24 insertions(+), 1 deletion(-) diff --git a/region_grower/context.py b/region_grower/context.py index 0e84427..e7493e0 100644 --- a/region_grower/context.py +++ b/region_grower/context.py @@ -216,6 +216,22 @@ def get_distance_to_mesh(mesh, ray_origin, ray_direction, mesh_type): boundary["mesh_name"] = Path(mesh_paths[mesh_id]).stem mesh = meshes[mesh_id] + if boundary.get("multimesh_mode", "closest") == "closest_y": + mesh_paths = list(Path(boundary["path"]).iterdir()) + meshes = [trimesh.load_mesh(mesh_path) for mesh_path in mesh_paths] + # this is assuming the original direction can be used + # as straight line, refinement would be to follow the vector field + distances = [ + np.linalg.norm( + np.cross(soma_position - mesh.center_mass, boundary["direction"]) + ) + for mesh in meshes + ] + + mesh_id = np.argmin(distances) + boundary["mesh_name"] = Path(mesh_paths[mesh_id]).stem + mesh = meshes[mesh_id] + if boundary.get("multimesh_mode", "closest") == "inside": mesh = None for mesh_path in Path(boundary["path"]).iterdir(): @@ -303,7 +319,7 @@ def lookup_target_reference_depths(self, depth: float) -> np.array: # we round to avoid being outside due to numerical precision if np.round(depth, 3) <= np.round(layer_depth, 3): return layer_depth, cortical_depth - + print(depth, self.layer_depths, self.cortical_depths) raise RegionGrowerError(f"Current depth ({depth}) is outside of circuit boundaries") def distance_to_constraint(self, depth: float, constraint: Dict) -> Optional[float]: @@ -384,6 +400,13 @@ def _correct_position_orientation_scaling(self, params_orig: Dict) -> Dict: PIA_DIRECTION ).tolist() + if self.context.boundaries is not None: + boundaries = json.loads(self.context.boundaries) + for boundary in boundaries: + if boundary.get("multimesh_mode", "closest") == "closest_y": + boundary["direction"] = self.cell.lookup_orientation(PIA_DIRECTION).tolist() + self.context.boundaries = json.dumps(boundaries) + for neurite_type in params["grow_types"]: if isinstance(params[neurite_type]["orientation"], dict): params["pia_direction"] = self.cell.lookup_orientation(PIA_DIRECTION).tolist() From 78912d608978f60b71430d5f253f95bf168a3c76 Mon Sep 17 00:00:00 2001 From: arnaudon Date: Fri, 2 Feb 2024 09:34:35 +0100 Subject: [PATCH 25/63] remove neuroc for now --- setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.py b/setup.py index e663045..f829367 100644 --- a/setup.py +++ b/setup.py @@ -20,7 +20,7 @@ "diameter-synthesis>=0.5.4,<1", "morphio>=3.3.6,<4", "morph-tool[nrn]>=2.9,<3", - "neuroc>=0.2.8,<1", + #"neuroc>=0.2.8,<1", "neurom>=3,<4", "neurots>=3.5,<4", "numpy>=1.22", From 1755b38ec076c134220e609ea52e816f67f0b563 Mon Sep 17 00:00:00 2001 From: arnaudon Date: Fri, 2 Feb 2024 09:40:39 +0100 Subject: [PATCH 26/63] revert --- setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.py b/setup.py index f829367..e663045 100644 --- a/setup.py +++ b/setup.py @@ -20,7 +20,7 @@ "diameter-synthesis>=0.5.4,<1", "morphio>=3.3.6,<4", "morph-tool[nrn]>=2.9,<3", - #"neuroc>=0.2.8,<1", + "neuroc>=0.2.8,<1", "neurom>=3,<4", "neurots>=3.5,<4", "numpy>=1.22", From 1cd25cf6492f508999d13fc04d6f471710cfa246 Mon Sep 17 00:00:00 2001 From: arnaudon Date: Fri, 1 Mar 2024 09:49:08 +0100 Subject: [PATCH 27/63] partial fix --- region_grower/context.py | 17 ++++++++++------- region_grower/synthesize_morphologies.py | 11 ++++++----- tests/test_synthesize_morphologies.py | 5 ++++- 3 files changed, 20 insertions(+), 13 deletions(-) diff --git a/region_grower/context.py b/region_grower/context.py index e7493e0..022d09d 100644 --- a/region_grower/context.py +++ b/region_grower/context.py @@ -6,16 +6,18 @@ TLDR: SpaceContext.synthesized() is being called by the placement_algorithm package to synthesize circuit morphologies. """ -import os + import json +import os from collections import defaultdict from copy import deepcopy -from pathlib import Path from functools import partial +from pathlib import Path from typing import Dict from typing import List from typing import Optional from typing import Union + import trimesh os.environ.setdefault("NEURON_MODULE_OPTIONS", "-nogui") # suppress NEURON warning @@ -111,15 +113,15 @@ def positions_to_indices(self, positions): return result def get_direction(self, mtype): + """Get direction data for the given mtype.""" directions = [] - for i, direction in enumerate(self.directions): + for direction in self.directions: mtypes = direction.pop("mtypes", None) if mtypes is not None and mtype not in mtypes: continue def section_prob(seg_direction, current_point, direction=None): """Probability function for the sections.""" - target_direction = direction["params"]["direction"] power = direction["params"].get("power", 1.0) mode = direction["params"].get("mode", "parallel") @@ -150,7 +152,9 @@ def section_prob(seg_direction, current_point, direction=None): directions.append(direction) return directions - def get_boundaries(self, mtype): + def get_boundaries( + self, mtype + ): # pylint: disable=too-many-locals,too-many-statements,too-many-branches """Returns a dict with boundaries data for NeuroTS.""" def get_distance_to_mesh(mesh, ray_origin, ray_direction, mesh_type): @@ -190,7 +194,7 @@ def get_distance_to_mesh(mesh, ray_origin, ray_direction, mesh_type): # add here necessary logic to convert raw config data to NeuroTS context data all_boundaries = json.loads(self.boundaries) boundaries = [] - for i, boundary in enumerate(all_boundaries): + for boundary in all_boundaries: mtypes = boundary.pop("mtypes", None) if mtypes is not None and mtype not in mtypes: continue @@ -319,7 +323,6 @@ def lookup_target_reference_depths(self, depth: float) -> np.array: # we round to avoid being outside due to numerical precision if np.round(depth, 3) <= np.round(layer_depth, 3): return layer_depth, cortical_depth - print(depth, self.layer_depths, self.cortical_depths) raise RegionGrowerError(f"Current depth ({depth}) is outside of circuit boundaries") def distance_to_constraint(self, depth: float, constraint: Dict) -> Optional[float]: diff --git a/region_grower/synthesize_morphologies.py b/region_grower/synthesize_morphologies.py index ae5993b..8ee8a79 100644 --- a/region_grower/synthesize_morphologies.py +++ b/region_grower/synthesize_morphologies.py @@ -11,6 +11,7 @@ import logging import os import subprocess +import time from pathlib import Path from shutil import which from typing import Optional @@ -305,7 +306,7 @@ def __init__( with open(tmd_distributions, "r", encoding="utf-8") as f: self.tmd_distributions = convert_from_legacy_neurite_type(json.load(f)) - for region, params in self.tmd_parameters.items(): + for params in self.tmd_parameters.values(): for param in params.values(): if synthesize_axons: if "axon" not in param["grow_types"]: @@ -486,11 +487,11 @@ def _init_parallel(self, mpi_only=False): def _close_parallel(self): if self._parallel_client is not None: LOGGER.debug("Closing the Dask client") - # self._parallel_client.retire_workers() - # time.sleep(1) - # self._parallel_client.shutdown() + self._parallel_client.retire_workers() + time.sleep(1) + self._parallel_client.shutdown() self._parallel_client.close() - # self._parallel_client = None + self._parallel_client = None def assign_atlas_data(self, min_depth=25, max_depth=5000): """Open an Atlas and compute depths and orientations according to the given positions.""" diff --git a/tests/test_synthesize_morphologies.py b/tests/test_synthesize_morphologies.py index e41fdd6..e506fcf 100644 --- a/tests/test_synthesize_morphologies.py +++ b/tests/test_synthesize_morphologies.py @@ -178,6 +178,7 @@ def test_synthesize_dask_config( None, None, 100, + DATA / "region_structure.yaml", ) custom_scratch_config = str(tmp_folder / "custom_scratch_config") @@ -299,7 +300,7 @@ def test_synthesize_boundary( mesh, with_sections, with_trunks, -): # pylint: disable=unused-argument +): # pylint: disable=unused-argument,too-many-locals """Test morphology synthesis but skip write step.""" with_axon = True with_NRN = True @@ -716,6 +717,7 @@ def test_inconsistent_params( axon_morph_tsv, "apical_NRN_sections.yaml", min_depth, + DATA / "region_structure.yaml", ) args["out_morph_ext"] = ["h5"] @@ -748,6 +750,7 @@ def test_inconsistent_context( axon_morph_tsv, "apical_NRN_sections.yaml", min_depth, + DATA / "region_structure.yaml", ) args["nb_processes"] = 0 From 202c268d3c06a2f4bc942d243bef1f917871b514 Mon Sep 17 00:00:00 2001 From: arnaudon Date: Wed, 6 Mar 2024 15:23:07 +0100 Subject: [PATCH 28/63] fix --- region_grower/synthesize_morphologies.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/region_grower/synthesize_morphologies.py b/region_grower/synthesize_morphologies.py index 8ee8a79..e4246c2 100644 --- a/region_grower/synthesize_morphologies.py +++ b/region_grower/synthesize_morphologies.py @@ -110,6 +110,7 @@ def inverse_mapper(self): def _parallel_wrapper( row, computation_parameters, + cortical_depths, rotational_jitter_std, scaling_jitter_std, min_hard_scale, @@ -126,7 +127,7 @@ def _parallel_wrapper( ) current_space_context = SpaceContext( layer_depths=row["layer_depths"], - cortical_depths=row["cortical_depths"], + cortical_depths=cortical_depths[row["synthesis_region"]], directions=row["directions"] if "directions" in row else None, boundaries=row["boundaries"] if "boundaries" in row else None, atlas_info=json.loads(row["atlas_info"]) if "atlas_info" in row else None, @@ -621,6 +622,7 @@ def compute(self): func_kwargs = { "computation_parameters": computation_parameters, "rotational_jitter_std": self.rotational_jitter_std, + "cortical_depths": self.cortical_depths, "scaling_jitter_std": self.scaling_jitter_std, "min_hard_scale": self.min_hard_scale, "tmd_parameters": self.tmd_parameters, From 3a159771b2a1e10641a170b08f7631ad7a1a755f Mon Sep 17 00:00:00 2001 From: arnaudon Date: Wed, 19 Apr 2023 14:30:00 +0200 Subject: [PATCH 29/63] squashed --- .pre-commit-config.yaml | 2 +- region_grower/atlas_helper.py | 2 +- region_grower/context.py | 100 ++++++++++++++++++++--- region_grower/synthesize_morphologies.py | 21 ++++- setup.py | 1 + tests/data_factories.py | 2 +- tests/test_context.py | 37 ++++----- tox.ini | 4 +- 8 files changed, 132 insertions(+), 37 deletions(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 6fa9e94..92d0ae0 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -1,5 +1,5 @@ default_language_version: - python: python3.8 + python: python3 repos: - repo: https://github.com/pre-commit/pre-commit-hooks rev: v4.4.0 diff --git a/region_grower/atlas_helper.py b/region_grower/atlas_helper.py index f41b8da..98c37fb 100644 --- a/region_grower/atlas_helper.py +++ b/region_grower/atlas_helper.py @@ -32,7 +32,7 @@ def __init__(self, atlas: Atlas, region_structure_path: str): with open(region_structure_path, "r", encoding="utf-8") as region_file: self.region_structure = yaml.safe_load(region_file) else: - raise ValueError("Please provide an existing region_structure file.") + raise ValueError(f"region_structure file at {region_structure_path} not found.") self.regions = list(self.region_structure.keys()) self.layers = {} diff --git a/region_grower/context.py b/region_grower/context.py index 5a46284..62e802c 100644 --- a/region_grower/context.py +++ b/region_grower/context.py @@ -7,12 +7,14 @@ the placement_algorithm package to synthesize circuit morphologies. """ import os +import json from collections import defaultdict from copy import deepcopy from typing import Dict from typing import List from typing import Optional from typing import Union +import trimesh os.environ.setdefault("NEURON_MODULE_OPTIONS", "-nogui") # suppress NEURON warning # This environment variable must be set before 'NEURON' is loaded in 'morph_tool.nrnhines'. @@ -83,6 +85,86 @@ class SpaceContext: layer_depths: List cortical_depths: List + boundaries: Dict = None + atlas_info: Dict = {} # voxel_dimnesions, offset and shape from atlas for indices conversion + + def indices_to_positions(self, indices): + """Local version of voxcel's function to prevent atlas loading.""" + return indices * np.array(self.atlas_info["voxel_dimensions"]) + np.array( + self.atlas_info["offset"] + ) + + def positions_to_indices(self, positions): + """Local and reduced version of voxcel's function to prevent atlas loading.""" + result = (positions - np.array(self.atlas_info["offset"])) / np.array( + self.atlas_info["voxel_dimensions"] + ) + result[np.abs(result) < 1e-7] = 0.0 # suppress rounding errors around 0 + result = np.floor(result).astype(int) + result[result < 0] = -1 + result[result >= self.atlas_info["shape"]] = -1 + return result + + def get_boundaries(self): + """Returns a dict with boundaries data for NeuroTS.""" + # add here necessary logic to convert raw config data to NeuroTS context data + self.boundaries = json.loads(self.boundaries) + mesh = trimesh.load_mesh(self.boundaries["path"]) + + # to save points for plotting while growing to debug the prob function + debug = True + + def get_distance_to_mesh(mesh, ray_origin, ray_direction): + """Compute distances from point/directions to a mesh.""" + vox_ray_origin = self.positions_to_indices(ray_origin) + locations = mesh.ray.intersects_location( + ray_origins=[vox_ray_origin], ray_directions=[ray_direction] + )[0] + if len(locations): + intersect = self.indices_to_positions(locations[0]) + dist = np.linalg.norm(intersect - ray_origin) + + if debug: + x, y, z = ray_origin + vx, vy, vz = ray_direction + with open("data.csv", "a", encoding="utf8") as file: + print(f"{x}, {y}, {z}, {vx}, {vy}, {vz}, {dist}", file=file) + return dist + if debug: + x, y, z = ray_origin + vx, vy, vz = ray_direction + with open("data.csv", "a", encoding="utf8") as file: + print(f"{x}, {y}, {z}, {vx}, {vy}, {vz}, {-100}", file=file) + return np.inf + + def _prob(dist, params): + """Probability function from distance to boundary.""" + p = (dist - params.get("d_min", 0)) / ( + params.get("d_max", 100) - params.get("d_min", 0) + ) ** params.get("power", 1.5) + + # TODO: other variants for this function coud be included here + return np.clip(p, 0, 1) + + if "params_section" in self.boundaries: + + def section_prob(direction, current_point): + """Probability function for the sections.""" + dist = get_distance_to_mesh(mesh, current_point, direction) + return _prob(dist, self.boundaries["params_section"]) + + self.boundaries["section_prob"] = section_prob + + if "params_trunk" in self.boundaries: + + def trunk_prob(direction, current_point): + """Probability function for the trunks.""" + dist = get_distance_to_mesh(mesh, current_point, direction) + return _prob(dist, self.boundaries["params_trunk"]) + + self.boundaries["trunk_prob"] = trunk_prob + + return self.boundaries def layer_fraction_to_position(self, layer: int, layer_fraction: float) -> float: """Returns an absolute position from a layer and a fraction of the layer. @@ -141,7 +223,6 @@ class SynthesisParameters: axon_morph_scale: Optional[float] = None rotational_jitter_std: float = None scaling_jitter_std: float = None - recenter: bool = True seed: int = 0 min_hard_scale: float = 0 @@ -313,14 +394,7 @@ def completion(self, synthesized): def _synthesize_once(self, rng) -> SynthesisResult: """One try to synthesize the cell.""" params = self._correct_position_orientation_scaling(self.params.tmd_parameters) - - # Today we don't use the atlas during the synthesis (we just use it to - # generate the parameters) so we can - # grow the cell as if it was in [0, 0, 0] - # But the day we use it during the actual growth, we will need to grow the cell at its - # absolute position and translate to [0, 0, 0] after the growth - if self.params.recenter: - params["origin"] = [0, 0, 0] + params["origin"] = self.cell.position.tolist() if self.params.tmd_parameters["diameter_params"]["method"] == "external": external_diametrizer = build_diameters.build @@ -332,16 +406,20 @@ def _synthesize_once(self, rng) -> SynthesisResult: else: axon_morph = None + context = {"debug_data": self.debug_infos["input_scaling"]} + if self.context.boundaries is not None: + context.update(self.context.get_boundaries()) + grower = NeuronGrower( input_parameters=params, input_distributions=self.params.tmd_distributions, external_diametrizer=external_diametrizer, skip_preprocessing=True, - context={"debug_data": self.debug_infos["input_scaling"]}, + context=context, rng_or_seed=rng, ) grower.grow() - + mt.translate(grower.neuron, -self.cell.position) self._post_growth_rescaling(grower, params) if axon_morph is not None: diff --git a/region_grower/synthesize_morphologies.py b/region_grower/synthesize_morphologies.py index 785c2a3..9720b27 100644 --- a/region_grower/synthesize_morphologies.py +++ b/region_grower/synthesize_morphologies.py @@ -126,7 +126,9 @@ def _parallel_wrapper( ) current_space_context = SpaceContext( layer_depths=row["layer_depths"], - cortical_depths=cortical_depths[row["synthesis_region"]], + cortical_depths=cortical_depths, + boundaries=row["boundaries"] if "boundaries" in row else None, + atlas_info=json.loads(row["atlas_info"]) if "atlas_info" in row else None, ) axon_scale = row.get("axon_scale", None) @@ -144,7 +146,6 @@ def _parallel_wrapper( axon_morph_scale=axon_scale, rotational_jitter_std=rotational_jitter_std, scaling_jitter_std=scaling_jitter_std, - recenter=True, seed=row["seed"], min_hard_scale=min_hard_scale, ) @@ -326,12 +327,28 @@ def __init__( self.cells_data["synthesis_region"] = self.cells_data["region"].apply( lambda region: self.region_mapper[region] ) + self.cells_data["current_depth"] = current_depths + self.cells_data["layer_depths"] = layer_depths.T.tolist() + + if "boundaries" in self.atlas.region_structure: + self.cells_data["boundaries"] = json.dumps(self.atlas.region_structure["boundaries"]) + self.cells_data["atlas_info"] = json.dumps( + { + "voxel_dimensions": self.atlas.depths.voxel_dimensions.tolist(), + "offset": self.atlas.depths.offset.tolist(), + "shape": self.atlas.depths.shape, + } + ) + + if not self.cells.orientations: # pragma: no cover + self.cells_data["orientation"] = orientations.tolist() LOGGER.info("Checking TMD parameters and distributions according to cells mtypes") self.verify() LOGGER.info("Fetching atlas data from %s", atlas) self.assign_atlas_data(min_depth, max_depth) + if morph_axon is not None: LOGGER.info("Loading axon morphologies from %s", morph_axon) self.axon_morph_list = load_morphology_list(morph_axon, self.task_ids) diff --git a/setup.py b/setup.py index e16c5e2..ed538bb 100644 --- a/setup.py +++ b/setup.py @@ -28,6 +28,7 @@ "tqdm>=4.28.1", "urllib3>=1.26,<2; python_version < '3.9'", "voxcell>=3.1.1,<4", + "trimesh", ] mpi_extras = [ diff --git a/tests/data_factories.py b/tests/data_factories.py index a496ba5..565ce2e 100644 --- a/tests/data_factories.py +++ b/tests/data_factories.py @@ -132,7 +132,7 @@ def get_tmd_distributions(filename): def get_cell_position(): """The cell position.""" - return [0, 500, 0] + return np.array([0, 500, 0]) def get_cell_mtype(): diff --git a/tests/test_context.py b/tests/test_context.py index bc1736e..934e585 100644 --- a/tests/test_context.py +++ b/tests/test_context.py @@ -116,19 +116,18 @@ def test_context( assert_array_almost_equal( synthesis_parameters.tmd_parameters["apical_dendrite"]["orientation"], [[0.0, 1.0, 0.0]] ) - assert_array_almost_equal( result.neuron.soma.points, np.array( [ - [-5.785500526428223, 4.9841227531433105, 0.0], - [-7.5740227699279785, -0.9734848141670227, 0.0], - [-1.966903805732727, -7.378671169281006, 0.0], - [3.565324068069458, -6.752922534942627, 0.0], - [7.266839027404785, -2.346604108810425, 0.0], - [7.384983062744141, 1.059786081314087, 0.0], - [6.818241119384766, 3.4387624263763428, 0.0], - [4.675901919111924e-16, 7.636327266693115, 0.0], + [-5.785500526428223, 4.984130859375, 0.0], + [-7.5740227699279785, -0.973480224609375, 0.0], + [-1.966903805732727, -7.378662109375, 0.0], + [3.565324068069458, -6.7529296875, 0.0], + [7.266839027404785, -2.34661865234375, 0.0], + [7.384983062744141, 1.059783935546875, 0.0], + [6.818241119384766, 3.438751220703125, 0.0], + [4.675901919111924e-16, 7.636322021484375, 0.0], ], dtype=np.float32, ), @@ -189,25 +188,25 @@ def test_context_external_diametrizer( assert_array_almost_equal( synthesis_parameters.tmd_parameters["apical_dendrite"]["orientation"], [[0.0, 1.0, 0.0]] ) - assert_array_almost_equal( result.neuron.soma.points, np.array( [ - [-5.785500526428223, 4.9841227531433105, 0.0], - [-7.5740227699279785, -0.9734848141670227, 0.0], - [-1.966903805732727, -7.378671169281006, 0.0], - [3.565324068069458, -6.752922534942627, 0.0], - [7.266839027404785, -2.346604108810425, 0.0], - [7.384983062744141, 1.059786081314087, 0.0], - [6.818241119384766, 3.4387624263763428, 0.0], - [4.675901919111924e-16, 7.636327266693115, 0.0], + [-5.785500526428223, 4.984130859375, 0.0], + [-7.5740227699279785, -0.973480224609375, 0.0], + [-1.966903805732727, -7.378662109375, 0.0], + [3.565324068069458, -6.7529296875, 0.0], + [7.266839027404785, -2.34661865234375, 0.0], + [7.384983062744141, 1.059783935546875, 0.0], + [6.818241119384766, 3.438751220703125, 0.0], + [4.675901919111924e-16, 7.636322021484375, 0.0], ], dtype=np.float32, ), ) assert len(result.neuron.root_sections) == 4 + assert_array_almost_equal( next(result.neuron.iter()).points, np.array( @@ -258,7 +257,7 @@ def test_scale(self, small_context_worker, tmd_parameters): assert_array_almost_equal( [ # Check only first and last points of neurites - np.around(np.array([neu.points[0], neu.points[-1]]), 6) + np.around(np.array([neu.points[0], neu.points[-1]]), 5) for neu in result.neuron.root_sections ], [ diff --git a/tox.ini b/tox.ini index 5227512..649008e 100644 --- a/tox.ini +++ b/tox.ini @@ -81,7 +81,7 @@ commands = basepython = python3.8 [testenv:lint] -basepython = python3.8 +basepython = python3 deps = pre-commit pylint @@ -90,7 +90,7 @@ commands = pylint -j {env:PYLINT_NPROCS:1} {[base]files} [testenv:format] -basepython = python3.8 +basepython = python3 skip_install = true deps = codespell From 93d6260037eff077b107e7ffb893c57517271b5c Mon Sep 17 00:00:00 2001 From: arnaudon Date: Mon, 15 May 2023 16:13:50 +0200 Subject: [PATCH 30/63] improve --- region_grower/synthesize_morphologies.py | 31 ++++++++++-------------- 1 file changed, 13 insertions(+), 18 deletions(-) diff --git a/region_grower/synthesize_morphologies.py b/region_grower/synthesize_morphologies.py index 9720b27..449431b 100644 --- a/region_grower/synthesize_morphologies.py +++ b/region_grower/synthesize_morphologies.py @@ -109,7 +109,6 @@ def inverse_mapper(self): def _parallel_wrapper( row, computation_parameters, - cortical_depths, rotational_jitter_std, scaling_jitter_std, min_hard_scale, @@ -126,7 +125,7 @@ def _parallel_wrapper( ) current_space_context = SpaceContext( layer_depths=row["layer_depths"], - cortical_depths=cortical_depths, + cortical_depths=row["cortical_depths"], boundaries=row["boundaries"] if "boundaries" in row else None, atlas_info=json.loads(row["atlas_info"]) if "atlas_info" in row else None, ) @@ -327,21 +326,6 @@ def __init__( self.cells_data["synthesis_region"] = self.cells_data["region"].apply( lambda region: self.region_mapper[region] ) - self.cells_data["current_depth"] = current_depths - self.cells_data["layer_depths"] = layer_depths.T.tolist() - - if "boundaries" in self.atlas.region_structure: - self.cells_data["boundaries"] = json.dumps(self.atlas.region_structure["boundaries"]) - self.cells_data["atlas_info"] = json.dumps( - { - "voxel_dimensions": self.atlas.depths.voxel_dimensions.tolist(), - "offset": self.atlas.depths.offset.tolist(), - "shape": self.atlas.depths.shape, - } - ) - - if not self.cells.orientations: # pragma: no cover - self.cells_data["orientation"] = orientations.tolist() LOGGER.info("Checking TMD parameters and distributions according to cells mtypes") self.verify() @@ -542,6 +526,18 @@ def assign_atlas_data(self, min_depth=25, max_depth=5000): dtype=object, ) + if "boundaries" in self.atlas.region_structure[_region]: + self.cells_data.loc[region_mask, "boundaries"] = json.dumps( + self.atlas.region_structure[_region]["boundaries"] + ) + self.cells_data.loc[region_mask, "atlas_info"] = json.dumps( + { + "voxel_dimensions": self.atlas.depths[_region].voxel_dimensions.tolist(), + "offset": self.atlas.depths[_region].offset.tolist(), + "shape": self.atlas.depths[_region].shape, + } + ) + @property def task_ids(self): """Task IDs (= CellCollection IDs).""" @@ -605,7 +601,6 @@ def compute(self): func_kwargs = { "computation_parameters": computation_parameters, - "cortical_depths": self.cortical_depths, "rotational_jitter_std": self.rotational_jitter_std, "scaling_jitter_std": self.scaling_jitter_std, "min_hard_scale": self.min_hard_scale, From c4dc5300529d87fc24cce914954a738cddcb5652 Mon Sep 17 00:00:00 2001 From: arnaudon Date: Mon, 15 May 2023 18:50:47 +0200 Subject: [PATCH 31/63] fix CI --- region_grower/context.py | 2 +- setup.py | 3 +- tests/conftest.py | 31 ++- tests/data/apical_boundary.yaml | 16 ++ tests/data_factories.py | 22 ++ tests/test_context.py | 14 +- .../test_synthesize_morphologies_boundary.py | 212 ++++++++++++++++++ tox.ini | 4 +- 8 files changed, 297 insertions(+), 7 deletions(-) create mode 100644 tests/data/apical_boundary.yaml create mode 100644 tests/test_synthesize_morphologies_boundary.py diff --git a/region_grower/context.py b/region_grower/context.py index 62e802c..cb25055 100644 --- a/region_grower/context.py +++ b/region_grower/context.py @@ -143,7 +143,7 @@ def _prob(dist, params): params.get("d_max", 100) - params.get("d_min", 0) ) ** params.get("power", 1.5) - # TODO: other variants for this function coud be included here + # TODO: other variants for this function could be included here return np.clip(p, 0, 1) if "params_section" in self.boundaries: diff --git a/setup.py b/setup.py index ed538bb..833bdde 100644 --- a/setup.py +++ b/setup.py @@ -28,7 +28,7 @@ "tqdm>=4.28.1", "urllib3>=1.26,<2; python_version < '3.9'", "voxcell>=3.1.1,<4", - "trimesh", + "trimesh>=3.21", ] mpi_extras = [ @@ -55,6 +55,7 @@ "pytest-html>=2", "pytest-mock>=3.5", "pytest-xdist>=3.0.2", + "neurocollage>=0.3.0", ] setup( diff --git a/tests/conftest.py b/tests/conftest.py index ad40e30..75a12e9 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -19,6 +19,8 @@ from .data_factories import generate_cell_collection from .data_factories import generate_cells_df from .data_factories import generate_input_cells +from .data_factories import generate_mesh +from .data_factories import generate_region_structure_boundary from .data_factories import generate_small_O1 from .data_factories import get_cell_mtype from .data_factories import get_cell_orientation @@ -37,6 +39,25 @@ def small_O1_path(tmpdir_factory): return atlas_directory +@pytest.fixture(scope="session") +def mesh(small_O1_path, tmpdir_factory): + """Generate mesh from atlas.""" + mesh_path = str(tmpdir_factory.mktemp("mesh") / "mesh.obj") + + atlas = {"atlas": small_O1_path, "structure": DATA / "region_structure.yaml"} + generate_mesh(atlas, mesh_path) + return mesh_path + + +@pytest.fixture(scope="session") +def region_structure_boundary(tmpdir_factory, mesh): + """Generate region_structure file with boundary entries.""" + out_path = str(tmpdir_factory.mktemp("structure") / "region_structure.yaml") + region_structure_path = DATA / "region_structure.yaml" + generate_region_structure_boundary(region_structure_path, out_path, mesh) + return out_path + + @pytest.fixture(scope="session") def small_O1(small_O1_path): """Open the atlas.""" @@ -113,7 +134,15 @@ def space_context(): """A space context object.""" layer_depth = [0.0, 200.0, 300.0, 400.0, 500.0, 600.0, 800.0] thicknesses = [165, 149, 353, 190, 525, 700] - return SpaceContext(layer_depths=layer_depth, cortical_depths=np.cumsum(thicknesses)) + atlas_info = { + "voxel_dimensions": [100.0, 100.0, 100.0], + "offset": [-1100.0, -100.0, -1000.0], + "shape": (22, 10, 20), + } + + return SpaceContext( + layer_depths=layer_depth, cortical_depths=np.cumsum(thicknesses), atlas_info=atlas_info + ) @pytest.fixture(scope="function") diff --git a/tests/data/apical_boundary.yaml b/tests/data/apical_boundary.yaml new file mode 100644 index 0000000..6663c73 --- /dev/null +++ b/tests/data/apical_boundary.yaml @@ -0,0 +1,16 @@ +14a03569d26b949692e5dfe8cb1855fe: +- 2.271697998046875 +- 114.72030639648438 +- -5.1029510498046875 +216363698b529b4a97b750923ceb3ffd: +- -3.509429931640625 +- 115.28182983398438 +- 5.8795013427734375 +4462ebfc5f915ef09cfbac6e7687a66e: +- 12.257843017578125 +- 85.43954467773438 +- -19.165298461914062 +e3e70682c2094cac629f6fbed82c07cd: +- -38.110015869140625 +- 152.74227905273438 +- -20.128509521484375 diff --git a/tests/data_factories.py b/tests/data_factories.py index 565ce2e..b2566e6 100644 --- a/tests/data_factories.py +++ b/tests/data_factories.py @@ -8,6 +8,8 @@ import numpy as np import pandas as pd +import yaml +from neurocollage.mesh_helper import MeshHelper from voxcell import CellCollection DF_SIZE = 12 @@ -33,6 +35,26 @@ def generate_small_O1(directory): return str(directory) +def generate_mesh(atlas, mesh_path): + """Generate a mesh from atlas to test boundary code.""" + mesh_helper = MeshHelper(atlas, "O0") + mesh = mesh_helper.get_boundary_mesh() + mesh.export(mesh_path) # pylint: disable=no-member + + +def generate_region_structure_boundary(region_structure_path, out_path, mesh): + """Generate region_structure file with boundary entries.""" + with open(region_structure_path, encoding="utf-8") as f: + structure = yaml.safe_load(f) + structure["O0"]["boundaries"] = { + "params_section": {"d_min": 5, "d_max": 50}, + "params_trunk": {"d_min": 5, "d_max": 500}, + "path": mesh, + } + with open(out_path, "w", encoding="utf-8") as f: + yaml.dump(structure, f) + + def generate_cells_df(): """Raw data for the cell collection.""" x = [200] * DF_SIZE diff --git a/tests/test_context.py b/tests/test_context.py index 934e585..f3f16b8 100644 --- a/tests/test_context.py +++ b/tests/test_context.py @@ -78,9 +78,17 @@ def test_distance_to_constraint(self, space_context): assert space_context.distance_to_constraint(-100, {"layer": 6, "fraction": 0}) == -900 assert space_context.distance_to_constraint(1000, {"layer": 6, "fraction": 0}) == 200 - for i in [np.nan, None, []]: - space_context.layer_depths = i - assert space_context.distance_to_constraint(1000, {"layer": 6, "fraction": 0}) is None + def test_indices_to_positions(self, space_context): + pos = space_context.indices_to_positions([0, 0, 0]) + assert_array_equal(pos, [-1100.0, -100.0, -1000.0]) + pos = space_context.indices_to_positions([100, 100, 100]) + assert_array_equal(pos, [8900.0, 9900.0, 9000.0]) + + def test_positions_to_indices(self, space_context): + indices = space_context.positions_to_indices([-1100.0, -100.0, -1000.0]) + assert_array_equal(indices, [0, 0, 0]) + indices = space_context.positions_to_indices([8900.0, 9900.0, -9000.0]) + assert_array_equal(indices, [-1, -1, -1]) class TestSpaceWorker: diff --git a/tests/test_synthesize_morphologies_boundary.py b/tests/test_synthesize_morphologies_boundary.py new file mode 100644 index 0000000..ef5b74b --- /dev/null +++ b/tests/test_synthesize_morphologies_boundary.py @@ -0,0 +1,212 @@ +"""Test the region_grower.synthesize_morphologies module.""" +# pylint: disable=missing-function-docstring +import logging +import shutil +from pathlib import Path +from uuid import uuid4 + +from morph_tool.utils import iter_morphology_files +from numpy.testing import assert_allclose + +from region_grower.synthesize_morphologies import SynthesizeMorphologies + +from .test_synthesize_morphologies import check_yaml + +DATA = Path(__file__).parent / "data" + + +def create_args( + with_mpi, + tmp_folder, + input_cells, + atlas_path, + axon_morph_tsv, + out_apical_NRN_sections, + min_depth, + region_structure_boundary, +): + """Create the arguments used for tests.""" + args = {} + + # Circuit + args["input_cells"] = input_cells + + # Atlas + args["atlas"] = atlas_path + + # Parameters + args["tmd_distributions"] = DATA / "distributions.json" + args["tmd_parameters"] = DATA / "parameters.json" + args["seed"] = 0 + args["min_depth"] = min_depth + + # Internals + args["overwrite"] = True + args["out_morph_ext"] = ["h5", "swc", "asc"] + args["out_morph_dir"] = tmp_folder + args["out_apical"] = tmp_folder / "apical.yaml" + args["out_cells"] = str(tmp_folder / "test_cells.mvd3") + if out_apical_NRN_sections: + args["out_apical_nrn_sections"] = tmp_folder / out_apical_NRN_sections + else: + args["out_apical_nrn_sections"] = None + if with_mpi: + args["with_mpi"] = with_mpi + else: + args["nb_processes"] = 2 + + # Axons + args["base_morph_dir"] = str(DATA / "input-cells") + args["morph_axon"] = axon_morph_tsv + args["max_drop_ratio"] = 0.5 + args["rotational_jitter_std"] = 10 + args["scaling_jitter_std"] = 0.5 + args["region_structure"] = region_structure_boundary + + return args + + +def test_synthesize( + tmpdir, + small_O1_path, + input_cells, + axon_morph_tsv, + mesh, + region_structure_boundary, +): # pylint: disable=unused-argument + """Test morphology synthesis but skip write step.""" + with_axon = True + with_NRN = True + min_depth = 25 + tmp_folder = Path(tmpdir) + + args = create_args( + False, + tmp_folder, + input_cells, + small_O1_path, + axon_morph_tsv if with_axon else None, + "apical_NRN_sections.yaml" if with_NRN else None, + min_depth, + region_structure_boundary, + ) + args["skip_write"] = True + + synthesizer = SynthesizeMorphologies(**args) + res = synthesizer.synthesize() + + assert (res["x"] == 200).all() + assert (res["y"] == 200).all() + assert (res["z"] == 200).all() + assert res["name"].tolist() == [ + "e3e70682c2094cac629f6fbed82c07cd", + None, + "216363698b529b4a97b750923ceb3ffd", + None, + "14a03569d26b949692e5dfe8cb1855fe", + None, + "4462ebfc5f915ef09cfbac6e7687a66e", + None, + ] + print(res["apical_points"]) + assert_allclose(res["apical_points"][0], [[-38.110016, 152.74228, -20.12851]]) + assert res["apical_points"][1] is None + + assert_allclose(res["apical_points"][3], [[-3.50943, 115.28183, 5.8795013]]) + assert res["apical_points"][4] is None + + assert_allclose(res["apical_points"][6], [[2.271698, 114.72031, -5.102951]]) + assert res["apical_points"][7] is None + + assert_allclose(res["apical_points"][9], [[12.257843, 85.439545, -19.165298]]) + assert res["apical_points"][10] is None + + # Check that the morphologies were not written + res_files = tmpdir.listdir() + assert len(res_files) == 5 + assert sorted(i.basename for i in res_files) == [ + "apical.yaml", + "apical_NRN_sections.yaml", + "axon_morphs.tsv", + "input_cells.mvd3", + "test_cells.mvd3", + ] + + +def run_with_mpi(): + """Test morphology synthesis with MPI.""" + # pylint: disable=import-outside-toplevel, too-many-locals, import-error + from data_factories import generate_axon_morph_tsv + from data_factories import generate_cell_collection + from data_factories import generate_cells_df + from data_factories import generate_input_cells + from data_factories import generate_mesh + from data_factories import generate_region_structure_boundary + from data_factories import generate_small_O1 + from data_factories import input_cells_path + from mpi4py import MPI + + from region_grower.utils import setup_logger + + COMM = MPI.COMM_WORLD # pylint: disable=c-extension-no-member + rank = COMM.Get_rank() + MASTER_RANK = 0 + is_master = rank == MASTER_RANK + + tmp_folder = Path("/tmp/test-run-synthesis_" + str(uuid4())) + tmp_folder = COMM.bcast(tmp_folder, root=MASTER_RANK) + input_cells = input_cells_path(tmp_folder) + small_O1_path = str(tmp_folder / "atlas") + region_structure_path = str(tmp_folder / "region_structure_boundary.yaml") + args = create_args( + True, + tmp_folder, + input_cells, + small_O1_path, + tmp_folder / "axon_morphs.tsv", + "apical_NRN_sections.yaml", + min_depth=25, + region_structure_boundary=region_structure_path, + ) + + setup_logger("info") + logging.getLogger("distributed").setLevel(logging.ERROR) + + if is_master: + tmp_folder.mkdir(exist_ok=True) + print(f"============= #{rank}: Create data") + cells_raw_data = generate_cells_df() + cell_collection = generate_cell_collection(cells_raw_data) + generate_input_cells(cell_collection, tmp_folder) + generate_small_O1(small_O1_path) + atlas = {"atlas": small_O1_path, "structure": DATA / "region_structure.yaml"} + generate_mesh(atlas, tmp_folder / "mesh.obj") + region_structure_path = generate_region_structure_boundary( + DATA / "region_structure.yaml", region_structure_path, str(tmp_folder / "mesh.obj") + ) + + generate_axon_morph_tsv(tmp_folder) + for dest in range(1, COMM.Get_size()): + req = COMM.isend("done", dest=dest) + else: + req = COMM.irecv(source=0) + req.wait() + + synthesizer = SynthesizeMorphologies(**args) + try: + print(f"============= #{rank}: Start synthesize") + synthesizer.synthesize() + + # Check results + print(f"============= #{rank}: Checking results") + expected_size = 12 + assert len(list(iter_morphology_files(tmp_folder))) == expected_size + check_yaml(DATA / "apical_boundary.yaml", args["out_apical"]) + finally: + # Clean the directory + print(f"============= #{rank}: Cleaning") + shutil.rmtree(tmp_folder) + + +if __name__ == "__main__": # pragma: no cover + run_with_mpi() diff --git a/tox.ini b/tox.ini index 649008e..5669993 100644 --- a/tox.ini +++ b/tox.ini @@ -49,7 +49,9 @@ deps = pytest allowlist_externals = mpiexec -commands = mpiexec -n 4 {envbindir}/python tests/test_synthesize_morphologies.py +commands = + mpiexec -n 4 {envbindir}/python tests/test_synthesize_morphologies.py + mpiexec -n 4 {envbindir}/python tests/test_synthesize_morphologies_boundary.py [testenv:coverage] skip_install = true From 36f69f0a6fd85271639ad4c13a81aa8f58b23737 Mon Sep 17 00:00:00 2001 From: arnaudon Date: Mon, 15 May 2023 21:09:10 +0200 Subject: [PATCH 32/63] set debug off --- region_grower/context.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/region_grower/context.py b/region_grower/context.py index cb25055..c970e93 100644 --- a/region_grower/context.py +++ b/region_grower/context.py @@ -112,7 +112,7 @@ def get_boundaries(self): mesh = trimesh.load_mesh(self.boundaries["path"]) # to save points for plotting while growing to debug the prob function - debug = True + debug = False def get_distance_to_mesh(mesh, ray_origin, ray_direction): """Compute distances from point/directions to a mesh.""" @@ -163,7 +163,6 @@ def trunk_prob(direction, current_point): return _prob(dist, self.boundaries["params_trunk"]) self.boundaries["trunk_prob"] = trunk_prob - return self.boundaries def layer_fraction_to_position(self, layer: int, layer_fraction: float) -> float: From 3da2b358766cbb0a4cdbb658ace4550a088da8c7 Mon Sep 17 00:00:00 2001 From: arnaudon Date: Tue, 16 May 2023 11:57:38 +0200 Subject: [PATCH 33/63] get cov to 100 --- region_grower/context.py | 7 +- setup.py | 1 + tests/conftest.py | 10 - tests/data/apical_boundary.yaml | 24 +- tests/data/apical_no_axon.yaml | 30 +-- tests/data/parameters.json | 58 +++-- tests/data_factories.py | 17 +- tests/test_context.py | 54 ++--- tests/test_synthesize_morphologies.py | 195 +++++++++++++++- .../test_synthesize_morphologies_boundary.py | 212 ------------------ tox.ini | 2 +- 11 files changed, 302 insertions(+), 308 deletions(-) delete mode 100644 tests/test_synthesize_morphologies_boundary.py diff --git a/region_grower/context.py b/region_grower/context.py index c970e93..eed8bf9 100644 --- a/region_grower/context.py +++ b/region_grower/context.py @@ -124,18 +124,18 @@ def get_distance_to_mesh(mesh, ray_origin, ray_direction): intersect = self.indices_to_positions(locations[0]) dist = np.linalg.norm(intersect - ray_origin) - if debug: + if debug: # pragma: no cover x, y, z = ray_origin vx, vy, vz = ray_direction with open("data.csv", "a", encoding="utf8") as file: print(f"{x}, {y}, {z}, {vx}, {vy}, {vz}, {dist}", file=file) return dist - if debug: + if debug: # pragma: no cover x, y, z = ray_origin vx, vy, vz = ray_direction with open("data.csv", "a", encoding="utf8") as file: print(f"{x}, {y}, {z}, {vx}, {vy}, {vz}, {-100}", file=file) - return np.inf + return np.inf # pragma: no cover def _prob(dist, params): """Probability function from distance to boundary.""" @@ -322,7 +322,6 @@ def _post_growth_rescaling(self, grower: NeuronGrower, params: Dict) -> None: target_min_length=target_min_length, target_max_length=target_max_length, ) - if scale > self.params.min_hard_scale: scale_section(root_section, ScaleParameters(mean=scale), recursive=True) is_deleted = False diff --git a/setup.py b/setup.py index 833bdde..3a8c636 100644 --- a/setup.py +++ b/setup.py @@ -29,6 +29,7 @@ "urllib3>=1.26,<2; python_version < '3.9'", "voxcell>=3.1.1,<4", "trimesh>=3.21", + "rtree>=1.0.1", ] mpi_extras = [ diff --git a/tests/conftest.py b/tests/conftest.py index 75a12e9..7eaf84e 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -20,7 +20,6 @@ from .data_factories import generate_cells_df from .data_factories import generate_input_cells from .data_factories import generate_mesh -from .data_factories import generate_region_structure_boundary from .data_factories import generate_small_O1 from .data_factories import get_cell_mtype from .data_factories import get_cell_orientation @@ -49,15 +48,6 @@ def mesh(small_O1_path, tmpdir_factory): return mesh_path -@pytest.fixture(scope="session") -def region_structure_boundary(tmpdir_factory, mesh): - """Generate region_structure file with boundary entries.""" - out_path = str(tmpdir_factory.mktemp("structure") / "region_structure.yaml") - region_structure_path = DATA / "region_structure.yaml" - generate_region_structure_boundary(region_structure_path, out_path, mesh) - return out_path - - @pytest.fixture(scope="session") def small_O1(small_O1_path): """Open the atlas.""" diff --git a/tests/data/apical_boundary.yaml b/tests/data/apical_boundary.yaml index 6663c73..fe76a95 100644 --- a/tests/data/apical_boundary.yaml +++ b/tests/data/apical_boundary.yaml @@ -1,16 +1,16 @@ 14a03569d26b949692e5dfe8cb1855fe: -- 2.271697998046875 -- 114.72030639648438 -- -5.1029510498046875 +- 4.369110107421875 +- 55.535736083984375 +- -6.2556915283203125 216363698b529b4a97b750923ceb3ffd: -- -3.509429931640625 -- 115.28182983398438 -- 5.8795013427734375 +- 5.216094970703125 +- 114.55657958984375 +- -12.4947509765625 4462ebfc5f915ef09cfbac6e7687a66e: -- 12.257843017578125 -- 85.43954467773438 -- -19.165298461914062 +- -11.084274291992188 +- 111.49124145507812 +- 0.980438232421875 e3e70682c2094cac629f6fbed82c07cd: -- -38.110015869140625 -- 152.74227905273438 -- -20.128509521484375 +- 16.126876831054688 +- 160.52789306640625 +- 22.863800048828125 diff --git a/tests/data/apical_no_axon.yaml b/tests/data/apical_no_axon.yaml index bfd6aa6..424b442 100644 --- a/tests/data/apical_no_axon.yaml +++ b/tests/data/apical_no_axon.yaml @@ -19,9 +19,9 @@ - 127.51496124267578 - 16.63180923461914 6513270e269e0d37f2a74de452e6b438: -- 1.8200105428695679 -- 113.30069732666016 -- 0.27308279275894165 +- -8.12664794921875 +- 114.60183715820312 +- 1.4770965576171875 7b89296c6dcbac5008577eb1924770d3: - -2.5671095848083496 - 58.11355972290039 @@ -31,21 +31,21 @@ - 371.3540344238281 - 6.843407154083252 b8a1abcd1a6916c74da4f9fc3c6da5d7: -- 2.9995338916778564 -- 114.8524169921875 -- 12.829068183898926 +- -13.87677001953125 +- 162.44619750976562 +- 1.412567138671875 cd613e30d8f16adf91b7584a2265b1f5: -- 14.439825057983398 -- 114.19844818115234 -- -7.658063888549805 +- 1.1701507568359375 +- 97.20172119140625 +- 1.061187744140625 d95bafc8f2a4d27bdcf4bb99f4bea973: -- 3.320021867752075 -- 59.75845718383789 -- -2.885958433151245 +- -2.2343902587890625 +- 59.40399169921875 +- -0.331756591796875 db5b5fab8f4d3e27dda1494c73cf256d: -- 10.5548677444458 -- 93.02537536621094 -- 0.9866926074028015 +- -1.606964111328125 +- 58.9842529296875 +- 1.571624755859375 e3e70682c2094cac629f6fbed82c07cd: - -4.8529953956604 - 158.9239959716797 diff --git a/tests/data/parameters.json b/tests/data/parameters.json index bb571ac..3a024fb 100644 --- a/tests/data/parameters.json +++ b/tests/data/parameters.json @@ -7,13 +7,16 @@ "has_apical_tuft": true, "metric": "path_distances", "modify": null, - "orientation": [ - [ - 0.0, - 1.0, - 0.0 - ] - ], + "orientation": { + "mode": "normal_pia_constraint", + "values": { + "direction": { + "mean": 0.0, + "std": 0.0 + } + } + }, + "radius": 0.3, "randomness": 0.15, "step_size": { "norm": { @@ -30,7 +33,17 @@ "growth_method": "tmd", "metric": "path_distances", "modify": null, - "orientation": null, + "orientation": { + "mode": "pia_constraint", + "values": { + "form": "step", + "params": [ + 0.44484279339093336, + 0.5461595884090537 + ] + } + }, + "radius": 0.3, "randomness": 0.15, "step_size": { "norm": { @@ -66,13 +79,16 @@ "has_apical_tuft": true, "metric": "path_distances", "modify": null, - "orientation": [ - [ - 0.0, - 1.0, - 0.0 - ] - ], + "orientation": { + "mode": "normal_pia_constraint", + "values": { + "direction": { + "mean": 0.0, + "std": 0.0 + } + } + }, + "radius": 0.3, "randomness": 0.15, "step_size": { "norm": { @@ -89,7 +105,17 @@ "growth_method": "tmd", "metric": "path_distances", "modify": null, - "orientation": null, + "orientation": { + "mode": "pia_constraint", + "values": { + "form": "step", + "params": [ + 0.44484279339093336, + 0.5461595884090537 + ] + } + }, + "radius": 0.3, "randomness": 0.15, "step_size": { "norm": { diff --git a/tests/data_factories.py b/tests/data_factories.py index b2566e6..457f65d 100644 --- a/tests/data_factories.py +++ b/tests/data_factories.py @@ -42,15 +42,20 @@ def generate_mesh(atlas, mesh_path): mesh.export(mesh_path) # pylint: disable=no-member -def generate_region_structure_boundary(region_structure_path, out_path, mesh): +def generate_region_structure_boundary( + region_structure_path, out_path, mesh, with_sections=True, with_trunks=True +): """Generate region_structure file with boundary entries.""" with open(region_structure_path, encoding="utf-8") as f: structure = yaml.safe_load(f) - structure["O0"]["boundaries"] = { - "params_section": {"d_min": 5, "d_max": 50}, - "params_trunk": {"d_min": 5, "d_max": 500}, - "path": mesh, - } + structure["O0"]["boundaries"] = {"path": mesh} + + if with_sections: + structure["O0"]["boundaries"]["params_section"] = {"d_min": 5, "d_max": 50} + + if with_trunks: + structure["O0"]["boundaries"]["params_trunk"] = {"d_min": 5, "d_max": 5000} + with open(out_path, "w", encoding="utf-8") as f: yaml.dump(structure, f) diff --git a/tests/test_context.py b/tests/test_context.py index f3f16b8..bc8237a 100644 --- a/tests/test_context.py +++ b/tests/test_context.py @@ -121,21 +121,30 @@ def test_context( ) # This tests that input orientations are not mutated by the synthesize() call - assert_array_almost_equal( - synthesis_parameters.tmd_parameters["apical_dendrite"]["orientation"], [[0.0, 1.0, 0.0]] + assert ( + list( + dictdiffer.diff( + synthesis_parameters.tmd_parameters["apical_dendrite"]["orientation"], + { + "mode": "normal_pia_constraint", + "values": {"direction": {"mean": 0.0, "std": 0.0}}, + }, + ) + ) + == [] ) assert_array_almost_equal( result.neuron.soma.points, np.array( [ - [-5.785500526428223, 4.984130859375, 0.0], - [-7.5740227699279785, -0.973480224609375, 0.0], - [-1.966903805732727, -7.378662109375, 0.0], - [3.565324068069458, -6.7529296875, 0.0], - [7.266839027404785, -2.34661865234375, 0.0], - [7.384983062744141, 1.059783935546875, 0.0], - [6.818241119384766, 3.438751220703125, 0.0], - [4.675901919111924e-16, 7.636322021484375, 0.0], + [7.0957971, -2.8218994, 0.0], + [0.0, 7.6363220, 0.0], + [-3.0748143, 6.9899292, 0.0], + [-4.8864875, 5.8681946, 0.0], + [-7.5141201, -1.3606873, 0.0], + [-6.1610136, -3.6210632, 0.0], + [-1.2113731, -7.5396423, 0.0], + [0.67872918, -7.6061096, 0.0], ], dtype=np.float32, ), @@ -447,6 +456,11 @@ def test_scale(self, small_context_worker, tmd_parameters): } } result = small_context_worker.synthesize() + assert [i.type for i in result.neuron.root_sections] == expected_types[1:] + + # test removing basals if scale is too small + small_context_worker.params.min_hard_scale = 2.2 + result = small_context_worker.synthesize() assert [i.type for i in result.neuron.root_sections] == [expected_types[-1]] params["basal_dendrite"]["orientation"] = {} @@ -491,18 +505,9 @@ def test_debug_scales(self, small_context_worker, tmd_parameters): "min_target_thickness": 1.0, }, "scaling": [ - { - "max_ph": 126.96580088057513, - "scaling_ratio": 0.9554140127388535, - }, - { - "max_ph": 139.93140451880743, - "scaling_ratio": 0.9554140127388535, - }, - { - "max_ph": 143.28779114733433, - "scaling_ratio": 0.9554140127388535, - }, + {"max_ph": 78.53590605172101, "scaling_ratio": 0.9554140127388535}, + {"max_ph": 126.96580088057513, "scaling_ratio": 0.9554140127388535}, + {"max_ph": 126.96580088057513, "scaling_ratio": 0.9554140127388535}, ], }, "target_func": { @@ -514,10 +519,7 @@ def test_debug_scales(self, small_context_worker, tmd_parameters): "min_target_path_distance": 1.0, }, "scaling": [ - { - "max_ph": 420.70512274498446, - "scaling_ratio": 0.18064909574697127, - } + {"max_ph": 255.45843100194077, "scaling_ratio": 0.2975043716581138} ], }, }, diff --git a/tests/test_synthesize_morphologies.py b/tests/test_synthesize_morphologies.py index d1d72e2..6e12098 100644 --- a/tests/test_synthesize_morphologies.py +++ b/tests/test_synthesize_morphologies.py @@ -1,6 +1,7 @@ """Test the region_grower.synthesize_morphologies module.""" # pylint: disable=missing-function-docstring import json +import sys import logging import os import shutil @@ -26,6 +27,7 @@ from region_grower.synthesize_morphologies import RegionMapper from region_grower.synthesize_morphologies import SynthesizeMorphologies + DATA = Path(__file__).parent / "data" @@ -53,6 +55,7 @@ def create_args( axon_morph_tsv, out_apical_NRN_sections, min_depth, + region_structure, ): """Create the arguments used for tests.""" args = {} @@ -90,7 +93,7 @@ def create_args( args["max_drop_ratio"] = 0.5 args["rotational_jitter_std"] = 10 args["scaling_jitter_std"] = 0.5 - args["region_structure"] = DATA / "region_structure.yaml" + args["region_structure"] = region_structure return args @@ -118,6 +121,7 @@ def test_synthesize( axon_morph_tsv if with_axon else None, "apical_NRN_sections.yaml" if with_NRN else None, min_depth, + DATA / "region_structure.yaml", ) synthesizer = SynthesizeMorphologies(**args) @@ -233,6 +237,7 @@ def test_synthesize_skip_write( axon_morph_tsv if with_axon else None, "apical_NRN_sections.yaml" if with_NRN else None, min_depth, + DATA / "region_structure.yaml", ) args["skip_write"] = True args["nb_processes"] = nb_processes @@ -285,6 +290,109 @@ def test_synthesize_skip_write( ] +@pytest.mark.parametrize("with_sections", [True, False]) +@pytest.mark.parametrize("with_trunks", [True, False]) +def test_synthesize_boundary( + tmpdir, + small_O1_path, + input_cells, + axon_morph_tsv, + mesh, + with_sections, + with_trunks, +): # pylint: disable=unused-argument + """Test morphology synthesis but skip write step.""" + with_axon = True + with_NRN = True + min_depth = 25 + tmp_folder = Path(tmpdir) + + # pylint: disable=import-outside-toplevel + from .data_factories import generate_region_structure_boundary + + region_structure = "region_structure.yaml" + generate_region_structure_boundary( + DATA / "region_structure.yaml", region_structure, mesh, with_sections, with_trunks + ) + args = create_args( + False, + tmp_folder, + input_cells, + small_O1_path, + axon_morph_tsv if with_axon else None, + "apical_NRN_sections.yaml" if with_NRN else None, + min_depth, + region_structure, + ) + args["skip_write"] = True + + synthesizer = SynthesizeMorphologies(**args) + res = synthesizer.synthesize() + + assert (res["x"] == 200).all() + assert (res["y"] == 200).all() + assert (res["z"] == 200).all() + assert res["name"].tolist() == [ + "e3e70682c2094cac629f6fbed82c07cd", + None, + "216363698b529b4a97b750923ceb3ffd", + None, + "14a03569d26b949692e5dfe8cb1855fe", + None, + "4462ebfc5f915ef09cfbac6e7687a66e", + None, + ] + if with_sections and with_trunks: + assert_allclose(res["apical_points"][0], [[16.126877, 160.5279, 22.8638]]) + assert res["apical_points"][1] is None + + assert_allclose(res["apical_points"][3], [[5.216095, 114.55658, -12.494751]]) + assert res["apical_points"][4] is None + + assert_allclose(res["apical_points"][6], [[4.36911, 55.535736, -6.2556915]]) + assert res["apical_points"][7] is None + + assert_allclose(res["apical_points"][9], [[-11.084274, 111.49124, 0.98043823]]) + assert res["apical_points"][10] is None + + if with_sections and not with_trunks: + assert_allclose(res["apical_points"][0], [[-0.29234314, 58.81488, 0.5384369]]) + assert res["apical_points"][1] is None + + assert_allclose(res["apical_points"][3], [[8.230438, 116.570435, -7.345169]]) + assert res["apical_points"][4] is None + + assert_allclose(res["apical_points"][6], [[11.181992, 96.11371, -12.706863]]) + assert res["apical_points"][7] is None + + assert_allclose(res["apical_points"][9], [[3.5267944, 164.42618, 1.2018433]]) + assert res["apical_points"][10] is None + + if not with_sections and with_trunks: + assert_allclose(res["apical_points"][0], [[8.792313, 131.91104, -4.2198334]]) + assert res["apical_points"][1] is None + + assert_allclose(res["apical_points"][3], [[2.0048676, 116.672, -5.334656]]) + assert res["apical_points"][4] is None + + assert_allclose(res["apical_points"][6], [[-2.347641, 58.229614, -0.56985474]]) + assert res["apical_points"][7] is None + + assert_allclose(res["apical_points"][9], [[6.861679, 112.31445, -13.528854]]) + assert res["apical_points"][10] is None + + # Check that the morphologies were not written + res_files = tmpdir.listdir() + assert len(res_files) == 5 + assert sorted(i.basename for i in res_files) == [ + "apical.yaml", + "apical_NRN_sections.yaml", + "axon_morphs.tsv", + "input_cells.mvd3", + "test_cells.mvd3", + ] + + def run_with_mpi(): """Test morphology synthesis with MPI.""" # pylint: disable=import-outside-toplevel, too-many-locals, import-error @@ -316,6 +424,7 @@ def run_with_mpi(): tmp_folder / "axon_morphs.tsv", "apical_NRN_sections.yaml", min_depth=25, + region_structure=DATA / "region_structure.yaml", ) setup_logger("debug", prefix=f"Rank = {rank} - ") @@ -346,12 +455,83 @@ def run_with_mpi(): print(f"============= #{rank}: Checking results") expected_size = 18 assert len(list(iter_morphology_files(tmp_folder))) == expected_size + check_yaml(DATA / "apical.yaml", args["out_apical"]) + check_yaml(DATA / "apical_NRN_sections.yaml", args["out_apical_nrn_sections"]) + finally: + # Clean the directory + print(f"============= #{rank}: Cleaning") + shutil.rmtree(tmp_folder) + + +def run_with_mpi_boundary(): + """Test morphology synthesis with MPI.""" + # pylint: disable=import-outside-toplevel, too-many-locals, import-error + from data_factories import generate_axon_morph_tsv + from data_factories import generate_cell_collection + from data_factories import generate_cells_df + from data_factories import generate_input_cells + from data_factories import generate_mesh + from data_factories import generate_region_structure_boundary + from data_factories import generate_small_O1 + from data_factories import input_cells_path + from mpi4py import MPI + + from region_grower.utils import setup_logger + + COMM = MPI.COMM_WORLD # pylint: disable=c-extension-no-member + rank = COMM.Get_rank() + MASTER_RANK = 0 + is_master = rank == MASTER_RANK + + tmp_folder = Path("/tmp/test-run-synthesis_" + str(uuid4())) + tmp_folder = COMM.bcast(tmp_folder, root=MASTER_RANK) + input_cells = input_cells_path(tmp_folder) + small_O1_path = str(tmp_folder / "atlas") + region_structure_path = str(tmp_folder / "region_structure_boundary.yaml") + args = create_args( + True, + tmp_folder, + input_cells, + small_O1_path, + tmp_folder / "axon_morphs.tsv", + "apical_NRN_sections.yaml", + min_depth=25, + region_structure=region_structure_path, + ) + + setup_logger("info") + logging.getLogger("distributed").setLevel(logging.ERROR) - check_yaml(DATA / ("apical.yaml"), args["out_apical"]) - check_yaml( - DATA / ("apical_NRN_sections.yaml"), - args["out_apical_nrn_sections"], + if is_master: + tmp_folder.mkdir(exist_ok=True) + print(f"============= #{rank}: Create data") + cells_raw_data = generate_cells_df() + cell_collection = generate_cell_collection(cells_raw_data) + generate_input_cells(cell_collection, tmp_folder) + generate_small_O1(small_O1_path) + atlas = {"atlas": small_O1_path, "structure": DATA / "region_structure.yaml"} + generate_mesh(atlas, tmp_folder / "mesh.obj") + region_structure_path = generate_region_structure_boundary( + DATA / "region_structure.yaml", region_structure_path, str(tmp_folder / "mesh.obj") ) + + generate_axon_morph_tsv(tmp_folder) + for dest in range(1, COMM.Get_size()): + req = COMM.isend("done", dest=dest) + else: + req = COMM.irecv(source=0) + req.wait() + + synthesizer = SynthesizeMorphologies(**args) + try: + print(f"============= #{rank}: Start synthesize") + synthesizer.synthesize() + + # Check results + print(f"============= #{rank}: Checking results") + expected_size = 12 + assert len(list(iter_morphology_files(tmp_folder))) == expected_size + check_yaml(DATA / "apical_boundary.yaml", args["out_apical"]) finally: # Clean the directory print(f"============= #{rank}: Cleaning") @@ -603,4 +783,7 @@ def test_inconsistent_context( if __name__ == "__main__": # pragma: no cover - run_with_mpi() + if sys.argv[-1] == "boundary": + run_with_mpi_boundary() + else: + run_with_mpi() diff --git a/tests/test_synthesize_morphologies_boundary.py b/tests/test_synthesize_morphologies_boundary.py deleted file mode 100644 index ef5b74b..0000000 --- a/tests/test_synthesize_morphologies_boundary.py +++ /dev/null @@ -1,212 +0,0 @@ -"""Test the region_grower.synthesize_morphologies module.""" -# pylint: disable=missing-function-docstring -import logging -import shutil -from pathlib import Path -from uuid import uuid4 - -from morph_tool.utils import iter_morphology_files -from numpy.testing import assert_allclose - -from region_grower.synthesize_morphologies import SynthesizeMorphologies - -from .test_synthesize_morphologies import check_yaml - -DATA = Path(__file__).parent / "data" - - -def create_args( - with_mpi, - tmp_folder, - input_cells, - atlas_path, - axon_morph_tsv, - out_apical_NRN_sections, - min_depth, - region_structure_boundary, -): - """Create the arguments used for tests.""" - args = {} - - # Circuit - args["input_cells"] = input_cells - - # Atlas - args["atlas"] = atlas_path - - # Parameters - args["tmd_distributions"] = DATA / "distributions.json" - args["tmd_parameters"] = DATA / "parameters.json" - args["seed"] = 0 - args["min_depth"] = min_depth - - # Internals - args["overwrite"] = True - args["out_morph_ext"] = ["h5", "swc", "asc"] - args["out_morph_dir"] = tmp_folder - args["out_apical"] = tmp_folder / "apical.yaml" - args["out_cells"] = str(tmp_folder / "test_cells.mvd3") - if out_apical_NRN_sections: - args["out_apical_nrn_sections"] = tmp_folder / out_apical_NRN_sections - else: - args["out_apical_nrn_sections"] = None - if with_mpi: - args["with_mpi"] = with_mpi - else: - args["nb_processes"] = 2 - - # Axons - args["base_morph_dir"] = str(DATA / "input-cells") - args["morph_axon"] = axon_morph_tsv - args["max_drop_ratio"] = 0.5 - args["rotational_jitter_std"] = 10 - args["scaling_jitter_std"] = 0.5 - args["region_structure"] = region_structure_boundary - - return args - - -def test_synthesize( - tmpdir, - small_O1_path, - input_cells, - axon_morph_tsv, - mesh, - region_structure_boundary, -): # pylint: disable=unused-argument - """Test morphology synthesis but skip write step.""" - with_axon = True - with_NRN = True - min_depth = 25 - tmp_folder = Path(tmpdir) - - args = create_args( - False, - tmp_folder, - input_cells, - small_O1_path, - axon_morph_tsv if with_axon else None, - "apical_NRN_sections.yaml" if with_NRN else None, - min_depth, - region_structure_boundary, - ) - args["skip_write"] = True - - synthesizer = SynthesizeMorphologies(**args) - res = synthesizer.synthesize() - - assert (res["x"] == 200).all() - assert (res["y"] == 200).all() - assert (res["z"] == 200).all() - assert res["name"].tolist() == [ - "e3e70682c2094cac629f6fbed82c07cd", - None, - "216363698b529b4a97b750923ceb3ffd", - None, - "14a03569d26b949692e5dfe8cb1855fe", - None, - "4462ebfc5f915ef09cfbac6e7687a66e", - None, - ] - print(res["apical_points"]) - assert_allclose(res["apical_points"][0], [[-38.110016, 152.74228, -20.12851]]) - assert res["apical_points"][1] is None - - assert_allclose(res["apical_points"][3], [[-3.50943, 115.28183, 5.8795013]]) - assert res["apical_points"][4] is None - - assert_allclose(res["apical_points"][6], [[2.271698, 114.72031, -5.102951]]) - assert res["apical_points"][7] is None - - assert_allclose(res["apical_points"][9], [[12.257843, 85.439545, -19.165298]]) - assert res["apical_points"][10] is None - - # Check that the morphologies were not written - res_files = tmpdir.listdir() - assert len(res_files) == 5 - assert sorted(i.basename for i in res_files) == [ - "apical.yaml", - "apical_NRN_sections.yaml", - "axon_morphs.tsv", - "input_cells.mvd3", - "test_cells.mvd3", - ] - - -def run_with_mpi(): - """Test morphology synthesis with MPI.""" - # pylint: disable=import-outside-toplevel, too-many-locals, import-error - from data_factories import generate_axon_morph_tsv - from data_factories import generate_cell_collection - from data_factories import generate_cells_df - from data_factories import generate_input_cells - from data_factories import generate_mesh - from data_factories import generate_region_structure_boundary - from data_factories import generate_small_O1 - from data_factories import input_cells_path - from mpi4py import MPI - - from region_grower.utils import setup_logger - - COMM = MPI.COMM_WORLD # pylint: disable=c-extension-no-member - rank = COMM.Get_rank() - MASTER_RANK = 0 - is_master = rank == MASTER_RANK - - tmp_folder = Path("/tmp/test-run-synthesis_" + str(uuid4())) - tmp_folder = COMM.bcast(tmp_folder, root=MASTER_RANK) - input_cells = input_cells_path(tmp_folder) - small_O1_path = str(tmp_folder / "atlas") - region_structure_path = str(tmp_folder / "region_structure_boundary.yaml") - args = create_args( - True, - tmp_folder, - input_cells, - small_O1_path, - tmp_folder / "axon_morphs.tsv", - "apical_NRN_sections.yaml", - min_depth=25, - region_structure_boundary=region_structure_path, - ) - - setup_logger("info") - logging.getLogger("distributed").setLevel(logging.ERROR) - - if is_master: - tmp_folder.mkdir(exist_ok=True) - print(f"============= #{rank}: Create data") - cells_raw_data = generate_cells_df() - cell_collection = generate_cell_collection(cells_raw_data) - generate_input_cells(cell_collection, tmp_folder) - generate_small_O1(small_O1_path) - atlas = {"atlas": small_O1_path, "structure": DATA / "region_structure.yaml"} - generate_mesh(atlas, tmp_folder / "mesh.obj") - region_structure_path = generate_region_structure_boundary( - DATA / "region_structure.yaml", region_structure_path, str(tmp_folder / "mesh.obj") - ) - - generate_axon_morph_tsv(tmp_folder) - for dest in range(1, COMM.Get_size()): - req = COMM.isend("done", dest=dest) - else: - req = COMM.irecv(source=0) - req.wait() - - synthesizer = SynthesizeMorphologies(**args) - try: - print(f"============= #{rank}: Start synthesize") - synthesizer.synthesize() - - # Check results - print(f"============= #{rank}: Checking results") - expected_size = 12 - assert len(list(iter_morphology_files(tmp_folder))) == expected_size - check_yaml(DATA / "apical_boundary.yaml", args["out_apical"]) - finally: - # Clean the directory - print(f"============= #{rank}: Cleaning") - shutil.rmtree(tmp_folder) - - -if __name__ == "__main__": # pragma: no cover - run_with_mpi() diff --git a/tox.ini b/tox.ini index 5669993..ced097a 100644 --- a/tox.ini +++ b/tox.ini @@ -51,7 +51,7 @@ allowlist_externals = mpiexec commands = mpiexec -n 4 {envbindir}/python tests/test_synthesize_morphologies.py - mpiexec -n 4 {envbindir}/python tests/test_synthesize_morphologies_boundary.py + mpiexec -n 4 {envbindir}/python tests/test_synthesize_morphologies.py boundary [testenv:coverage] skip_install = true From b491ff4c68a4aec2d94288f07e28224f064b0c7f Mon Sep 17 00:00:00 2001 From: Adrien Berchet Date: Mon, 15 May 2023 19:26:02 +0000 Subject: [PATCH 34/63] Apply 1 suggestion(s) to 1 file(s) --- region_grower/context.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/region_grower/context.py b/region_grower/context.py index eed8bf9..ae75502 100644 --- a/region_grower/context.py +++ b/region_grower/context.py @@ -86,7 +86,7 @@ class SpaceContext: layer_depths: List cortical_depths: List boundaries: Dict = None - atlas_info: Dict = {} # voxel_dimnesions, offset and shape from atlas for indices conversion + atlas_info: Dict = {} # voxel_dimensions, offset and shape from atlas for indices conversion def indices_to_positions(self, indices): """Local version of voxcel's function to prevent atlas loading.""" From b7f033f59dffeb0ee9433bb26a55eac1ccc89e7e Mon Sep 17 00:00:00 2001 From: Adrien Berchet Date: Mon, 15 May 2023 19:28:15 +0000 Subject: [PATCH 35/63] Apply 1 suggestion(s) to 1 file(s) --- region_grower/atlas_helper.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/region_grower/atlas_helper.py b/region_grower/atlas_helper.py index 98c37fb..cdac2af 100644 --- a/region_grower/atlas_helper.py +++ b/region_grower/atlas_helper.py @@ -32,7 +32,7 @@ def __init__(self, atlas: Atlas, region_structure_path: str): with open(region_structure_path, "r", encoding="utf-8") as region_file: self.region_structure = yaml.safe_load(region_file) else: - raise ValueError(f"region_structure file at {region_structure_path} not found.") + raise ValueError(f"region_structure file not found at {region_structure_path}.") self.regions = list(self.region_structure.keys()) self.layers = {} From e4219290bf926c1cd5937fc57af9179049820b77 Mon Sep 17 00:00:00 2001 From: arnaudon Date: Tue, 16 May 2023 13:34:39 +0200 Subject: [PATCH 36/63] lint --- tests/test_synthesize_morphologies.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/tests/test_synthesize_morphologies.py b/tests/test_synthesize_morphologies.py index 6e12098..e41fdd6 100644 --- a/tests/test_synthesize_morphologies.py +++ b/tests/test_synthesize_morphologies.py @@ -1,10 +1,10 @@ """Test the region_grower.synthesize_morphologies module.""" # pylint: disable=missing-function-docstring import json -import sys import logging import os import shutil +import sys from copy import deepcopy from itertools import combinations from pathlib import Path @@ -27,7 +27,6 @@ from region_grower.synthesize_morphologies import RegionMapper from region_grower.synthesize_morphologies import SynthesizeMorphologies - DATA = Path(__file__).parent / "data" From e4a41661feaae0b6988e6c6b59e21347924954ff Mon Sep 17 00:00:00 2001 From: arnaudon Date: Thu, 22 Jun 2023 16:32:13 +0200 Subject: [PATCH 37/63] more --- region_grower/atlas_helper.py | 1 + region_grower/context.py | 52 ++++++++++++++---------- region_grower/synthesize_morphologies.py | 10 +++-- 3 files changed, 38 insertions(+), 25 deletions(-) diff --git a/region_grower/atlas_helper.py b/region_grower/atlas_helper.py index cdac2af..63ddbb1 100644 --- a/region_grower/atlas_helper.py +++ b/region_grower/atlas_helper.py @@ -31,6 +31,7 @@ def __init__(self, atlas: Atlas, region_structure_path: str): if region_structure_path is not None and Path(region_structure_path).exists(): with open(region_structure_path, "r", encoding="utf-8") as region_file: self.region_structure = yaml.safe_load(region_file) + self.region_structure_base_path = Path(region_structure_path).parent else: raise ValueError(f"region_structure file not found at {region_structure_path}.") diff --git a/region_grower/context.py b/region_grower/context.py index ae75502..825a928 100644 --- a/region_grower/context.py +++ b/region_grower/context.py @@ -85,7 +85,7 @@ class SpaceContext: layer_depths: List cortical_depths: List - boundaries: Dict = None + boundaries: List = None atlas_info: Dict = {} # voxel_dimensions, offset and shape from atlas for indices conversion def indices_to_positions(self, indices): @@ -105,16 +105,10 @@ def positions_to_indices(self, positions): result[result >= self.atlas_info["shape"]] = -1 return result - def get_boundaries(self): + def get_boundaries(self, mtype, neurite_types): """Returns a dict with boundaries data for NeuroTS.""" - # add here necessary logic to convert raw config data to NeuroTS context data - self.boundaries = json.loads(self.boundaries) - mesh = trimesh.load_mesh(self.boundaries["path"]) - - # to save points for plotting while growing to debug the prob function - debug = False - def get_distance_to_mesh(mesh, ray_origin, ray_direction): + def get_distance_to_mesh(mesh, ray_origin, ray_direction, debug=False): """Compute distances from point/directions to a mesh.""" vox_ray_origin = self.positions_to_indices(ray_origin) locations = mesh.ray.intersects_location( @@ -146,23 +140,35 @@ def _prob(dist, params): # TODO: other variants for this function could be included here return np.clip(p, 0, 1) - if "params_section" in self.boundaries: + # add here necessary logic to convert raw config data to NeuroTS context data + self.boundaries = json.loads(self.boundaries) + for i, boundary in enumerate(self.boundaries): + + mtypes = boundary.pop('mtypes', None) + if mtypes is not None and mtype not in mtypes: + continue + # TODO: add check a similar check on neurite_types + + mesh = trimesh.load_mesh(boundary["path"]) + + if "params_section" in boundary: - def section_prob(direction, current_point): - """Probability function for the sections.""" - dist = get_distance_to_mesh(mesh, current_point, direction) - return _prob(dist, self.boundaries["params_section"]) + def section_prob(direction, current_point): + """Probability function for the sections.""" + dist = get_distance_to_mesh(mesh, current_point, direction) + return _prob(dist, boundary["params_section"]) - self.boundaries["section_prob"] = section_prob + self.boundaries[i]["section_prob"] = section_prob - if "params_trunk" in self.boundaries: + if "params_trunk" in boundary: - def trunk_prob(direction, current_point): - """Probability function for the trunks.""" - dist = get_distance_to_mesh(mesh, current_point, direction) - return _prob(dist, self.boundaries["params_trunk"]) + def trunk_prob(direction, current_point): + """Probability function for the trunks.""" + dist = get_distance_to_mesh(mesh, current_point, direction) + return _prob(dist, boundary["params_trunk"]) + + self.boundaries[i]["trunk_prob"] = trunk_prob - self.boundaries["trunk_prob"] = trunk_prob return self.boundaries def layer_fraction_to_position(self, layer: int, layer_fraction: float) -> float: @@ -406,7 +412,9 @@ def _synthesize_once(self, rng) -> SynthesisResult: context = {"debug_data": self.debug_infos["input_scaling"]} if self.context.boundaries is not None: - context.update(self.context.get_boundaries()) + context["boundaries"] = self.context.get_boundaries( + mtype=self.cell.mtype, neurite_types=params["grow_types"] + ) grower = NeuronGrower( input_parameters=params, diff --git a/region_grower/synthesize_morphologies.py b/region_grower/synthesize_morphologies.py index 449431b..4225b12 100644 --- a/region_grower/synthesize_morphologies.py +++ b/region_grower/synthesize_morphologies.py @@ -527,9 +527,13 @@ def assign_atlas_data(self, min_depth=25, max_depth=5000): ) if "boundaries" in self.atlas.region_structure[_region]: - self.cells_data.loc[region_mask, "boundaries"] = json.dumps( - self.atlas.region_structure[_region]["boundaries"] - ) + boundaries = self.atlas.region_structure[_region]["boundaries"] + for boundary in boundaries: + if not Path(boundary["path"]).is_absolute(): + boundary["path"] = str( + self.atlas.region_structure_base_path / boundary["path"] + ) + self.cells_data.loc[region_mask, "boundaries"] = json.dumps(boundaries) self.cells_data.loc[region_mask, "atlas_info"] = json.dumps( { "voxel_dimensions": self.atlas.depths[_region].voxel_dimensions.tolist(), From bff52c8f6f5a7984a56e3c2d4269b59a4138de6f Mon Sep 17 00:00:00 2001 From: arnaudon Date: Thu, 22 Jun 2023 21:27:11 +0200 Subject: [PATCH 38/63] small fix --- region_grower/synthesize_morphologies.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/region_grower/synthesize_morphologies.py b/region_grower/synthesize_morphologies.py index 4225b12..0b82b89 100644 --- a/region_grower/synthesize_morphologies.py +++ b/region_grower/synthesize_morphologies.py @@ -531,7 +531,7 @@ def assign_atlas_data(self, min_depth=25, max_depth=5000): for boundary in boundaries: if not Path(boundary["path"]).is_absolute(): boundary["path"] = str( - self.atlas.region_structure_base_path / boundary["path"] + (self.atlas.region_structure_base_path / boundary["path"]).absolute() ) self.cells_data.loc[region_mask, "boundaries"] = json.dumps(boundaries) self.cells_data.loc[region_mask, "atlas_info"] = json.dumps( From b4d1b08d1977355076e31d93f392399fe23d3b1b Mon Sep 17 00:00:00 2001 From: arnaudon Date: Fri, 30 Jun 2023 11:51:33 +0200 Subject: [PATCH 39/63] more --- region_grower/context.py | 29 +++++++++++++++++++++-------- 1 file changed, 21 insertions(+), 8 deletions(-) diff --git a/region_grower/context.py b/region_grower/context.py index 825a928..a436685 100644 --- a/region_grower/context.py +++ b/region_grower/context.py @@ -108,14 +108,22 @@ def positions_to_indices(self, positions): def get_boundaries(self, mtype, neurite_types): """Returns a dict with boundaries data for NeuroTS.""" - def get_distance_to_mesh(mesh, ray_origin, ray_direction, debug=False): + def get_distance_to_mesh(mesh, ray_origin, ray_direction, mesh_type): """Compute distances from point/directions to a mesh.""" - vox_ray_origin = self.positions_to_indices(ray_origin) + debug = os.environ.get("REGION_GROWER_BOUNDARY_DEBUG", 0) + + if mesh_type == "voxel": + ray_origin = self.positions_to_indices(ray_origin) + locations = mesh.ray.intersects_location( - ray_origins=[vox_ray_origin], ray_directions=[ray_direction] + ray_origins=[ray_origin], ray_directions=[ray_direction] )[0] + if len(locations): - intersect = self.indices_to_positions(locations[0]) + intersect = locations[0] + if mesh_type == "voxel": + intersect = self.indices_to_positions(intersect) + dist = np.linalg.norm(intersect - ray_origin) if debug: # pragma: no cover @@ -124,11 +132,13 @@ def get_distance_to_mesh(mesh, ray_origin, ray_direction, debug=False): with open("data.csv", "a", encoding="utf8") as file: print(f"{x}, {y}, {z}, {vx}, {vy}, {vz}, {dist}", file=file) return dist + if debug: # pragma: no cover x, y, z = ray_origin vx, vy, vz = ray_direction with open("data.csv", "a", encoding="utf8") as file: print(f"{x}, {y}, {z}, {vx}, {vy}, {vz}, {-100}", file=file) + return np.inf # pragma: no cover def _prob(dist, params): @@ -143,8 +153,7 @@ def _prob(dist, params): # add here necessary logic to convert raw config data to NeuroTS context data self.boundaries = json.loads(self.boundaries) for i, boundary in enumerate(self.boundaries): - - mtypes = boundary.pop('mtypes', None) + mtypes = boundary.pop("mtypes", None) if mtypes is not None and mtype not in mtypes: continue # TODO: add check a similar check on neurite_types @@ -155,7 +164,9 @@ def _prob(dist, params): def section_prob(direction, current_point): """Probability function for the sections.""" - dist = get_distance_to_mesh(mesh, current_point, direction) + dist = get_distance_to_mesh( + mesh, current_point, direction, mesh_type=boundary.get("mesh_type", "voxel") + ) return _prob(dist, boundary["params_section"]) self.boundaries[i]["section_prob"] = section_prob @@ -164,7 +175,9 @@ def section_prob(direction, current_point): def trunk_prob(direction, current_point): """Probability function for the trunks.""" - dist = get_distance_to_mesh(mesh, current_point, direction) + dist = get_distance_to_mesh( + mesh, current_point, direction, mesh_type=boundary.get("mesh_type", "voxel") + ) return _prob(dist, boundary["params_trunk"]) self.boundaries[i]["trunk_prob"] = trunk_prob From 5765db532cf8c35b22d802e29bb8282982a3cda0 Mon Sep 17 00:00:00 2001 From: arnaudon Date: Tue, 5 Sep 2023 16:14:52 +0200 Subject: [PATCH 40/63] more --- region_grower/context.py | 48 +++++++++++++++++++----- region_grower/synthesize_morphologies.py | 5 +++ 2 files changed, 43 insertions(+), 10 deletions(-) diff --git a/region_grower/context.py b/region_grower/context.py index a436685..42fafe3 100644 --- a/region_grower/context.py +++ b/region_grower/context.py @@ -86,6 +86,7 @@ class SpaceContext: layer_depths: List cortical_depths: List boundaries: List = None + directions: Dict = None atlas_info: Dict = {} # voxel_dimensions, offset and shape from atlas for indices conversion def indices_to_positions(self, indices): @@ -105,7 +106,25 @@ def positions_to_indices(self, positions): result[result >= self.atlas_info["shape"]] = -1 return result - def get_boundaries(self, mtype, neurite_types): + def get_direction(self, mtype): + for i, direction in enumerate(self.directions): + target_direction = direction["params"]["direction"] + strength = direction["params"]["strength"] + + def section_prob(seg_direction, current_point): + """Probability function for the sections.""" + p = np.clip( + 1.0 - np.arccos(seg_direction.dot(target_direction)) * strength / np.pi, + 0, + 1, + ) + return p + + self.directions[i]["section_prob"] = section_prob + + return self.directions + + def get_boundaries(self, mtype): """Returns a dict with boundaries data for NeuroTS.""" def get_distance_to_mesh(mesh, ray_origin, ray_direction, mesh_type): @@ -151,12 +170,11 @@ def _prob(dist, params): return np.clip(p, 0, 1) # add here necessary logic to convert raw config data to NeuroTS context data - self.boundaries = json.loads(self.boundaries) - for i, boundary in enumerate(self.boundaries): + boundaries = json.loads(self.boundaries) + for i, boundary in enumerate(boundaries): mtypes = boundary.pop("mtypes", None) if mtypes is not None and mtype not in mtypes: continue - # TODO: add check a similar check on neurite_types mesh = trimesh.load_mesh(boundary["path"]) @@ -169,7 +187,7 @@ def section_prob(direction, current_point): ) return _prob(dist, boundary["params_section"]) - self.boundaries[i]["section_prob"] = section_prob + boundaries[i]["section_prob"] = section_prob if "params_trunk" in boundary: @@ -180,9 +198,9 @@ def trunk_prob(direction, current_point): ) return _prob(dist, boundary["params_trunk"]) - self.boundaries[i]["trunk_prob"] = trunk_prob + boundaries[i]["trunk_prob"] = trunk_prob - return self.boundaries + return boundaries def layer_fraction_to_position(self, layer: int, layer_fraction: float) -> float: """Returns an absolute position from a layer and a fraction of the layer. @@ -282,6 +300,12 @@ def __init__( def _correct_position_orientation_scaling(self, params_orig: Dict) -> Dict: """Return a copy of the parameter with the correct orientation and scaling.""" params = deepcopy(params_orig) + self.context.directions = json.loads(self.context.directions) + if self.context.directions is not None: + for i, direction in enumerate(self.context.directions): + self.context.directions[i]["params"]["direction"] = self.cell.lookup_orientation( + direction["params"]["direction"] + ) for neurite_type in params["grow_types"]: if isinstance(params[neurite_type]["orientation"], dict): @@ -424,10 +448,14 @@ def _synthesize_once(self, rng) -> SynthesisResult: axon_morph = None context = {"debug_data": self.debug_infos["input_scaling"]} + if self.context.directions is not None: + if "constraints" not in context: + context["constraints"] = [] + context["constraints"] += self.context.get_direction(mtype=self.cell.mtype) if self.context.boundaries is not None: - context["boundaries"] = self.context.get_boundaries( - mtype=self.cell.mtype, neurite_types=params["grow_types"] - ) + if "constraints" not in context: + context["constraints"] = [] + context["constraints"] += self.context.get_boundaries(mtype=self.cell.mtype) grower = NeuronGrower( input_parameters=params, diff --git a/region_grower/synthesize_morphologies.py b/region_grower/synthesize_morphologies.py index 0b82b89..79a7bea 100644 --- a/region_grower/synthesize_morphologies.py +++ b/region_grower/synthesize_morphologies.py @@ -126,6 +126,7 @@ def _parallel_wrapper( current_space_context = SpaceContext( layer_depths=row["layer_depths"], cortical_depths=row["cortical_depths"], + directions=row["directions"] if "directions" in row else None, boundaries=row["boundaries"] if "boundaries" in row else None, atlas_info=json.loads(row["atlas_info"]) if "atlas_info" in row else None, ) @@ -526,6 +527,10 @@ def assign_atlas_data(self, min_depth=25, max_depth=5000): dtype=object, ) + if "directions" in self.atlas.region_structure[_region]: + self.cells_data.loc[region_mask, "directions"] = json.dumps( + self.atlas.region_structure[_region]["directions"] + ) if "boundaries" in self.atlas.region_structure[_region]: boundaries = self.atlas.region_structure[_region]["boundaries"] for boundary in boundaries: From 2209b3f0f598d735acf2772c0ef189814a5cb06a Mon Sep 17 00:00:00 2001 From: arnaudon Date: Thu, 21 Sep 2023 15:23:17 +0200 Subject: [PATCH 41/63] more --- region_grower/context.py | 3 ++- region_grower/synthesize_morphologies.py | 6 +++--- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/region_grower/context.py b/region_grower/context.py index 42fafe3..14ffb97 100644 --- a/region_grower/context.py +++ b/region_grower/context.py @@ -142,6 +142,7 @@ def get_distance_to_mesh(mesh, ray_origin, ray_direction, mesh_type): intersect = locations[0] if mesh_type == "voxel": intersect = self.indices_to_positions(intersect) + ray_origin = self.indices_to_positions(ray_origin) dist = np.linalg.norm(intersect - ray_origin) @@ -300,8 +301,8 @@ def __init__( def _correct_position_orientation_scaling(self, params_orig: Dict) -> Dict: """Return a copy of the parameter with the correct orientation and scaling.""" params = deepcopy(params_orig) - self.context.directions = json.loads(self.context.directions) if self.context.directions is not None: + self.context.directions = json.loads(self.context.directions) for i, direction in enumerate(self.context.directions): self.context.directions[i]["params"]["direction"] = self.cell.lookup_orientation( direction["params"]["direction"] diff --git a/region_grower/synthesize_morphologies.py b/region_grower/synthesize_morphologies.py index 79a7bea..4f45a03 100644 --- a/region_grower/synthesize_morphologies.py +++ b/region_grower/synthesize_morphologies.py @@ -541,9 +541,9 @@ def assign_atlas_data(self, min_depth=25, max_depth=5000): self.cells_data.loc[region_mask, "boundaries"] = json.dumps(boundaries) self.cells_data.loc[region_mask, "atlas_info"] = json.dumps( { - "voxel_dimensions": self.atlas.depths[_region].voxel_dimensions.tolist(), - "offset": self.atlas.depths[_region].offset.tolist(), - "shape": self.atlas.depths[_region].shape, + "voxel_dimensions": self.atlas.brain_regions.voxel_dimensions.tolist(), + "offset": self.atlas.brain_regions.offset.tolist(), + "shape": self.atlas.brain_regions.shape, } ) From f4fd1a85983debb88ff3d0fa4a351e7f833206e5 Mon Sep 17 00:00:00 2001 From: arnaudon Date: Mon, 25 Sep 2023 15:32:20 +0200 Subject: [PATCH 42/63] more --- region_grower/context.py | 15 +++++++++------ region_grower/modify.py | 2 +- 2 files changed, 10 insertions(+), 7 deletions(-) diff --git a/region_grower/context.py b/region_grower/context.py index 14ffb97..6f3418e 100644 --- a/region_grower/context.py +++ b/region_grower/context.py @@ -110,15 +110,17 @@ def get_direction(self, mtype): for i, direction in enumerate(self.directions): target_direction = direction["params"]["direction"] strength = direction["params"]["strength"] + mode = direction["params"]["mode"] def section_prob(seg_direction, current_point): """Probability function for the sections.""" - p = np.clip( - 1.0 - np.arccos(seg_direction.dot(target_direction)) * strength / np.pi, - 0, - 1, - ) - return p + if mode == "parallel": + p = 1.0 - np.arccos(seg_direction.dot(target_direction)) * strength / np.pi + if mode == "perpendicular": + p = 1.0 - strength * abs( + np.arccos(seg_direction.dot(target_direction)) / np.pi * 2.0 - 1.0 + ) + return np.clip(p, 0, 1) self.directions[i]["section_prob"] = section_prob @@ -402,6 +404,7 @@ def synthesize(self) -> SynthesisResult: """ rng = np.random.default_rng(self.params.seed) + return self._synthesize_once(rng) for _ in range(self.internals.retries): try: return self._synthesize_once(rng) diff --git a/region_grower/modify.py b/region_grower/modify.py index 4bf42a0..618c786 100644 --- a/region_grower/modify.py +++ b/region_grower/modify.py @@ -134,7 +134,7 @@ def input_scaling( "with_debug_info": debug_info is not None, }, } - else: + elif neurite_type != "axon": if debug_info is not None: debug_info.update( { From 257503c9ca76b4a5c09a4fafc29803719acea0d7 Mon Sep 17 00:00:00 2001 From: arnaudon Date: Wed, 27 Sep 2023 15:51:00 +0200 Subject: [PATCH 43/63] more --- region_grower/context.py | 60 ++++++++++++++++-------- region_grower/synthesize_morphologies.py | 2 + 2 files changed, 43 insertions(+), 19 deletions(-) diff --git a/region_grower/context.py b/region_grower/context.py index 6f3418e..563ddc6 100644 --- a/region_grower/context.py +++ b/region_grower/context.py @@ -10,6 +10,7 @@ import json from collections import defaultdict from copy import deepcopy +from functools import partial from typing import Dict from typing import List from typing import Optional @@ -85,6 +86,8 @@ class SpaceContext: layer_depths: List cortical_depths: List + soma_position: float + soma_depth: float boundaries: List = None directions: Dict = None atlas_info: Dict = {} # voxel_dimensions, offset and shape from atlas for indices conversion @@ -107,13 +110,30 @@ def positions_to_indices(self, positions): return result def get_direction(self, mtype): + directions = [] for i, direction in enumerate(self.directions): - target_direction = direction["params"]["direction"] - strength = direction["params"]["strength"] - mode = direction["params"]["mode"] + mtypes = direction.pop("mtypes", None) + if mtypes is not None and mtype not in mtypes: + continue - def section_prob(seg_direction, current_point): + def section_prob(seg_direction, current_point, direction=None): """Probability function for the sections.""" + + target_direction = direction["params"]["direction"] + strength = direction["params"]["strength"] + mode = direction["params"]["mode"] + layers = direction["params"].get("layers", []) + if layers: + depth = self.soma_depth + (self.soma_position - current_point).dot( + direction["pia_direction"] + ) + in_layer = False + for layer in layers: + if self.layer_depths[layer - 1] < depth < self.layer_depths[layer]: + in_layer = True + if not in_layer: + return 1.0 + if mode == "parallel": p = 1.0 - np.arccos(seg_direction.dot(target_direction)) * strength / np.pi if mode == "perpendicular": @@ -122,9 +142,9 @@ def section_prob(seg_direction, current_point): ) return np.clip(p, 0, 1) - self.directions[i]["section_prob"] = section_prob - - return self.directions + direction["section_prob"] = partial(section_prob, direction=direction) + directions.append(direction) + return directions def get_boundaries(self, mtype): """Returns a dict with boundaries data for NeuroTS.""" @@ -173,35 +193,34 @@ def _prob(dist, params): return np.clip(p, 0, 1) # add here necessary logic to convert raw config data to NeuroTS context data - boundaries = json.loads(self.boundaries) - for i, boundary in enumerate(boundaries): + all_boundaries = json.loads(self.boundaries) + boundaries = [] + for i, boundary in enumerate(all_boundaries): mtypes = boundary.pop("mtypes", None) if mtypes is not None and mtype not in mtypes: continue mesh = trimesh.load_mesh(boundary["path"]) + mesh_type = boundary.get("mesh_type", "voxel") if "params_section" in boundary: - def section_prob(direction, current_point): + def section_prob(direction, current_point, mesh=None, mesh_type=None): """Probability function for the sections.""" - dist = get_distance_to_mesh( - mesh, current_point, direction, mesh_type=boundary.get("mesh_type", "voxel") - ) + dist = get_distance_to_mesh(mesh, current_point, direction, mesh_type=mesh_type) return _prob(dist, boundary["params_section"]) - boundaries[i]["section_prob"] = section_prob + boundary["section_prob"] = partial(section_prob, mesh=mesh, mesh_type=mesh_type) if "params_trunk" in boundary: - def trunk_prob(direction, current_point): + def trunk_prob(direction, current_point, mesh=None, mesh_type=None): """Probability function for the trunks.""" - dist = get_distance_to_mesh( - mesh, current_point, direction, mesh_type=boundary.get("mesh_type", "voxel") - ) + dist = get_distance_to_mesh(mesh, current_point, direction, mesh_type) return _prob(dist, boundary["params_trunk"]) - boundaries[i]["trunk_prob"] = trunk_prob + boundary["trunk_prob"] = partial(trunk_prob, mesh=mesh, mesh_type=mesh_type) + boundaries.append(boundary) return boundaries @@ -309,6 +328,9 @@ def _correct_position_orientation_scaling(self, params_orig: Dict) -> Dict: self.context.directions[i]["params"]["direction"] = self.cell.lookup_orientation( direction["params"]["direction"] ) + self.context.directions[i]["pia_direction"] = self.cell.lookup_orientation( + PIA_DIRECTION + ).tolist() for neurite_type in params["grow_types"]: if isinstance(params[neurite_type]["orientation"], dict): diff --git a/region_grower/synthesize_morphologies.py b/region_grower/synthesize_morphologies.py index 4f45a03..a872ed6 100644 --- a/region_grower/synthesize_morphologies.py +++ b/region_grower/synthesize_morphologies.py @@ -129,6 +129,8 @@ def _parallel_wrapper( directions=row["directions"] if "directions" in row else None, boundaries=row["boundaries"] if "boundaries" in row else None, atlas_info=json.loads(row["atlas_info"]) if "atlas_info" in row else None, + soma_position=current_cell.position, + soma_depth=row["current_depth"], ) axon_scale = row.get("axon_scale", None) From cf0d128d7c840ae2bae6f299d0c983a27d2a0aff Mon Sep 17 00:00:00 2001 From: arnaudon Date: Tue, 3 Oct 2023 14:49:00 +0200 Subject: [PATCH 44/63] more --- region_grower/context.py | 82 ++++++++++++++++++++++++++++------------ 1 file changed, 57 insertions(+), 25 deletions(-) diff --git a/region_grower/context.py b/region_grower/context.py index 563ddc6..2f5c04b 100644 --- a/region_grower/context.py +++ b/region_grower/context.py @@ -10,6 +10,7 @@ import json from collections import defaultdict from copy import deepcopy +from pathlib import Path from functools import partial from typing import Dict from typing import List @@ -120,7 +121,7 @@ def section_prob(seg_direction, current_point, direction=None): """Probability function for the sections.""" target_direction = direction["params"]["direction"] - strength = direction["params"]["strength"] + power = direction["params"].get("power", 1.0) mode = direction["params"]["mode"] layers = direction["params"].get("layers", []) if layers: @@ -135,11 +136,12 @@ def section_prob(seg_direction, current_point, direction=None): return 1.0 if mode == "parallel": - p = 1.0 - np.arccos(seg_direction.dot(target_direction)) * strength / np.pi + p = (1.0 - np.arccos(seg_direction.dot(target_direction)) / np.pi) ** power if mode == "perpendicular": - p = 1.0 - strength * abs( - np.arccos(seg_direction.dot(target_direction)) / np.pi * 2.0 - 1.0 - ) + p = ( + 1.0 + - abs(np.arccos(seg_direction.dot(target_direction)) / np.pi * 2.0 - 1.0) + ) ** power return np.clip(p, 0, 1) direction["section_prob"] = partial(section_prob, direction=direction) @@ -183,15 +185,6 @@ def get_distance_to_mesh(mesh, ray_origin, ray_direction, mesh_type): return np.inf # pragma: no cover - def _prob(dist, params): - """Probability function from distance to boundary.""" - p = (dist - params.get("d_min", 0)) / ( - params.get("d_max", 100) - params.get("d_min", 0) - ) ** params.get("power", 1.5) - - # TODO: other variants for this function could be included here - return np.clip(p, 0, 1) - # add here necessary logic to convert raw config data to NeuroTS context data all_boundaries = json.loads(self.boundaries) boundaries = [] @@ -200,26 +193,65 @@ def _prob(dist, params): if mtypes is not None and mtype not in mtypes: continue - mesh = trimesh.load_mesh(boundary["path"]) + meshes = [] + if Path(boundary["path"]).is_dir(): + for mesh_path in Path(boundary["path"]).iterdir(): + meshes.append(trimesh.load_mesh(mesh_path)) + distances = [ + trimesh.proximity.closest_point(mesh, [self.soma_position])[1][0] + for mesh in meshes + ] + mesh = meshes[np.argmin(distances)] + + else: + mesh = trimesh.load_mesh(boundary["path"]) + mesh_type = boundary.get("mesh_type", "voxel") + mode = boundary.get("mode", "repulsive") - if "params_section" in boundary: + if mode == "repulsive": - def section_prob(direction, current_point, mesh=None, mesh_type=None): + def prob(direction, current_point, mesh=None, mesh_type=None, params=None): """Probability function for the sections.""" dist = get_distance_to_mesh(mesh, current_point, direction, mesh_type=mesh_type) - return _prob(dist, boundary["params_section"]) + p = (dist - params.get("d_min", 0)) / ( + params.get("d_max", 100) - params.get("d_min", 0) + ) ** params.get("power", 1.5) + return np.clip(p, 0, 1) - boundary["section_prob"] = partial(section_prob, mesh=mesh, mesh_type=mesh_type) + elif mode == "attractive": - if "params_trunk" in boundary: + def prob(direction, current_point, mesh=None, mesh_type=None, params=None): + """Probability function for the sections.""" + closest_point = trimesh.proximity.closest_point(mesh, [current_point])[0][0] + + if mesh_type == "voxel": + closest_point = self.indices_to_positions(closest_point) + + mesh_direction = closest_point - current_point + distance = np.linalg.norm(mesh_direction) + if params.get("d_min", 0) < distance < params.get("d_max", 100): + if mesh_type == "voxel": + current_point = self.positions_to_indices(current_point) + p = ( + 1 - np.arccos(direction.dot(mesh_direction) / distance) / np.pi + ) ** params.get("power", 1.0) + return np.clip(p, 0, 1) + return 1.0 + + else: + raise ValueError(f"boundary mode {mode} not understood!") + + if "params_section" in boundary: + boundary["section_prob"] = partial( + prob, mesh=mesh, mesh_type=mesh_type, params=boundary["params_section"] + ) - def trunk_prob(direction, current_point, mesh=None, mesh_type=None): - """Probability function for the trunks.""" - dist = get_distance_to_mesh(mesh, current_point, direction, mesh_type) - return _prob(dist, boundary["params_trunk"]) + if "params_trunk" in boundary: + boundary["trunk_prob"] = partial( + prob, mesh=mesh, mesh_type=mesh_type, params=boundary["params_ trunk"] + ) - boundary["trunk_prob"] = partial(trunk_prob, mesh=mesh, mesh_type=mesh_type) boundaries.append(boundary) return boundaries From 20e5aa429cdf3369a56e037c79bff9e75ecc78a6 Mon Sep 17 00:00:00 2001 From: arnaudon Date: Tue, 10 Oct 2023 11:38:00 +0200 Subject: [PATCH 45/63] option to synthesize axons --- region_grower/cli.py | 6 ++++++ region_grower/synthesize_morphologies.py | 11 ++++++++++- 2 files changed, 16 insertions(+), 1 deletion(-) diff --git a/region_grower/cli.py b/region_grower/cli.py index 3c0f93f..b7842fe 100644 --- a/region_grower/cli.py +++ b/region_grower/cli.py @@ -284,6 +284,12 @@ def generate_distributions( is_flag=True, default=False, ) +@click.option( + "--synthesize-axons", + help="Display the versions of all the accessible modules in a logger entry", + is_flag=True, + default=False, +) def synthesize_morphologies(**kwargs): # pylint: disable=too-many-arguments, too-many-locals """Synthesize morphologies.""" if mpi_enabled and kwargs.get("with_mpi", False): # pragma: no cover diff --git a/region_grower/synthesize_morphologies.py b/region_grower/synthesize_morphologies.py index a872ed6..2c0d47f 100644 --- a/region_grower/synthesize_morphologies.py +++ b/region_grower/synthesize_morphologies.py @@ -201,6 +201,7 @@ class SynthesizeMorphologies: should be used to graft the axon on each synthesized morphology. base_morph_dir: the path containing the morphologies listed in the TSV file given in ``morph_axon``. + synthesize_axons: set to True to synthesize axons instead of grafting atlas_cache: the path to the directory used for the atlas cache. seed: the starting seed to use (note that the GID of each cell is added to this seed to ensure all cells has different seeds). @@ -254,6 +255,7 @@ def __init__( hide_progress_bar=False, dask_config=None, chunksize=None, + synthesize_axons=False, ): # pylint: disable=too-many-arguments, too-many-locals, too-many-statements self.seed = seed self.scaling_jitter_std = scaling_jitter_std @@ -303,6 +305,13 @@ def __init__( with open(tmd_distributions, "r", encoding="utf-8") as f: self.tmd_distributions = convert_from_legacy_neurite_type(json.load(f)) + for region, params in self.tmd_parameters.items(): + for param in params.values(): + if synthesize_axons: + assert "axon" in param["grow_types"] + elif "axon" in param["grow_types"]: + param["grow_types"].remove("axon") + # Set default values to tmd_parameters and tmd_distributions self.set_default_params_and_distrs() @@ -336,7 +345,7 @@ def __init__( LOGGER.info("Fetching atlas data from %s", atlas) self.assign_atlas_data(min_depth, max_depth) - if morph_axon is not None: + if morph_axon is not None and not synthesize_axons: LOGGER.info("Loading axon morphologies from %s", morph_axon) self.axon_morph_list = load_morphology_list(morph_axon, self.task_ids) check_na_morphologies( From a654639d7a5c89977a5a9bec476e2638ba1b8520 Mon Sep 17 00:00:00 2001 From: arnaudon Date: Wed, 18 Oct 2023 15:44:54 +0200 Subject: [PATCH 46/63] add option for barrel --- region_grower/context.py | 54 +++++++++++++++++++++++++--------------- 1 file changed, 34 insertions(+), 20 deletions(-) diff --git a/region_grower/context.py b/region_grower/context.py index 2f5c04b..2291b07 100644 --- a/region_grower/context.py +++ b/region_grower/context.py @@ -193,26 +193,39 @@ def get_distance_to_mesh(mesh, ray_origin, ray_direction, mesh_type): if mtypes is not None and mtype not in mtypes: continue + mesh_type = boundary.get("mesh_type", "voxel") + mode = boundary.get("mode", "repulsive") + meshes = [] if Path(boundary["path"]).is_dir(): - for mesh_path in Path(boundary["path"]).iterdir(): - meshes.append(trimesh.load_mesh(mesh_path)) - distances = [ - trimesh.proximity.closest_point(mesh, [self.soma_position])[1][0] - for mesh in meshes - ] - mesh = meshes[np.argmin(distances)] + soma_position = self.soma_position + if mesh_type == "voxel": + soma_position = self.positions_to_indices(soma_position) + + if boundary.get("multimesh_mode", "closest") == "closest": + for mesh_path in Path(boundary["path"]).iterdir(): + meshes.append(trimesh.load_mesh(mesh_path)) + distances = [ + trimesh.proximity.closest_point(mesh, [soma_position])[1][0] + for mesh in meshes + ] + mesh = meshes[np.argmin(distances)] + + if boundary.get("multimesh_mode", "closest") == "inside": + mesh = None + for mesh_path in Path(boundary["path"]).iterdir(): + _mesh = trimesh.load_mesh(mesh_path) + if _mesh.contains([soma_position])[0]: + mesh = _mesh + break else: mesh = trimesh.load_mesh(boundary["path"]) - mesh_type = boundary.get("mesh_type", "voxel") - mode = boundary.get("mode", "repulsive") - if mode == "repulsive": def prob(direction, current_point, mesh=None, mesh_type=None, params=None): - """Probability function for the sections.""" + """Probability function for repulsive mode.""" dist = get_distance_to_mesh(mesh, current_point, direction, mesh_type=mesh_type) p = (dist - params.get("d_min", 0)) / ( params.get("d_max", 100) - params.get("d_min", 0) @@ -222,7 +235,7 @@ def prob(direction, current_point, mesh=None, mesh_type=None, params=None): elif mode == "attractive": def prob(direction, current_point, mesh=None, mesh_type=None, params=None): - """Probability function for the sections.""" + """Probability function for attractive mode.""" closest_point = trimesh.proximity.closest_point(mesh, [current_point])[0][0] if mesh_type == "voxel": @@ -242,15 +255,16 @@ def prob(direction, current_point, mesh=None, mesh_type=None, params=None): else: raise ValueError(f"boundary mode {mode} not understood!") - if "params_section" in boundary: - boundary["section_prob"] = partial( - prob, mesh=mesh, mesh_type=mesh_type, params=boundary["params_section"] - ) + if mesh is not None: + if "params_section" in boundary: + boundary["section_prob"] = partial( + prob, mesh=mesh, mesh_type=mesh_type, params=boundary["params_section"] + ) - if "params_trunk" in boundary: - boundary["trunk_prob"] = partial( - prob, mesh=mesh, mesh_type=mesh_type, params=boundary["params_ trunk"] - ) + if "params_trunk" in boundary: + boundary["trunk_prob"] = partial( + prob, mesh=mesh, mesh_type=mesh_type, params=boundary["params_trunk"] + ) boundaries.append(boundary) From 6367e59e60c2c033fba3e90187c7d9f1c6c20ba1 Mon Sep 17 00:00:00 2001 From: arnaudon Date: Thu, 9 Nov 2023 15:24:41 +0100 Subject: [PATCH 47/63] more --- region_grower/context.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/region_grower/context.py b/region_grower/context.py index 2291b07..1f8a171 100644 --- a/region_grower/context.py +++ b/region_grower/context.py @@ -137,11 +137,13 @@ def section_prob(seg_direction, current_point, direction=None): if mode == "parallel": p = (1.0 - np.arccos(seg_direction.dot(target_direction)) / np.pi) ** power - if mode == "perpendicular": + elif mode == "perpendicular": p = ( 1.0 - abs(np.arccos(seg_direction.dot(target_direction)) / np.pi * 2.0 - 1.0) ) ** power + else: + raise ValueError(f"mode {mode} not understood for direction probability") return np.clip(p, 0, 1) direction["section_prob"] = partial(section_prob, direction=direction) @@ -228,7 +230,7 @@ def prob(direction, current_point, mesh=None, mesh_type=None, params=None): """Probability function for repulsive mode.""" dist = get_distance_to_mesh(mesh, current_point, direction, mesh_type=mesh_type) p = (dist - params.get("d_min", 0)) / ( - params.get("d_max", 100) - params.get("d_min", 0) + params.get("d_max", 100) - max(0, params.get("d_min", 0)) ) ** params.get("power", 1.5) return np.clip(p, 0, 1) @@ -246,6 +248,7 @@ def prob(direction, current_point, mesh=None, mesh_type=None, params=None): if params.get("d_min", 0) < distance < params.get("d_max", 100): if mesh_type == "voxel": current_point = self.positions_to_indices(current_point) + p = ( 1 - np.arccos(direction.dot(mesh_direction) / distance) / np.pi ) ** params.get("power", 1.0) From 8eb757d710c21844261affbc8a64b1159f85a4f4 Mon Sep 17 00:00:00 2001 From: arnaudon Date: Tue, 21 Nov 2023 13:25:41 +0100 Subject: [PATCH 48/63] more --- region_grower/context.py | 1 - region_grower/synthesize_morphologies.py | 14 ++++++++------ 2 files changed, 8 insertions(+), 7 deletions(-) diff --git a/region_grower/context.py b/region_grower/context.py index 1f8a171..028e177 100644 --- a/region_grower/context.py +++ b/region_grower/context.py @@ -475,7 +475,6 @@ def synthesize(self) -> SynthesisResult: """ rng = np.random.default_rng(self.params.seed) - return self._synthesize_once(rng) for _ in range(self.internals.retries): try: return self._synthesize_once(rng) diff --git a/region_grower/synthesize_morphologies.py b/region_grower/synthesize_morphologies.py index 2c0d47f..c0bd83f 100644 --- a/region_grower/synthesize_morphologies.py +++ b/region_grower/synthesize_morphologies.py @@ -308,7 +308,10 @@ def __init__( for region, params in self.tmd_parameters.items(): for param in params.values(): if synthesize_axons: - assert "axon" in param["grow_types"] + if "axon" not in param["grow_types"]: + LOGGER.warning( + "No axon data, but axon synthesis requested, you will not have axons" + ) elif "axon" in param["grow_types"]: param["grow_types"].remove("axon") @@ -368,7 +371,7 @@ def set_cortical_depths(self): for region in self.regions: if ( region not in self.atlas.region_structure - or self.atlas.region_structure[region]["thicknesses"] is None + or self.atlas.region_structure[region].get("thicknesses") is None ): # pragma: no cover self.cortical_depths[region] = self.cortical_depths["default"] else: @@ -507,8 +510,8 @@ def assign_atlas_data(self, min_depth=25, max_depth=5000): LOGGER.debug("Extract atlas data for %s region", _region) if ( _region in self.atlas.regions - and self.atlas.region_structure[_region].get("thicknesses", None) is not None - and self.atlas.region_structure[_region].get("layers", None) is not None + and self.atlas.region_structure[_region].get("thicknesses") is not None + and self.atlas.region_structure[_region].get("layers") is not None ): layers = self.atlas.layers[_region] thicknesses = [self.atlas.layer_thickness(layer) for layer in layers] @@ -519,8 +522,7 @@ def assign_atlas_data(self, min_depth=25, max_depth=5000): current_depths = np.clip(depths.lookup(positions), min_depth, max_depth) else: LOGGER.warning( - "We are not able to synthesize the region %s, we fallback to 'default' region", - _region, + "We are not able to synthesize the region %s with thicknesses", _region ) layer_depths = None current_depths = None From 047d1821ce94e72c770a145ace85d2be51da8696 Mon Sep 17 00:00:00 2001 From: arnaudon Date: Thu, 21 Dec 2023 14:55:45 +0100 Subject: [PATCH 49/63] more --- region_grower/context.py | 2 +- region_grower/synthesize_morphologies.py | 8 ++++---- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/region_grower/context.py b/region_grower/context.py index 028e177..b132d00 100644 --- a/region_grower/context.py +++ b/region_grower/context.py @@ -122,7 +122,7 @@ def section_prob(seg_direction, current_point, direction=None): target_direction = direction["params"]["direction"] power = direction["params"].get("power", 1.0) - mode = direction["params"]["mode"] + mode = direction["params"].get("mode", "parallel") layers = direction["params"].get("layers", []) if layers: depth = self.soma_depth + (self.soma_position - current_point).dot( diff --git a/region_grower/synthesize_morphologies.py b/region_grower/synthesize_morphologies.py index c0bd83f..6983798 100644 --- a/region_grower/synthesize_morphologies.py +++ b/region_grower/synthesize_morphologies.py @@ -486,11 +486,11 @@ def _init_parallel(self, mpi_only=False): def _close_parallel(self): if self._parallel_client is not None: LOGGER.debug("Closing the Dask client") - self._parallel_client.retire_workers() - time.sleep(1) - self._parallel_client.shutdown() + #self._parallel_client.retire_workers() + #time.sleep(1) + #self._parallel_client.shutdown() self._parallel_client.close() - self._parallel_client = None + #self._parallel_client = None def assign_atlas_data(self, min_depth=25, max_depth=5000): """Open an Atlas and compute depths and orientations according to the given positions.""" From 1e7d5dc624b79c9016120996b7b18761c44c794a Mon Sep 17 00:00:00 2001 From: arnaudon Date: Thu, 11 Jan 2024 14:40:35 +0100 Subject: [PATCH 50/63] save mesh_name --- region_grower/context.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/region_grower/context.py b/region_grower/context.py index b132d00..0e84427 100644 --- a/region_grower/context.py +++ b/region_grower/context.py @@ -205,13 +205,16 @@ def get_distance_to_mesh(mesh, ray_origin, ray_direction, mesh_type): soma_position = self.positions_to_indices(soma_position) if boundary.get("multimesh_mode", "closest") == "closest": - for mesh_path in Path(boundary["path"]).iterdir(): + mesh_paths = list(Path(boundary["path"]).iterdir()) + for mesh_path in mesh_paths: meshes.append(trimesh.load_mesh(mesh_path)) distances = [ trimesh.proximity.closest_point(mesh, [soma_position])[1][0] for mesh in meshes ] - mesh = meshes[np.argmin(distances)] + mesh_id = np.argmin(distances) + boundary["mesh_name"] = Path(mesh_paths[mesh_id]).stem + mesh = meshes[mesh_id] if boundary.get("multimesh_mode", "closest") == "inside": mesh = None @@ -530,6 +533,8 @@ def _synthesize_once(self, rng) -> SynthesisResult: if "constraints" not in context: context["constraints"] = [] context["constraints"] += self.context.get_boundaries(mtype=self.cell.mtype) + if context["constraints"] and "mesh_name" in context["constraints"][-1]: + self.debug_infos["mesh_name"] = context["constraints"][-1]["mesh_name"] grower = NeuronGrower( input_parameters=params, From f3c18751d9e790397f8d1ed70ae08cc5d9c72931 Mon Sep 17 00:00:00 2001 From: arnaudon Date: Fri, 12 Jan 2024 11:47:18 +0100 Subject: [PATCH 51/63] shuffle rows --- region_grower/synthesize_morphologies.py | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/region_grower/synthesize_morphologies.py b/region_grower/synthesize_morphologies.py index 6983798..9593665 100644 --- a/region_grower/synthesize_morphologies.py +++ b/region_grower/synthesize_morphologies.py @@ -6,11 +6,11 @@ - assign identity cell rotations to MVD3/sonata - optional axon grafting "on-the-fly" """ + import json import logging import os import subprocess -import time from pathlib import Path from shutil import which from typing import Optional @@ -486,11 +486,11 @@ def _init_parallel(self, mpi_only=False): def _close_parallel(self): if self._parallel_client is not None: LOGGER.debug("Closing the Dask client") - #self._parallel_client.retire_workers() - #time.sleep(1) - #self._parallel_client.shutdown() + # self._parallel_client.retire_workers() + # time.sleep(1) + # self._parallel_client.shutdown() self._parallel_client.close() - #self._parallel_client = None + # self._parallel_client = None def assign_atlas_data(self, min_depth=25, max_depth=5000): """Open an Atlas and compute depths and orientations according to the given positions.""" @@ -643,6 +643,9 @@ def compute(self): } ) + # shuffle rows to get more even loads on tasks + self.cells_data = self.cells_data.sample(frac=1.0).reset_index() + if self.nb_processes == 0: LOGGER.info("Start computation") computed = self.cells_data.apply( @@ -662,8 +665,7 @@ def compute(self): computed = future.compute() LOGGER.info("Format results") - res = self.cells_data.join(computed) - return res + return self.cells_data.join(computed).sort_values(by="index").set_index("index") def finalize(self, result: pd.DataFrame): """Finalize master work. From 28649f829a928c6e5bd63f6e552d78fc89106146 Mon Sep 17 00:00:00 2001 From: arnaudon Date: Thu, 1 Feb 2024 10:21:24 +0100 Subject: [PATCH 52/63] closest_y option --- region_grower/context.py | 25 ++++++++++++++++++++++++- 1 file changed, 24 insertions(+), 1 deletion(-) diff --git a/region_grower/context.py b/region_grower/context.py index 0e84427..e7493e0 100644 --- a/region_grower/context.py +++ b/region_grower/context.py @@ -216,6 +216,22 @@ def get_distance_to_mesh(mesh, ray_origin, ray_direction, mesh_type): boundary["mesh_name"] = Path(mesh_paths[mesh_id]).stem mesh = meshes[mesh_id] + if boundary.get("multimesh_mode", "closest") == "closest_y": + mesh_paths = list(Path(boundary["path"]).iterdir()) + meshes = [trimesh.load_mesh(mesh_path) for mesh_path in mesh_paths] + # this is assuming the original direction can be used + # as straight line, refinement would be to follow the vector field + distances = [ + np.linalg.norm( + np.cross(soma_position - mesh.center_mass, boundary["direction"]) + ) + for mesh in meshes + ] + + mesh_id = np.argmin(distances) + boundary["mesh_name"] = Path(mesh_paths[mesh_id]).stem + mesh = meshes[mesh_id] + if boundary.get("multimesh_mode", "closest") == "inside": mesh = None for mesh_path in Path(boundary["path"]).iterdir(): @@ -303,7 +319,7 @@ def lookup_target_reference_depths(self, depth: float) -> np.array: # we round to avoid being outside due to numerical precision if np.round(depth, 3) <= np.round(layer_depth, 3): return layer_depth, cortical_depth - + print(depth, self.layer_depths, self.cortical_depths) raise RegionGrowerError(f"Current depth ({depth}) is outside of circuit boundaries") def distance_to_constraint(self, depth: float, constraint: Dict) -> Optional[float]: @@ -384,6 +400,13 @@ def _correct_position_orientation_scaling(self, params_orig: Dict) -> Dict: PIA_DIRECTION ).tolist() + if self.context.boundaries is not None: + boundaries = json.loads(self.context.boundaries) + for boundary in boundaries: + if boundary.get("multimesh_mode", "closest") == "closest_y": + boundary["direction"] = self.cell.lookup_orientation(PIA_DIRECTION).tolist() + self.context.boundaries = json.dumps(boundaries) + for neurite_type in params["grow_types"]: if isinstance(params[neurite_type]["orientation"], dict): params["pia_direction"] = self.cell.lookup_orientation(PIA_DIRECTION).tolist() From 502f403d2b6da66248e12532d4be1cfe73bd59e7 Mon Sep 17 00:00:00 2001 From: arnaudon Date: Fri, 2 Feb 2024 09:34:35 +0100 Subject: [PATCH 53/63] remove neuroc for now --- setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.py b/setup.py index 3a8c636..c15ae74 100644 --- a/setup.py +++ b/setup.py @@ -20,7 +20,7 @@ "diameter-synthesis>=0.5.4,<1", "morphio>=3.3.6,<4", "morph-tool[nrn]>=2.9,<3", - "neuroc>=0.2.8,<1", + #"neuroc>=0.2.8,<1", "neurom>=3,<4", "neurots>=3.5,<4", "numpy>=1.22", From 556d6885305f671f0a0ea09b34c2e853e2cf8261 Mon Sep 17 00:00:00 2001 From: arnaudon Date: Fri, 2 Feb 2024 09:40:39 +0100 Subject: [PATCH 54/63] revert --- setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.py b/setup.py index c15ae74..3a8c636 100644 --- a/setup.py +++ b/setup.py @@ -20,7 +20,7 @@ "diameter-synthesis>=0.5.4,<1", "morphio>=3.3.6,<4", "morph-tool[nrn]>=2.9,<3", - #"neuroc>=0.2.8,<1", + "neuroc>=0.2.8,<1", "neurom>=3,<4", "neurots>=3.5,<4", "numpy>=1.22", From 83ac062a02f3f1f17eefd563eccab42cb463cd3c Mon Sep 17 00:00:00 2001 From: arnaudon Date: Fri, 1 Mar 2024 09:49:08 +0100 Subject: [PATCH 55/63] partial fix --- region_grower/context.py | 17 ++++++++++------- region_grower/synthesize_morphologies.py | 11 ++++++----- tests/test_synthesize_morphologies.py | 5 ++++- 3 files changed, 20 insertions(+), 13 deletions(-) diff --git a/region_grower/context.py b/region_grower/context.py index e7493e0..022d09d 100644 --- a/region_grower/context.py +++ b/region_grower/context.py @@ -6,16 +6,18 @@ TLDR: SpaceContext.synthesized() is being called by the placement_algorithm package to synthesize circuit morphologies. """ -import os + import json +import os from collections import defaultdict from copy import deepcopy -from pathlib import Path from functools import partial +from pathlib import Path from typing import Dict from typing import List from typing import Optional from typing import Union + import trimesh os.environ.setdefault("NEURON_MODULE_OPTIONS", "-nogui") # suppress NEURON warning @@ -111,15 +113,15 @@ def positions_to_indices(self, positions): return result def get_direction(self, mtype): + """Get direction data for the given mtype.""" directions = [] - for i, direction in enumerate(self.directions): + for direction in self.directions: mtypes = direction.pop("mtypes", None) if mtypes is not None and mtype not in mtypes: continue def section_prob(seg_direction, current_point, direction=None): """Probability function for the sections.""" - target_direction = direction["params"]["direction"] power = direction["params"].get("power", 1.0) mode = direction["params"].get("mode", "parallel") @@ -150,7 +152,9 @@ def section_prob(seg_direction, current_point, direction=None): directions.append(direction) return directions - def get_boundaries(self, mtype): + def get_boundaries( + self, mtype + ): # pylint: disable=too-many-locals,too-many-statements,too-many-branches """Returns a dict with boundaries data for NeuroTS.""" def get_distance_to_mesh(mesh, ray_origin, ray_direction, mesh_type): @@ -190,7 +194,7 @@ def get_distance_to_mesh(mesh, ray_origin, ray_direction, mesh_type): # add here necessary logic to convert raw config data to NeuroTS context data all_boundaries = json.loads(self.boundaries) boundaries = [] - for i, boundary in enumerate(all_boundaries): + for boundary in all_boundaries: mtypes = boundary.pop("mtypes", None) if mtypes is not None and mtype not in mtypes: continue @@ -319,7 +323,6 @@ def lookup_target_reference_depths(self, depth: float) -> np.array: # we round to avoid being outside due to numerical precision if np.round(depth, 3) <= np.round(layer_depth, 3): return layer_depth, cortical_depth - print(depth, self.layer_depths, self.cortical_depths) raise RegionGrowerError(f"Current depth ({depth}) is outside of circuit boundaries") def distance_to_constraint(self, depth: float, constraint: Dict) -> Optional[float]: diff --git a/region_grower/synthesize_morphologies.py b/region_grower/synthesize_morphologies.py index 9593665..96c8320 100644 --- a/region_grower/synthesize_morphologies.py +++ b/region_grower/synthesize_morphologies.py @@ -11,6 +11,7 @@ import logging import os import subprocess +import time from pathlib import Path from shutil import which from typing import Optional @@ -305,7 +306,7 @@ def __init__( with open(tmd_distributions, "r", encoding="utf-8") as f: self.tmd_distributions = convert_from_legacy_neurite_type(json.load(f)) - for region, params in self.tmd_parameters.items(): + for params in self.tmd_parameters.values(): for param in params.values(): if synthesize_axons: if "axon" not in param["grow_types"]: @@ -486,11 +487,11 @@ def _init_parallel(self, mpi_only=False): def _close_parallel(self): if self._parallel_client is not None: LOGGER.debug("Closing the Dask client") - # self._parallel_client.retire_workers() - # time.sleep(1) - # self._parallel_client.shutdown() + self._parallel_client.retire_workers() + time.sleep(1) + self._parallel_client.shutdown() self._parallel_client.close() - # self._parallel_client = None + self._parallel_client = None def assign_atlas_data(self, min_depth=25, max_depth=5000): """Open an Atlas and compute depths and orientations according to the given positions.""" diff --git a/tests/test_synthesize_morphologies.py b/tests/test_synthesize_morphologies.py index e41fdd6..e506fcf 100644 --- a/tests/test_synthesize_morphologies.py +++ b/tests/test_synthesize_morphologies.py @@ -178,6 +178,7 @@ def test_synthesize_dask_config( None, None, 100, + DATA / "region_structure.yaml", ) custom_scratch_config = str(tmp_folder / "custom_scratch_config") @@ -299,7 +300,7 @@ def test_synthesize_boundary( mesh, with_sections, with_trunks, -): # pylint: disable=unused-argument +): # pylint: disable=unused-argument,too-many-locals """Test morphology synthesis but skip write step.""" with_axon = True with_NRN = True @@ -716,6 +717,7 @@ def test_inconsistent_params( axon_morph_tsv, "apical_NRN_sections.yaml", min_depth, + DATA / "region_structure.yaml", ) args["out_morph_ext"] = ["h5"] @@ -748,6 +750,7 @@ def test_inconsistent_context( axon_morph_tsv, "apical_NRN_sections.yaml", min_depth, + DATA / "region_structure.yaml", ) args["nb_processes"] = 0 From 5fa9707367a776e9b0a7a83bf014bf25489c9ea1 Mon Sep 17 00:00:00 2001 From: arnaudon Date: Wed, 6 Mar 2024 15:23:07 +0100 Subject: [PATCH 56/63] fix --- region_grower/synthesize_morphologies.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/region_grower/synthesize_morphologies.py b/region_grower/synthesize_morphologies.py index 96c8320..9eb25e6 100644 --- a/region_grower/synthesize_morphologies.py +++ b/region_grower/synthesize_morphologies.py @@ -110,6 +110,7 @@ def inverse_mapper(self): def _parallel_wrapper( row, computation_parameters, + cortical_depths, rotational_jitter_std, scaling_jitter_std, min_hard_scale, @@ -126,7 +127,7 @@ def _parallel_wrapper( ) current_space_context = SpaceContext( layer_depths=row["layer_depths"], - cortical_depths=row["cortical_depths"], + cortical_depths=cortical_depths[row["synthesis_region"]], directions=row["directions"] if "directions" in row else None, boundaries=row["boundaries"] if "boundaries" in row else None, atlas_info=json.loads(row["atlas_info"]) if "atlas_info" in row else None, @@ -625,6 +626,7 @@ def compute(self): func_kwargs = { "computation_parameters": computation_parameters, "rotational_jitter_std": self.rotational_jitter_std, + "cortical_depths": self.cortical_depths, "scaling_jitter_std": self.scaling_jitter_std, "min_hard_scale": self.min_hard_scale, "tmd_parameters": self.tmd_parameters, From 0d4fcf6d570366517ba83d646746c64b32472126 Mon Sep 17 00:00:00 2001 From: arnaudon Date: Wed, 12 Jun 2024 10:44:39 +0200 Subject: [PATCH 57/63] resume --- region_grower/atlas_helper.py | 12 ++++++++++-- region_grower/cli.py | 6 ++++++ region_grower/morph_io.py | 21 +++++++++++++++++++-- region_grower/synthesize_morphologies.py | 22 +++++++++++++++++++--- 4 files changed, 54 insertions(+), 7 deletions(-) diff --git a/region_grower/atlas_helper.py b/region_grower/atlas_helper.py index 63ddbb1..e7d27d5 100644 --- a/region_grower/atlas_helper.py +++ b/region_grower/atlas_helper.py @@ -2,6 +2,7 @@ This helper allows simple lookups without having to reason in term of [PH][1-6] and [PH]y. """ + import operator from pathlib import Path from typing import List @@ -39,9 +40,16 @@ def __init__(self, atlas: Atlas, region_structure_path: str): self.layers = {} for region in self.regions: self.layers[region] = self.region_structure[region]["layers"] - self.y = atlas.load_data("[PH]y") + try: + self.y = atlas.load_data("[PH]y") + except: + self.y = None + self.brain_regions = atlas.load_data("brain_regions") - self.orientations = atlas.load_data("orientation", cls=OrientationField) + try: + self.orientations = atlas.load_data("orientation", cls=OrientationField) + except: + self.orientations = None def layer_thickness(self, layer: Union[int, str]) -> Atlas: """Returns an atlas of the layer thickness.""" diff --git a/region_grower/cli.py b/region_grower/cli.py index b7842fe..1bf414e 100644 --- a/region_grower/cli.py +++ b/region_grower/cli.py @@ -195,6 +195,12 @@ def generate_distributions( is_flag=True, default=False, ) +@click.option( + "--resume", + help="Resume synthesis, do not synthesize if file exists (default: False)", + is_flag=True, + default=False, +) @click.option( "--max-drop-ratio", help="Max drop ratio for any mtype (default: 0)", diff --git a/region_grower/morph_io.py b/region_grower/morph_io.py index 9106ee1..33f5545 100644 --- a/region_grower/morph_io.py +++ b/region_grower/morph_io.py @@ -1,4 +1,5 @@ """Module to load and write morphologies.""" + import os import random import uuid @@ -39,11 +40,12 @@ def get(self, name, options=None): class MorphWriter: """Helper class for writing morphologies.""" - def __init__(self, output_dir, file_exts, skip_write=False): + def __init__(self, output_dir, file_exts, skip_write=False, resume=False): self.output_dir = os.path.realpath(output_dir) self.file_exts = file_exts self._dir_depth = None self.skip_write = skip_write + self.resume = resume @staticmethod def _calc_dir_depth(num_files, max_files_per_dir=None): @@ -75,6 +77,9 @@ def prepare(self, num_files, max_files_per_dir=None, overwrite=False): - ensure it either does not exist, or is empty - if it does not exist, create an empty one """ + if self.resume: + overwrite = True + if self.skip_write: return self._dir_depth = MorphWriter._calc_dir_depth( @@ -117,5 +122,17 @@ def __call__(self, morph, seed): if not self.skip_write: for path in ext_paths: - morph.write(path) + if not (self.resume and path.exists()): + morph.write(path) return str(full_stem), ext_paths + + def check_resume(self, seed): + """Check if morphology exists if we are in resume mode, and return its path.""" + if not self.resume: + return False + + morph_name, subdirs = self.generate_name(seed) + full_stem = Path(subdirs, morph_name) + for ext_path in self.filepaths(full_stem): + if ext_path.exists(): + return ext_path diff --git a/region_grower/synthesize_morphologies.py b/region_grower/synthesize_morphologies.py index 9eb25e6..2021ce7 100644 --- a/region_grower/synthesize_morphologies.py +++ b/region_grower/synthesize_morphologies.py @@ -42,6 +42,7 @@ from region_grower.context import SpaceContext from region_grower.context import SpaceWorker from region_grower.context import SynthesisParameters +from region_grower.context import SynthesisResult from region_grower.morph_io import MorphLoader from region_grower.morph_io import MorphWriter from region_grower.utils import assign_morphologies @@ -159,7 +160,15 @@ def _parallel_wrapper( current_synthesis_parameters, computation_parameters, ) - new_cell = space_worker.synthesize() + resume_path = space_worker.internals.morph_writer.check_resume(row["seed"]) + if resume_path: + try: + new_cell = SynthesisResult(morphio.mut.Morphology(resume_path), [], []) + except morphio.RawDataError: + # in case it created an invalid morphology, just re-synthesize + new_cell = space_worker.synthesize() + else: + new_cell = space_worker.synthesize() res = space_worker.completion(new_cell) res["debug_infos"] = dict(space_worker.debug_infos) @@ -212,6 +221,7 @@ class SynthesizeMorphologies: max_files_per_dir: the maximum number of file in each directory (will create subdirectories if needed). overwrite: if set to False, the directory given to ``out_morph_dir`` must be empty. + resume: set to True to resume synthesis, do not synthesize if file exists. max_drop_ratio: the maximum ratio that scaling_jitter_std: the std of the scaling jitter. rotational_jitter_std: the std of the rotational jitter. @@ -242,6 +252,7 @@ def __init__( out_apical_nrn_sections=None, max_files_per_dir=None, overwrite=False, + resume=False, max_drop_ratio=0, scaling_jitter_std=None, rotational_jitter_std=None, @@ -324,7 +335,9 @@ def __init__( self.set_cortical_depths() LOGGER.info("Preparing morphology output folder in %s", out_morph_dir) - self.morph_writer = MorphWriter(out_morph_dir, out_morph_ext or ["swc"], skip_write) + self.morph_writer = MorphWriter( + out_morph_dir, out_morph_ext or ["swc"], skip_write, resume=resume + ) self.morph_writer.prepare( num_files=len(self.cells.positions), max_files_per_dir=max_files_per_dir, @@ -535,7 +548,10 @@ def assign_atlas_data(self, min_depth=25, max_depth=5000): ) LOGGER.debug("Extract orientations data for %s region", _region) - orientations = self.atlas.orientations.lookup(positions) + if self.atlas.orientations is not None: + orientations = self.atlas.orientations.lookup(positions) + else: + orientations = np.array(len(positions) * [np.eye(3)]) self.cells_data.loc[region_mask, "orientation"] = pd.Series( data=orientations.tolist(), index=self.cells_data.loc[region_mask].index, From d4ce5c2480ec0612e80bd720122554b940ae5591 Mon Sep 17 00:00:00 2001 From: arnaudon Date: Wed, 31 Jul 2024 15:23:14 +0200 Subject: [PATCH 58/63] added mclass --- region_grower/context.py | 21 +++++++++++++++------ region_grower/synthesize_morphologies.py | 9 +++++++++ 2 files changed, 24 insertions(+), 6 deletions(-) diff --git a/region_grower/context.py b/region_grower/context.py index 022d09d..bd81b59 100644 --- a/region_grower/context.py +++ b/region_grower/context.py @@ -74,6 +74,7 @@ class CellState: orientation: Matrix mtype: str depth: float + other_parameters: Dict def lookup_orientation(self, vector: Optional[Point] = None) -> np.array: """Returns the looked-up orientation for the given position. @@ -112,12 +113,16 @@ def positions_to_indices(self, positions): result[result >= self.atlas_info["shape"]] = -1 return result - def get_direction(self, mtype): + def get_direction(self, cell): """Get direction data for the given mtype.""" directions = [] for direction in self.directions: mtypes = direction.pop("mtypes", None) - if mtypes is not None and mtype not in mtypes: + if mtypes is not None and cell.mtype not in mtypes: + continue + # specific for barrel cortex + mclasses = direction.pop("mclass", None) + if mclasses is not None and cell.other_parameters["mclass"] not in mclasses: continue def section_prob(seg_direction, current_point, direction=None): @@ -153,7 +158,7 @@ def section_prob(seg_direction, current_point, direction=None): return directions def get_boundaries( - self, mtype + self, cell ): # pylint: disable=too-many-locals,too-many-statements,too-many-branches """Returns a dict with boundaries data for NeuroTS.""" @@ -196,7 +201,11 @@ def get_distance_to_mesh(mesh, ray_origin, ray_direction, mesh_type): boundaries = [] for boundary in all_boundaries: mtypes = boundary.pop("mtypes", None) - if mtypes is not None and mtype not in mtypes: + if mtypes is not None and cell.mtype not in mtypes: + continue + # specific for barrel cortex + mclasses = boundary.pop("mclass", None) + if mclasses is not None and cell.other_parameters["mclass"] not in mclasses: continue mesh_type = boundary.get("mesh_type", "voxel") @@ -554,11 +563,11 @@ def _synthesize_once(self, rng) -> SynthesisResult: if self.context.directions is not None: if "constraints" not in context: context["constraints"] = [] - context["constraints"] += self.context.get_direction(mtype=self.cell.mtype) + context["constraints"] += self.context.get_direction(self.cell) if self.context.boundaries is not None: if "constraints" not in context: context["constraints"] = [] - context["constraints"] += self.context.get_boundaries(mtype=self.cell.mtype) + context["constraints"] += self.context.get_boundaries(self.cell) if context["constraints"] and "mesh_name" in context["constraints"][-1]: self.debug_infos["mesh_name"] = context["constraints"][-1]["mesh_name"] diff --git a/region_grower/synthesize_morphologies.py b/region_grower/synthesize_morphologies.py index 2021ce7..ed7725f 100644 --- a/region_grower/synthesize_morphologies.py +++ b/region_grower/synthesize_morphologies.py @@ -125,6 +125,11 @@ def _parallel_wrapper( orientation=np.array([row["orientation"]]), mtype=row["mtype"], depth=row["current_depth"], + other_parameters={ + p: row[p] + for p in row.keys() + if p not in ["x", "y", "z", "orientation", "mtype", "current_depth"] + }, ) current_space_context = SpaceContext( layer_depths=row["layer_depths"], @@ -309,6 +314,10 @@ def __init__( LOGGER.info("Loading CellCollection from %s", input_cells) self.cells = CellCollection.load(input_cells) + df = self.cells.as_dataframe() + df = df[df.mtype == "L4_SSC"].reset_index() + df.index += 1 + self.cells = CellCollection.from_dataframe(df) LOGGER.info("Loading TMD parameters from %s", tmd_parameters) with open(tmd_parameters, "r", encoding="utf-8") as f: From 2537acdd30b833f5b209bfbdd8535d817ca974ac Mon Sep 17 00:00:00 2001 From: arnaudon Date: Fri, 2 Aug 2024 11:45:32 +0200 Subject: [PATCH 59/63] fix --- region_grower/context.py | 4 ++-- region_grower/synthesize_morphologies.py | 4 ---- 2 files changed, 2 insertions(+), 6 deletions(-) diff --git a/region_grower/context.py b/region_grower/context.py index bd81b59..3f3d799 100644 --- a/region_grower/context.py +++ b/region_grower/context.py @@ -122,7 +122,7 @@ def get_direction(self, cell): continue # specific for barrel cortex mclasses = direction.pop("mclass", None) - if mclasses is not None and cell.other_parameters["mclass"] not in mclasses: + if mclasses is not None and cell.other_parameters.get("mclass", None) not in mclasses: continue def section_prob(seg_direction, current_point, direction=None): @@ -205,7 +205,7 @@ def get_distance_to_mesh(mesh, ray_origin, ray_direction, mesh_type): continue # specific for barrel cortex mclasses = boundary.pop("mclass", None) - if mclasses is not None and cell.other_parameters["mclass"] not in mclasses: + if mclasses is not None and cell.other_parameters.get("mclass", None) not in mclasses: continue mesh_type = boundary.get("mesh_type", "voxel") diff --git a/region_grower/synthesize_morphologies.py b/region_grower/synthesize_morphologies.py index ed7725f..c0116c3 100644 --- a/region_grower/synthesize_morphologies.py +++ b/region_grower/synthesize_morphologies.py @@ -314,10 +314,6 @@ def __init__( LOGGER.info("Loading CellCollection from %s", input_cells) self.cells = CellCollection.load(input_cells) - df = self.cells.as_dataframe() - df = df[df.mtype == "L4_SSC"].reset_index() - df.index += 1 - self.cells = CellCollection.from_dataframe(df) LOGGER.info("Loading TMD parameters from %s", tmd_parameters) with open(tmd_parameters, "r", encoding="utf-8") as f: From 3c5ce3be8e82afee4b11b794f8b3cdafdde88f1f Mon Sep 17 00:00:00 2001 From: arnaudon Date: Fri, 13 Sep 2024 08:26:43 +0200 Subject: [PATCH 60/63] mob update --- region_grower/context.py | 43 +++++++++++++++++++----- region_grower/synthesize_morphologies.py | 29 ++++++---------- 2 files changed, 45 insertions(+), 27 deletions(-) diff --git a/region_grower/context.py b/region_grower/context.py index 3f3d799..0872e8e 100644 --- a/region_grower/context.py +++ b/region_grower/context.py @@ -40,6 +40,7 @@ from neurots import NeuroTSError # noqa: E402 ; pylint: disable=C0413 from neurots.utils import PIA_DIRECTION # noqa: E402 ; pylint: disable=C0413 from voxcell.cell_collection import CellCollection # noqa: E402 ; pylint: disable=C0413 +from voxcell.voxel_data import OrientationField from region_grower import RegionGrowerError # noqa: E402 ; pylint: disable=C0413 from region_grower import SkipSynthesisError # noqa: E402 ; pylint: disable=C0413 @@ -95,6 +96,7 @@ class SpaceContext: boundaries: List = None directions: Dict = None atlas_info: Dict = {} # voxel_dimensions, offset and shape from atlas for indices conversion + _orientations = None def indices_to_positions(self, indices): """Local version of voxcel's function to prevent atlas loading.""" @@ -127,7 +129,17 @@ def get_direction(self, cell): def section_prob(seg_direction, current_point, direction=None): """Probability function for the sections.""" - target_direction = direction["params"]["direction"] + if self._orientations is None: + self._orientations = OrientationField.load_nrrd( + self.atlas_info["direction_nrrd_path"] + ) + try: + target_direction = np.dot( + self._orientations.lookup(current_point), direction["params"]["direction"] + )[0] + except: + # if we are outside, we do nothing + return 1.0 power = direction["params"].get("power", 1.0) mode = direction["params"].get("mode", "parallel") layers = direction["params"].get("layers", []) @@ -149,8 +161,13 @@ def section_prob(seg_direction, current_point, direction=None): 1.0 - abs(np.arccos(seg_direction.dot(target_direction)) / np.pi * 2.0 - 1.0) ) ** power + else: raise ValueError(f"mode {mode} not understood for direction probability") + + if np.isnan(p): + # sometimes we get nan from this computation, to check why + return 1.0 return np.clip(p, 0, 1) direction["section_prob"] = partial(section_prob, direction=direction) @@ -253,6 +270,12 @@ def get_distance_to_mesh(mesh, ray_origin, ray_direction, mesh_type): mesh = _mesh break + if boundary.get("multimesh_mode", "closest") == "territories": + mesh = trimesh.load_mesh( + Path(boundary["path"]) + / f"glomeruli_{int(cell.other_parameters['glomerulus_id']):03d}.obj" + ) + else: mesh = trimesh.load_mesh(boundary["path"]) @@ -270,17 +293,22 @@ def prob(direction, current_point, mesh=None, mesh_type=None, params=None): def prob(direction, current_point, mesh=None, mesh_type=None, params=None): """Probability function for attractive mode.""" - closest_point = trimesh.proximity.closest_point(mesh, [current_point])[0][0] + if mesh_type == "voxel": + current_point_id = self.positions_to_indices(current_point) + else: + current_point_id = current_point + + if mesh.contains([current_point_id])[0]: + # when we get inside the mesh, we are free + return 1.0 + closest_point = trimesh.proximity.closest_point(mesh, [current_point_id])[0][0] if mesh_type == "voxel": closest_point = self.indices_to_positions(closest_point) mesh_direction = closest_point - current_point distance = np.linalg.norm(mesh_direction) if params.get("d_min", 0) < distance < params.get("d_max", 100): - if mesh_type == "voxel": - current_point = self.positions_to_indices(current_point) - p = ( 1 - np.arccos(direction.dot(mesh_direction) / distance) / np.pi ) ** params.get("power", 1.0) @@ -302,7 +330,6 @@ def prob(direction, current_point, mesh=None, mesh_type=None, params=None): ) boundaries.append(boundary) - return boundaries def layer_fraction_to_position(self, layer: int, layer_fraction: float) -> float: @@ -405,9 +432,6 @@ def _correct_position_orientation_scaling(self, params_orig: Dict) -> Dict: if self.context.directions is not None: self.context.directions = json.loads(self.context.directions) for i, direction in enumerate(self.context.directions): - self.context.directions[i]["params"]["direction"] = self.cell.lookup_orientation( - direction["params"]["direction"] - ) self.context.directions[i]["pia_direction"] = self.cell.lookup_orientation( PIA_DIRECTION ).tolist() @@ -417,6 +441,7 @@ def _correct_position_orientation_scaling(self, params_orig: Dict) -> Dict: for boundary in boundaries: if boundary.get("multimesh_mode", "closest") == "closest_y": boundary["direction"] = self.cell.lookup_orientation(PIA_DIRECTION).tolist() + self.context.boundaries = json.dumps(boundaries) for neurite_type in params["grow_types"]: diff --git a/region_grower/synthesize_morphologies.py b/region_grower/synthesize_morphologies.py index c0116c3..c9f83d3 100644 --- a/region_grower/synthesize_morphologies.py +++ b/region_grower/synthesize_morphologies.py @@ -528,24 +528,11 @@ def assign_atlas_data(self, min_depth=25, max_depth=5000): positions = self.cells.positions[region_mask] LOGGER.debug("Extract atlas data for %s region", _region) - if ( - _region in self.atlas.regions - and self.atlas.region_structure[_region].get("thicknesses") is not None - and self.atlas.region_structure[_region].get("layers") is not None - ): - layers = self.atlas.layers[_region] - thicknesses = [self.atlas.layer_thickness(layer) for layer in layers] - depths = self.atlas.compute_region_depth(_region) - layer_depths = self.atlas.get_layer_boundary_depths( - positions, thicknesses - ).T.tolist() - current_depths = np.clip(depths.lookup(positions), min_depth, max_depth) - else: - LOGGER.warning( - "We are not able to synthesize the region %s with thicknesses", _region - ) - layer_depths = None - current_depths = None + layers = self.atlas.layers[_region] + thicknesses = [self.atlas.layer_thickness(layer) for layer in layers] + depths = self.atlas.compute_region_depth(_region) + layer_depths = self.atlas.get_layer_boundary_depths(positions, thicknesses).T.tolist() + current_depths = np.clip(depths.lookup(positions), min_depth, max_depth) self.cells_data.loc[region_mask, "current_depth"] = current_depths self.cells_data.loc[region_mask, "layer_depths"] = pd.Series( @@ -574,12 +561,18 @@ def assign_atlas_data(self, min_depth=25, max_depth=5000): boundary["path"] = str( (self.atlas.region_structure_base_path / boundary["path"]).absolute() ) + if boundary.get("multimesh_mode", "closest") == "territories": + territories = self.atlas.atlas.load_data("glomerular_territories") + pos = self.cells_data.loc[region_mask, ["x", "y", "z"]].to_numpy() + self.cells_data.loc[region_mask, "glomerulus_id"] = territories.lookup(pos) + self.cells_data.loc[region_mask, "boundaries"] = json.dumps(boundaries) self.cells_data.loc[region_mask, "atlas_info"] = json.dumps( { "voxel_dimensions": self.atlas.brain_regions.voxel_dimensions.tolist(), "offset": self.atlas.brain_regions.offset.tolist(), "shape": self.atlas.brain_regions.shape, + "direction_nrrd_path": self.atlas.atlas.fetch_data("orientation"), } ) From d3fac7f97c3dc1bc8b7be5dab6d73c3317b7fa99 Mon Sep 17 00:00:00 2001 From: arnaudon Date: Fri, 13 Sep 2024 11:22:20 +0200 Subject: [PATCH 61/63] more --- region_grower/context.py | 13 ++++++++----- region_grower/synthesize_morphologies.py | 4 ++++ 2 files changed, 12 insertions(+), 5 deletions(-) diff --git a/region_grower/context.py b/region_grower/context.py index 0872e8e..ec3fe1a 100644 --- a/region_grower/context.py +++ b/region_grower/context.py @@ -41,6 +41,7 @@ from neurots.utils import PIA_DIRECTION # noqa: E402 ; pylint: disable=C0413 from voxcell.cell_collection import CellCollection # noqa: E402 ; pylint: disable=C0413 from voxcell.voxel_data import OrientationField +from voxcell import VoxcellError from region_grower import RegionGrowerError # noqa: E402 ; pylint: disable=C0413 from region_grower import SkipSynthesisError # noqa: E402 ; pylint: disable=C0413 @@ -96,7 +97,7 @@ class SpaceContext: boundaries: List = None directions: Dict = None atlas_info: Dict = {} # voxel_dimensions, offset and shape from atlas for indices conversion - _orientations = None + orientations = None def indices_to_positions(self, indices): """Local version of voxcel's function to prevent atlas loading.""" @@ -129,21 +130,22 @@ def get_direction(self, cell): def section_prob(seg_direction, current_point, direction=None): """Probability function for the sections.""" - if self._orientations is None: - self._orientations = OrientationField.load_nrrd( + if self.orientations is None: + self.orientations = OrientationField.load_nrrd( self.atlas_info["direction_nrrd_path"] ) try: target_direction = np.dot( - self._orientations.lookup(current_point), direction["params"]["direction"] + self.orientations.lookup(current_point), direction["params"]["direction"] )[0] - except: + except VoxcellError: # if we are outside, we do nothing return 1.0 power = direction["params"].get("power", 1.0) mode = direction["params"].get("mode", "parallel") layers = direction["params"].get("layers", []) if layers: + # this does not work well in curved atlas depth = self.soma_depth + (self.soma_position - current_point).dot( direction["pia_direction"] ) @@ -168,6 +170,7 @@ def section_prob(seg_direction, current_point, direction=None): if np.isnan(p): # sometimes we get nan from this computation, to check why return 1.0 + return np.clip(p, 0, 1) direction["section_prob"] = partial(section_prob, direction=direction) diff --git a/region_grower/synthesize_morphologies.py b/region_grower/synthesize_morphologies.py index c9f83d3..cd7cd84 100644 --- a/region_grower/synthesize_morphologies.py +++ b/region_grower/synthesize_morphologies.py @@ -567,6 +567,10 @@ def assign_atlas_data(self, min_depth=25, max_depth=5000): self.cells_data.loc[region_mask, "glomerulus_id"] = territories.lookup(pos) self.cells_data.loc[region_mask, "boundaries"] = json.dumps(boundaries) + if ( + "directions" in self.atlas.region_structure[_region] + or "boundaries" in self.atlas.region_structure[_region] + ): self.cells_data.loc[region_mask, "atlas_info"] = json.dumps( { "voxel_dimensions": self.atlas.brain_regions.voxel_dimensions.tolist(), From d9a3b5713e61fff0f6d62ff34cb170ff7522b2a3 Mon Sep 17 00:00:00 2001 From: arnaudon Date: Fri, 13 Sep 2024 14:21:58 +0200 Subject: [PATCH 62/63] more --- region_grower/context.py | 33 +++++++++++++++++------- region_grower/synthesize_morphologies.py | 5 ++-- 2 files changed, 26 insertions(+), 12 deletions(-) diff --git a/region_grower/context.py b/region_grower/context.py index ec3fe1a..2415a1e 100644 --- a/region_grower/context.py +++ b/region_grower/context.py @@ -120,11 +120,12 @@ def get_direction(self, cell): """Get direction data for the given mtype.""" directions = [] for direction in self.directions: - mtypes = direction.pop("mtypes", None) + mtypes = direction.get("mtypes", None) if mtypes is not None and cell.mtype not in mtypes: continue + # specific for barrel cortex - mclasses = direction.pop("mclass", None) + mclasses = direction.get("mclass", None) if mclasses is not None and cell.other_parameters.get("mclass", None) not in mclasses: continue @@ -220,11 +221,12 @@ def get_distance_to_mesh(mesh, ray_origin, ray_direction, mesh_type): all_boundaries = json.loads(self.boundaries) boundaries = [] for boundary in all_boundaries: - mtypes = boundary.pop("mtypes", None) + mtypes = boundary.get("mtypes", None) if mtypes is not None and cell.mtype not in mtypes: continue + # specific for barrel cortex - mclasses = boundary.pop("mclass", None) + mclasses = boundary.get("mclass", None) if mclasses is not None and cell.other_parameters.get("mclass", None) not in mclasses: continue @@ -274,10 +276,12 @@ def get_distance_to_mesh(mesh, ray_origin, ray_direction, mesh_type): break if boundary.get("multimesh_mode", "closest") == "territories": - mesh = trimesh.load_mesh( - Path(boundary["path"]) - / f"glomeruli_{int(cell.other_parameters['glomerulus_id']):03d}.obj" - ) + glom_id = int(cell.other_parameters["glomerulus_id"]) + mesh = None + if glom_id >= 0: + mesh = trimesh.load_mesh( + Path(boundary["path"]) / f"glomeruli_{glom_id:03d}.obj" + ) else: mesh = trimesh.load_mesh(boundary["path"]) @@ -433,7 +437,8 @@ def _correct_position_orientation_scaling(self, params_orig: Dict) -> Dict: """Return a copy of the parameter with the correct orientation and scaling.""" params = deepcopy(params_orig) if self.context.directions is not None: - self.context.directions = json.loads(self.context.directions) + if isinstance(self.context.directions, str): + self.context.directions = json.loads(self.context.directions) for i, direction in enumerate(self.context.directions): self.context.directions[i]["pia_direction"] = self.cell.lookup_orientation( PIA_DIRECTION @@ -544,7 +549,10 @@ def synthesize(self) -> SynthesisResult: for _ in range(self.internals.retries): try: return self._synthesize_once(rng) - except NeuroTSError: + except ( + NeuroTSError, + ValueError, + ): # value error is raised when nan in orientation results in soma generation error pass raise SkipSynthesisError( @@ -565,6 +573,11 @@ def completion(self, synthesized): synthesized.apical_points, morph_path ) + # enforce freeing memory of orientation nrrd data + if hasattr(self.context, "orientations"): + if self.context.orientations is not None: + del self.context.orientations + return { "name": morph_name, "apical_points": synthesized.apical_points, diff --git a/region_grower/synthesize_morphologies.py b/region_grower/synthesize_morphologies.py index cd7cd84..a773188 100644 --- a/region_grower/synthesize_morphologies.py +++ b/region_grower/synthesize_morphologies.py @@ -175,7 +175,6 @@ def _parallel_wrapper( else: new_cell = space_worker.synthesize() res = space_worker.completion(new_cell) - res["debug_infos"] = dict(space_worker.debug_infos) except (SkipSynthesisError, RegionGrowerError, NoDendriteException) as exc: # pragma: no cover LOGGER.error("Skip %s because of the following error: %s", row.name, exc) @@ -564,7 +563,9 @@ def assign_atlas_data(self, min_depth=25, max_depth=5000): if boundary.get("multimesh_mode", "closest") == "territories": territories = self.atlas.atlas.load_data("glomerular_territories") pos = self.cells_data.loc[region_mask, ["x", "y", "z"]].to_numpy() - self.cells_data.loc[region_mask, "glomerulus_id"] = territories.lookup(pos) + self.cells_data.loc[region_mask, "glomerulus_id"] = territories.lookup( + pos, outer_value=-1 + ) self.cells_data.loc[region_mask, "boundaries"] = json.dumps(boundaries) if ( From 9b3812c134a01324967c6127903525564dc845fd Mon Sep 17 00:00:00 2001 From: arnaudon Date: Tue, 24 Sep 2024 11:44:16 +0200 Subject: [PATCH 63/63] keep track of meshes --- region_grower/context.py | 33 +++++++++++++++++++++------------ 1 file changed, 21 insertions(+), 12 deletions(-) diff --git a/region_grower/context.py b/region_grower/context.py index 2415a1e..07f7d40 100644 --- a/region_grower/context.py +++ b/region_grower/context.py @@ -98,6 +98,7 @@ class SpaceContext: directions: Dict = None atlas_info: Dict = {} # voxel_dimensions, offset and shape from atlas for indices conversion orientations = None + mesh = None def indices_to_positions(self, indices): """Local version of voxcel's function to prevent atlas loading.""" @@ -207,6 +208,7 @@ def get_distance_to_mesh(mesh, ray_origin, ray_direction, mesh_type): vx, vy, vz = ray_direction with open("data.csv", "a", encoding="utf8") as file: print(f"{x}, {y}, {z}, {vx}, {vy}, {vz}, {dist}", file=file) + return dist if debug: # pragma: no cover @@ -249,7 +251,7 @@ def get_distance_to_mesh(mesh, ray_origin, ray_direction, mesh_type): ] mesh_id = np.argmin(distances) boundary["mesh_name"] = Path(mesh_paths[mesh_id]).stem - mesh = meshes[mesh_id] + self.mesh = meshes[mesh_id] if boundary.get("multimesh_mode", "closest") == "closest_y": mesh_paths = list(Path(boundary["path"]).iterdir()) @@ -265,26 +267,24 @@ def get_distance_to_mesh(mesh, ray_origin, ray_direction, mesh_type): mesh_id = np.argmin(distances) boundary["mesh_name"] = Path(mesh_paths[mesh_id]).stem - mesh = meshes[mesh_id] + self.mesh = meshes[mesh_id] if boundary.get("multimesh_mode", "closest") == "inside": - mesh = None for mesh_path in Path(boundary["path"]).iterdir(): _mesh = trimesh.load_mesh(mesh_path) if _mesh.contains([soma_position])[0]: - mesh = _mesh + self.mesh = _mesh break if boundary.get("multimesh_mode", "closest") == "territories": glom_id = int(cell.other_parameters["glomerulus_id"]) - mesh = None if glom_id >= 0: - mesh = trimesh.load_mesh( + self.mesh = trimesh.load_mesh( Path(boundary["path"]) / f"glomeruli_{glom_id:03d}.obj" ) else: - mesh = trimesh.load_mesh(boundary["path"]) + self.mesh = trimesh.load_mesh(boundary["path"]) if mode == "repulsive": @@ -325,15 +325,15 @@ def prob(direction, current_point, mesh=None, mesh_type=None, params=None): else: raise ValueError(f"boundary mode {mode} not understood!") - if mesh is not None: + if self.mesh is not None: if "params_section" in boundary: boundary["section_prob"] = partial( - prob, mesh=mesh, mesh_type=mesh_type, params=boundary["params_section"] + prob, mesh=self.mesh, mesh_type=mesh_type, params=boundary["params_section"] ) if "params_trunk" in boundary: boundary["trunk_prob"] = partial( - prob, mesh=mesh, mesh_type=mesh_type, params=boundary["params_trunk"] + prob, mesh=self.mesh, mesh_type=mesh_type, params=boundary["params_trunk"] ) boundaries.append(boundary) @@ -350,6 +350,9 @@ def layer_fraction_to_position(self, layer: int, layer_fraction: float) -> float Returns: target depth """ layer_thickness = self.layer_depths[layer] - self.layer_depths[layer - 1] + if layer_thickness <= 0: + raise SkipSynthesisError(f"Layer {layer} has no/negative thickness") + return self.layer_depths[layer - 1] + (1.0 - layer_fraction) * layer_thickness def lookup_target_reference_depths(self, depth: float) -> np.array: @@ -382,7 +385,10 @@ def distance_to_constraint(self, depth: float, constraint: Dict) -> Optional[flo constraint_position = self.layer_fraction_to_position( constraint["layer"], constraint["fraction"] ) - return depth - constraint_position + out = depth - constraint_position + if np.isnan(out): + raise SkipSynthesisError("Nan in some PH values") + return out @attr.s(auto_attribs=True) @@ -573,10 +579,13 @@ def completion(self, synthesized): synthesized.apical_points, morph_path ) - # enforce freeing memory of orientation nrrd data + # enforce freeing memory of orientation nrrd data and meshes if hasattr(self.context, "orientations"): if self.context.orientations is not None: del self.context.orientations + if hasattr(self.context, "mesh"): + if self.context.mesh is not None: + del self.context.mesh return { "name": morph_name,