From f86b8950c130d4cb7d5b6a14fc2e85f830de930d Mon Sep 17 00:00:00 2001 From: Alexis Arnaudon Date: Wed, 25 Sep 2024 15:03:06 +0200 Subject: [PATCH] Minor updates (#51) --- emodel_generalisation/adaptation.py | 29 +++++--- .../bluecellulab_evaluator.py | 72 +++++-------------- emodel_generalisation/cli.py | 51 ++++++++----- emodel_generalisation/mcmc.py | 22 ++++-- emodel_generalisation/model/bpopt.py | 17 ++--- emodel_generalisation/model/ecodes.py | 2 +- emodel_generalisation/model/evaluation.py | 4 +- emodel_generalisation/utils.py | 69 ++++++++++++++---- tests/data/adapt_df.csv | 6 +- tests/data/adapted_evaluation_df.csv | 10 +-- tests/data/evaluation_df.csv | 10 +-- tests/test_cli.py | 10 +-- tests/test_evaluation.py | 8 +-- 13 files changed, 178 insertions(+), 132 deletions(-) diff --git a/emodel_generalisation/adaptation.py b/emodel_generalisation/adaptation.py index 0e8f179..301d6dd 100644 --- a/emodel_generalisation/adaptation.py +++ b/emodel_generalisation/adaptation.py @@ -24,7 +24,6 @@ from functools import partial from pathlib import Path -import matplotlib import matplotlib.pyplot as plt import numpy as np import pandas as pd @@ -44,8 +43,7 @@ from emodel_generalisation.model.modifiers import remove_soma from emodel_generalisation.utils import FEATURE_FILTER from emodel_generalisation.utils import get_scores - -matplotlib.use("Agg") +from emodel_generalisation.utils import isolate def get_scales(scales_params, with_unity=False): @@ -85,7 +83,6 @@ def build_resistance_models( # it seems dask does not quite work on this (to investigate, but multiprocessing is fast enough) df = func(df, access_point, parallel_factory="multiprocessing") - models = {} for emodel in emodels: _df = df[df.emodel == emodel] @@ -94,13 +91,16 @@ def build_resistance_models( if len(rin[rin < 0]) == 0: try: coeffs, extra = Polynomial.fit(np.log10(scaler), np.log10(rin), 3, full=True) - if extra[0] < rcond_min: - models[emodel] = { - "resistance": {"polyfit_params": coeffs.convert().coef.tolist()}, - "shape": exemplar_data[key], - } + if extra[0] > rcond_min: + print(f"resistance fit for {key} of {emodel} is not so good") + models[emodel] = { + "resistance": {"polyfit_params": coeffs.convert().coef.tolist()}, + "shape": exemplar_data[key], + } except (np.linalg.LinAlgError, TypeError): print(f"fail to fit emodel {emodel}") + else: + print(f"resistance fit for {key} of {emodel} has negative rin") return df[df.emodel.isin(models)], models @@ -145,7 +145,7 @@ def _adapt_combo(combo, models, rhos, key="soma", min_scale=0.01, max_scale=10): return {f"{key}_scaler": np.clip(scale, min_scale, max_scale)} -def _adapt_single_soma_ais( +def __adapt_single_soma_ais( combo, access_point=None, models=None, @@ -199,6 +199,15 @@ def _adapt_single_soma_ais( return {k: combo[k] for k in ["soma_scaler", "ais_scaler"]} +def _adapt_single_soma_ais(*args, **kwargs): + timeout = kwargs.pop("timeout", 30 * 60) + res = isolate(__adapt_single_soma_ais, timeout=timeout)(*args, **kwargs) + if res is None: + print("timeout", args, kwargs) + return {k: 1.0 for k in ["soma_scaler", "ais_scaler"]} + return res + + def make_evaluation_df(combos_df, emodels, exemplar_data, rhos=None): """Make a df to be evaluated.""" df = pd.DataFrame() diff --git a/emodel_generalisation/bluecellulab_evaluator.py b/emodel_generalisation/bluecellulab_evaluator.py index 6f2e8c3..9884765 100644 --- a/emodel_generalisation/bluecellulab_evaluator.py +++ b/emodel_generalisation/bluecellulab_evaluator.py @@ -3,60 +3,19 @@ import logging import os from copy import copy -from multiprocessing.context import TimeoutError # pylint: disable=redefined-builtin from pathlib import Path import bluecellulab import efel from bluecellulab.simulation.neuron_globals import NeuronGlobals from bluepyparallel.evaluator import evaluate -from bluepyparallel.parallel import NestedPool + +from emodel_generalisation.utils import isolate logger = logging.getLogger(__name__) AXON_LOC = "self.axonal[1](0.5)._ref_v" -def isolate(func, timeout=None): - """Isolate a generic function for independent NEURON instances. - - It must be used in conjunction with NestedPool. - - Example: - - .. code-block:: python - - def _to_be_isolated(morphology_path, point): - cell = nrnhines.get_NRN_cell(morphology_path) - return nrnhines.point_to_section_end(cell.icell.all, point) - - def _isolated(morph_data): - return nrnhines.isolate(_to_be_isolated)(*morph_data) - - with nrnhines.NestedPool(processes=n_workers) as pool: - result = pool.imap_unordered(_isolated, data) - - - Args: - func (function): function to isolate - - Returns: - the isolated function - - Note: it does not work as decorator. - """ - - def func_isolated(*args, **kwargs): - with NestedPool(1, maxtasksperchild=1) as pool: - res = pool.apply_async(func, args, kwargs) - try: - out = res.get(timeout=timeout) - except TimeoutError: # pragma: no cover - out = None - return out - - return func_isolated - - def calculate_threshold_current(cell, config, holding_current): """Calculate threshold current""" min_current_spike_count = run_spike_sim( @@ -84,7 +43,7 @@ def calculate_threshold_current(cell, config, holding_current): if max_current_spike_count < 1: logger.debug("Cell is not firing at max current, we multiply by 2") config["min_threshold_current"] = copy(config["max_threshold_current"]) - config["max_threshold_current"] *= 2.0 + config["max_threshold_current"] *= 1.2 return calculate_threshold_current(cell, config, holding_current) return binsearch_threshold_current( @@ -99,8 +58,7 @@ def calculate_threshold_current(cell, config, holding_current): def binsearch_threshold_current(cell, config, holding_current, min_current, max_current): """Binary search for threshold currents""" mid_current = (min_current + max_current) / 2 - - if abs(max_current - min_current) < config["threshold_current_precision"]: + if abs(max_current - min_current) < config.get("threshold_current_precision", 0.001): spike_count = run_spike_sim( cell, config, @@ -218,13 +176,21 @@ def calculate_rmp_and_rin(cell, config): "stim_end": [config["rin"]["step_stop"]], "stimulus_current": [config["rin"]["step_amp"]], } - features = efel.getFeatureValues([trace], ["voltage_base", "ohmic_input_resistance_vb_ssse"])[0] - rmp = None + features = efel.getFeatureValues( + [trace], ["spike_count", "voltage_base", "ohmic_input_resistance_vb_ssse"] + )[0] + + rmp = 0 if features["voltage_base"] is not None: rmp = features["voltage_base"][0] - rin = None + + rin = -1.0 if features["ohmic_input_resistance_vb_ssse"] is not None: rin = features["ohmic_input_resistance_vb_ssse"][0] + + if features["spike_count"] > 0: + logger.warning("SPIKES! %s, %s", rmp, rin) + return 0.0, -1.0 return rmp, rin @@ -289,12 +255,12 @@ def _isolated_current_evaluation(*args, **kwargs): res = isolate(_current_evaluation, timeout=timeout)(*args, **kwargs) if res is None: res = { - "resting_potential": None, - "input_resistance": None, + "resting_potential": 0.0, + "input_resistance": -1.0, } if not kwargs.get("only_rin", False): - res["holding_current"] = None - res["threshold_current"] = None + res["holding_current"] = 0.0 + res["threshold_current"] = 0.0 return res diff --git a/emodel_generalisation/cli.py b/emodel_generalisation/cli.py index 3c439d9..5f1584d 100644 --- a/emodel_generalisation/cli.py +++ b/emodel_generalisation/cli.py @@ -164,7 +164,7 @@ def compute_currents( "step_stop": 2000.0, "threshold_current_precision": 0.001, "min_threshold_current": 0.0, - "max_threshold_current": 0.2, + "max_threshold_current": 0.1, "spike_at_ais": False, # does not work with placeholder "deterministic": True, "celsius": 34.0, @@ -206,12 +206,15 @@ def compute_currents( ) failed_cells = unique_cells_df[ - unique_cells_df["input_resistance"].isna() | (unique_cells_df["input_resistance"] < 0) + unique_cells_df["input_resistance"].isna() | (unique_cells_df["input_resistance"] <= 0) ].index if len(failed_cells) > 0: - L.info("still %s failed cells (we drop):", len(failed_cells)) + L.info("still %s failed cells (we set default values):", len(failed_cells)) L.info(unique_cells_df.loc[failed_cells]) - unique_cells_df.loc[failed_cells, "mtype"] = None + unique_cells_df.loc[failed_cells, "@dynamics:holding_current"] = 0.0 + unique_cells_df.loc[failed_cells, "@dynamics:threshold_current"] = 0.0 + unique_cells_df.loc[failed_cells, "@dynamics:input_resistance"] = 0.0 + unique_cells_df.loc[failed_cells, "@dynamics:resting_potential"] = -80.0 cols = ["resting_potential", "input_resistance", "exception"] if not only_rin: @@ -427,7 +430,7 @@ def evaluate( "name": pd.Series(dtype=str), } ) - for gid, emodel in enumerate(cells_df.emodel.unique): + for gid, emodel in enumerate(cells_df.emodel.unique()): morph = access_point.get_morphologies(emodel) exemplar_df.loc[gid, "emodel"] = emodel exemplar_df.loc[gid, "path"] = morph["path"] @@ -755,8 +758,13 @@ def _get_resistance_models(exemplar_df, exemplar_data, scales_params): """We fit the scale/Rin relation for AIS and soma.""" models = {} for emodel in exemplar_df.emodel: + Path(f"local/{emodel}").mkdir(parents=True, exist_ok=True) models[emodel] = build_all_resistance_models( - access_point, [emodel], exemplar_data[emodel], scales_params + access_point, + [emodel], + exemplar_data[emodel], + scales_params, + fig_path=Path(f"local/{emodel}"), ) return models @@ -765,6 +773,11 @@ def _get_resistance_models(exemplar_df, exemplar_data, scales_params): _get_resistance_models, exemplar_df.loc[placeholder_mask], exemplar_data, scales_params ) + # needed for internal reuse + for col in cells_df.columns: + if cells_df[col].dtype == "category": + cells_df[col] = cells_df[col].astype("object") + L.info("Adapting AIS and soma of all cells..") cells_df["ais_scaler"] = 0.0 cells_df["soma_scaler"] = 0.0 @@ -777,7 +790,7 @@ def _adapt(): mask = cells_df["emodel"] == emodel if emodel in exemplar_data and not exemplar_data[emodel]["placeholder"]: - L.info("Adapting a non placeholder model...") + L.info("Adapting the non placeholder model %s...", emodel) if len(Morphology(cells_df[mask].head(1)["path"].tolist()[0]).root_sections) == 1: raise ValueError( @@ -792,16 +805,20 @@ def _adapt(): .set_index("emodel") .to_dict() ) - cells_df.loc[mask] = adapt_soma_ais( - cells_df.loc[mask], - access_point, - resistance_models[emodel], - rhos, - parallel_factory=parallel_factory, - min_scale=min_scale, - max_scale=max_scale, - n_steps=2, - ) + with Reuse( + local_dir / f"adapt_df_{emodel}.csv", disable=no_reuse, index_col=0 + ) as reuse: + cells_df.loc[mask] = reuse( + adapt_soma_ais, + cells_df.loc[mask], + access_point, + resistance_models[emodel], + rhos, + parallel_factory=parallel_factory, + min_scale=min_scale, + max_scale=max_scale, + n_steps=2, + ).drop(columns=["exception"]) else: if len(Morphology(cells_df[mask].head(1)["path"].tolist()[0]).root_sections) > 1: diff --git a/emodel_generalisation/mcmc.py b/emodel_generalisation/mcmc.py index 73d77ba..7e6f926 100644 --- a/emodel_generalisation/mcmc.py +++ b/emodel_generalisation/mcmc.py @@ -26,7 +26,7 @@ from functools import partial from pathlib import Path -import matplotlib +import matplotlib as mpl import matplotlib.pyplot as plt import numpy as np import pandas as pd @@ -52,8 +52,6 @@ # pylint: disable=too-many-lines,too-many-locals -matplotlib.use("Agg") - class MarkovChain: """Class to setup and run a markov chain on emodel parameter space.""" @@ -1060,6 +1058,8 @@ def plot_corner( cmap="gnuplot", normalize=False, highlights=None, + sort_params=True, + with_pearson=False, ): """Make a corner plot which consists of scatter plots of all pairs. @@ -1067,9 +1067,10 @@ def plot_corner( feature (str): name of feature for coloring heatmap filename (str): name of figure for corner plot """ - params = np.array(sorted(df.normalized_parameters.columns.to_list())) + params = np.array(df.normalized_parameters.columns.to_list()) _params = np.array([PARAM_LABELS.get(p, p) for p in params]) - params = params[np.argsort(_params)] + if sort_params: + params = params[np.argsort(_params)] n_params = len(params) # get feature data @@ -1096,6 +1097,10 @@ def plot_corner( fig = plt.figure(figsize=(5 + 0.5 * n_params, 5 + 0.5 * n_params)) gs = fig.add_gridspec(n_params, n_params, hspace=0.1, wspace=0.1) im = None + if with_pearson: + _cmap = plt.get_cmap("coolwarm") + norm = mpl.colors.Normalize(vmin=-0.8, vmax=0.8) + pearson_colors = mpl.cm.ScalarMappable(norm=norm, cmap=_cmap) # pylint: disable=too-many-nested-blocks for i, param1 in enumerate(params): _param1 = PARAM_LABELS.get(param1, param1) @@ -1112,6 +1117,13 @@ def plot_corner( ax.set_frame_on(False) elif j < i: ax.set_frame_on(True) + if with_pearson: + pearson = pearsonr( + df[("normalized_parameters", param1)].to_numpy(), + df[("normalized_parameters", param2)].to_numpy(), + ) + ax.spines[:].set_color(pearson_colors.to_rgba(pearson[0])) + ax.spines[:].set(lw=3.0) im = _plot_2d_data( ax, m[i][j], vmin, vmax, rev=feature is not None, cmap=cmap, normalize=normalize ) diff --git a/emodel_generalisation/model/bpopt.py b/emodel_generalisation/model/bpopt.py index 93a9513..9b23981 100644 --- a/emodel_generalisation/model/bpopt.py +++ b/emodel_generalisation/model/bpopt.py @@ -113,8 +113,6 @@ def calculate_bpo_score(self, responses): if feature_value is None: return self.max_score - if self.efel_feature_name == "peak_voltage": - print(feature_value, self.exp_mean, self.exp_std) return abs(feature_value - self.exp_mean) / self.exp_std def _construct_efel_trace(self, responses): @@ -617,11 +615,11 @@ def get_holding_current_from_voltage( target_voltage=target_voltage, stimulus_duration=2000.0, max_depth=max_depth, + strict_bounds=False, ) if self.stimulus.holding_current is not None: hold_prot.stimulus.delay = 1000 hold_prot.stimulus.holding_current = self.stimulus.holding_current - return hold_prot.run(cell_model, param_values, sim, isolate, timeout, responses)[ "bpo_holding_current" ] @@ -632,6 +630,7 @@ def run( """Run protocol""" if self.stimulus.holding_voltage is not None: + self.stimulus.holding_current = None self.stimulus.holding_current = self.get_holding_current_from_voltage( self.stimulus.holding_voltage, cell_model, @@ -833,7 +832,7 @@ def __init__( location, target_voltage=None, voltage_precision=0.1, - stimulus_duration=500.0, + stimulus_duration=1000.0, upper_bound=0.2, lower_bound=-0.2, strict_bounds=True, @@ -900,11 +899,11 @@ def __init__( # start/end are not used by this feature self.spike_feature = ephys.efeatures.eFELFeature( - name=f"{name}.Spikecount", - efel_feature_name="Spikecount", + name=f"{name}.spike_count_stimint", + efel_feature_name="spike_count_stimint", recording_names={"": f"{name}.{location.name}.v"}, - stim_start=0, - stim_end=1000, + stim_start=500, # discard beginning to allow for a couple of initial APs + stim_end=2000, exp_mean=1, exp_std=0.1, ) @@ -913,7 +912,6 @@ def get_voltage_base( self, holding_current, cell_model, param_values, sim, isolate, timeout=None ): """Calculate voltage base for a certain holding current""" - self.stimuli[0].amp = holding_current response = BPEMProtocol.run( self, cell_model, param_values, sim=sim, isolate=isolate, timeout=timeout @@ -1103,7 +1101,6 @@ def __init__( "totduration": stimulus_totduration, "holding_current": None, } - self.recording_name = f"{name}.{location.name}.v" stimulus = eCodes["step"](location=location, **stimulus_definition) recordings = [ diff --git a/emodel_generalisation/model/ecodes.py b/emodel_generalisation/model/ecodes.py index d5396a6..c25e161 100644 --- a/emodel_generalisation/model/ecodes.py +++ b/emodel_generalisation/model/ecodes.py @@ -133,7 +133,7 @@ def __init__(self, location, **kwargs): location(Location): location of stimulus """ - self.amp = kwargs.get("amp", None) + self.amp = kwargs.get("amp", 0.0) self.amp_rel = kwargs.get("thresh_perc", 200.0) self.amp_voltage = kwargs.get("amp_voltage", None) diff --git a/emodel_generalisation/model/evaluation.py b/emodel_generalisation/model/evaluation.py index ed32b0d..e5a7138 100644 --- a/emodel_generalisation/model/evaluation.py +++ b/emodel_generalisation/model/evaluation.py @@ -451,7 +451,7 @@ def _define_Rin_protocol( def _define_holding_protocol( - efeatures, strict_bounds=False, ais_recording=False, max_depth=7, stimulus_duration=500.0 + efeatures, strict_bounds=False, ais_recording=False, max_depth=7, stimulus_duration=1000.0 ): """Define the search holding current protocol""" target_voltage = None @@ -1013,7 +1013,7 @@ def __init__( self.rin_step_duration = 500.0 self.rin_step_amp = -0.02 self.rin_totduration = self.rin_step_delay + self.rin_step_duration - self.search_holding_duration = 500.0 + self.search_holding_duration = 1000.0 self.search_threshold_step_delay = 500.0 self.search_threshold_step_duration = 2000.0 self.search_threshold_totduration = 3000.0 diff --git a/emodel_generalisation/utils.py b/emodel_generalisation/utils.py index 4233013..4305910 100644 --- a/emodel_generalisation/utils.py +++ b/emodel_generalisation/utils.py @@ -25,12 +25,14 @@ from functools import partial from hashlib import sha256 from itertools import cycle +from multiprocessing.context import TimeoutError # pylint: disable=redefined-builtin from pathlib import Path import matplotlib.pyplot as plt import numpy as np import pandas as pd from bluepyopt.ephys.responses import TimeVoltageResponse +from bluepyparallel.parallel import NestedPool from matplotlib.backends.backend_pdf import PdfPages from scipy.cluster.hierarchy import dendrogram from scipy.cluster.hierarchy import linkage @@ -94,14 +96,14 @@ def get_scores(morphs_combos_df, features_to_ignore=None, features_to_keep=None, morphs_combos_df.apply(filter_features, axis=1) morphs_combos_df.loc[:, "median_score"] = morphs_combos_df["scores"].apply( - lambda score: np.clip(np.median(list(score.values())), 0, clip) - if isinstance(score, dict) - else np.nan + lambda score: ( + np.clip(np.median(list(score.values())), 0, clip) if isinstance(score, dict) else np.nan + ) ) morphs_combos_df.loc[:, "max_score"] = morphs_combos_df["scores"].apply( - lambda score: np.clip(np.max(list(score.values())), 0, clip) - if isinstance(score, dict) - else np.nan + lambda score: ( + np.clip(np.max(list(score.values())), 0, clip) if isinstance(score, dict) else np.nan + ) ) morphs_combos_df.loc[:, "cost"] = morphs_combos_df["scores"].apply( lambda score: np.max(list(score.values())) if isinstance(score, dict) else np.nan @@ -112,9 +114,9 @@ def get_scores(morphs_combos_df, features_to_ignore=None, features_to_keep=None, def get_feature_df(df, filters=None): """Get feature df from complete df.""" feature_df = df["features"].apply( - lambda json_str: pd.Series(json.loads(json_str)) - if isinstance(json_str, str) - else pd.Series(dtype=float) + lambda json_str: ( + pd.Series(json.loads(json_str)) if isinstance(json_str, str) else pd.Series(dtype=float) + ) ) if filters is not None: feature_df = feature_df.drop( @@ -129,9 +131,11 @@ def get_feature_df(df, filters=None): def get_score_df(df, filters=None): """Get score df from complete df.""" score_df = df["scores"].apply( - lambda json_str: pd.Series(json.loads(json_str)) - if isinstance(json_str, str) - else lambda json_str: pd.Series(json_str, dtype=float) + lambda json_str: ( + pd.Series(json.loads(json_str)) + if isinstance(json_str, str) + else lambda json_str: pd.Series(json_str, dtype=float) + ) ) if filters is not None: score_df = score_df.drop( @@ -214,3 +218,44 @@ def load_json(filepath, **kwargs): def write_json(filepath, data, **kwargs): """Write json file.""" Path(filepath).write_text(json.dumps(data, **kwargs), encoding="utf-8") + + +def isolate(func, timeout=None): + """Isolate a generic function for independent NEURON instances. + + It must be used in conjunction with NestedPool. + + Example: + + .. code-block:: python + + def _to_be_isolated(morphology_path, point): + cell = nrnhines.get_NRN_cell(morphology_path) + return nrnhines.point_to_section_end(cell.icell.all, point) + + def _isolated(morph_data): + return nrnhines.isolate(_to_be_isolated)(*morph_data) + + with nrnhines.NestedPool(processes=n_workers) as pool: + result = pool.imap_unordered(_isolated, data) + + + Args: + func (function): function to isolate + + Returns: + the isolated function + + Note: it does not work as decorator. + """ + + def func_isolated(*args, **kwargs): + with NestedPool(1, maxtasksperchild=1) as pool: + res = pool.apply_async(func, args, kwargs) + try: + out = res.get(timeout=timeout) + except TimeoutError: # pragma: no cover + out = None + return out + + return func_isolated diff --git a/tests/data/adapt_df.csv b/tests/data/adapt_df.csv index e9130cd..e463da9 100644 --- a/tests/data/adapt_df.csv +++ b/tests/data/adapt_df.csv @@ -1,7 +1,7 @@ -Unnamed: 0,morph_class,mtype,region,hemisphere,model_template,synapse_class,etype,layer,morphology,x,y,z,orientation,emodel,ais_scaler,soma_scaler,ais_model,soma_model -1,PYR,L5_TPC:A,MOs@right,right,hoc:generic_model,EXC,generic,5,hashed/be/fb/befbcea89657fdd07ae181b977b9e320,3887.062211566772,1924.0204607290975,6775.602532718369,"[[1. 0. 0.] +Unnamed: 0,morph_class,mtype,morphology,model_template,region,synapse_class,layer,hemisphere,etype,x,y,z,orientation,emodel,ais_scaler,soma_scaler,ais_model,soma_model +1,PYR,L5_TPC:A,hashed/be/fb/befbcea89657fdd07ae181b977b9e320,hoc:generic_model,MOs@right,EXC,5,right,generic,3887.062211566772,1924.0204607290975,6775.602532718369,"[[1. 0. 0.] [0. 1. 0.] [0. 0. 1.]]",generic_model,1.1,1.1,"{""popt"": [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]}","{""soma_surface"": 1167.714655051592, ""soma_radius"": 10.0}" -2,PYR,L5_TPC:A,MOp@right,right,hoc:generic_model,EXC,generic,5,hashed/d6/a9/d6a9976c7588ae0111f17f28e39a1e6d,6291.256509715739,1139.1852640664094,6812.970413610967,"[[1. 0. 0.] +2,PYR,L5_TPC:A,hashed/d6/a9/d6a9976c7588ae0111f17f28e39a1e6d,hoc:generic_model,MOp@right,EXC,5,right,generic,6291.256509715739,1139.1852640664094,6812.970413610967,"[[1. 0. 0.] [0. 1. 0.] [0. 0. 1.]]",generic_model,1.1,1.1,"{""popt"": [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]}","{""soma_surface"": 1167.714655051592, ""soma_radius"": 10.0}" diff --git a/tests/data/adapted_evaluation_df.csv b/tests/data/adapted_evaluation_df.csv index 015d3a7..33ff02d 100644 --- a/tests/data/adapted_evaluation_df.csv +++ b/tests/data/adapted_evaluation_df.csv @@ -1,7 +1,7 @@ -morph_class,mtype,region,hemisphere,model_template,synapse_class,etype,layer,morphology,ais_scaler,soma_scaler,x,y,z,orientation,emodel,ais_model,soma_model,exception,features,scores,trace_data -PYR,L5_TPC:A,MOs@right,right,hoc:generic_model,EXC,generic,5,hashed/be/fb/befbcea89657fdd07ae181b977b9e320,1.1,1.1,3887.062211566772,1924.0204607290975,6775.602532718369,"[[1. 0. 0.] +morph_class,mtype,morphology,model_template,region,synapse_class,layer,hemisphere,etype,ais_scaler,soma_scaler,x,y,z,orientation,emodel,ais_model,soma_model,exception,features,scores,trace_data +PYR,L5_TPC:A,hashed/be/fb/befbcea89657fdd07ae181b977b9e320,hoc:generic_model,MOs@right,EXC,5,right,generic,1.1,1.1,3887.062211566772,1924.0204607290975,6775.602532718369,"[[1. 0. 0.] [0. 1. 0.] - [0. 0. 1.]]",generic_model,"{""popt"": [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]}","{""soma_radius"": 10.0, ""soma_surface"": 1167.714655051592}",,"{""IDRest_150.soma.v.mean_frequency"": 17.570535563419252, ""IDRest_150.soma.v.AP_amplitude"": 55.39737419147699, ""IDRest_200.soma.v.mean_frequency"": 24.347355606626557, ""IDRest_200.soma.v.AP_amplitude"": 55.783184825944716, ""IDRest_280.soma.v.mean_frequency"": 32.55285600443829, ""IDRest_280.soma.v.AP_amplitude"": 55.9149320756509, ""RMPProtocol.soma.v.voltage_base"": -70.28549140927481, ""RMPProtocol.soma.v.Spikecount"": 5.0, ""SearchHoldingCurrent.soma.v.bpo_holding_current"": -0.053125000000000006, ""SearchHoldingCurrent.soma.v.voltage_base"": -83.00675420359404, ""RinProtocol.soma.v.ohmic_input_resistance_vb_ssse"": 169.91847638730917, ""SearchThresholdCurrent.soma.v.bpo_threshold_current"": 0.046635434877685934}","{""IDRest_150.soma.v.mean_frequency"": 1.285267781709626, ""IDRest_150.soma.v.AP_amplitude"": 7.301312904261504, ""IDRest_200.soma.v.mean_frequency"": 2.1736778033132786, ""IDRest_200.soma.v.AP_amplitude"": 7.108407587027645, ""IDRest_280.soma.v.mean_frequency"": 3.776428002219145, ""IDRest_280.soma.v.AP_amplitude"": 7.0425339621745495, ""RMPProtocol.soma.v.voltage_base"": 1.3429017181450376, ""RMPProtocol.soma.v.Spikecount"": 250.0, ""SearchHoldingCurrent.soma.v.bpo_holding_current"": 0.09375, ""SearchHoldingCurrent.soma.v.voltage_base"": 0.0067542035940419964, ""RinProtocol.soma.v.ohmic_input_resistance_vb_ssse"": 3.3008152361269083, ""SearchThresholdCurrent.soma.v.bpo_threshold_current"": 0.7327086975537186}", -PYR,L5_TPC:A,MOp@right,right,hoc:generic_model,EXC,generic,5,hashed/d6/a9/d6a9976c7588ae0111f17f28e39a1e6d,1.1,1.1,6291.256509715739,1139.1852640664094,6812.970413610967,"[[1. 0. 0.] + [0. 0. 1.]]",generic_model,"{""popt"": [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]}","{""soma_radius"": 10.0, ""soma_surface"": 1167.714655051592}",,"{""IDRest_150.soma.v.mean_frequency"": 17.570535563419252, ""IDRest_150.soma.v.AP_amplitude"": 55.39737419147699, ""IDRest_200.soma.v.mean_frequency"": 24.347355606626557, ""IDRest_200.soma.v.AP_amplitude"": 55.783184825944716, ""IDRest_280.soma.v.mean_frequency"": 32.55285600443829, ""IDRest_280.soma.v.AP_amplitude"": 55.9149320756509, ""RMPProtocol.soma.v.voltage_base"": -70.28549140927481, ""RMPProtocol.soma.v.Spikecount"": 5.0, ""SearchHoldingCurrent.soma.v.bpo_holding_current"": -0.053125000000000006, ""SearchHoldingCurrent.soma.v.voltage_base"": -83.00777720496266, ""RinProtocol.soma.v.ohmic_input_resistance_vb_ssse"": 169.91847638730917, ""SearchThresholdCurrent.soma.v.bpo_threshold_current"": 0.046635434877685934}","{""IDRest_150.soma.v.mean_frequency"": 1.285267781709626, ""IDRest_150.soma.v.AP_amplitude"": 7.301312904261504, ""IDRest_200.soma.v.mean_frequency"": 2.1736778033132786, ""IDRest_200.soma.v.AP_amplitude"": 7.108407587027645, ""IDRest_280.soma.v.mean_frequency"": 3.776428002219145, ""IDRest_280.soma.v.AP_amplitude"": 7.0425339621745495, ""RMPProtocol.soma.v.voltage_base"": 1.3429017181450376, ""RMPProtocol.soma.v.Spikecount"": 250.0, ""SearchHoldingCurrent.soma.v.bpo_holding_current"": 0.09375, ""SearchHoldingCurrent.soma.v.voltage_base"": 0.007777204962664541, ""RinProtocol.soma.v.ohmic_input_resistance_vb_ssse"": 3.3008152361269083, ""SearchThresholdCurrent.soma.v.bpo_threshold_current"": 0.7327086975537186}", +PYR,L5_TPC:A,hashed/d6/a9/d6a9976c7588ae0111f17f28e39a1e6d,hoc:generic_model,MOp@right,EXC,5,right,generic,1.1,1.1,6291.256509715739,1139.1852640664094,6812.970413610967,"[[1. 0. 0.] [0. 1. 0.] - [0. 0. 1.]]",generic_model,"{""popt"": [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]}","{""soma_radius"": 10.0, ""soma_surface"": 1167.714655051592}",,"{""IDRest_150.soma.v.mean_frequency"": 19.51974691636182, ""IDRest_150.soma.v.AP_amplitude"": 53.64663901095507, ""IDRest_200.soma.v.mean_frequency"": 26.99191929412969, ""IDRest_200.soma.v.AP_amplitude"": 53.955235597646244, ""IDRest_280.soma.v.mean_frequency"": 36.01800900445903, ""IDRest_280.soma.v.AP_amplitude"": 54.06750345591785, ""RMPProtocol.soma.v.voltage_base"": -78.16946608707919, ""RMPProtocol.soma.v.Spikecount"": 6.0, ""SearchHoldingCurrent.soma.v.bpo_holding_current"": -0.06875, ""SearchHoldingCurrent.soma.v.voltage_base"": -83.09765493391816, ""RinProtocol.soma.v.ohmic_input_resistance_vb_ssse"": 147.51245669415312, ""SearchThresholdCurrent.soma.v.bpo_threshold_current"": 0.060956196156172524}","{""IDRest_150.soma.v.mean_frequency"": 2.259873458180911, ""IDRest_150.soma.v.AP_amplitude"": 8.176680494522463, ""IDRest_200.soma.v.mean_frequency"": 3.4959596470648453, ""IDRest_200.soma.v.AP_amplitude"": 8.022382201176882, ""IDRest_280.soma.v.mean_frequency"": 5.509004502229516, ""IDRest_280.soma.v.AP_amplitude"": 7.96624827204107, ""RMPProtocol.soma.v.voltage_base"": 0.2338932174158373, ""RMPProtocol.soma.v.Spikecount"": 250.0, ""SearchHoldingCurrent.soma.v.bpo_holding_current"": 0.0625, ""SearchHoldingCurrent.soma.v.voltage_base"": 0.09765493391816449, ""RinProtocol.soma.v.ohmic_input_resistance_vb_ssse"": 3.524875433058469, ""SearchThresholdCurrent.soma.v.bpo_threshold_current"": 1.0191239231234503}", + [0. 0. 1.]]",generic_model,"{""popt"": [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]}","{""soma_radius"": 10.0, ""soma_surface"": 1167.714655051592}",,"{""IDRest_150.soma.v.mean_frequency"": 19.51974691636182, ""IDRest_150.soma.v.AP_amplitude"": 53.64663901095507, ""IDRest_200.soma.v.mean_frequency"": 26.99191929412969, ""IDRest_200.soma.v.AP_amplitude"": 53.955235597646244, ""IDRest_280.soma.v.mean_frequency"": 36.01800900445903, ""IDRest_280.soma.v.AP_amplitude"": 54.06750345591785, ""RMPProtocol.soma.v.voltage_base"": -78.16946608707919, ""RMPProtocol.soma.v.Spikecount"": 6.0, ""SearchHoldingCurrent.soma.v.bpo_holding_current"": -0.06875, ""SearchHoldingCurrent.soma.v.voltage_base"": -83.09902934622046, ""RinProtocol.soma.v.ohmic_input_resistance_vb_ssse"": 147.51245669415312, ""SearchThresholdCurrent.soma.v.bpo_threshold_current"": 0.060956196156172524}","{""IDRest_150.soma.v.mean_frequency"": 2.259873458180911, ""IDRest_150.soma.v.AP_amplitude"": 8.176680494522463, ""IDRest_200.soma.v.mean_frequency"": 3.4959596470648453, ""IDRest_200.soma.v.AP_amplitude"": 8.022382201176882, ""IDRest_280.soma.v.mean_frequency"": 5.509004502229516, ""IDRest_280.soma.v.AP_amplitude"": 7.96624827204107, ""RMPProtocol.soma.v.voltage_base"": 0.2338932174158373, ""RMPProtocol.soma.v.Spikecount"": 250.0, ""SearchHoldingCurrent.soma.v.bpo_holding_current"": 0.0625, ""SearchHoldingCurrent.soma.v.voltage_base"": 0.09902934622046189, ""RinProtocol.soma.v.ohmic_input_resistance_vb_ssse"": 3.524875433058469, ""SearchThresholdCurrent.soma.v.bpo_threshold_current"": 1.0191239231234503}", diff --git a/tests/data/evaluation_df.csv b/tests/data/evaluation_df.csv index ee07c68..8e93642 100644 --- a/tests/data/evaluation_df.csv +++ b/tests/data/evaluation_df.csv @@ -1,7 +1,7 @@ -mtype,etype,model_template,morph_class,synapse_class,layer,region,hemisphere,morphology,x,y,z,orientation,emodel,ais_model,soma_model,exception,features,scores,trace_data -L5_TPC:A,generic,hoc:generic_model,PYR,EXC,5,MOs@right,right,hashed/be/fb/befbcea89657fdd07ae181b977b9e320,3887.062211566772,1924.0204607290975,6775.602532718369,"[[1. 0. 0.] +morph_class,mtype,morphology,model_template,region,synapse_class,layer,hemisphere,etype,x,y,z,orientation,emodel,exception,features,scores,trace_data +PYR,L5_TPC:A,hashed/be/fb/befbcea89657fdd07ae181b977b9e320,hoc:generic_model,MOs@right,EXC,5,right,generic,3887.062211566772,1924.0204607290975,6775.602532718369,"[[1. 0. 0.] [0. 1. 0.] - [0. 0. 1.]]",generic_model,"{""popt"": [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]}","{""soma_radius"": 10.0, ""soma_surface"": 1167.714655051592}",,"{""IDRest_150.soma.v.mean_frequency"": 17.01809022989483, ""IDRest_150.soma.v.AP_amplitude"": 57.74731715766996, ""IDRest_200.soma.v.mean_frequency"": 24.038862828210206, ""IDRest_200.soma.v.AP_amplitude"": 58.078713635979184, ""IDRest_280.soma.v.mean_frequency"": 32.58766421123901, ""IDRest_280.soma.v.AP_amplitude"": 58.245488451959865, ""RMPProtocol.soma.v.voltage_base"": -66.35232862203746, ""RMPProtocol.soma.v.Spikecount"": 5.0, ""SearchHoldingCurrent.soma.v.bpo_holding_current"": -0.051562500000000004, ""SearchHoldingCurrent.soma.v.voltage_base"": -82.94181734045071, ""RinProtocol.soma.v.ohmic_input_resistance_vb_ssse"": 179.45614761441178, ""SearchThresholdCurrent.soma.v.bpo_threshold_current"": 0.043736971847982836}","{""IDRest_150.soma.v.mean_frequency"": 1.0090451149474156, ""IDRest_150.soma.v.AP_amplitude"": 6.126341421165019, ""IDRest_200.soma.v.mean_frequency"": 2.019431414105103, ""IDRest_200.soma.v.AP_amplitude"": 5.960643182010407, ""IDRest_280.soma.v.mean_frequency"": 3.7938321056195043, ""IDRest_280.soma.v.AP_amplitude"": 5.877255774020064, ""RMPProtocol.soma.v.voltage_base"": 2.1295342755925075, ""RMPProtocol.soma.v.Spikecount"": 250.0, ""SearchHoldingCurrent.soma.v.bpo_holding_current"": 0.096875, ""SearchHoldingCurrent.soma.v.voltage_base"": 0.058182659549288473, ""RinProtocol.soma.v.ohmic_input_resistance_vb_ssse"": 3.2054385238558827, ""SearchThresholdCurrent.soma.v.bpo_threshold_current"": 0.6747394369596567}", -L5_TPC:A,generic,hoc:generic_model,PYR,EXC,5,MOp@right,right,hashed/d6/a9/d6a9976c7588ae0111f17f28e39a1e6d,6291.256509715739,1139.1852640664094,6812.970413610967,"[[1. 0. 0.] + [0. 0. 1.]]",generic_model,,"{""IDRest_150.soma.v.mean_frequency"": 17.01809022989483, ""IDRest_150.soma.v.AP_amplitude"": 57.74731715766996, ""IDRest_200.soma.v.mean_frequency"": 24.038862828210206, ""IDRest_200.soma.v.AP_amplitude"": 58.078713635979184, ""IDRest_280.soma.v.mean_frequency"": 32.58766421123901, ""IDRest_280.soma.v.AP_amplitude"": 58.245488451959865, ""RMPProtocol.soma.v.voltage_base"": -66.35232862203746, ""RMPProtocol.soma.v.Spikecount"": 5.0, ""SearchHoldingCurrent.soma.v.bpo_holding_current"": -0.051562500000000004, ""SearchHoldingCurrent.soma.v.voltage_base"": -82.94169679712304, ""RinProtocol.soma.v.ohmic_input_resistance_vb_ssse"": 179.45614761441178, ""SearchThresholdCurrent.soma.v.bpo_threshold_current"": 0.043736971847982836}","{""IDRest_150.soma.v.mean_frequency"": 1.0090451149474156, ""IDRest_150.soma.v.AP_amplitude"": 6.126341421165019, ""IDRest_200.soma.v.mean_frequency"": 2.019431414105103, ""IDRest_200.soma.v.AP_amplitude"": 5.960643182010407, ""IDRest_280.soma.v.mean_frequency"": 3.7938321056195043, ""IDRest_280.soma.v.AP_amplitude"": 5.877255774020064, ""RMPProtocol.soma.v.voltage_base"": 2.1295342755925075, ""RMPProtocol.soma.v.Spikecount"": 250.0, ""SearchHoldingCurrent.soma.v.bpo_holding_current"": 0.096875, ""SearchHoldingCurrent.soma.v.voltage_base"": 0.05830320287695656, ""RinProtocol.soma.v.ohmic_input_resistance_vb_ssse"": 3.2054385238558827, ""SearchThresholdCurrent.soma.v.bpo_threshold_current"": 0.6747394369596567}", +PYR,L5_TPC:A,hashed/d6/a9/d6a9976c7588ae0111f17f28e39a1e6d,hoc:generic_model,MOp@right,EXC,5,right,generic,6291.256509715739,1139.1852640664094,6812.970413610967,"[[1. 0. 0.] [0. 1. 0.] - [0. 0. 1.]]",generic_model,"{""popt"": [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]}","{""soma_radius"": 10.0, ""soma_surface"": 1167.714655051592}",,"{""IDRest_150.soma.v.mean_frequency"": 21.307672480429662, ""IDRest_150.soma.v.AP_amplitude"": 51.276480916739004, ""IDRest_200.soma.v.mean_frequency"": 30.406432649716162, ""IDRest_200.soma.v.AP_amplitude"": 51.25679986852141, ""IDRest_280.soma.v.mean_frequency"": 41.679271284949614, ""IDRest_280.soma.v.AP_amplitude"": 50.94403679026918, ""RMPProtocol.soma.v.voltage_base"": -72.26510227533872, ""RMPProtocol.soma.v.Spikecount"": 6.0, ""SearchHoldingCurrent.soma.v.bpo_holding_current"": -0.06718750000000001, ""SearchHoldingCurrent.soma.v.voltage_base"": -83.08857871429298, ""RinProtocol.soma.v.ohmic_input_resistance_vb_ssse"": 155.80389971422335, ""SearchThresholdCurrent.soma.v.bpo_threshold_current"": 0.059734491885581026}","{""IDRest_150.soma.v.mean_frequency"": 3.153836240214831, ""IDRest_150.soma.v.AP_amplitude"": 9.361759541630494, ""IDRest_200.soma.v.mean_frequency"": 5.203216324858081, ""IDRest_200.soma.v.AP_amplitude"": 9.37160006573929, ""IDRest_280.soma.v.mean_frequency"": 8.339635642474807, ""IDRest_280.soma.v.AP_amplitude"": 9.52798160486541, ""RMPProtocol.soma.v.voltage_base"": 0.9469795449322561, ""RMPProtocol.soma.v.Spikecount"": 250.0, ""SearchHoldingCurrent.soma.v.bpo_holding_current"": 0.06562499999999999, ""SearchHoldingCurrent.soma.v.voltage_base"": 0.08857871429297859, ""RinProtocol.soma.v.ohmic_input_resistance_vb_ssse"": 3.4419610028577665, ""SearchThresholdCurrent.soma.v.bpo_threshold_current"": 0.9946898377116205}", + [0. 0. 1.]]",generic_model,,"{""IDRest_150.soma.v.mean_frequency"": 21.307672480429662, ""IDRest_150.soma.v.AP_amplitude"": 51.276480916739004, ""IDRest_200.soma.v.mean_frequency"": 30.406432649716162, ""IDRest_200.soma.v.AP_amplitude"": 51.25679986852141, ""IDRest_280.soma.v.mean_frequency"": 41.679271284949614, ""IDRest_280.soma.v.AP_amplitude"": 50.94403679026918, ""RMPProtocol.soma.v.voltage_base"": -72.26510227533872, ""RMPProtocol.soma.v.Spikecount"": 6.0, ""SearchHoldingCurrent.soma.v.bpo_holding_current"": -0.06718750000000001, ""SearchHoldingCurrent.soma.v.voltage_base"": -83.08653117477492, ""RinProtocol.soma.v.ohmic_input_resistance_vb_ssse"": 155.80389971422335, ""SearchThresholdCurrent.soma.v.bpo_threshold_current"": 0.059734491885581026}","{""IDRest_150.soma.v.mean_frequency"": 3.153836240214831, ""IDRest_150.soma.v.AP_amplitude"": 9.361759541630494, ""IDRest_200.soma.v.mean_frequency"": 5.203216324858081, ""IDRest_200.soma.v.AP_amplitude"": 9.37160006573929, ""IDRest_280.soma.v.mean_frequency"": 8.339635642474807, ""IDRest_280.soma.v.AP_amplitude"": 9.52798160486541, ""RMPProtocol.soma.v.voltage_base"": 0.9469795449322561, ""RMPProtocol.soma.v.Spikecount"": 250.0, ""SearchHoldingCurrent.soma.v.bpo_holding_current"": 0.06562499999999999, ""SearchHoldingCurrent.soma.v.voltage_base"": 0.0865311747749189, ""RinProtocol.soma.v.ohmic_input_resistance_vb_ssse"": 3.4419610028577665, ""SearchThresholdCurrent.soma.v.bpo_threshold_current"": 0.9946898377116205}", diff --git a/tests/test_cli.py b/tests/test_cli.py index c81609f..1b33d5a 100644 --- a/tests/test_cli.py +++ b/tests/test_cli.py @@ -1,4 +1,5 @@ """Test cli module.""" + import json from pathlib import Path @@ -66,7 +67,7 @@ def test_compute_currents(cli_runner, tmpdir): def test_evaluate(cli_runner, tmpdir): - """Tetst cli evaluate.""" + """Test cli evaluate.""" # fmt: off response = cli_runner.invoke( tested.cli, @@ -195,18 +196,17 @@ def test_adapt(cli_runner, tmpdir): df = CellCollection().load_sonata(tmpdir / "sonata_currents_adapted.h5").as_dataframe() npt.assert_allclose( df["@dynamics:resting_potential"].to_list(), - [-72.841806, -71.32893], + [0.0, 0.0], # we get spikes at rest rtol=1e-5, ) - df.to_csv("test.csv") npt.assert_allclose( df["@dynamics:input_resistance"].to_list(), - [105.342194, 1863.809101], + [-1.0, -1.0], # we get spikes at rest rtol=1e-5, ) npt.assert_allclose( df["@dynamics:holding_current"].to_list(), - [-0.06431251604510635, -0.08120532523605561], + [-0.064316, -0.081205], rtol=1e-5, ) npt.assert_allclose( diff --git a/tests/test_evaluation.py b/tests/test_evaluation.py index 7dc152e..bbe1f2f 100644 --- a/tests/test_evaluation.py +++ b/tests/test_evaluation.py @@ -39,8 +39,8 @@ def test_feature_evaluation(morphs_combos_df, access_point): "IDRest_280.soma.v.AP_amplitude": 0.4711507560288358, "RMPProtocol.soma.v.voltage_base": 0.33131119446020707, "RMPProtocol.soma.v.Spikecount": 0.0, - "SearchHoldingCurrent.soma.v.bpo_holding_current": 0.18281250000000002, - "SearchHoldingCurrent.soma.v.voltage_base": 0.043970045446585004, + "SearchHoldingCurrent.soma.v.bpo_holding_current": 0.1828125, + "SearchHoldingCurrent.soma.v.voltage_base": 0.038499, "RinProtocol.soma.v.ohmic_input_resistance_vb_ssse": 0.17312024013731503, "SearchThresholdCurrent.soma.v.bpo_threshold_current": 0.10935385147499851, } @@ -84,7 +84,7 @@ def test_evaluate_soma_rin(morphs_combos_df, access_point): def test_evaluate_ais_rin(morphs_combos_df, access_point): """ """ df = evaluation.evaluate_ais_rin(morphs_combos_df, access_point) - assert_allclose(df.loc[0, "rin_ais"], 131668.9953893715, rtol=1e-5) + assert_allclose(df.loc[0, "rin_ais"], 131668.995389, rtol=1e-5) def test_evaluate_somadend_rin(morphs_combos_df, access_point): @@ -102,4 +102,4 @@ def test_evaluate_rho(morphs_combos_df, access_point): def test_evaluate_rho_axon(morphs_combos_df, access_point): """ """ df = evaluation.evaluate_rho_axon(morphs_combos_df, access_point) - assert_allclose(df.loc[0, "rho_axon"], 244.2217237122562, rtol=1e-5) + assert_allclose(df.loc[0, "rho_axon"], 244.221723, rtol=1e-5)