From cb18a104ff0fa6b86e2a1ca24233cb3499ae790b Mon Sep 17 00:00:00 2001 From: Eleftherios Zisis Date: Fri, 26 Jan 2024 16:08:45 +0100 Subject: [PATCH] Add testing of example scripts (#1041) (#1096) --- examples/boxplot.py | 26 +++- examples/density_plot.py | 31 ++++- examples/end_to_end_distance.py | 15 ++- examples/extract_distribution.py | 38 +++--- examples/features_graph_table.py | 48 +++---- examples/get_features.py | 13 +- examples/histogram.py | 30 ++++- examples/iteration_analysis.py | 12 +- examples/nl_fst_compat.py | 57 ++++---- examples/plot_features.py | 202 ----------------------------- examples/plot_somas.py | 18 ++- examples/radius_of_gyration.py | 17 ++- examples/section_ids.py | 12 +- examples/soma_radius_fit.py | 30 ++--- examples/synthesis_json.py | 216 ------------------------------- tests/test_examples.py | 27 ++++ 16 files changed, 245 insertions(+), 547 deletions(-) delete mode 100755 examples/plot_features.py delete mode 100755 examples/synthesis_json.py create mode 100644 tests/test_examples.py diff --git a/examples/boxplot.py b/examples/boxplot.py index 6474bf46..b3bcab3d 100644 --- a/examples/boxplot.py +++ b/examples/boxplot.py @@ -27,11 +27,17 @@ # SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. """Box Plot function for multiple morphs.""" +from pathlib import Path +from neurom import load_morphologies from neurom.view import matplotlib_utils +from neurom.features import get -def boxplot(neurons, feature, new_fig=True, subplot=False): +PACKAGE_DIR = Path(__file__).resolve().parent.parent + + +def boxplot(neurons, feature, new_fig=True, subplot=111): """Plot a histogram of the selected feature for the population of morphologies. Plots x-axis versus y-axis on a scatter|histogram|binned values plot. @@ -52,12 +58,26 @@ def boxplot(neurons, feature, new_fig=True, subplot=False): Default is False, which returns a matplotlib figure object. If True, returns a matplotlib axis object, for use as a subplot. """ - feature_values = [getattr(neu, 'get_' + feature)() for neu in neurons] + feature_values = [get(feature, neuron) for neuron in neurons] _, ax = matplotlib_utils.get_figure(new_fig=new_fig, subplot=subplot) ax.boxplot(feature_values) - x_labels = ['neuron_id' for _ in neurons] + x_labels = [neuron.name for neuron in neurons] ax.set_xticklabels(x_labels) + + # uncomment below to show image + # pylab.show() + + +def main(): + + morphology_directory = Path(PACKAGE_DIR, "tests/data/valid_set") + neurons = load_morphologies(morphology_directory) + boxplot(neurons, "section_lengths") + + +if __name__ == "__main__": + main() diff --git a/examples/density_plot.py b/examples/density_plot.py index 6e008eaa..e7fd4ad2 100644 --- a/examples/density_plot.py +++ b/examples/density_plot.py @@ -27,20 +27,28 @@ # SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. """Example for generating density plots.""" +from pathlib import Path import pylab as plt +import matplotlib as mpl import numpy as np from neurom import get as get_feat from neurom.view import matplotlib_utils, matplotlib_impl from neurom.core.types import NeuriteType +from neurom import load_morphologies + +PACKAGE_DIR = Path(__file__).resolve().parent.parent + def extract_density(population, plane='xy', bins=100, neurite_type=NeuriteType.basal_dendrite): """Extracts the 2d histogram of the center coordinates of segments in the selected plane. """ - segment_midpoints = get_feat('segment_midpoints', population, neurite_type=neurite_type) + segment_midpoints = np.array( + get_feat('segment_midpoints', population, neurite_type=neurite_type) + ) horiz = segment_midpoints[:, 'xyz'.index(plane[0])] vert = segment_midpoints[:, 'xyz'.index(plane[1])] return np.histogram2d(np.array(horiz), np.array(vert), bins=(bins, bins)) @@ -62,12 +70,13 @@ def plot_density(population, # pylint: disable=too-many-arguments, too-many-loc mask = H1 < threshold # mask = H1==0 H2 = np.ma.masked_array(H1, mask) - getattr(plt.cm, color_map).set_bad(color='white', alpha=None) + colormap = mpl.cm.get_cmap(color_map).copy() + colormap.set_bad(color='white', alpha=None) plots = ax.contourf((xedges1[:-1] + xedges1[1:]) / 2, (yedges1[:-1] + yedges1[1:]) / 2, np.transpose(H2), # / np.max(H2), - cmap=getattr(plt.cm, color_map), levels=levels) + cmap=colormap, levels=levels) if not no_colorbar: cbar = plt.colorbar(plots) @@ -91,9 +100,23 @@ def plot_neuron_on_density(population, # pylint: disable=too-many-arguments """ _, ax = matplotlib_utils.get_figure(new_fig=new_fig) - matplotlib_impl.plot_tree(population.neurites[0], ax) + ref_neuron = population[0] + matplotlib_impl.plot_tree(ref_neuron.neurites[0], ax) return plot_density(population, plane=plane, bins=bins, new_fig=False, subplot=subplot, colorlabel=colorlabel, labelfontsize=labelfontsize, levels=levels, color_map=color_map, no_colorbar=no_colorbar, threshold=threshold, neurite_type=neurite_type, **kwargs) + + +def main(): + + morphology_directory = Path(PACKAGE_DIR, "tests/data/valid_set") + neurons = load_morphologies(morphology_directory) + + plot_density(neurons) + plot_neuron_on_density(neurons) + + +if __name__ == "__main__": + main() diff --git a/examples/end_to_end_distance.py b/examples/end_to_end_distance.py index 070a02e0..4e2adeb3 100755 --- a/examples/end_to_end_distance.py +++ b/examples/end_to_end_distance.py @@ -28,6 +28,7 @@ # SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. """Calculate and plot end-to-end distance of neurites.""" +from pathlib import Path import neurom as nm from neurom import morphmath @@ -35,6 +36,9 @@ import matplotlib.pyplot as plt +PACKAGE_DIR = Path(__file__).resolve().parent.parent + + def path_end_to_end_distance(neurite): """Calculate and return end-to-end-distance of a given neurite.""" trunk = neurite.root_node.points[0] @@ -54,7 +58,8 @@ def make_end_to_end_distance_plot(nb_segments, end_to_end_distance, neurite_type plt.title(neurite_type) plt.xlabel('Number of segments') plt.ylabel('End-to-end distance') - plt.show() + # uncomment to show + #plt.show() def calculate_and_plot_end_to_end_distance(neurite): @@ -71,9 +76,9 @@ def _dist(seg): end_to_end_distance, neurite.type) -if __name__ == '__main__': +def main(): # load a neuron from an SWC file - filename = 'tests/data/swc/Neuron_3_random_walker_branches.swc' + filename = Path(PACKAGE_DIR, 'tests/data/swc/Neuron_3_random_walker_branches.swc') m = nm.load_morphology(filename) # print mean end-to-end distance per neurite type @@ -92,3 +97,7 @@ def _dist(seg): # print (number of segments, end-to-end distance, neurite type) print(sum(len(s.points) - 1 for s in nrte.root_node.ipreorder()), path_end_to_end_distance(nrte), nrte.type) + + +if __name__ == '__main__': + main() diff --git a/examples/extract_distribution.py b/examples/extract_distribution.py index 128957bc..1288e5e5 100755 --- a/examples/extract_distribution.py +++ b/examples/extract_distribution.py @@ -31,6 +31,7 @@ """Extract a distribution for the selected feature of the population of morphologies among the exponential, normal and uniform distribution, according to the minimum ks distance. """ +from pathlib import Path from itertools import chain import argparse @@ -38,31 +39,19 @@ import neurom as nm from neurom import stats +from neurom.utils import NeuromJSON -def parse_args(): - """Parse command line arguments.""" - parser = argparse.ArgumentParser( - description='Morphology fit distribution extractor', - epilog='Note: Outputs json of the optimal distribution \ - and corresponding parameters.') +PACKAGE_DIR = Path(__file__).resolve().parent.parent - parser.add_argument('datapath', - help='Path to morphology data directory') - parser.add_argument('feature', - help='Feature to be extracted with neurom.get') - - return parser.parse_args() - - -def extract_data(data_path, feature): +def find_optimal_distribution(population_directory, feature): """Loads a list of morphologies, extracts feature and transforms the fitted distribution in the correct format. Returns the optimal distribution, corresponding parameters, minimun and maximum values. """ - population = nm.load_morphologies(data_path) + population = nm.load_morphologies(population_directory) feature_data = [nm.get(feature, n) for n in population] feature_data = list(chain(*feature_data)) @@ -70,13 +59,18 @@ def extract_data(data_path, feature): return stats.optimal_distribution(feature_data) -if __name__ == '__main__': - args = parse_args() +def main(): + + population_directory = Path(PACKAGE_DIR, "tests/data/valid_set") - d_path = args.datapath + result = stats.fit_results_to_dict( + find_optimal_distribution(population_directory, "section_lengths") + ) - feat = args.feature + print(json.dumps( + result, indent=2, separators=(',', ': '), cls=NeuromJSON + )) - _result = stats.fit_results_to_dict(extract_data(d_path, feat)) - print(json.dumps(_result, indent=2, separators=(',', ': '))) +if __name__ == '__main__': + main() diff --git a/examples/features_graph_table.py b/examples/features_graph_table.py index 95a98ea4..8d934a31 100755 --- a/examples/features_graph_table.py +++ b/examples/features_graph_table.py @@ -28,7 +28,6 @@ # SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. """Example for comparison of the same feature of multiple cells.""" -import argparse from pathlib import Path import pylab as pl @@ -36,25 +35,7 @@ from neurom.io.utils import get_morph_files -def parse_args(): - """Parse command line arguments.""" - parser = argparse.ArgumentParser(description='Feature Comparison Between Different Cells') - - parser.add_argument('-d', - '--datapath', - help='Data directory') - - parser.add_argument('-o', - '--odir', - default='.', - help='Output path') - - parser.add_argument('-f', - '--features', - nargs='+', - help='List features separated by spaces') - - return parser.parse_args() +PACKAGE_DIR = Path(__file__).resolve().parent.parent def stylize(ax, name, feature): @@ -88,7 +69,7 @@ def histogram(neuron, feature, ax, bins=15, normed=True, cumulative=False): feature_values = nm.get(feature, neuron) # generate histogram - ax.hist(feature_values, bins=bins, cumulative=cumulative, normed=normed) + ax.hist(feature_values, bins=bins, cumulative=cumulative, density=normed) def plot_feature(feature, cell): @@ -106,14 +87,25 @@ def plot_feature(feature, cell): return fig -if __name__ == '__main__': - args = parse_args() +def create_feature_plots(morphologies_dir, feature_list, output_dir): - for morph_file in get_morph_files(args.datapath): + for morph_file in get_morph_files(morphologies_dir): m = nm.load_morphology(morph_file) - for _feature in args.features: - f = plot_feature(_feature, m) - figname = "{0}_{1}.eps".format(_feature, m.name) - f.savefig(Path(args.odir, figname)) + for feature_name in feature_list: + f = plot_feature(feature_name, m) + figname = f"{feature_name}_{m.name}.eps" + f.savefig(Path(output_dir, figname)) pl.close(f) + + +def main(): + create_feature_plots( + morphologies_dir=Path(PACKAGE_DIR, "tests/data/valid_set"), + feature_list=["section_lengths"], + output_dir=".", + ) + + +if __name__ == '__main__': + main() diff --git a/examples/get_features.py b/examples/get_features.py index 1d223884..1fee0f16 100755 --- a/examples/get_features.py +++ b/examples/get_features.py @@ -33,13 +33,16 @@ morphometrics functionality. """ +from pathlib import Path -from __future__ import print_function from pprint import pprint import numpy as np import neurom as nm +PACKAGE_DIR = Path(__file__).resolve().parent.parent + + def stats(data): """Dictionary with summary stats for data @@ -60,9 +63,9 @@ def pprint_stats(data): pprint(stats(data)) -if __name__ == '__main__': +def main(): - filename = 'tests/data/swc/Neuron.swc' + filename = Path(PACKAGE_DIR, 'tests/data/swc/Neuron.swc') # load a neuron from an SWC file m = nm.load_morphology(filename) @@ -152,3 +155,7 @@ def pprint_stats(data): rem_bifangles = nm.get('remote_bifurcation_angles', m, neurite_type=ttype) print('Local bifurcation angles (', ttype, '):', sep='') pprint_stats(rem_bifangles) + + +if __name__ == '__main__': + main() diff --git a/examples/histogram.py b/examples/histogram.py index bd336637..0d55670b 100644 --- a/examples/histogram.py +++ b/examples/histogram.py @@ -27,13 +27,20 @@ # SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. """Simple Histogram function for multiple morphologies.""" +from pathlib import Path + from itertools import chain import numpy as np +import neurom.features from neurom.view import matplotlib_utils +from neurom import load_morphologies + +PACKAGE_DIR = Path(__file__).resolve().parent.parent -def histogram(neurons, feature, new_fig=True, subplot=False, normed=False, **kwargs): + +def histogram(neurons, feature, new_fig=True, subplot=111, normed=False, **kwargs): """ Plot a histogram of the selected feature for the population of morphologies. Plots x-axis versus y-axis on a scatter|histogram|binned values plot. @@ -80,11 +87,11 @@ def histogram(neurons, feature, new_fig=True, subplot=False, normed=False, **kwa kwargs['title'] = kwargs.get('title', feature + ' histogram') - feature_values = [getattr(neu, 'get_' + feature)() for neu in neurons] + feature_values = [neurom.features.get(feature, neu) for neu in neurons] neu_labels = [neu.name for neu in neurons] - ax.hist(feature_values, bins=bins, cumulative=cumulative, label=neu_labels, normed=normed) + ax.hist(feature_values, bins=bins, cumulative=cumulative, label=neu_labels, density=normed) kwargs['no_legend'] = len(neu_labels) == 1 @@ -98,7 +105,7 @@ def population_feature_values(pops, feature): for pop in pops: - feature_values = [getattr(neu, 'get_' + feature)() for neu in pop.morphologies] + feature_values = [neurom.features.get(feature, neu) for neu in pop] # ugly hack to chain in case of list of lists if any([isinstance(p, (list, np.ndarray)) for p in feature_values]): @@ -110,7 +117,7 @@ def population_feature_values(pops, feature): return pops_feature_values -def population_histogram(pops, feature, new_fig=True, normed=False, subplot=False, **kwargs): +def population_histogram(pops, feature, new_fig=True, normed=False, subplot=111, **kwargs): """ Plot a histogram of the selected feature for the population of morphologies. Plots x-axis versus y-axis on a scatter|histogram|binned values plot. @@ -160,8 +167,19 @@ def population_histogram(pops, feature, new_fig=True, normed=False, subplot=Fals pops_labels = [pop.name for pop in pops] - ax.hist(pops_feature_values, bins=bins, cumulative=cumulative, label=pops_labels, normed=normed) + ax.hist(pops_feature_values, bins=bins, cumulative=cumulative, label=pops_labels, density=normed) kwargs['no_legend'] = len(pops_labels) == 1 return matplotlib_utils.plot_style(fig=fig, ax=ax, **kwargs) + + +def main(): + + pop1 = load_morphologies(Path(PACKAGE_DIR, "tests/data/valid_set")) + pop2 = load_morphologies(Path(PACKAGE_DIR, "tests/data/valid_set")) + population_histogram([pop1, pop2], "section_lengths") + + +if __name__ == "__main__": + main() diff --git a/examples/iteration_analysis.py b/examples/iteration_analysis.py index a3a712e9..ffc6a6b1 100755 --- a/examples/iteration_analysis.py +++ b/examples/iteration_analysis.py @@ -33,8 +33,8 @@ morphometrics functionality using iterators. """ +from pathlib import Path -from __future__ import print_function from neurom.core.dataformat import COLS import neurom as nm from neurom import geom @@ -44,10 +44,12 @@ from neurom import morphmath as mm import numpy as np +PACKAGE_DIR = Path(__file__).resolve().parent.parent -if __name__ == '__main__': - filename = 'tests/data/swc/Neuron.swc' +def main(): + + filename = Path(PACKAGE_DIR, 'tests/data/swc/Neuron.swc') # load a neuron from an SWC file m = nm.load_morphology(filename) @@ -142,3 +144,7 @@ def n_points(sec): # Morphology's bounding box # Note: does not account for soma radius print('Bounding box ((min x, y, z), (max x, y, z))', geom.bounding_box(m)) + + +if __name__ == '__main__': + main() diff --git a/examples/nl_fst_compat.py b/examples/nl_fst_compat.py index b1607714..38671fc2 100755 --- a/examples/nl_fst_compat.py +++ b/examples/nl_fst_compat.py @@ -28,32 +28,39 @@ # SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. """Compatibility between NL and H5 files.""" # pylint: disable=protected-access +from pathlib import Path import numpy as np import neurom as nm -from neurom.features import neurite as _nf - -m_h5 = nm.load_morphology('tests/data/h5/v1/bio_neuron-001.h5') -m_asc = nm.load_morphology('tests/data/neurolucida/bio_neuron-001.asc') - -print('h5 number of sections: %s' % nm.get('number_of_sections', m_h5)[0]) -print('nl number of sections: %s\n' % nm.get('number_of_sections', m_asc)[0]) -print('h5 number of segments: %s' % nm.get('number_of_segments', m_h5)[0]) -print('nl number of segments: %s\n' % nm.get('number_of_segments', m_asc)[0]) -print('h5 total neurite length: %s' % - np.sum(nm.get('section_lengths', m_h5))) -print('nl total neurite length: %s\n' % - np.sum(nm.get('section_lengths', m_asc))) -print('h5 principal direction extents: %s' % - nm.get('principal_direction_extents', m_h5)) -print('nl principal direction extents: %s' % - nm.get('principal_direction_extents', m_asc)) - -print('\nNumber of neurites:') -for nt in iter(nm.NeuriteType): - print(nt, _nf.n_neurites(m_h5, neurite_type=nt), _nf.n_neurites(m_asc, neurite_type=nt)) - -print('\nNumber of segments:') -for nt in iter(nm.NeuriteType): - print(nt, _nf.n_segments(m_h5, neurite_type=nt), _nf.n_segments(m_asc, neurite_type=nt)) +from neurom.features import neurite as nf +from neurom.features import morphology as mf + +PACKAGE_DIR = Path(__file__).resolve().parent.parent + + +def main(): + + m_h5 = nm.load_morphology(Path(PACKAGE_DIR, 'tests/data/h5/v1/bio_neuron-001.h5')) + m_asc = nm.load_morphology(Path(PACKAGE_DIR, 'tests/data/neurolucida/bio_neuron-001.asc')) + + print('h5 number of sections:', nm.get('number_of_sections', m_h5)) + print('nl number of sections:', nm.get('number_of_sections', m_asc)) + print('h5 number of segments:', nm.get('number_of_segments', m_h5)) + print('nl number of segments:', nm.get('number_of_segments', m_asc)) + print('h5 total neurite length:', np.sum(nm.get('section_lengths', m_h5))) + print('nl total neurite length:', np.sum(nm.get('section_lengths', m_asc))) + print('h5 principal direction extents:', nm.get('principal_direction_extents', m_h5)) + print('nl principal direction extents:', nm.get('principal_direction_extents', m_asc)) + + print('\nNumber of neurites:') + for nt in iter(nm.NeuriteType): + print(nt, mf.number_of_neurites(m_h5, neurite_type=nt), mf.number_of_neurites(m_asc, neurite_type=nt)) + + print('\nNumber of segments:') + for nt in iter(nm.NeuriteType): + print(nt, nf.number_of_segments(m_h5.neurites[0]), nf.number_of_segments(m_asc.neurites[0])) + + +if __name__ == "__main__": + main() diff --git a/examples/plot_features.py b/examples/plot_features.py deleted file mode 100755 index b5bb0f2c..00000000 --- a/examples/plot_features.py +++ /dev/null @@ -1,202 +0,0 @@ -#!/usr/bin/env python -# Copyright (c) 2015, Ecole Polytechnique Federale de Lausanne, Blue Brain Project -# All rights reserved. -# -# This file is part of NeuroM -# -# Redistribution and use in source and binary forms, with or without -# modification, are permitted provided that the following conditions are met: -# -# 1. Redistributions of source code must retain the above copyright -# notice, this list of conditions and the following disclaimer. -# 2. Redistributions in binary form must reproduce the above copyright -# notice, this list of conditions and the following disclaimer in the -# documentation and/or other materials provided with the distribution. -# 3. Neither the name of the copyright holder nor the names of -# its contributors may be used to endorse or promote products -# derived from this software without specific prior written permission. -# -# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND -# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED -# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE -# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY -# DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES -# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; -# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND -# ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT -# (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS -# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -"""Plot a selection of features from a morphology population.""" - -from collections import defaultdict -from collections import namedtuple -import sys -import json -import argparse - -import numpy as np -import neurom as nm -from neurom.view import matplotlib_utils -import scipy.stats as _st -from matplotlib.backends.backend_pdf import PdfPages - - -DISTS = { - 'normal': lambda p, bins: _st.norm.pdf(bins, p['mu'], p['sigma']), - 'uniform': lambda p, bins: _st.uniform.pdf(bins, p['min'], p['max'] - p['min']), - 'constant': lambda p, bins: None -} - - -def bin_centers(bin_edges): - """Return array of bin centers given an array of bin edges""" - return (bin_edges[1:] + bin_edges[:-1]) / 2.0 - - -def bin_widths(bin_edges): - """Return array of bin widths given an array of bin edges""" - return bin_edges[1:] - bin_edges[:-1] - - -def histo_entries(histo): - """Calculate the number of entries in a histogram - - This is the sum of bin height * bin width - """ - bw = bin_widths(histo[1]) - return np.sum(histo[0] * bw) - - -def dist_points(bin_edges, d): - """Return an array of values according to a distribution - - Points are calculated at the center of each bin - """ - bc = bin_centers(bin_edges) - if d is not None: - d = DISTS[d['type']](d, bc) - return d, bc - - -def calc_limits(data, dist=None, padding=0.25): - """Calculate a suitable range for a histogram - - Returns: - tuple of (min, max) - """ - dmin = sys.float_info.max if dist is None else dist.get('min', - sys.float_info.max) - dmax = sys.float_info.min if dist is None else dist.get('max', - sys.float_info.min) - _min = min(min(data), dmin) - _max = max(max(data), dmax) - - padding = padding * (_max - _min) - return _min - padding, _max + padding - - -# Neurite types of interest -NEURITES_ = (nm.NeuriteType.axon, - nm.NeuriteType.apical_dendrite, - nm.NeuriteType.basal_dendrite,) - - -# Features of interest -FEATURES = ('segment_lengths', - 'section_lengths', - 'section_path_distances', - 'section_radial_distances', - 'trunk_origin_radii') - - -def load_neurite_features(filepath): - """Unpack relevant data into megadict.""" - stuff = defaultdict(lambda: defaultdict(list)) - morphs = nm.load_morphologies(filepath) - # unpack data into arrays - for m in morphs: - for t in NEURITES_: - for feat in FEATURES: - stuff[feat][str(t).split('.')[1]].extend( - nm.get(feat, m, neurite_type=t) - ) - return stuff - - -Plot = namedtuple('Plot', 'fig, ax') - - -def parse_args(): - """Parse command line arguments.""" - parser = argparse.ArgumentParser( - description='Morphology feature plotter', - epilog='Note: Makes plots of various features and superimposes\ - input distributions. Plots are saved to PDF file.', - formatter_class=argparse.ArgumentDefaultsHelpFormatter) - - parser.add_argument('datapath', - help='Morphology data directory path') - - parser.add_argument('--mtypeconfig', - required=True, - help='Get mtype JSON configuration file') - - parser.add_argument('--output', - default='plots.pdf', - help='Output PDF file name') - return parser.parse_args() - - -def main(data_dir, mtype_file): # pylint: disable=too-many-locals - """Run the stuff.""" - # data structure to store results - stuff = load_neurite_features(data_dir) - sim_params = json.load(open(mtype_file)) - - # load histograms, distribution parameter sets and figures into arrays. - # To plot figures, do - # plots[i].fig.show() - # To modify an axis, do - # plots[i].ax.something() - - _plots = [] - - for feat, d in stuff.items(): - for typ, data in d.items(): - dist = sim_params['components'][typ].get(feat, None) - print('Type = %s, Feature = %s, Distribution = %s' % (typ, feat, dist)) - # if no data available, skip this feature - if not data: - print("No data found for feature %s (%s)" % (feat, typ)) - continue - # print 'DATA', data - num_bins = 100 - limits = calc_limits(data, dist) - bin_edges = np.linspace(limits[0], limits[1], num_bins + 1) - histo = np.histogram(data, bin_edges, normed=True) - print('PLOT LIMITS:', limits) - # print 'DATA:', data - # print 'BIN HEIGHT', histo[0] - plot = Plot(*matplotlib_utils.get_figure(new_fig=True, subplot=111)) - plot.ax.set_xlim(*limits) - plot.ax.bar(histo[1][:-1], histo[0], width=bin_widths(histo[1])) - dp, bc = dist_points(histo[1], dist) - # print 'BIN CENTERS:', bc, len(bc) - if dp is not None: - # print 'DIST POINTS:', dp, len(dp) - plot.ax.plot(bc, dp, 'r*') - plot.ax.set_title('%s (%s)' % (feat, typ)) - _plots.append(plot) - - return _plots - - -if __name__ == '__main__': - args = parse_args() - print('MTYPE FILE:', args.mtypeconfig) - plots = main(args.datapath, args.mtypeconfig) - - pp = PdfPages(args.output) - for p in plots: - pp.savefig(p.fig) - pp.close() diff --git a/examples/plot_somas.py b/examples/plot_somas.py index 0c7838ff..46457888 100755 --- a/examples/plot_somas.py +++ b/examples/plot_somas.py @@ -36,30 +36,34 @@ import matplotlib.pyplot as plt import numpy as np -DATA_PATH = Path(__file__).parent.parent / 'tests/data' -SWC_PATH = Path(DATA_PATH, 'swc') +DATA_PATH = Path(__file__).resolve().parent.parent / 'tests/data/swc' def random_color(): """Random color generation.""" - return np.random.rand(3, 1) + return np.random.rand(4) def plot_somas(somas): """Plot set of somas on same figure as spheres, each with different color.""" _, ax = matplotlib_utils.get_figure(new_fig=True, subplot=111, - params={'projection': '3d', 'aspect': 'equal'}) + params={'projection': '3d', 'aspect': 'auto'}) for s in somas: matplotlib_utils.plot_sphere(ax, s.center, s.radius, color=random_color(), alpha=1) - plt.show() + # uncomment to show + # plt.show() -if __name__ == '__main__': + +def main(): # define set of files containing relevant morphs - file_nms = [Path(SWC_PATH, file_nm) for file_nm in ['Soma_origin.swc', + file_nms = [Path(DATA_PATH, file_nm) for file_nm in ['Soma_origin.swc', 'Soma_translated_1.swc', 'Soma_translated_2.swc']] # load from file and plot sms = [load_morphology(file_nm).soma for file_nm in file_nms] plot_somas(sms) + +if __name__ == '__main__': + main() diff --git a/examples/radius_of_gyration.py b/examples/radius_of_gyration.py index 991dd967..0bcd4306 100755 --- a/examples/radius_of_gyration.py +++ b/examples/radius_of_gyration.py @@ -28,6 +28,7 @@ # SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. """Calculate radius of gyration of neurites.""" +from pathlib import Path import neurom as nm from neurom import morphmath as mm @@ -35,6 +36,9 @@ import numpy as np +PACKAGE_DIR = Path(__file__).resolve().parent.parent + + def segment_centre_of_mass(seg): """Calculate and return centre of mass of a segment. @@ -53,8 +57,8 @@ def neurite_centre_of_mass(neurite): centre_of_mass = np.zeros(3) total_volume = 0 - seg_vol = np.array(map(mm.segment_volume, nm.iter_segments(neurite))) - seg_centre_of_mass = np.array(map(segment_centre_of_mass, nm.iter_segments(neurite))) + seg_vol = np.array(list(map(mm.segment_volume, nm.iter_segments(neurite)))) + seg_centre_of_mass = np.array(list(map(segment_centre_of_mass, nm.iter_segments(neurite)))) # multiply array of scalars with array of arrays # http://stackoverflow.com/questions/5795700/multiply-numpy-array-of-scalars-by-array-of-vectors @@ -87,9 +91,10 @@ def mean_rad_of_gyration(neurites): return np.mean([radius_of_gyration(n) for n in neurites]) -if __name__ == '__main__': +def main(): + # load a neuron from an SWC file - filename = 'tests/data/swc/Neuron.swc' + filename = Path(PACKAGE_DIR, 'tests/data/swc/Neuron.swc') m = nm.load_morphology(filename) # for every neurite, print (number of segments, radius of gyration, neurite type) @@ -104,3 +109,7 @@ def mean_rad_of_gyration(neurites): print('Mean radius of gyration for apical dendrites: ', mean_rad_of_gyration(n for n in m.neurites if n.type == nm.APICAL_DENDRITE)) + + +if __name__ == '__main__': + main() diff --git a/examples/section_ids.py b/examples/section_ids.py index 754fac32..00cac5c8 100755 --- a/examples/section_ids.py +++ b/examples/section_ids.py @@ -28,12 +28,16 @@ # SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. """Get sections and segments by ID.""" +from pathlib import Path import neurom as nm from neurom import morphmath as mm from neurom.core.dataformat import COLS +PACKAGE_DIR = Path(__file__).resolve().parent.parent + + def get_segment(neuron, section_id, segment_id): """Get a segment given a section and segment id @@ -44,11 +48,15 @@ def get_segment(neuron, section_id, segment_id): return sec.points[segment_id:segment_id + 2][:, COLS.XYZR] -if __name__ == '__main__': +def main(): - m = nm.load_morphology('tests/data/h5/v1/Neuron.h5') + m = nm.load_morphology(Path(PACKAGE_DIR, 'tests/data/h5/v1/Neuron.h5')) seg = get_segment(m, 3, 2) print('Segment:\n', seg) print('Mid-point (x, y, z):\n', mm.linear_interpolate(seg[0], seg[1], 0.5)) print('Mid-point R:\n', mm.interpolate_radius(seg[0][COLS.R], seg[1][COLS.R], 0.5)) + + +if __name__ == '__main__': + main() diff --git a/examples/soma_radius_fit.py b/examples/soma_radius_fit.py index 3861cf08..91e06db3 100755 --- a/examples/soma_radius_fit.py +++ b/examples/soma_radius_fit.py @@ -31,23 +31,13 @@ """Extract a distribution for the soma radii of the population (list) of morphologies. for the soma radii of the population (list) of morphologies. """ - -import argparse +from pathlib import Path import neurom as nm from neurom import stats as st -def parse_args(): - """Parse command line arguments.""" - parser = argparse.ArgumentParser( - description='Morphology fit distribution extractor', - epilog='Note: Prints the optimal distribution and corresponding parameters.') - - parser.add_argument('datapath', - help='Path to morphology data file or directory') - - return parser.parse_args() +PACKAGE_DIR = Path(__file__).resolve().parent.parent def test_multiple_distr(filepath): @@ -61,17 +51,19 @@ def test_multiple_distr(filepath): distr_to_check = ('norm', 'expon', 'uniform') # Get the soma radii of a population of morphs - soma_size = nm.get('soma_radii', population) + soma_size = nm.get('soma_radius', population) # Find the best fit distribution return st.optimal_distribution(soma_size, distr_to_check) -if __name__ == '__main__': - args = parse_args() +def main(): - data_path = args.datapath + morphology_path = Path(PACKAGE_DIR, "tests/data/swc/Neuron.swc") - result = test_multiple_distr(data_path) - print("Optimal distribution fit for soma radius is: %s with parameters %s" % - (result.type, result.params)) + result = test_multiple_distr(morphology_path) + print(f"Optimal distribution fit for soma radius is: {result.type} with parameters {result.params}") + + +if __name__ == '__main__': + main() diff --git a/examples/synthesis_json.py b/examples/synthesis_json.py deleted file mode 100755 index 4d67feb0..00000000 --- a/examples/synthesis_json.py +++ /dev/null @@ -1,216 +0,0 @@ -#!/usr/bin/env python - -# Copyright (c) 2015, Ecole Polytechnique Federale de Lausanne, Blue Brain Project -# All rights reserved. -# -# This file is part of NeuroM -# -# Redistribution and use in source and binary forms, with or without -# modification, are permitted provided that the following conditions are met: -# -# 1. Redistributions of source code must retain the above copyright -# notice, this list of conditions and the following disclaimer. -# 2. Redistributions in binary form must reproduce the above copyright -# notice, this list of conditions and the following disclaimer in the -# documentation and/or other materials provided with the distribution. -# 3. Neither the name of the copyright holder nor the names of -# its contributors may be used to endorse or promote products -# derived from this software without specific prior written permission. -# -# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND -# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED -# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE -# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY -# DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES -# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; -# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND -# ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT -# (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS -# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. - -"""Extract the optimal distributions for the following features of the population of morphologies: - soma: radius - basal dendrites: n_neurites - apical dendrites: n_neurites - axons: n_neurites -""" - -import argparse -from collections import OrderedDict -from collections import defaultdict -from collections import namedtuple -from itertools import chain -import os -import json -from json import encoder - -from neurom import get, load_morphologies, NeuriteType, stats -from neurom.io.utils import get_morph_files - -encoder.FLOAT_REPR = lambda o: format(o, '.2f') - -Limits = namedtuple('Limits', 'min, max') -Feature = namedtuple('Feature', 'name, limits') -Component = namedtuple('Component', 'name, features') - - -# Map feature names to data getter namer -FEATURE_MAP = { - 'radius': 'soma_radii', - 'number': 'number_of_neurites', - 'segment_length': 'segment_lengths', - 'initial_radius': 'trunk_origin_radii', - # 'taper_rate': lambda n, kwargs: np.array(list(n.iter_segments(mm.segment_taper_rate, - # **kwargs))), -} - -# Map component name to filtering parameters where applicable -PARAM_MAP = { - 'basal_dendrite': {'neurite_type': NeuriteType.basal_dendrite}, - 'apical_dendrite': {'neurite_type': NeuriteType.apical_dendrite}, - 'axon': {'neurite_type': NeuriteType.axon}, - 'soma': None -} - - -def extract_data(neurons, feature, params=None): - """Extracts feature from a list of morphologies - and transforms the fitted distribution in the correct format. - Returns the optimal distribution and corresponding parameters. - Normal distribution params (mean, std) - Exponential distribution params (loc, scale) - Uniform distribution params (min, range) - """ - if params is None: - params = {} - - feature_data = [get(FEATURE_MAP[feature], n, **params) for n in neurons] - feature_data = list(chain(*feature_data)) - return stats.optimal_distribution(feature_data) - - -def transform_header(mtype_name): - """Add header to json output to wrap around distribution data. - """ - head_dict = OrderedDict() - - head_dict["m-type"] = mtype_name - head_dict["components"] = defaultdict(OrderedDict) - - return head_dict - - -def transform_package(mtype, files, components): - """Put together header and list of data into one json output. - feature_list contains all the information about the data to be extracted: - features, feature_names, feature_components, feature_min, feature_max - """ - data_dict = transform_header(mtype) - neurons = load_morphologies(files) - - for comp in components: - params = PARAM_MAP[comp.name] - for feature in comp.features: - result = stats.fit_results_to_dict( - extract_data(neurons, feature.name, params), - feature.limits.min, feature.limits.max - ) - - # When the distribution is normal with sigma = 0 - # it will be replaced with constant - if result['type'] == 'normal' and result['sigma'] == 0.0: - replace_result = OrderedDict((('type', 'constant'), - ('val', result['mu']))) - result = replace_result - - data_dict["components"][comp.name].update({feature.name: result}) - - return data_dict - - -def get_mtype_from_filename(filename, sep='_'): - """Get mtype of a morphology file from file name - - Assumes file name has structure 'a/b/c/d/mtype_xyx.abc' - """ - return os.path.basename(filename).split(sep)[0] - - -def get_mtype_from_directory(filename): - """Get mtype of a morphology file from file's parent directory name - - Assumes file name has structure 'a/b/c/mtype/xyx.abc' - """ - return os.path.split(os.path.dirname(filename))[-1] - - -MTYPE_GETTERS = { - 'filename': get_mtype_from_filename, - 'directory': get_mtype_from_directory -} - - -def parse_args(): - """Parse command line arguments.""" - parser = argparse.ArgumentParser( - description='Morphology fit distribution extractor', - epilog='Note: Outputs json of the optimal distribution \ - and corresponding parameters.') - - parser.add_argument('datapaths', - nargs='+', - help='Morphology data directory paths') - - parser.add_argument('--mtype', choices=MTYPE_GETTERS.keys(), - help='Get mtype from filename or parent directory') - - return parser.parse_args() - - -# Note: new features require a data getter function in -# FEATURE_MAP -# TODO: Set this config from a JSON or YAML file -COMPONENTS = [ - Component('soma', - [Feature('radius', Limits(None, None))]), - Component('basal_dendrite', - [ - Feature('number', Limits(0, None)), - Feature('segment_length', Limits(0, None)), - Feature('initial_radius', Limits(0, None)), - # Feature('taper_rate', Limits(0, None)) - ]), - Component('apical_dendrite', - [ - Feature('number', Limits(0, None)), - Feature('segment_length', Limits(0, None)), - Feature('initial_radius', Limits(0, None)), - # Feature('taper_rate', Limits(0, None)) - ]), - Component('axon', - [ - Feature('number', Limits(0, None)), - Feature('segment_length', Limits(0, None)), - Feature('initial_radius', Limits(0, None)), - # Feature('taper_rate', Limits(0, None)) - ]), -] - - -if __name__ == '__main__': - args = parse_args() - - data_dirs = args.datapaths - - mtype_getter = MTYPE_GETTERS.get(args.mtype, lambda f: 'UNKNOWN') - - for d in data_dirs: - mtype_files = defaultdict(list) - for f in get_morph_files(d): - mtype_files[mtype_getter(f)].append(f) - - _results = [transform_package(mtype_, files_, COMPONENTS) - for mtype_, files_ in mtype_files.items()] - - for res in _results: - print(json.dumps(res, indent=2), '\n') diff --git a/tests/test_examples.py b/tests/test_examples.py new file mode 100644 index 00000000..0b3c2693 --- /dev/null +++ b/tests/test_examples.py @@ -0,0 +1,27 @@ +import os +import tempfile +import importlib.util +from pathlib import Path + +import pytest + +TESTS_DIR = Path(__file__).resolve().parent + +EXAMPLES_DIR = TESTS_DIR.parent / "examples" +print(EXAMPLES_DIR) + +@pytest.mark.parametrize("filepath", EXAMPLES_DIR.glob("*.py")) +def test_example(filepath): + + spec = importlib.util.spec_from_file_location(filepath.stem, filepath) + module = spec.loader.load_module() + + with tempfile.TemporaryDirectory() as tempdir: + + # change directory to avoid creating files in the root folder + try: + cwd = os.getcwd() + os.chdir(tempdir) + module.main() + finally: + os.chdir(cwd)