Feat: Use 3d_angles and simpler diametrizer by default

.pre-commit-config.yaml

@@ -1,5 +1,5 @@
 default_language_version:
-  python: python3.8
+  python: python3.9
 repos:
   - repo: https://github.com/pre-commit/pre-commit-hooks
     rev: v4.4.0

examples/trunk_validation/luigi.cfg

@@ -1,66 +1,32 @@
-# TODO: Update this file according to new target classes
-
 # luigi parameters
-# [core]
-# logging_conf_file = logging.conf
+[core]
+logging_conf_file = logging.conf
 
-# global parameters
 [SynthesisConfig]
-tmd_parameters_path = out_vacuum/tmd_parameters.json
-tmd_distributions_path = out_vacuum/tmd_distributions.json
-mtypes = ["L5_TPC:A"]
-
-[PathConfig]
-morphs_df_path = morphs_df.csv
-synth_morphs_df_path = out_vacuum/synth_morphs_df.csv
-synth_output_path = out_vacuum/synthesized_morphologies
-substituted_morphs_df_path = out_vacuum/substituted_morphs_df.csv
-morphology_path = repaired_morphology_path
+# uncomment below to select specific mtypes
+mtypes = ["L5_TPC:A", "L6_BPC"]
 
 # prepare step
 [GetSynthesisInputs]
-url = git@bbpgitlab.epfl.ch:circuits/thncx.git
-git_synthesis_input_path = entities/synthesis_workflow/rat/vacuum/synthesis_input
-local_synthesis_input_path = synthesis_input
-
-[BuildMorphsDF]
-neurondb_path = /gpfs/bbp.cscs.ch/project/proj83/home/gevaert/morph-release/morph_release_old_code-2020-07-27/output/06_RepairUnravel-asc/neuronDB.xml
-morphology_dirs = {"repaired_morphology_path": "/gpfs/bbp.cscs.ch/project/proj83/home/gevaert/morph-release/morph_release_old_code-2020-07-27/output/06_RepairUnravel-asc", "repaired_morphology_path_h5": "/gpfs/bbp.cscs.ch/project/proj83/home/gevaert/morph-release/morph_release_old_code-2020-07-27/output/06_RepairUnravel-h5"}
-
-
-# synthesis setup
-[ApplySubstitutionRules]
-substitution_rules_path = substitution_rules.yaml
+url = git@bbpgitlab.epfl.ch:neuromath/synthdb.git
+git_synthesis_input_path = synthdb/insitu_synthesis_inputs/mouse_isocortex
 
-[BuildSynthesisParameters]
-tmd_parameters_path = out_vacuum/tmd_parameters_no_scaling.json
-#input_tmd_parameters_path = tmd_specific_parameters.json
-morphology_path = repaired_morphology_path_h5
+[GetDefaultParameters]
+trunk_method = 3d_angles
 
-[BuildSynthesisDistributions]
-morphology_path = repaired_morphology_path_h5
+[BuildMorphsDF]
+neurondb_path = /gpfs/bbp.cscs.ch/project/proj82/home/gevaert/morphology_release/mouse-scaled/scaled_output/06_RepairUnravel-asc/neuronDB.xml
+morphology_dirs = {"morphology_path": "/gpfs/bbp.cscs.ch/project/proj82/home/gevaert/morphology_release/mouse-scaled/scaled_output/06_RepairUnravel-asc", "morphology_path_h5": "gpfs/bbp.cscs.ch/project/proj82/home/gevaert/morphology_release/mouse-scaled/scaled_output/06_RepairUnravel-h5"}
 
-# synthesize in vacuum
 [VacuumSynthesize]
-vacuum_synth_morphology_path = vacuum_synth_morphologies
-vacuum_synth_morphs_df_path = out_vacuum/vacuum_synth_morphs_df.csv
-n_cells = 100
+n_cells = 200
 
 # validation plots
 [ValidateVacuumSynthesis]
+with_trunk_validation=True
 with_morphometrics = False
 with_density_profiles = False
 with_vacuum_morphologies= False
 with_score_matrix_reports=False
-with_trunk_validation=True
 with_morphology_validation_reports=False
 with_scale_statistics=False
-
-[PlotVacuumMorphologies]
-pdf_filename = figures_vacuum/vacuum_morphologies.pdf
-
-[PlotMorphometrics]
-morphometrics_path = figures_vacuum/morphometrics
-
-[PlotDensityProfiles]
-density_profiles_path = figures_vacuum/density_profiles.pdf

requirements/base.pip

@@ -4,7 +4,7 @@ bluepy-configfile>=0.1.19
 bluepymm>=0.8.5
 bluepyparallel>=0.0.8
 brainbuilder>=0.18.3
-diameter_synthesis>=0.5.2
+diameter_synthesis>=0.5.3
 dictdiffer>=0.9
 gitpython>=3.1.30
 jinja2>=3
@@ -18,12 +18,12 @@ morphio>=3.3.4,<4
 neuroc>=0.2.8,<1
 neurocollage>=0.2.1
 neurom>=3.2.2,<4
-neurots>=3.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.2
+region_grower>=0.4.3
 scipy>=1.10
 seaborn>=0.12.2
 tmd>=2.2

src/synthesis_workflow/synthesis.py

@@ -81,7 +81,6 @@ def _build_distributions_single_mtype(
     diameter_model_function=None,
     config=None,
     morphology_path=None,
-    trunk_method="simple",
 ):
     """Internal function for multiprocessing of tmd_distribution building."""
     data = {}
@@ -98,9 +97,6 @@ def _build_distributions_single_mtype(
             "diameter_input_morph": morphology_paths,
             "diameter_model": partial(diameter_model_function, config=config),
         }
-        if trunk_method != "simple":
-            kwargs["trunk_method"] = trunk_method
-
         _data = extract_input.distributions(morphology_paths, **kwargs)
         data[neurite_type] = _data[neurite_type]
         data["diameter"] = _data["diameter"]
@@ -118,7 +114,6 @@ def build_distributions(
     region_structure_path,
     nb_jobs=-1,
     joblib_verbose=10,
-    trunk_method="simple",
 ):
     """Build tmd_distribution dictionary for synthesis.
 
@@ -130,7 +125,6 @@ def build_distributions(
         region_structure_path (str): path to region_structure yaml file with thicknesses
         nb_jobs (int): number of jobs to run in parallal with joblib
         joblib_verbose (int): verbose level of joblib
-        trunk_method (str): method to set trunk on soma for synthesis (simple|3d_angle)
 
     Returns:
         dict: dict to save to tmd_distribution.json
@@ -142,7 +136,6 @@ def build_distributions(
         diameter_model_function=diameter_model_function,
         config=config,
         morphology_path=morphology_path,
-        trunk_method=trunk_method,
     )
     thicknesses = []
     if Path(region_structure_path).exists():

src/synthesis_workflow/tasks/synthesis.py

@@ -19,6 +19,8 @@
 from neurom import load_morphology
 from neurom.check.morphology_checks import has_apical_dendrite
 from neurots import extract_input
+from neurots.generate.orientations import fit_3d_angles
+from neurots.validator import validate_neuron_distribs
 from neurots.validator import validate_neuron_params
 from placement_algorithm.app.compact_annotations import _collect_annotations
 from region_grower.synthesize_morphologies import SynthesizeMorphologies
@@ -139,7 +141,6 @@ class GetDefaultParameters(WorkflowTask):
         default="neurots_input/tmd_parameters_default.json",
         description=":str: Path to default tmd_parameters.json.",
     )
-    trunk_method = luigi.ChoiceParameter(default="simple", choices=["simple", "3d_angle"])
 
     def requires(self):
         """Required input tasks."""
@@ -170,8 +171,6 @@ def run(self):
             config = DiametrizerConfig().config_diametrizer
             config["neurite_types"] = neurite_types
             kwargs = {"neurite_types": neurite_types, "diameter_parameters": config}
-            if self.trunk_method != "simple":
-                kwargs["trunk_method"] = self.trunk_method
             tmd_parameters[mtype] = extract_input.parameters(**kwargs)
 
         with self.output().open("w") as f:
@@ -188,7 +187,7 @@ class BuildSynthesisParameters(WorkflowTask):
 
     def requires(self):
         """Required input tasks."""
-        return {"tmd_parameters": OverwriteCustomParameters()}
+        return {"tmd_parameters": AddTrunkFitToParameters()}
 
     def run(self):
         """Actual process of the task."""
@@ -217,8 +216,6 @@ class BuildSynthesisDistributions(WorkflowTask):
         nb_jobs (int): Number of workers.
     """
 
-    trunk_method = luigi.ChoiceParameter(default="simple", choices=["simple", "3d_angle"])
-
     def requires(self):
         """Required input tasks."""
         return {"rules": ApplySubstitutionRules(), "synthesis": GetSynthesisInputs()}
@@ -245,9 +242,11 @@ def run(self):
             self.morphology_path,
             self.input()["synthesis"].pathlib_path / CircuitConfig().region_structure_path,
             nb_jobs=self.nb_jobs,
-            trunk_method=self.trunk_method,
         )
 
+        for distr in tmd_distributions["mtypes"].values():
+            validate_neuron_distribs(distr)
+
         with self.output().open("w") as f:
             json.dump(tmd_distributions, f, cls=NumpyEncoder, indent=4, sort_keys=True)
 
@@ -601,7 +600,6 @@ class AddScalingRulesToParameters(WorkflowTask):
     Attributes:
         morphology_path (str): Column name to use in the DF to compute
             axon_morphs_base_dir if it is not provided.
-        tmd_parameters_path (str): The path to the TMD parameters.
         nb_jobs (int): Number of threads used for synthesis.
     """
 
@@ -653,6 +651,56 @@ def output(self):
         return SynthesisLocalTarget(self.scaling_tmd_parameters_path)
 
 
+@copy_params(
+    morphology_path=ParamRef(PathConfig),
+    nb_jobs=ParamRef(RunnerConfig),
+)
+class AddTrunkFitToParameters(WorkflowTask):
+    """Add fits to trunk angles to tmd_parameter.json.
+
+    Attributes:
+        morphology_path (str): Column name to use in the DF to compute
+            axon_morphs_base_dir if it is not provided.
+        nb_jobs (int): Number of threads used for synthesis.
+    """
+
+    trunk_tmd_parameters_path = luigi.PathParameter(
+        default="neurots_input/tmd_parameters_trunk.json",
+        description=":str: Path to tmd_parameters.json with trunk fit added.",
+    )
+
+    scaling_rules_path = luigi.Parameter(
+        default="scaling_rules.yaml",
+        description=":str: Path to the file containing the scaling rules.",
+    )
+
+    def requires(self):
+        """Required input tasks."""
+        return {
+            "tmd_parameters": OverwriteCustomParameters(),
+            "tmd_distributions": BuildSynthesisDistributions(),
+        }
+
+    def run(self):
+        """Actual process of the task."""
+        with self.input()["tmd_parameters"].open("r") as f:
+            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]
+            )
+
+        with self.output().open("w") as f:
+            json.dump(tmd_parameters, f, cls=NumpyEncoder, indent=4, sort_keys=True)
+
+    def output(self):
+        """Outputs of the task."""
+        return SynthesisLocalTarget(self.trunk_tmd_parameters_path)
+
+
 class OverwriteCustomParameters(WorkflowTask):
     """Overwrite parameters with custom parameters.
 

src/synthesis_workflow/tasks/validation.py

@@ -32,10 +32,10 @@
 from synthesis_workflow.tasks.config import RunnerConfig
 from synthesis_workflow.tasks.config import SynthesisConfig
 from synthesis_workflow.tasks.config import ValidationLocalTarget
-from synthesis_workflow.tasks.synthesis import AddScalingRulesToParameters
 from synthesis_workflow.tasks.synthesis import ApplySubstitutionRules
 from synthesis_workflow.tasks.synthesis import BuildMorphsDF
 from synthesis_workflow.tasks.synthesis import BuildSynthesisDistributions
+from synthesis_workflow.tasks.synthesis import BuildSynthesisParameters
 from synthesis_workflow.tasks.synthesis import Synthesize
 from synthesis_workflow.tasks.vacuum_synthesis import VacuumSynthesize
 from synthesis_workflow.vacuum_synthesis import VACUUM_SYNTH_MORPHOLOGY_PATH
@@ -467,7 +467,7 @@ class PlotPathDistanceFits(WorkflowTask):
     def requires(self):
         """Required input tasks."""
         return {
-            "scaling_rules": AddScalingRulesToParameters(),
+            "parameters": BuildSynthesisParameters(),
             "rescaled": ApplySubstitutionRules(),
             "distributions": BuildSynthesisDistributions(),
         }
@@ -476,7 +476,7 @@ def run(self):
         """Actual process of the task."""
         L.debug("output_path = %s", self.output().path)
         plot_path_distance_fits(
-            self.input()["scaling_rules"].path,
+            self.input()["parameters"].path,
             self.input()["distributions"].path,
             self.input()["rescaled"].path,
             self.morphology_path,
@@ -669,7 +669,7 @@ class TrunkValidation(WorkflowTask):
         description=":str: Path to output directory (relative from ``PathConfig.result_path``).",
     )
     base_key = luigi.Parameter(
-        default="repaired_morphology_path",
+        default="morphology_path",
         description=":str: Base key to use in the morphology DataFrame.",
     )
     comp_key = luigi.Parameter(
@@ -683,10 +683,15 @@ class TrunkValidation(WorkflowTask):
 
     def requires(self):
         """Required input tasks."""
+        tasks = {
+            "distributions": BuildSynthesisDistributions(),
+            "parameters": BuildSynthesisParameters(),
+        }
         if self.in_atlas:
-            return {"morphs": BuildMorphsDF(), "circuit": ConvertCircuit()}
+            tasks.update({"morphs": BuildMorphsDF(), "circuit": ConvertCircuit()})
         else:
-            return {"vacuum": VacuumSynthesize(), "morphs": ApplySubstitutionRules()}
+            tasks.update({"vacuum": VacuumSynthesize(), "morphs": ApplySubstitutionRules()})
+        return tasks
 
     def run(self):
         """Actual process of the task."""
@@ -702,10 +707,12 @@ def run(self):
             morphs_df,
             synth_morphs_df,
             self.output().pathlib_path,
-            base_key=self.base_key,
-            comp_key=comp_key,
-            base_label=self.base_label,
-            comp_label=self.comp_label,
+            self.base_key,
+            comp_key,
+            self.base_label,
+            self.comp_label,
+            self.input()["parameters"].path,
+            self.input()["distributions"].path,
         )
 
     def output(self):

src/synthesis_workflow/utils.py

@@ -105,10 +105,23 @@ def create_parameter_diff(param, param_spec):
 def _create_entry(param, entry):
     """Create a dict entry if it does not exist."""
     if entry[0].endswith("]"):
-        param[entry[0].split("[")[0]].append([None])
+        # add list for specific orientations
+        if entry[0].split("[")[0] not in param:
+            param[entry[0].split("[")[0]] = [None]
+        else:
+            param[entry[0].split("[")[0]].append([None])
+    elif entry[0] in param and not isinstance(param[entry[0]], dict):
+        # erase anything that is not a dict to be able to go deeper
+        param[entry[0]] = {}
     elif entry[0] not in param:
-        param[entry[0]] = None
+        # add empty entries depending on depths
+        if len(entry) > 1:
+            param[entry[0]] = {}
+        else:
+            param[entry[0]] = None
+
     if len(entry) > 1:
+        # go deeper if needed
         _create_entry(param[entry[0]], entry[1:])
 
 

src/synthesis_workflow/vacuum_synthesis.py

@@ -38,6 +38,7 @@ def _grow_morphology(
         input_parameters=tmd_parameters,
         input_distributions=tmd_distributions,
         external_diametrizer=external_diametrizer,
+        skip_preprocessing=True,
     )
     grower.grow()
     grower.neuron.write(morphology_path)