Improve test coverage and minor fixes

Change-Id: I496bf2e519961dc5f7cf9459749d87eb2e02772f

.gitignore

@@ -80,6 +80,7 @@ nosetests.xml
 coverage.xml
 *.cover
 .hypothesis/
+reports/
 
 # Translations
 *.mo

.pylintrc

@@ -67,5 +67,6 @@ ignore-docstrings=yes
 # (useful for modules/projects where namespaces are manipulated during runtime
 # and thus existing member attributes cannot be deduced by static analysis. It
 # supports qualified module names, as well as Unix pattern matching.
-ignored-modules=numpy,numpy.*,scipy.stats,scipy.spatial,morphio,neurom,seaborn
+ignored-modules=numpy,numpy.*,scipy.stats,scipy.spatial,morphio,neurom,seaborn,morphio,morphio.mut
+
 extension-pkg-whitelist=morphio

diameter_synthesis/build_diameters.py

@@ -1,6 +1,5 @@
 """Build neurite diameters from a pre-generated model. TO REVIEW!."""
 import logging
-import random
 from collections import deque
 from functools import partial
 
@@ -71,6 +70,7 @@ def _sample_sibling_ratio(
         if asymmetry_value > asymmetry_threshold:
             return 0.0
         return sample_distribution(params["sibling_ratios"][neurite_type])
+    # This case should never happen since the mode is already checked in `_select_model`
     raise DiameterSynthesisError("mode not understood {}".format(mode))
 
 
@@ -96,7 +96,9 @@ def _sample_diameter_power_relation(
             return 1.0
         return sample_distribution(params["diameter_power_relation"][neurite_type])
     if mode == "exact":
+        # This case should never happen since this mode is not known by `_select_model`
         return 1.0
+    # This case should never happen since the mode is already checked in `_select_model`
     raise DiameterSynthesisError("mode not understood {}".format(mode))
 
 
@@ -198,12 +200,15 @@ def _sample_daughter_diameters(section, params, params_tree):
 
     diams = [diam_1] + (len(section.children) - 1) * [diam_2]
     if params_tree["with_asymmetry"]:
+        # This case should always happen since the `with_asymmetry` attribute is always set to True
+        # in `_select_model`
+
         # returns child diameters sorted by child length (major/secondary for apical tree)
         child_sort = np.argsort(
             [morph_funcs.n_children_downstream(child) for child in section.children]
         )[::-1]
         return list(np.array(diams)[child_sort])
-    random.shuffle(diams)
+    np.random.shuffle(diams)
     return diams
 
 

diameter_synthesis/cli.py

@@ -48,7 +48,8 @@ def run_diameters(config_file, models_params_file):
 @click.argument("diametrized_folder", type=click.Path(exists=True))
 @click.argument("plot_folder", type=click.Path())
 @click.option("-n", "--ncells", help="max number of cells to plot")
-def plot_diff(original_folder, diametrized_folder, plot_folder, ncells=None):
+@click.option("-e", "--ext", help="figures extention")
+def plot_diff(original_folder, diametrized_folder, plot_folder, ncells=None, ext=".png"):
     """Plot original and new neurons as well as their differences."""
     from .plotting import plot_diameter_diff
 
@@ -76,7 +77,7 @@ def plot_diff(original_folder, diametrized_folder, plot_folder, ncells=None):
             neuron_new = nm.load_neuron(Path(diametrized_folder) / neuron.name)
 
             plot_diameter_diff(
-                neuron, neuron_new, neurite_types, plot_folder_mtype, ext=".png",
+                neuron, neuron_new, neurite_types, plot_folder_mtype, ext=ext,
             )
 
 
@@ -90,10 +91,13 @@ def plot_diff(original_folder, diametrized_folder, plot_folder, ncells=None):
 @click.option("--cumulative", help="Cumulative distribution plots", is_flag=True)
 @click.option("--individual", help="Output a plot for each neuron", is_flag=True)
 @click.option("--violin", help="Violin distribution plots", is_flag=True)
+@click.option("-e", "--ext", help="Figures extention")
 def run_analysis(
-    orig_path, diam_path, out_dir, cumulative, individual, violin, mtypes_file=None
+    orig_path, diam_path, out_dir, cumulative, individual, violin, mtypes_file=None, ext=".png"
 ):
     """Produce figures for validation/analysis."""
+    if not cumulative and not violin:
+        raise ValueError("Should at least set one of '--cumulative' or '--violin' to True.")
     if cumulative:
         from .plotting import cumulative_analysis
 
@@ -104,6 +108,7 @@ def run_analysis(
             individual,
             mtypes_file=mtypes_file,
             neurite_types=["basal"],
+            ext=ext,
         )
         cumulative_analysis(
             orig_path,
@@ -112,6 +117,7 @@ def run_analysis(
             individual,
             mtypes_file=mtypes_file,
             neurite_types=["axon"],
+            ext=ext,
         )
         cumulative_analysis(
             orig_path,
@@ -120,6 +126,7 @@ def run_analysis(
             individual,
             mtypes_file=mtypes_file,
             neurite_types=["apical"],
+            ext=ext,
         )
     if violin:
         from .plotting import violin_analysis

diameter_synthesis/distribution_fitting.py

@@ -15,15 +15,24 @@
 PERCENTILE = 5
 MIN_SAMPLE_NUM = 10
 
+MAX_TRUNCATE_TRIES = 100
+
 L = logging.getLogger(__name__)
 np.seterr(invalid="ignore", divide="ignore")
 
 
 def _truncate(sample_func, min_value, max_value):
     """Ensure sample is within bounds."""
     sample = sample_func()
+    n_tries = 0
     while sample > max_value or sample < min_value:
         sample = sample_func()
+        n_tries += 1
+        if n_tries >= MAX_TRUNCATE_TRIES:
+            raise DiameterSynthesisError(
+                f"Could not truncate the sample between {min_value} and {max_value} "
+                f"(the last value was {sample})"
+            )
     return sample
 
 

diameter_synthesis/morph_functions.py

@@ -170,7 +170,7 @@ def trunk_diameter(neurite, attribute_name=None, bounds=None, method="last"):
     )
 
 
-def taper(neurite, params=None, attribute_name=None):
+def taper(neurite, params, attribute_name=None):
     """Get the taper rates (neurom only)."""
     if params is None:
         raise DiameterSynthesisError(
@@ -228,7 +228,7 @@ def get_additional_attribute(
             return nm.features.bifurcationfunc.sibling_ratio(section, method="mean")
 
     raise DiameterSynthesisError(
-        "Please provide either a neurite or a section, not both"
+        "Please provide a valid attribute_name and provide either a neurite or a section, not both"
     )
 
 

diameter_synthesis/plotting.py

@@ -19,7 +19,6 @@
 
 from diameter_synthesis import utils
 from diameter_synthesis.distribution_fitting import evaluate_distribution
-from diameter_synthesis.build_diameters import STR_TO_TYPES
 
 # pylint: disable=too-many-statements,too-many-locals,too-many-arguments
 
@@ -80,6 +79,12 @@
     "Remote bif angles",
 ]
 
+NEURITE_STR_TO_TYPES = {
+    "basal": nm.NeuriteType.basal_dendrite,
+    "apical": nm.NeuriteType.apical_dendrite,
+    "axon": nm.NeuriteType.axon,
+}
+
 
 def _compute_neurite_diff(
     neuron_orig, neuron_new, neuron_diff_pos, neuron_diff_neg, neurite_types
@@ -89,22 +94,22 @@ def _compute_neurite_diff(
         neurites_orig = [
             neurite
             for neurite in neuron_orig.neurites
-            if neurite.type == STR_TO_TYPES[neurite_type]
+            if neurite.type == NEURITE_STR_TO_TYPES[neurite_type]
         ]
         neurites_new = [
             neurite
             for neurite in neuron_new.neurites
-            if neurite.type == STR_TO_TYPES[neurite_type]
+            if neurite.type == NEURITE_STR_TO_TYPES[neurite_type]
         ]
         neurites_diff_neg = [
             neurite
             for neurite in neuron_diff_neg.neurites
-            if neurite.type == STR_TO_TYPES[neurite_type]
+            if neurite.type == NEURITE_STR_TO_TYPES[neurite_type]
         ]
         neurites_diff_pos = [
             neurite
             for neurite in neuron_diff_pos.neurites
-            if neurite.type == STR_TO_TYPES[neurite_type]
+            if neurite.type == NEURITE_STR_TO_TYPES[neurite_type]
         ]
 
         for neurite_orig, neurite_new, neurite_diff_pos, neurite_diff_neg in zip(
@@ -122,13 +127,13 @@ def _compute_neurite_diff(
                 diff = diam_new[j] - diam_orig[j]
                 diff_pos = diff.copy()
                 diff_pos[diff_pos < 0] = 0
-                utils._set_diameters(section, diff_pos)
+                section.points[:, nm.COLS.R] = diff_pos
 
             for j, section in enumerate(iter_sections(neurite_diff_neg)):
                 diff = diam_new[j] - diam_orig[j]
                 diff_neg = -diff.copy()
                 diff_neg[diff_neg < 0] = 0
-                utils._set_diameters(section, diff_neg)
+                section.points[:, nm.COLS.R] = diff_neg
 
 
 def plot_diameter_diff(neuron_name, neuron_new, neurite_types, folder, ext=".png"):
@@ -263,7 +268,7 @@ def _create_data(
     feature1, feature2, original_cells, diametrized_cells, step_size, neurite_types
 ):  # noqa, pylint: disable=too-many-locals,too-many-arguments
     def feature_data(cell, neurite_type):
-        nm_neurite_type = STR_TO_TYPES[neurite_type]
+        nm_neurite_type = NEURITE_STR_TO_TYPES[neurite_type]
         return [
             get(feat, cell, neurite_type=nm_neurite_type)
             for feat in (feature1, feature2)
@@ -480,8 +485,12 @@ def make_cumulative_figures(
     out_dir,
     individual=False,
     figname_prefix="",
+    ext=".png",
 ):
     """Make plots for cumulative distributions for a pair of features."""
+    out_dir = Path(out_dir)
+    out_dir.mkdir(parents=True, exist_ok=True)
+
     prefix1, basename1 = _split_prefix(feature1)
     prefix2, basename2 = _split_prefix(feature2)
 
@@ -495,12 +504,12 @@ def make_cumulative_figures(
         prefix1, basename1, basename2
     )
 
-    fig.savefig(Path(out_dir) / (figure_name + ".png"), bbox_inches="tight")
+    fig.savefig(out_dir / (figure_name + ext), bbox_inches="tight")
     plt.close(fig)
 
     if individual:
-        if not (Path(out_dir) / (figure_name + "_individual")).exists():
-            os.mkdir(Path(out_dir) / (figure_name + "_individual"))
+        if not (out_dir / (figure_name + "_individual")).exists():
+            os.mkdir(out_dir / (figure_name + "_individual"))
 
         for i, (original_cell, diametrized_cell) in enumerate(
             zip(original_cells, diametrized_cells)
@@ -513,9 +522,9 @@ def make_cumulative_figures(
                 neurite_types,
                 auto_limit=False,
             )
-            fname = "{}_{}.png".format(figure_name, original_cell.name)
+            fname = "{}_{}{}".format(figure_name, original_cell.name, ext)
             f.savefig(
-                Path(out_dir) / (figure_name + "_individual") / (str(i) + "_" + fname),
+                out_dir / (figure_name + "_individual") / (str(i) + "_" + fname),
                 bbox_inches="tight",
             )
             plt.close(f)
@@ -535,13 +544,14 @@ def cumulative_analysis(
     original_path,
     diametrized_path,
     out_dir,
-    individual,
+    individual=False,
     mtypes_file=None,
     neurite_types=None,
+    ext=".png",
 ):
     """Make plots for cumulative distributions."""
-    if not Path(out_dir).exists():
-        os.mkdir(out_dir)
+    out_dir = Path(out_dir)
+    out_dir.mkdir(parents=True, exist_ok=True)
 
     all_original_cells = _load_morphologies(original_path, mtypes_file=mtypes_file)
     all_diametrized_cells = _load_morphologies(
@@ -560,6 +570,7 @@ def cumulative_analysis(
                 out_dir,
                 individual=individual,
                 figname_prefix=mtype,
+                ext=ext,
             )
 
 
@@ -630,11 +641,16 @@ def plot_violins(data, x="Morphological features", y="Values", hues="Data", ax=N
 
 
 def violin_analysis(
-    original_path, diametrized_path, out_dir, mtypes_file=None, max_cells=200
+    original_path,
+    diametrized_path,
+    out_dir,
+    mtypes_file=None,
+    max_cells=200,
+    with_axon=False,
 ):
     """Plot violin distributions."""
-    if not Path(out_dir).exists():
-        os.mkdir(out_dir)
+    out_dir = Path(out_dir)
+    out_dir.mkdir(parents=True, exist_ok=True)
 
     orig_morphologies_dict = utils.create_morphologies_dict(
         original_path, mtypes_file=mtypes_file
@@ -647,27 +663,28 @@ def violin_analysis(
         [orig_morphologies_dict[mtype], diametrized_morphologies_dict[mtype], mtype]
         for mtype in orig_morphologies_dict
     ]
-    analyze_from_dict = partial(_analyze_from_dict, out_dir, max_cells)
-    with multiprocessing.Pool() as pool:
+    analyze_from_dict = partial(_analyze_from_dict, max_cells, with_axon=with_axon)
+
+    pool = multiprocessing.Pool()
+    try:
         figs = list(
             tqdm(
                 pool.imap_unordered(analyze_from_dict, cells_data),
                 total=len(cells_data),
             )
         )
+    finally:
+        pool.close()
+        pool.join()
 
-    with PdfPages("morphometrics.pdf") as pdf:
+    with PdfPages(out_dir / "morphometrics.pdf") as pdf:
         for mtype, fig in figs:
             if mtype is not None:
-                # fig.suptitle(mtype)
                 pdf.savefig(fig)
 
 
-def _analyze_from_dict(out_dir, max_cells, cells, with_axon=False):
+def _analyze_from_dict(max_cells, cells, with_axon=False):
     cell_orig, cell_diametrized, mtype = cells
-    print(mtype)
-    if mtype != "L6_TPC:C":
-        return None, None
     cell_diametrized = cell_diametrized[:max_cells]
     cell_orig = cell_orig[:max_cells]
     original_cells = nm.load_neurons(cell_orig)
@@ -692,9 +709,11 @@ def _analyze_from_dict(out_dir, max_cells, cells, with_axon=False):
     axes[0].set_ylim(-3, 5)
     axes[0].title.set_text("basal dendrites")
 
-    if nm.check.neuron_checks.has_apical_dendrite(
-        original_cells[0]
-    ) and nm.check.neuron_checks.has_apical_dendrite(diametrized_cells[0]):
+    if any(
+        [i.type == nm.NeuriteType.apical_dendrite for i in original_cells[0].neurites]
+    ) and any(
+        [i.type == nm.NeuriteType.apical_dendrite for i in diametrized_cells[0].neurites]
+    ):
         data = get_features_all(
             original_cells,
             diametrized_cells,
@@ -707,8 +726,6 @@ def _analyze_from_dict(out_dir, max_cells, cells, with_axon=False):
         axes[1].set_ylim(-3, 5)
         axes[1].title.set_text("apical dendrites")
 
-    fig.tight_layout()
-    fig.savefig(Path(out_dir) / ("violin_" + mtype + ".png"))
     if with_axon:
         data = get_features_all(
             original_cells,
@@ -719,4 +736,5 @@ def _analyze_from_dict(out_dir, max_cells, cells, with_axon=False):
         data_frame = transform2DataFrame(data, pop_names, flist=VIOLIN_FEATURES_NAME)
         plot_violins(data_frame.replace([np.inf, -np.inf], np.nan).dropna(), ax=axes[2])
         axes[2].set_ylim(-3, 5)
+
     return mtype, fig