Boundary examples

.pylintrc

@@ -15,7 +15,7 @@ max-line-length=100
 
 [DESIGN]
 # Maximum number of arguments for function / method
-max-args=10
+max-args=12
 # Argument names that match this expression will be ignored. Default to name
 # with leading underscore
 ignored-argument-names=_.*

README.rst

@@ -4,8 +4,8 @@ Synthesis Workflow
 This project contains several workflows used for neuron synthesis and the validation of this process.
 It is divided into two packages:
 
-* **synthesis-workflow**, which contains the workflow tasks and tools.
-* **MorphVal**, which is a library used for morphology validation and can be used as a standalone.
+* **synthesis-workflow**, the main package, which contains the workflow tasks and tools.
+* **MorphVal**, which is a legacy library used for morphology validation.
 
 The complete documentation of this project is available here: `<https://bbpteam.epfl.ch/documentation/projects/synthesis-workflow/latest/index.html>`_
 

examples/O1_example/luigi.cfg

@@ -1,24 +1,61 @@
 # run syntheis in an O1 atlas
 
 [core]
-workers = 1
+workers = 10
 logging_conf_file = logging.conf
 
-[SynthesisConfig]
-mtypes = ["L5_TPC:A"]
-
 [RunnerConfig]
 nb_jobs = 5
 
+[SynthesisConfig]
+axon_method = synthesis
+
+# L1 IN
+#mtypes = ["L1_DAC", "L1_HAC", "L1_NGC-SA", "L1_NGC-DA", "L1_LAC", "L1_SAC"]
+
+# L23 IN
+# mtypes = ["L23_MC", "L23_ChC", "L23_NBC", "L23_DBC", "L23_LBC", "L23_BTC", "L23_NGC", "L23_BP"]
+
+# L23 PC
+# mtypes = ["L2_TPC:A", "L2_TPC:B", "L2_IPC", "L3_TPC:A", "L3_TPC:C"]
+
+# L4 IN
+# mtypes = ["L4_SSC", "L4_ChC", "L4_BP", "L4_BTC", "L4_NBC", "L4_NGC","L4_SBC", "L4_DBC", "L4_LBC", "L4_MC"]
+
+# L4 PC
+# mtypes = ["L4_TPC", "L4_UPC"]
+
+# L5 IN
+# mtypes = ["L5_ChC", "L5_BP", "L5_BTC", "L5_NBC", "L5_NGC","L5_SBC", "L5_DBC", "L5_LBC", "L5_MC"]
+
+# L5 PC
+# mtypes = ["L5_TPC:A", "L5_TPC:B", "L5_TPC:C", "L5_UPC"]
+
+# L6 PC
+#mtypes = ["L6_TPC:A", "L6_TPC:C", "L6_HPC", "L6_IPC", "L6_BPC", "L6_UPC"]
+
+# L6 IN
+#mtypes = ["L6_ChC", "L6_BP", "L6_BTC", "L6_NBC", "L6_NGC", "L6_SBC", "L6_DBC", "L6_LBC", "L6_MC"]
+
+# all
+mtypes = ["L1_DAC", "L1_HAC", "L1_NGC-SA", "L1_NGC-DA", "L1_LAC", "L1_SAC", "L23_MC", "L23_ChC", "L23_NBC", "L23_DBC", "L23_LBC", "L23_BTC", "L23_NGC", "L23_BP", "L2_TPC:A", "L2_TPC:B", "L2_IPC", "L3_TPC:A", "L3_TPC:C", "L4_SSC", "L4_ChC", "L4_BP", "L4_BTC", "L4_NBC", "L4_NGC","L4_SBC", "L4_DBC", "L4_LBC", "L4_MC", "L4_TPC", "L4_UPC", "L5_ChC", "L5_BP", "L5_BTC", "L5_NBC", "L5_NGC","L5_SBC", "L5_DBC", "L5_LBC", "L5_MC", "L5_TPC:A", "L5_TPC:B", "L5_TPC:C", "L5_UPC", "L6_TPC:A", "L6_TPC:C", "L6_HPC", "L6_IPC", "L6_BPC", "L6_UPC", "L6_ChC", "L6_BP", "L6_BTC", "L6_NBC", "L6_NGC", "L6_SBC", "L6_DBC", "L6_LBC", "L6_MC"]
+
+#mtypes = ["L4_MC"]
 [CircuitConfig]
 atlas_path = atlas
 region = O0
 
 [BuildCircuit]
-density_factor = 1
+density_factor = 1.0
+
+[CreateBoundaryMask]
+boundary_thickness = 2
+
+[CreateBoundaryMask]
+boundary_thickness = 2
 
 [SliceCircuit]
-n_cells = 5
+n_cells = 10
 
 [CreateAtlasPlanes]
 plane_type = centerline_straight
@@ -27,18 +64,19 @@ slice_thickness = 50
 
 [GetSynthesisInputs]
 url = [email protected]:neuromath/synthdb.git
-git_synthesis_input_path = synthdb/insitu_synthesis_inputs/rat_O1
-branch = main
+git_synthesis_input_path = synthdb/insitu_synthesis_inputs/rat_O0
+branch = o0o1
 
 [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 = {"morphology_path": "/gpfs/bbp.cscs.ch/project/proj83/home/gevaert/morph-release/morph_release_old_code-2020-07-27/output/06_RepairUnravel-asc"}
+morphology_dirs = {"path": "/gpfs/bbp.cscs.ch/project/proj83/home/gevaert/morph-release/morph_release_old_code-2020-07-27/output/06_RepairUnravel-asc"}
 
 [ValidateSynthesis]
 with_collage = True
-with_path_distance_fits = False
-with_morphometrics = False
-with_density_profiles = False
-with_scale_statistics = False
-with_score_matrix_reports=False
-with_morphology_validation_reports = False
+with_trunk_validation = True
+with_path_distance_fits = True
+with_morphometrics = True
+with_density_profiles = True
+with_scale_statistics = True
+with_score_matrix_reports = True
+with_morphology_validation_reports = True

examples/O1_example/run.sh

@@ -1,9 +1,9 @@
 #!/bin/bash -l
 
-export OMP_NUM_THREADS=1
+#export OMP_NUM_THREADS=1
 
-rm -rf atlas
-brainbuilder atlases -n 6,5,4,3,2,1 -t 700,525,190,353,149,165 -d 10 -o atlas column -a 1000
+#rm -rf atlas
+#brainbuilder atlases -n 6,5,4,3,2,1 -t 700,525,190,353,149,165 -d 10 -o atlas column -a 1000
 
 python -m luigi --module synthesis_workflow.tasks.workflows ValidateSynthesis \
     --local-scheduler \

examples/boundary_example/logging.conf

@@ -0,0 +1,38 @@
+[loggers]
+keys=root,luigi,luigi_interface
+
+[handlers]
+keys=consoleHandler,fileHandler
+
+[formatters]
+keys=PrettyFormatter
+
+[logger_root]
+level=INFO
+handlers=consoleHandler,fileHandler
+
+[logger_luigi]
+level=INFO
+handlers=consoleHandler,fileHandler
+qualname=luigi
+propagate=0
+
+[logger_luigi_interface]
+level=INFO
+handlers=consoleHandler,fileHandler
+qualname=luigi-interface
+propagate=0
+
+[handler_consoleHandler]
+class=StreamHandler
+formatter=PrettyFormatter
+args=(sys.stdout,)
+
+[handler_fileHandler]
+class=FileHandler
+formatter=PrettyFormatter
+args=('synthesis_workflow.log',)
+
+[formatter_PrettyFormatter]
+format=%(asctime)s - %(name)s - %(levelname)s - %(message)s
+datefmt=%Y-%m-%d %H:%M:%S

examples/boundary_example/luigi.cfg

@@ -0,0 +1,51 @@
+# run syntheis in an O1 atlas
+
+[core]
+workers = 10
+logging_conf_file = logging.conf
+
+[RunnerConfig]
+nb_jobs = 5
+
+[SynthesisConfig]
+axon_method = synthesis
+
+#mtypes = ["L1_open", "L1_hard", "L1_leaky", "L1_hard_trunk", "L1_attractive"]
+mtypes = ["L1_open", "L1_left", "L1_down"]
+
+[CircuitConfig]
+atlas_path = atlas
+region = O0
+
+[BuildCircuit]
+density_factor = 1.0
+
+[CreateBoundaryMask]
+boundary_thickness = 1
+
+[SliceCircuit]
+n_cells = 10
+
+[CreateAtlasPlanes]
+plane_type = centerline_straight
+plane_count = 1
+slice_thickness = 20
+
+[GetSynthesisInputs]
+url = [email protected]:neuromath/synthdb.git
+git_synthesis_input_path = synthdb/insitu_synthesis_inputs/boundary_example
+branch = boundary_example
+
+[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 = {"path": "/gpfs/bbp.cscs.ch/project/proj83/home/gevaert/morph-release/morph_release_old_code-2020-07-27/output/06_RepairUnravel-asc"}
+
+[ValidateSynthesis]
+with_collage = True
+with_trunk_validation = True
+with_path_distance_fits = False
+with_morphometrics = False
+with_density_profiles = False
+with_scale_statistics = False
+with_score_matrix_reports=False
+with_morphology_validation_reports = False

examples/boundary_example/run.sh

@@ -0,0 +1,12 @@
+#!/bin/bash -l
+
+export OMP_NUM_THREADS=1
+export REGION_GROWER_BOUNDARY_DEBUG=1
+
+rm -rf atlas
+brainbuilder atlases -n 2,1 -t 100,100 -d 10 -o atlas column -a 200
+
+
+python -m luigi --module synthesis_workflow.tasks.workflows ValidateSynthesis \
+    --local-scheduler \
+    --log-level INFO \

src/synthesis_workflow/__init__.py

@@ -1,6 +1,5 @@
 """Workflow for neuronal synthesis validation."""
-
-import pkg_resources
+import importlib.metadata
 from morphio import SectionType  # noqa
 
-__version__ = pkg_resources.get_distribution("synthesis_workflow").version
+__version__ = importlib.metadata.version("synthesis-workflow")

src/synthesis_workflow/circuit.py

@@ -197,3 +197,43 @@ def get_layer_tags(atlas_dir, region_structure_path, region=None):
             L.warning("No voxel found for layer %s.", layer)
     brain_regions.raw = layers_data
     return brain_regions, layer_mapping
+
+
+def add_mclass(cells_path, inh_probmap):
+    """Add mclass column to circuit for barrel cortex."""
+    cells = CellCollection.load(cells_path)
+    cells_df = cells.as_dataframe()
+    cells_df["metype"] = [
+        f"{mtype}_{etype}" for mtype, etype in cells_df[["mtype", "etype"]].values
+    ]
+
+    exc_df = cells_df[cells_df["synapse_class"] == "EXC"]
+
+    def generate_exc_marker_labels(mtype):
+        layer = mtype.split("_")[0]
+        if layer == "L2" or layer == "L3":
+            layer = "L23"
+        return f"{layer}_EXC"
+
+    exc_df["marker"] = exc_df.apply(lambda x: generate_exc_marker_labels(x["mtype"]), axis=1)
+
+    group_by_marker = exc_df.groupby("metype")
+    metype_to_mclass = group_by_marker["marker"].value_counts() / group_by_marker["metype"].count()
+    metype_to_mclass.name = "probability"
+
+    exc_probmap = (
+        metype_to_mclass.reset_index()
+        .pivot_table(index="marker", columns="metype", values="probability")
+        .T
+    )
+    prob_map = pd.concat([inh_probmap, exc_probmap]).fillna(0)
+
+    for metype, mapping in prob_map.iterrows():
+        cells_metype = cells_df[cells_df["metype"] == metype]
+        subtype_mapping = mapping[mapping > 0]
+        cells_df.loc[cells_metype.index, "mclass"] = np.random.choice(
+            subtype_mapping.index.values, size=len(cells_metype), p=subtype_mapping.values
+        )
+
+    cells.add_properties({"mclass": cells_df.mclass.values})
+    return cells

src/synthesis_workflow/synthesis.py

@@ -44,12 +44,11 @@ def get_neurite_types(morphs_df):
         morph_class = list(set(_df["morph_class"]))
 
         if len(morph_class) > 1:
-            raise ValueError("f{mtype} has not consistent morph_class, we stop here")
+            L.warning(f"{mtype} has not consistent morph_class, we take PC if apical, see {_df}")
 
-        if morph_class[0] == "IN":
-            neurite_types[mtype] = ["basal_dendrite"]
-        if morph_class[0] == "PC":
-            neurite_types[mtype] = ["basal_dendrite", "apical_dendrite"]
+        neurite_types[mtype] = ["basal_dendrite"]
+        if "PC" in morph_class:
+            neurite_types[mtype] += ["apical_dendrite"]
     return neurite_types
 
 
@@ -86,24 +85,32 @@ def _build_distributions_single_mtype(
     """Internal function for multiprocessing of tmd_distribution building."""
     data = {}
     for neurite_type in neurite_types[mtype]:
-        if "use_axon" in morphs_df.columns:
-            morphology_paths = morphs_df.loc[
-                (morphs_df.mtype == mtype) & morphs_df.use_axon, morphology_path
-            ].to_list()
-        else:
-            morphology_paths = morphs_df.loc[morphs_df.mtype == mtype, morphology_path].to_list()
+        _morphs_df = morphs_df[morphs_df.mtype == mtype]
+        if neurite_type == "axon" and "use_axon" in morphs_df.columns:
+            _morphs_df = _morphs_df.loc[morphs_df.use_axon]
+        if neurite_type == "apical_dendrite" and "use_apical" in morphs_df.columns:
+            _morphs_df = _morphs_df.loc[morphs_df.use_apical]
+        if neurite_type == "basal_dendrite" and "use_basal" in morphs_df.columns:
+            _morphs_df = _morphs_df.loc[morphs_df.use_basal]
+
+        morphology_paths = _morphs_df[morphology_path].to_list()
+
         config["neurite_types"] = neurite_types[mtype]
         kwargs = {
             "neurite_types": neurite_types[mtype],
             "diameter_input_morph": morphology_paths,
+            "threshold_sec": 0,
+            "diameter_model": partial(diameter_model_function, config=config),
         }
-        if config["models"][0] != "simpler":
-            config["diameter_model"] = partial(diameter_model_function, config=config)
-        _data = extract_input.distributions(morphology_paths, **kwargs)
+        try:
+            _data = extract_input.distributions(morphology_paths, **kwargs)
+
+            data[neurite_type] = _data[neurite_type]
+            data["diameter"] = _data["diameter"]
+            data["soma"] = _data["soma"]
+        except Exception as exc:
+            print(mtype, neurite_type, exc)
 
-        data[neurite_type] = _data[neurite_type]
-        data["diameter"] = _data["diameter"]
-        data["soma"] = _data["soma"]
     return mtype, data