Expose subtree processing from the morph_stats api (#1034)

BlueBrain · May 12, 2022 · 5901eb2 · 5901eb2
1 parent 006cafc
commit 5901eb2
Show file tree

Hide file tree

Showing 5 changed files with 194 additions and 16 deletions.
diff --git a/CHANGELOG.rst b/CHANGELOG.rst
@@ -4,6 +4,7 @@ Changelog
 Version 4.0.0
 -------------
 
+- Mixed subtree processing can be used in morph_stats app via the use_subtrees flag.
 - ``neurom.view.[plot_tree|plot_tree3d|plot_soma|plot_soma3D]`` were hidden from the
   neurom.view module. They can still be imported from neurom.view.matplotlib_impl.
 - Deprecated modules and classes were removed.

diff --git a/neurom/apps/cli.py b/neurom/apps/cli.py
@@ -95,9 +95,13 @@ def view(input_file, is_3d, plane, backend, realistic_diameters):
               help='If enabled the directory is treated as a population')
 @click.option('-I', '--ignored-exceptions', help='Exception to ignore',
               type=click.Choice(morph_stats.IGNORABLE_EXCEPTIONS.keys()))
-def stats(datapath, config, output, full_config, as_population, ignored_exceptions):
+@click.option('--use-subtrees', is_flag=True, show_default=True, default=False,
+              help="Enable mixed subtree processing.")
+def stats(datapath, config, output, full_config, as_population, ignored_exceptions, use_subtrees):
     """Cli for apps/morph_stats."""
-    morph_stats.main(datapath, config, output, full_config, as_population, ignored_exceptions)
+    morph_stats.main(
+        datapath, config, output, full_config, as_population, ignored_exceptions, use_subtrees
+    )
 
 
 @cli.command(short_help='Perform checks on morphologies, more details at'

diff --git a/neurom/apps/morph_stats.py b/neurom/apps/morph_stats.py
@@ -59,14 +59,14 @@
 IGNORABLE_EXCEPTIONS = {'SomaError': SomaError}
 
 
-def _run_extract_stats(morph, config):
+def _run_extract_stats(morph, config, use_subtrees=False):
     """The function to be called by multiprocessing.Pool.imap_unordered."""
     if not isinstance(morph, Morphology):
         morph = nm.load_morphology(morph)
-    return morph.name, extract_stats(morph, config)
+    return morph.name, extract_stats(morph, config, use_subtrees=use_subtrees)
 
 
-def extract_dataframe(morphs, config, n_workers=1):
+def extract_dataframe(morphs, config, n_workers=1, use_subtrees=False):
     """Extract stats grouped by neurite type from morphs.
 
     Arguments:
@@ -83,6 +83,7 @@ def extract_dataframe(morphs, config, n_workers=1):
             - morphology: same as neurite entry, but it will not be run on each neurite_type,
               but only once on the whole morphology.
         n_workers (int): number of workers for multiprocessing (on collection of morphs)
+        use_subtrees (bool): Enable of heterogeneous subtree processing.
 
     Returns:
         The extracted statistics
@@ -94,7 +95,7 @@ def extract_dataframe(morphs, config, n_workers=1):
     if isinstance(morphs, Morphology):
         morphs = [morphs]
 
-    func = partial(_run_extract_stats, config=config)
+    func = partial(_run_extract_stats, config=config, use_subtrees=use_subtrees)
     if n_workers == 1:
         stats = list(map(func, morphs))
     else:
@@ -114,12 +115,12 @@ def extract_dataframe(morphs, config, n_workers=1):
 extract_dataframe.__doc__ += str(EXAMPLE_CONFIG)
 
 
-def _get_feature_stats(feature_name, morphs, modes, kwargs):
+def _get_feature_stats(feature_name, morphs, modes, use_subtrees=False, **kwargs):
     """Insert the stat data in the dict.
 
     If the feature is 2-dimensional, the feature is flattened on its last axis
     """
-    def stat_name_format(mode, feature_name, kwargs):
+    def stat_name_format(mode, feature_name, **kwargs):
         """Returns the key name for the data dictionary.
 
         The key is a combination of the mode, feature_name and an optional suffix of all the extra
@@ -135,14 +136,16 @@ def stat_name_format(mode, feature_name, kwargs):
         return f"{mode}_{feature_name}"
 
     data = {}
-    value, func = _get_feature_value_and_func(feature_name, morphs, **kwargs)
+    value, func = _get_feature_value_and_func(
+        feature_name, morphs, use_subtrees=use_subtrees, **kwargs
+    )
     shape = func.shape
     if len(shape) > 2:
         raise ValueError(f'Len of "{feature_name}" feature shape must be <= 2')  # pragma: no cover
 
     for mode in modes:
 
-        stat_name = stat_name_format(mode, feature_name, kwargs)
+        stat_name = stat_name_format(mode, feature_name, **kwargs)
 
         stat = value
         if isinstance(value, Sized):
@@ -161,7 +164,7 @@ def stat_name_format(mode, feature_name, kwargs):
     return data
 
 
-def extract_stats(morphs, config):
+def extract_stats(morphs, config, use_subtrees=False):
     """Extract stats from morphs.
 
     Arguments:
@@ -180,6 +183,7 @@ def extract_stats(morphs, config):
                   ['min', 'max', 'median', 'mean', 'std', 'raw', 'sum']
             - morphology: same as neurite entry, but it will not be run on each neurite_type,
               but only once on the whole morphology.
+        use_subtrees (bool): Enable of heterogeneous subtree processing.
 
     Returns:
         The extracted statistics
@@ -215,12 +219,18 @@ def extract_stats(morphs, config):
                     for neurite_type in types:
                         feature_kwargs["neurite_type"] = neurite_type
                         stats[neurite_type.name].update(
-                            _get_feature_stats(feature_name, morphs, modes, feature_kwargs)
+                            _get_feature_stats(
+                                feature_name, morphs, modes,
+                                use_subtrees=use_subtrees, **feature_kwargs
+                            )
                         )
 
                 else:
                     stats[category].update(
-                        _get_feature_stats(feature_name, morphs, modes, feature_kwargs)
+                        _get_feature_stats(
+                            feature_name, morphs, modes,
+                            use_subtrees=use_subtrees, **feature_kwargs
+                        )
                     )
 
     return dict(stats)
@@ -347,7 +357,15 @@ def _sanitize_config(config):
     return config
 
 
-def main(datapath, config, output_file, is_full_config, as_population, ignored_exceptions):
+def main(
+    datapath,
+    config,
+    output_file,
+    is_full_config,
+    as_population,
+    ignored_exceptions,
+    use_subtrees=False
+):
     """Main function that get statistics for morphologies.
 
     Args:
@@ -357,6 +375,7 @@ def main(datapath, config, output_file, is_full_config, as_population, ignored_e
         is_full_config (bool): should be statistics made over all possible features, modes, neurites
         as_population (bool): treat ``datapath`` as directory of morphologies population
         ignored_exceptions (list|tuple|None): exceptions to ignore when loading a morphology
+        use_subtrees (bool): Enable of heterogeneous subtree processing
     """
     config = full_config() if is_full_config else get_config(config, EXAMPLE_CONFIG)
 
@@ -374,9 +393,9 @@ def main(datapath, config, output_file, is_full_config, as_population, ignored_e
     )
 
     if as_population:
-        results = {datapath: extract_stats(morphs, config)}
+        results = {datapath: extract_stats(morphs, config, use_subtrees=use_subtrees)}
     else:
-        results = {m.name: extract_stats(m, config) for m in morphs}
+        results = {m.name: extract_stats(m, config, use_subtrees=use_subtrees) for m in morphs}
 
     if not output_file:
         print(json.dumps(results, indent=2, separators=(',', ':'), cls=NeuromJSON))

diff --git a/tests/apps/test_cli.py b/tests/apps/test_cli.py
@@ -83,6 +83,16 @@ def test_morph_stat_full_config():
         assert not df.empty
 
 
+def test_morph_stat_full_config__subtrees():
+    runner = CliRunner()
+    filename = DATA / 'h5/v1/Neuron.h5'
+    with tempfile.NamedTemporaryFile() as f:
+        result = runner.invoke(cli, ['stats', str(filename), '--full-config', '--use-subtrees', '--output', f.name])
+        assert result.exit_code == 0
+        df = pd.read_csv(f)
+        assert not df.empty
+
+
 def test_morph_stat_invalid_config():
     runner = CliRunner()
     with tempfile.NamedTemporaryFile('w') as config_f:

diff --git a/tests/test_mixed.py b/tests/test_mixed.py
@@ -3,6 +3,7 @@
 import pytest
 import neurom
 import numpy as np
+import pandas as pd
 import numpy.testing as npt
 from neurom import NeuriteType
 from neurom.features import get
@@ -14,6 +15,7 @@
 
 import neurom.core.morphology
 import neurom.features.neurite
+import neurom.apps.morph_stats
 
 
 @pytest.fixture
@@ -240,6 +242,148 @@ def assert_sections(neurite, section_type, iterator_type, expected_section_ids):
     )
 
 
+def test_mixed_morph_stats(mixed_morph):
+
+    def assert_stats_equal(actual_dict, expected_dict):
+        assert actual_dict.keys() == expected_dict.keys()
+        for (key, value) in actual_dict.items():
+            expected_value = expected_dict[key]
+            if value is None or expected_value is None:
+                assert expected_value is value
+            else:
+                npt.assert_almost_equal(value, expected_value, decimal=3, err_msg=f"\nKey: {key}")
+
+    cfg = {
+        'neurite': {
+            'max_radial_distance': ['mean'],
+            'number_of_sections': ['min'],
+            'number_of_bifurcations': ['max'],
+            'number_of_leaves': ['median'],
+            'total_length': ['min'],
+            'total_area': ['max'],
+            'total_volume': ['median'],
+            'section_lengths': ['mean'],
+            'section_term_lengths': ['mean'],
+            'section_bif_lengths': ['mean'],
+            'section_branch_orders': ['mean'],
+            'section_bif_branch_orders': ['mean'],
+            'section_term_branch_orders': ['mean'],
+            'section_path_distances': ['mean'],
+            'section_taper_rates': ['median'],
+            'local_bifurcation_angles': ['mean'],
+            'remote_bifurcation_angles': ['mean'],
+            'partition_asymmetry': ['mean'],
+            'partition_asymmetry_length': ['mean'],
+            'sibling_ratios': ['mean'],
+            'diameter_power_relations': ['median'],
+            'section_radial_distances': ['mean'],
+            'section_term_radial_distances': ['mean'],
+            'section_bif_radial_distances': ['mean'],
+            'terminal_path_lengths': ['mean'],
+            'section_volumes': ['min'],
+            'section_areas': ['mean'],
+            'section_tortuosity': ['mean'],
+            'section_strahler_orders': ['min']
+        },
+        'morphology': {
+            'soma_surface_area': ['mean'],
+            'soma_radius': ['max'],
+            'max_radial_distance': ['mean'],
+            'number_of_sections_per_neurite': ['median'],
+            'total_length_per_neurite': ['mean'],
+            'total_area_per_neurite': ['mean'],
+            'total_volume_per_neurite': ['mean'],
+            'number_of_neurites': ['median']
+        },
+        'neurite_type': ['AXON', 'BASAL_DENDRITE', 'APICAL_DENDRITE']
+    }
+
+    res = neurom.apps.morph_stats.extract_stats(mixed_morph, cfg, use_subtrees=False)
+
+    expected_axon_wout_subtrees = {
+        'max_number_of_bifurcations': 0,
+        'max_total_area': 0,
+        'mean_local_bifurcation_angles': None,
+        'mean_max_radial_distance': 0.0,
+        'mean_partition_asymmetry': None,
+        'mean_partition_asymmetry_length': None,
+        'mean_remote_bifurcation_angles': None,
+        'mean_section_areas': None,
+        'mean_section_bif_branch_orders': None,
+        'mean_section_bif_lengths': None,
+        'mean_section_bif_radial_distances': None,
+        'mean_section_branch_orders': None,
+        'mean_section_lengths': None,
+        'mean_section_path_distances': None,
+        'mean_section_radial_distances': None,
+        'mean_section_term_branch_orders': None,
+        'mean_section_term_lengths': None,
+        'mean_section_term_radial_distances': None,
+        'mean_section_tortuosity': None,
+        'mean_sibling_ratios': None,
+        'mean_terminal_path_lengths': None,
+        'median_diameter_power_relations': None,
+        'median_number_of_leaves': 0,
+        'median_section_taper_rates': None,
+        'median_total_volume': 0,
+        'min_number_of_sections': 0,
+        'min_section_strahler_orders': None,
+        'min_section_volumes': None,
+        'min_total_length': 0
+    }
+
+    assert_stats_equal(res["axon"], expected_axon_wout_subtrees)
+
+    res_df = neurom.apps.morph_stats.extract_dataframe(mixed_morph, cfg, use_subtrees=False)
+
+    # get axon column and tranform it to look like the expected values above
+    values = res_df.loc[pd.IndexSlice[:, "axon"]].iloc[0, :].to_dict()
+    assert_stats_equal(values, expected_axon_wout_subtrees)
+
+
+    res = neurom.apps.morph_stats.extract_stats(mixed_morph, cfg, use_subtrees=True)
+
+    expected_axon_with_subtrees = {
+        'max_number_of_bifurcations': 2,
+        'max_total_area': 3.4018507611950346,
+        'mean_local_bifurcation_angles': 2.356194490192345,
+        'mean_max_radial_distance': 4.472136,
+        'mean_partition_asymmetry': 0.25,
+        'mean_partition_asymmetry_length': 0.1846990320847273,
+        'mean_remote_bifurcation_angles': 2.356194490192345,
+        'mean_section_areas': 0.6803701522390069,
+        'mean_section_bif_branch_orders': 1.5,
+        'mean_section_bif_lengths': 1.2071068,
+        'mean_section_bif_radial_distances': 3.9240959,
+        'mean_section_branch_orders': 2.2,
+        'mean_section_lengths': 1.0828427,
+        'mean_section_path_distances': 2.614213538169861,
+        'mean_section_radial_distances': 4.207625,
+        'mean_section_term_branch_orders': 2.6666666666666665,
+        'mean_section_term_lengths': 1.0,
+        'mean_section_term_radial_distances': 4.396645,
+        'mean_section_tortuosity': 1.0,
+        'mean_sibling_ratios': 1.0,
+        'mean_terminal_path_lengths': 3.0808802048365274,
+        'median_diameter_power_relations': 2.0,
+        'median_number_of_leaves': 3,
+        'median_section_taper_rates': 8.6268466e-17,
+        'median_total_volume': 0.17009254152367845,
+        'min_number_of_sections': 5,
+        'min_section_strahler_orders': 1,
+        'min_section_volumes': 0.03141592778425469,
+        'min_total_length': 5.414213538169861
+    }
+
+    assert_stats_equal(res["axon"], expected_axon_with_subtrees)
+
+    res_df = neurom.apps.morph_stats.extract_dataframe(mixed_morph, cfg, use_subtrees=True)
+
+    # get axon column and tranform it to look like the expected values above
+    values = res_df.loc[pd.IndexSlice[:, "axon"]].iloc[0, :].to_dict()
+    assert_stats_equal(values, expected_axon_with_subtrees)
+
+
 @pytest.fixture
 def population(mixed_morph):
     return Population([mixed_morph, mixed_morph])