Feat: Add region support

requirements/base.pip

@@ -16,14 +16,14 @@ matplotlib>=3.6.3
 morph_tool>=2.9.1,<3
 morphio>=3.3.4,<4
 neuroc>=0.2.8,<1
-neurocollage>=0.2.1
+neurocollage>=0.3.0
 neurom>=3.2.2,<4
 neurots>=3.3.1,<4
 numpy>=1.24.1
 pandas>=1.5.3
 placement_algorithm>=2.3.1
 PyYAML>=6
-region_grower>=0.4.3,<1
+region_grower>=1.0.0,<2
 scipy>=1.10
 seaborn>=0.12.2
 tmd>=2.2

src/synthesis_workflow/circuit.py

@@ -86,6 +86,7 @@ def build_circuit(
     """
     if seed is not None:
         np.random.seed(seed)
+
     return place(
         composition_path=cell_composition_path,
         mtype_taxonomy_path=mtype_taxonomy_path,

src/synthesis_workflow/extra/insitu_validation.py

@@ -257,7 +257,7 @@ def plot_layer_collage(
                 _points = points + pos[np.newaxis]
                 layers = layer_annotation.lookup(_points, outer_value=0)
 
-                depth = atlas.depths.lookup(pos)
+                depth = atlas.depths[region].lookup(pos)
 
                 rotation = atlas.orientations.lookup(pos)[0]
 

src/synthesis_workflow/synthesis.py

@@ -8,7 +8,6 @@
 import matplotlib
 import numpy as np
 import pandas as pd
-import yaml
 from joblib import Parallel
 from joblib import delayed
 from morphio.mut import Morphology
@@ -111,7 +110,7 @@ def build_distributions(
     diameter_model_function,
     config,
     morphology_path,
-    region_structure_path,
+    region,
     nb_jobs=-1,
     joblib_verbose=10,
 ):
@@ -122,7 +121,7 @@ def build_distributions(
         morphs_df (DataFrame): morphology dataframe with reconstruction to use
         diameter_model_function (function): diametrizer function (from diameter-synthesis)
         morphology_path (str): name of the column in morphs_df to use for paths to morphologies
-        region_structure_path (str): path to region_structure yaml file with thicknesses
+        region (str): region we are building
         nb_jobs (int): number of jobs to run in parallal with joblib
         joblib_verbose (int): verbose level of joblib
 
@@ -137,19 +136,11 @@ def build_distributions(
         config=config,
         morphology_path=morphology_path,
     )
-    thicknesses = []
-    if Path(region_structure_path).exists():
-        with open(region_structure_path, "r", encoding="utf-8") as r_p:
-            region_structure = yaml.safe_load(r_p)
-        thicknesses = [
-            region_structure["thicknesses"][layer] for layer in region_structure["layers"]
-        ]
-
-    tmd_distributions = {"mtypes": {}, "metadata": {"cortical_thickness": thicknesses}}
+    tmd_distributions = {region: {}}
     for mtype, distribution in Parallel(nb_jobs, verbose=joblib_verbose)(
         delayed(build_distributions_single_mtype)(mtype) for mtype in mtypes
     ):
-        tmd_distributions["mtypes"][mtype] = distribution
+        tmd_distributions[region][mtype] = distribution
     return tmd_distributions
 
 

src/synthesis_workflow/tasks/circuit.py

@@ -217,6 +217,9 @@ def run(self):
             seed=self.seed,
             region=CircuitConfig().region,
         )
+        # new version of brainbuilder only assigns subregion, not sure why, to investigate
+        if "region" not in cells.properties:
+            cells.properties["region"] = cells.properties["subregion"]
         cells.save(self.output().path)
 
     def output(self):

src/synthesis_workflow/tasks/config.py

@@ -197,7 +197,7 @@ class CircuitConfig(luigi.Config):
         default="region_structure.yaml",
         description=":str: Path to the file containing the layer structure data.",
     )
-    region = luigi.OptionalParameter(default=None)
+    region = luigi.OptionalParameter(default="default")
     hemisphere = OptionalChoiceParameter(
         default=None,
         choices=["left", "right"],

src/synthesis_workflow/tasks/synthesis.py

@@ -155,13 +155,13 @@ def run(self):
         if SynthesisConfig().axon_method != "no_axon":
             for neurite_type in neurite_types.values():
                 neurite_type.append("axon")
-
-        tmd_parameters = {}
+        region = CircuitConfig().region
+        tmd_parameters = {region: {}}
         for mtype in tqdm(mtypes):
             config = DiametrizerConfig().config_diametrizer
             config["neurite_types"] = neurite_types[mtype]
             kwargs = {"neurite_types": neurite_types[mtype], "diameter_parameters": config}
-            tmd_parameters[mtype] = extract_input.parameters(**kwargs)
+            tmd_parameters[region][mtype] = extract_input.parameters(**kwargs)
 
         with self.output().open("w") as f:
             json.dump(tmd_parameters, f, cls=NumpyEncoder, indent=4, sort_keys=True)
@@ -183,8 +183,8 @@ def run(self):
         """Actual process of the task."""
         # possibly other fine tuning here or validate intermediate steps
         tmd_parameters = json.load(self.input()["tmd_parameters"].open("r"))
-        for mtype in tmd_parameters:
-            validate_neuron_params(tmd_parameters[mtype])
+        for mtype in tmd_parameters[CircuitConfig().region]:
+            validate_neuron_params(tmd_parameters[CircuitConfig().region][mtype])
 
         with self.output().open("w") as f:
             json.dump(tmd_parameters, f, cls=NumpyEncoder, indent=4, sort_keys=True)
@@ -230,11 +230,11 @@ def run(self):
             build_diameter_models,
             DiametrizerConfig().config_model,
             self.morphology_path,
-            self.input()["synthesis"].pathlib_path / CircuitConfig().region_structure_path,
+            region=CircuitConfig().region,
             nb_jobs=self.nb_jobs,
         )
 
-        for distr in tmd_distributions["mtypes"].values():
+        for distr in tmd_distributions[CircuitConfig().region].values():
             validate_neuron_distribs(distr)
 
         with self.output().open("w") as f:
@@ -621,7 +621,7 @@ def run(self):
                     scaling_rules = yaml.full_load(f)
 
         add_scaling_rules_to_parameters(
-            tmd_parameters,
+            tmd_parameters[CircuitConfig().region],
             self.input()["morphologies"].path,
             self.morphology_path,
             scaling_rules,
@@ -672,10 +672,10 @@ def run(self):
             tmd_parameters = json.load(f)
         with self.input()["tmd_distributions"].open("r") as f:
             tmd_distributions = json.load(f)
-
-        for mtype in tmd_parameters:
-            tmd_parameters[mtype] = fit_3d_angles(
-                tmd_parameters[mtype], tmd_distributions["mtypes"][mtype]
+        region = CircuitConfig().region
+        for mtype in tmd_parameters[region]:
+            tmd_parameters[region][mtype] = fit_3d_angles(
+                tmd_parameters[region][mtype], tmd_distributions[region][mtype]
             )
 
         with self.output().open("w") as f:
@@ -715,12 +715,12 @@ def run(self):
         custom_path = self.input()["synthesis_input"].pathlib_path / self.custom_parameters_path
         if custom_path.exists():
             custom_parameters = pd.read_csv(custom_path)
-            apply_parameter_diff(tmd_parameters, custom_parameters)
+            apply_parameter_diff(tmd_parameters[CircuitConfig().region], custom_parameters)
 
         # if we are no_axon, ensure tmd_parameters has no axon data, or json schema may crash
         if SynthesisConfig().axon_method == "no_axon":
-            for mtype in tmd_parameters:
-                tmd_parameters[mtype]["axon"] = {}
+            for mtype in tmd_parameters[CircuitConfig().region]:
+                tmd_parameters[CircuitConfig().region][mtype]["axon"] = {}
 
         with self.output().open("w") as f:
             json.dump(tmd_parameters, f, cls=NumpyEncoder, indent=4, sort_keys=True)

src/synthesis_workflow/tasks/vacuum_synthesis.py

@@ -10,6 +10,7 @@
 from luigi_tools.task import WorkflowTask
 from luigi_tools.task import copy_params
 
+from synthesis_workflow.tasks.config import CircuitConfig
 from synthesis_workflow.tasks.config import MorphsDfLocalTarget
 from synthesis_workflow.tasks.config import RunnerConfig
 from synthesis_workflow.tasks.config import SynthesisConfig
@@ -68,7 +69,7 @@ def run(self):
         tmd_distributions = json.load(self.input()["tmd_distributions"].open())
 
         if self.mtypes is None:
-            mtypes = list(tmd_parameters.keys())
+            mtypes = list(tmd_parameters[CircuitConfig().region].keys())
         else:
             mtypes = self.mtypes
 
@@ -80,6 +81,7 @@ def run(self):
             tmd_parameters,
             tmd_distributions,
             morphology_base_path.absolute(),
+            CircuitConfig().region,
             vacuum_morphology_path=self.vacuum_synth_morphology_path,
             diametrizer=self.diametrizer,
             joblib_verbose=self.joblib_verbose,

src/synthesis_workflow/tasks/validation.py

@@ -441,6 +441,7 @@ def output(self):
 
 @copy_params(
     mtypes=ParamRef(SynthesisConfig),
+    region=ParamRef(CircuitConfig),
     morphology_path=ParamRef(PathConfig),
     nb_jobs=ParamRef(RunnerConfig),
 )
@@ -482,6 +483,7 @@ def run(self):
             self.morphology_path,
             self.output().path,
             self.mtypes,
+            self.region,
             self.outlier_percentage,
             self.nb_jobs,
         )
@@ -713,6 +715,7 @@ def run(self):
             self.comp_label,
             self.input()["parameters"].path,
             self.input()["distributions"].path,
+            CircuitConfig().region,
         )
 
     def output(self):

src/synthesis_workflow/vacuum_synthesis.py

@@ -56,6 +56,7 @@ def grow_vacuum_morphologies(
     tmd_parameters,
     tmd_distributions,
     morphology_base_path,
+    region,
     vacuum_morphology_path=VACUUM_SYNTH_MORPHOLOGY_PATH,
     diametrizer="external",
     joblib_verbose=0,
@@ -69,8 +70,8 @@ def grow_vacuum_morphologies(
     global_gid = 0
     rows = []
     for mtype in tqdm(mtypes):
-        tmd_parameters[mtype]["diameter_params"]["method"] = diametrizer
-        tmd_distributions["mtypes"][mtype]["diameter"]["method"] = diametrizer
+        tmd_parameters[region][mtype]["diameter_params"]["method"] = diametrizer
+        tmd_distributions[region][mtype]["diameter"]["method"] = diametrizer
 
         if diametrizer == "external":
             external_diametrizer = build_diameters.build
@@ -88,8 +89,8 @@ def grow_vacuum_morphologies(
             delayed(_grow_morphology)(
                 gid,
                 mtype,
-                tmd_parameters[mtype],
-                tmd_distributions["mtypes"][mtype],
+                tmd_parameters[region][mtype],
+                tmd_distributions[region][mtype],
                 morphology_base_path,
                 vacuum_morphology_path=vacuum_morphology_path,
                 external_diametrizer=external_diametrizer,

src/synthesis_workflow/validation.py

@@ -642,6 +642,7 @@ def plot_path_distance_fits(
     morphology_path,
     output_path,
     mtypes=None,
+    region=None,
     outlier_percentage=90,
     nb_jobs=-1,
 ):
@@ -662,7 +663,8 @@ def plot_path_distance_fits(
             [
                 mtype
                 for mtype in morphs_df.mtype.unique().tolist()
-                if tmd_parameters.get(mtype, {})
+                if tmd_parameters[region]
+                .get(mtype, {})
                 .get("context_constraints", {})
                 .get("apical_dendrite", {})
                 .get("extent_to_target")
@@ -677,16 +679,15 @@ def plot_path_distance_fits(
     ]
 
     L.debug("Number of files: %s", [(t, len(f)) for t, f in file_lists])
-
     ensure_dir(output_path)
     with PdfPages(output_path) as pdf:
         for mtype, x, y, x_clean, y_clean, x_synth, y_synth, msg in Parallel(nb_jobs)(
             delayed(_get_fit_population)(
                 mtype,
                 files,
                 outlier_percentage,
-                tmd_parameters[mtype],
-                tmd_distributions["mtypes"][mtype],
+                tmd_parameters[region][mtype],
+                tmd_distributions[region][mtype],
             )
             for mtype, files in file_lists
         ):
@@ -705,10 +706,10 @@ def plot_path_distance_fits(
                 plt.plot(
                     [
                         get_path_distance_from_extent(
-                            tmd_parameters[mtype]["context_constraints"]["apical_dendrite"][
+                            tmd_parameters[region][mtype]["context_constraints"]["apical_dendrite"][
                                 "extent_to_target"
                             ]["slope"],
-                            tmd_parameters[mtype]["context_constraints"]["apical_dendrite"][
+                            tmd_parameters[region][mtype]["context_constraints"]["apical_dendrite"][
                                 "extent_to_target"
                             ]["intercept"],
                             i,
@@ -1100,6 +1101,7 @@ def trunk_validation(
     comp_label,
     tmd_parameters_path,
     tmd_distributions_path,
+    region,
 ):
     """Create plots to validate trunk angles."""
     with open(tmd_parameters_path, "r", encoding="utf-8") as f:
@@ -1120,24 +1122,27 @@ def trunk_validation(
                 fit = None
                 if (
                     "params"
-                    in tmd_parameters[mtype].get(t2, {}).get("orientation", {}).get("values", {})
-                    and tmd_parameters[mtype][t2]["orientation"]["mode"].split("_")[0]
+                    in tmd_parameters[region][mtype]
+                    .get(t2, {})
+                    .get("orientation", {})
+                    .get("values", {})
+                    and tmd_parameters[region][mtype][t2]["orientation"]["mode"].split("_")[0]
                     == t1.split("_")[0]
                 ):
                     fit = get_probability_function(
-                        form=tmd_parameters[mtype][t2]["orientation"]["values"]["form"],
+                        form=tmd_parameters[region][mtype][t2]["orientation"]["values"]["form"],
                         with_density=True,
                     )(
                         angles,
-                        *tmd_parameters[mtype][t2]["orientation"]["values"]["params"],
+                        *tmd_parameters[region][mtype][t2]["orientation"]["values"]["params"],
                     )
                 plt.figure(figsize=(6, 4))
                 plt.axvline(0, c="k")
                 plt.axvline(np.pi, c="k")
                 plt.plot(*_get_hist(data_bio[data_type]), label=base_label)
                 key = t1.split("_")[0] + "_3d_angles"
-                if key in tmd_distributions["mtypes"][mtype].get(t2, {}).get("trunk", {}):
-                    bio_data = tmd_distributions["mtypes"][mtype][t2]["trunk"][key]["data"]
+                if key in tmd_distributions[region][mtype].get(t2, {}).get("trunk", {}):
+                    bio_data = tmd_distributions[region][mtype][t2]["trunk"][key]["data"]
                     plt.plot(
                         bio_data["bins"],
                         np.array(bio_data["weights"]) / max(bio_data["weights"]),

src/version.py

@@ -1,2 +1,2 @@
 """Package version."""
-VERSION = "0.1.3"
+VERSION = "1.0.0.dev0"

tests/conftest.py

@@ -66,7 +66,8 @@ def small_O1(tmp_path):
         {
             "atlas": str(atlas_dir),
             "structure": str(DATA / "synthesis_input" / "region_structure.yaml"),
-        }
+        },
+        "O0",
     )["annotation"]
 
     layer_tags = br.raw.copy()

tests/data/in_small_O1/luigi.cfg

@@ -6,6 +6,7 @@ nb_jobs = 1
 
 [CircuitConfig]
 circuit_somata_path = circuit_somata.mvd3
+region = O0
 atlas_path = <path_to_atlas_directory>
 
 [PathConfig]