Reorganised the code ; Merged PlotMorphometrics and PlotVacuumMorphom…

…etrics tasks

Change-Id: Ibb174bcd4f5fd0cdaa3bcd668421f1b9c26ddd9c

synthesis_workflow/__init__.py

@@ -1 +1,9 @@
 """Workflow for neuronal synthesis validation."""
+from morphio import SectionType
+
+
+STR_TO_TYPES = {
+    "basal": SectionType.basal_dendrite,
+    "apical": SectionType.apical_dendrite,
+    "axon": SectionType.axon,
+}

synthesis_workflow/synthesis.py

@@ -9,40 +9,29 @@
 from pathlib import Path
 
 import matplotlib
-import matplotlib.pyplot as plt
 import numpy as np
 import pandas as pd
 import yaml
 from joblib import delayed
 from joblib import Parallel
-from matplotlib.backends.backend_pdf import PdfPages
 from tqdm import tqdm
 
-from morphio import SectionType
 from morphio.mut import Morphology
 from neuroc.scale import scale_section
 from neuroc.scale import ScaleParameters
-from neurom import load_neuron
-from neurom import viewer
 from neurom.core.dataformat import COLS
 from placement_algorithm.app.mpi_app import run_master
 from placement_algorithm.app.synthesize_morphologies import Master
 from tns import extract_input
-from tns import NeuronGrower
 from voxcell import CellCollection
 
+from . import STR_TO_TYPES
+
 
 L = logging.getLogger(__name__)
 matplotlib.use("Agg")
 
 
-STR_TO_TYPES = {
-    "basal": SectionType.basal_dendrite,
-    "apical": SectionType.apical_dendrite,
-    "axon": SectionType.axon,
-}
-
-
 def get_neurite_types(pc_in_types_path, mtypes):
     """Get the neurite types to consider for PC or IN cells."""
     with open(pc_in_types_path, "rb") as pc_in_file:
@@ -243,98 +232,6 @@ def _percentile(q, a, *args, **kwargs):
     return {mtype: float(f(lengths)) for mtype, lengths in apical_lengths.items()}
 
 
-def _grow_morphology(
-    gid, mtype, tmd_parameters, tmd_distributions, morphology_base_path
-):
-    """Grow single morphology for parallel computations."""
-
-    name = f"vacuum_{gid}.asc"
-    morphology_path = morphology_base_path / name
-    vacuum_synth_morphs_df = pd.DataFrame()
-    np.random.seed(gid)
-
-    grower = NeuronGrower(
-        input_parameters=tmd_parameters,
-        input_distributions=tmd_distributions,
-        external_diametrizer=None,
-    )
-    grower.grow()
-    grower.neuron.write(morphology_path)
-
-    vacuum_synth_morphs_df.loc[gid, "name"] = name
-    vacuum_synth_morphs_df.loc[gid, "mtype"] = mtype
-    vacuum_synth_morphs_df.loc[gid, "vacuum_morphology_path"] = morphology_path
-    # vacuum_synth_morphs_df.loc[gid, 'apical_point'] = grower.apical_points
-
-    return vacuum_synth_morphs_df
-
-
-def grow_vacuum_morphologies(
-    mtypes,
-    n_cells,
-    tmd_parameters,
-    tmd_distributions,
-    morphology_base_path,
-    joblib_verbose=0,
-):
-    """Grow morphologies in vacuum."""
-
-    global_gid = 0
-    vacuum_synth_morphs_df = pd.DataFrame()
-    for mtype in tqdm(mtypes):
-        # no need to realistic diameters here, using internal TNS diametrizer
-        tmd_parameters[mtype]["diameter_params"]["method"] = "M1"
-        tmd_distributions["mtypes"][mtype]["diameter"]["method"] = "M1"
-
-        gids = range(global_gid, global_gid + n_cells)
-        global_gid += n_cells
-        for row in Parallel(-1, verbose=joblib_verbose)(
-            delayed(_grow_morphology)(
-                gid,
-                mtype,
-                tmd_parameters[mtype],
-                tmd_distributions["mtypes"][mtype],
-                morphology_base_path,
-            )
-            for gid in gids
-        ):
-            vacuum_synth_morphs_df = vacuum_synth_morphs_df.append(row)
-    return vacuum_synth_morphs_df
-
-
-def plot_vacuum_morphologies(
-    vacuum_synth_morphs_df, pdf_filename, morphology_path, mean_lengths
-):
-    """Plot synthesized morphologies on top of each others."""
-    colors = {"apical": "m", "basal": "r", "axon": "b"}
-
-    with PdfPages(pdf_filename) as pdf:
-        for mtype in tqdm(sorted(vacuum_synth_morphs_df.mtype.unique())):
-            plt.figure()
-            ax = plt.gca()
-            for gid in vacuum_synth_morphs_df[
-                vacuum_synth_morphs_df.mtype == mtype
-            ].index:
-                morphology = load_neuron(
-                    vacuum_synth_morphs_df.loc[gid, morphology_path]
-                )
-                viewer.plot_neuron(ax, morphology, plane="zy")
-                morph = Morphology(vacuum_synth_morphs_df.loc[gid, morphology_path])
-                for neurite in morph.root_sections:
-                    if neurite.type == STR_TO_TYPES["apical"]:
-                        max_len = get_max_len(neurite)
-                        ax.axhline(max_len, c="0.5", lw=0.5)
-            for neurite_type in mean_lengths:
-                if mtype in mean_lengths[neurite_type]:
-                    ax.axhline(
-                        mean_lengths[neurite_type][mtype], c=colors[neurite_type]
-                    )
-            ax.set_title(mtype)
-            ax.set_rasterized(True)
-            plt.axis([-800, 800, -800, 2000])
-            pdf.savefig()
-
-
 def get_target_length(soma_layer, target_layer, cortical_thicknesses):
     """Compute the target length of a neurite from soma and target layer."""
     cortical_depths = np.insert(np.cumsum(cortical_thicknesses), 0, 0.0)

synthesis_workflow/tasks/circuit.py

@@ -0,0 +1,134 @@
+"""Luigi tasks for morphology synthesis."""
+from functools import partial
+from pathlib import Path
+
+import luigi
+import numpy as np
+
+from atlas_analysis.planes.planes import _create_planes
+from atlas_analysis.planes.planes import create_centerline_planes
+from atlas_analysis.planes.planes import save_planes_centerline
+from voxcell import VoxelData
+
+from ..circuit_slicing import generic_slicer
+from ..circuit_slicing import halve_atlas
+from ..circuit_slicing import slice_circuit
+from ..tools import ensure_dir
+from .config import circuitconfigs
+from .utils import BaseTask
+
+
+class CreateAtlasPlanes(BaseTask):
+    """Crerate plane cuts of an atlas."""
+
+    plane_type = luigi.Parameter(default="centerline")
+    plane_count = luigi.IntParameter(default=10)
+    inter_plane_count = luigi.IntParameter(default=200)
+
+    use_half = luigi.Parameter(default=True)
+    half_axis = luigi.IntParameter(default=0)
+    half_side = luigi.IntParameter(default=0)
+
+    tmp_nrrd_path = luigi.Parameter(default="layer_annotation.nrrd")
+
+    centerline_first_bound = luigi.ListParameter(default=[126, 181, 220])
+    centerline_last_bound = luigi.ListParameter(default=[407, 110, 66])
+
+    centerline_axis = luigi.IntParameter(default=0)
+    centerline_start = luigi.FloatParameter(3000)
+    centerline_end = luigi.FloatParameter(10000)
+
+    atlas_planes_path = luigi.Parameter(default="atlas_planes")
+
+    def run(self):
+        """"""
+        layer_annotation_path = Path(circuitconfigs().atlas_path) / "layers.nrrd"
+        layer_annotation = VoxelData.load_nrrd(layer_annotation_path)
+        if self.use_half:
+            layer_annotation.raw = halve_atlas(
+                layer_annotation.raw, axis=self.half_axis, side=self.half_side
+            )
+        ensure_dir(self.tmp_nrrd_path)
+        layer_annotation.save_nrrd(self.tmp_nrrd_path)
+
+        if self.plane_type == "centerline":
+            bounds = [self.centerline_first_bound, self.centerline_last_bound]
+            create_centerline_planes(
+                self.tmp_nrrd_path,
+                self.output().path,
+                bounds,
+                plane_count=self.inter_plane_count,
+            )
+
+        if self.plane_type == "aligned":
+            centerline = np.zeros([100, 3])
+            centerline[:, self.centerline_axis] = np.linspace(
+                self.centerline_start, self.centerline_end, 100
+            )
+            planes = _create_planes(centerline, plane_count=self.inter_plane_count)
+            save_planes_centerline(self.output().path, planes, centerline)
+
+        all_planes = np.load(self.output().path)
+        selected_planes = []
+        di = int(self.inter_plane_count / self.plane_count)
+        for i in range(self.plane_count):
+            selected_planes.append(all_planes["planes"][di * i])
+            selected_planes.append(all_planes["planes"][di * i + 1])
+        np.savez(
+            self.output().path,
+            planes=np.array(selected_planes),
+            centerline=all_planes["centerline"],
+            plane_format=all_planes["plane_format"],
+        )
+
+    def output(self):
+        """"""
+        return luigi.LocalTarget(self.atlas_planes_path + ".npz")
+
+
+class SliceCircuit(BaseTask):
+    """Create a smaller circuit .mvd3 file for subsampling.
+
+    Args:
+        sliced_circuit_path (str): path to save sliced circuit somata mvd3
+        mtypes (list): list of mtypes to consider
+        n_cells (int): number of cells per mtype to consider
+        hemisphere (str): 'left' or 'right'
+    """
+
+    sliced_circuit_path = luigi.Parameter(default="sliced_circuit_somata.mvd3")
+    mtypes = luigi.ListParameter(default=None)
+    n_cells = luigi.IntParameter(default=10)
+    hemisphere = luigi.Parameter(default=None)
+
+    def run(self):
+        """"""
+        mtypes = self.mtypes
+        if "all" in mtypes:  # pylint: disable=unsupported-membership-test
+            mtypes = None
+
+        if self.hemisphere is not None:
+            atlas_planes = yield CreateAtlasPlanes()
+            planes = np.load(atlas_planes.path)["planes"]
+        else:
+            planes = None
+
+        slicer = partial(
+            generic_slicer,
+            n_cells=self.n_cells,
+            mtypes=mtypes,
+            planes=planes,
+            hemisphere=self.hemisphere,
+        )
+
+        ensure_dir(self.output().path)
+        cells = slice_circuit(
+            circuitconfigs().circuit_somata_path, self.output().path, slicer
+        )
+
+        if len(cells.index) == 0:
+            raise Exception("No cells will be synthtesized, better stop here.")
+
+    def output(self):
+        """"""
+        return luigi.LocalTarget(self.sliced_circuit_path)

synthesis_workflow/utils_tasks.py → synthesis_workflow/tasks/config.py

@@ -1,10 +1,8 @@
-"""utils functions for luigi tasks."""
+"""config functions for luigi tasks."""
 import logging
 import warnings
-from pathlib import Path
 
 import luigi
-import pandas as pd
 
 
 logger = logging.getLogger("luigi-interface")
@@ -16,50 +14,6 @@
 warnings.filterwarnings("ignore", module="scipy")
 
 
-def load_circuit(
-    path_to_mvd3=None,
-    path_to_morphologies=None,
-    path_to_atlas=None,
-    circuit_config=None,
-):
-    """Load a circuit with bluepy.v2."""
-    from bluepy.v2 import Circuit
-
-    if circuit_config:
-        return Circuit(circuit_config)
-    return Circuit(
-        {
-            "cells": path_to_mvd3,
-            "morphologies": path_to_morphologies,
-            "atlas": path_to_atlas,
-        }
-    )
-
-
-def ensure_dir(file_path):
-    """Create directory to save file."""
-    Path(file_path).parent.mkdir(parents=True, exist_ok=True)
-
-
-def get_morphs_df(
-    morphs_df_path, to_use_flag="all", morphology_path="morphology_path", h5_path=None,
-):
-    """Get valid morphs_df for diametrizer (select using flag + remove duplicates)."""
-    morphs_df = pd.read_csv(morphs_df_path)
-    if to_use_flag != "all":
-        morphs_df = morphs_df.loc[morphs_df[to_use_flag]]
-    if h5_path is not None:
-        morphs_df[morphology_path] = morphs_df[
-            "_".join(morphology_path.split("_")[:-1])
-        ].apply(lambda path: str((Path(h5_path) / Path(path).stem).with_suffix(".h5")))
-    return morphs_df[["name", "mtype", morphology_path]].drop_duplicates()
-
-
-def update_morphs_df(morphs_df_path, new_morphs_df):
-    """Update a morphs_df with new entries to preserve duplicates."""
-    return pd.read_csv(morphs_df_path).merge(new_morphs_df, how="left")
-
-
 class diametrizerconfigs(luigi.Config):
     """Diametrizer configuration."""
 
@@ -126,42 +80,3 @@ class pathconfigs(luigi.Config):
     synth_morphs_df_path = luigi.Parameter(default="synth_morphs_df.csv")
     synth_output_path = luigi.Parameter(default="synthesized_morphologies")
     substituted_morphs_df_path = luigi.Parameter(default="substituted_morphs_df.csv")
-
-
-@luigi.Task.event_handler(luigi.Event.START)
-def log_parameters(task):
-    class_name = task.__class__.__name__
-    logger.debug("Attributes of {} class after global processing:".format(class_name))
-    for name in task.get_param_names():
-        try:
-            logger.debug("Atribute: {} == {}".format(name, getattr(task, name)))
-        except Exception:
-            logger.debug("Can't print '{}' attribute for unknown reason".format(name))
-
-
-class BaseTask(luigi.Task):
-    """Base class used to add customisable global parameters"""
-    _global_configs = [diametrizerconfigs, synthesisconfigs, circuitconfigs, pathconfigs]
-
-    def __getattribute__(self, name):
-        tmp = super(BaseTask, self).__getattribute__(name)
-        if tmp is not None:
-            return tmp
-        for conf in self._global_configs:
-            tmp_conf = conf()
-            if hasattr(tmp_conf, name):
-                return getattr(tmp_conf, name)
-        return tmp
-
-    def __setattr__(self, name, value):
-        if value is None and name in self.get_param_names():
-            msg = (
-                "The Parameter '{}' of the task '{}' is set to None, thus the global "
-                "value will be taken frow now on"
-            ).format(name, self.__class__.__name__)
-            warnings.warn(msg)
-        return super(BaseTask, self).__setattr__(name, value)
-
-
-class BaseWrapperTask(BaseTask, luigi.WrapperTask):
-    pass