[MRG] Validate param files and network configs against model_variant (#14)

* Validate param files and network configs against model_variant

Introduce a 'model_variant' entry in the param .json files so that each
network model can check that it was given parameters written for it.

- preserve 'model_variant' in Params, which otherwise drops any key that
  is not part of the legacy defaults. This also fixes default values
  being silently substituted for the parameters of default_duecker_ET.json,
  whose 'sim_prefix' no longer marked it as a duecker_ET_model param file
- store the variant on the Network and write it to the network config
- accept abbreviated variant names, e.g. 'neymotin' or 'duecker_ET'
- check that the params define the cell types of the network being built:
  pyramidal and basket cells for neymotin_2020_model (and the law_2021_model
  and calcium_model networks built on it), pyramidal cells and interneurons
  for duecker_ET_model
- check the cell types of a network config against its variant when reading
  it back in, pointing out when a neymotin_2020_model network is being read
  as a duecker_ET_model

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

* test pass

* add model variant check method

* remove sim_prefix

* update model_variant bsl_cor check

* [MRG] MAINT cap ruff version due to breaking changes (#1333)

* MAINT cap ruff version due to breaking changes

* FIX forgot to update standalone ruff install

* removed model variant bsl_cor checks because tests failing

* remove model validation cheks

* remove stale import

* remove sim_prefix

* trivial change to re-run test

* ruff

* remove final sim_prefix

* update test and validation function

* DOC some docstring and comment changes

* FIX remove str type check on params_fname is none passed

---------

Co-authored-by: Claude Opus 5 <noreply@anthropic.com>
Co-authored-by: Austin E. Soplata <asoplata@users.noreply.github.com>
Co-authored-by: Austin E. Soplata <me@asoplata.com>

Author: 3 people

hnn_core/dipole.py

@@ -116,8 +116,8 @@ def simulate_dipole(
 
     if bsl_cor is None:
         bsl_cor = "jones"
-    elif bsl_cor not in {"jones", "duecker", "none"}:
-        raise ValueError("'bsl_cor' must be 'jones', 'duecker' or 'none'")
+    elif bsl_cor not in {"jones", "neymotin", "duecker", "none"}:
+        raise ValueError("'bsl_cor' must be 'jones', 'neymotin', 'duecker' or 'none'")
 
     net._instantiate_drives(n_trials=n_trials, tstop=tstop)
     net._reset_rec_arrays()
@@ -363,72 +363,8 @@ def _rmse(dpl, exp_dpl, tstart=0.0, tstop=0.0, weights=None):
     return np.sqrt((weights * ((dpl1 - dpl2) ** 2)).sum() / weights.sum())
 
 
-# # KDTODO very different _rmse, copied here
-# def _rmse(dpl, exp_dpl, tstart=0.0, tstop=0.0, weights=None):
-#     """Calculates RMSE between data in dpl and exp_dpl
-#     Parameters
-#     ----------
-#     dpl : instance of Dipole
-#         A dipole object with simulated data
-#     exp_dpl : instance of Dipole
-#         A dipole object with experimental data
-#     tstart : None | float
-#         Time at beginning of range over which to calculate RMSE
-#     tstop : None | float
-#         Time at end of range over which to calculate RMSE
-#     weights : None | array
-#         An array of weights to be applied to each point in
-#         simulated dpl. Must have length >= dpl.data
-#         If None, weights will be replaced with 1's for typical RMSE
-#         calculation.
-
-#     Returns
-#     -------
-#     err : float
-#         Weighted RMSE between data in dpl and exp_dpl
-#     """
-#     from scipy import signal
-
-#     exp_times = exp_dpl.times
-#     sim_times = dpl.times
-
-#     # do tstart and tstop fall within both datasets?
-#     # if not, use the closest data point as the new tstop/tstart
-#     for tseries in [exp_times, sim_times]:
-#         if tstart < tseries[0]:
-#             tstart = tseries[0]
-#         if tstop > tseries[-1]:
-#             tstop = tseries[-1]
-
-#     # make sure start and end times are valid for both dipoles
-#     exp_start_index = (np.abs(exp_times - tstart)).argmin()
-#     exp_end_index = (np.abs(exp_times - tstop)).argmin()
-#     exp_length = exp_end_index - exp_start_index
-
-#     sim_start_index = (np.abs(sim_times - tstart)).argmin()
-#     sim_end_index = (np.abs(sim_times - tstop)).argmin()
-#     sim_length = sim_end_index - sim_start_index
-
-#     if weights is None:
-#         # weighted RMSE with weights of all 1's is equivalent to
-#         # normal RMSE
-#         weights = np.ones(len(sim_times[0:sim_end_index]))
-#     weights = weights[sim_start_index:sim_end_index]
-
-#     dpl1 = dpl.data["agg"][sim_start_index:sim_end_index]
-#     dpl2 = exp_dpl.data["agg"][exp_start_index:exp_end_index]
-#     if sim_length > exp_length:
-#         # downsample simulation timeseries to match exp data
-#         dpl1 = signal.resample(dpl1, exp_length)
-#         weights = signal.resample(weights, exp_length)
-#         indices = np.where(weights < 1e-4)
-#         weights[indices] = 0
-#     elif sim_length < exp_length:
-#         # downsample exp timeseries to match simulation data
-#         dpl2 = signal.resample(dpl2, sim_length)
-
-#     return np.sqrt((weights * ((dpl1 - dpl2) ** 2)).sum() / weights.sum())
-# # end of KDTODO
+def exp_decay(t, A, C, b):
+    return ((C - A) * np.exp(-b * (t))) + A
 
 
 def _anticorr(dpl, exp_dpl, tstart=0.0, tstop=0.0, weights=None):

hnn_core/hnn_io.py

@@ -381,6 +381,7 @@ def network_to_dict(net, write_output=False):
 
     net_data = {
         "object_type": "Network",
+        "model_variant": net._model_variant,
         "legacy_mode": net._legacy_mode,
         "N_pyr_x": net._N_pyr_x,
         "N_pyr_y": net._N_pyr_y,
@@ -505,7 +506,7 @@ def dict_to_network(net_data, read_drives=True, read_external_biases=True):
     params = dict()
     params["celsius"] = net_data["celsius"]
     params["threshold"] = net_data["threshold"]
-
+    params["model_variant"] = net_data.get("model_variant", None)
     mesh_shape = (net_data["N_pyr_x"], net_data["N_pyr_y"])
 
     # Instantiating network
@@ -581,6 +582,36 @@ def read_network_configuration(fname, read_drives=True, read_external_biases=Tru
             "type %s" % (net_data.get("object_type"))
         )
 
+    # ensure the cell types match the model variant
+    check_var = net_data.get("model_variant", None)
+    if check_var is not None and "duecker_ET_model".startswith(check_var):
+        missing_cells = [
+            cell_name
+            for cell_name in [
+                "L2_pyramidal",
+                "L5_pyramidal",
+                "L2_inhibitory",
+                "L5_inhibitory",
+            ]
+            if cell_name not in net_data["cell_types"]
+        ]
+        if missing_cells:
+            hint = ""
+            if all(
+                cell_name in net_data["cell_types"]
+                for cell_name in ["L2_basket", "L5_basket"]
+            ):
+                hint = (
+                    " The network has basket cells instead, so you are likely"
+                    " trying to create a duecker_ET_model with"
+                    " neymotin_2020_model cell types."
+                )
+            raise ValueError(
+                f"The cell types of the network do not match "
+                f"model_variant duecker_ET_model: no "
+                f"{', '.join(missing_cells)} found.{hint}"
+            )
+
     net = dict_to_network(net_data, read_drives, read_external_biases)
     _check_global_synaptic_gains_uniformity(net)
 

hnn_core/network.py

@@ -450,6 +450,7 @@ def __init__(
         # Save the parameters used to create the Network
         _validate_type(params, dict, "params")
         self._params = params
+        self._model_variant = params.get("model_variant", None)
         # Initialise a dictionary of cell ID's, which get used when the
         # network is constructed ('built') in NetworkBuilder
         # We want it to remain in each Network object, so that the user can

hnn_core/network_models.py

@@ -75,6 +75,92 @@
 }
 
 
+def _validate_params_for_model(
+    net,
+    params,
+    model_variant,
+    alt_variants=[],
+    require_variant=False,
+    excluded_cells=[],
+):
+    """Check that a param file matches the network model it is used for.
+
+    Parameters
+    ----------
+    net : Instance of Network object
+        The network the parameters are used for.
+    params : dict
+        The parameters the network was built from.
+    model_variant : str
+        Name of the network model, e.g. 'duecker_ET_model'. The
+        'model_variant' entry of `params` must be this name (or of one of
+        `alt_variants`).
+    alt_variants : list of str, default=[]
+        Further model names that are accepted in the 'model_variant' entry of `params`,
+        e.g. the deprecated name of a model. If `params` defines a local 'model_variant'
+        that is in 'alt_variants', the returned value will be the value of
+        'model_variant' that is passed to this function instead.
+    require_variant : bool, default=False
+        If True, raise if `params` does not define 'model_variant'. Used for
+        models that share no parameters with the default model, and would
+        otherwise silently fall back to default values.
+    excluded_cells : list of str, default=[]
+        Short names of cells that are *not* part of this network, e.g.
+        ('L2Basket', 'L5Basket') for a model in which basket cells are
+        replaced. Parameters for these cells are rejected.
+
+    Returns
+    -------
+    model_variant : str
+        The official model variant name.
+    """
+    check_var = params.get("model_variant", None)
+    if check_var is None:
+        if require_variant:
+            raise ValueError(
+                f"'model_variant' is required for simulations with "
+                f"{model_variant}. If you are sure that you are using the "
+                f"correct parameters, add 'model_variant': '{model_variant}', "
+                "to the first line of the param .json file."
+            )
+    elif check_var not in [model_variant] + alt_variants:
+        raise ValueError(
+            f"Parameters for {check_var} used for {model_variant}."
+            " Ensure that your param .json file matches the network "
+            f"and that your model variant is one of {[model_variant] + alt_variants}. "
+        )
+
+    # check that the params define the cell types of this network
+    missing_cells = [
+        cell_name
+        for cell_name in net.cell_types
+        if not any(_short_name(cell_name) in key for key in params)
+    ]
+    if missing_cells:
+        raise ValueError(
+            f"No parameters found for {', '.join(missing_cells)}."
+            " Ensure that your param .json file matches the network. "
+            " Reach out to us if this doesn't solve your problem.  "
+            " https://github.com/jonescompneurolab/hnn-core/discussions"
+        )
+
+    # check that the params don't define cell types this network replaced
+    unexpected_cells = [
+        cell_name
+        for cell_name in excluded_cells
+        if any(cell_name in key for key in params)
+    ]
+    if unexpected_cells:
+        raise ValueError(
+            f"Parameters found for {', '.join(unexpected_cells)}, which"
+            f" are not part of {model_variant}. Ensure that your"
+            " param .json file matches the network."
+            " Reach out to us if this doesn't solve your problem.  "
+            " https://github.com/jonescompneurolab/hnn-core/discussions"
+        )
+    return model_variant
+
+
 def neymotin_2020_model(
     params=None,
     add_drives_from_params=False,
@@ -131,9 +217,8 @@ def neymotin_2020_model(
     """
     hnn_core_root = Path(hnn_core.__file__).parent
     if params is None:
-        params = hnn_core_root / "param" / "default.json"
-    if isinstance(params, (str, Path)):
-        params = read_params(params)
+        params_fname = hnn_core_root / "param" / "default.json"
+        params = read_params(params_fname)
 
     # Define cell types for Jones 2009 model
     # data is here in metaData format
@@ -185,6 +270,9 @@ def neymotin_2020_model(
 
     delay = net.delay
 
+    # Ensure model_variant and params' cell types match current model
+    net._model_variant = _validate_params_for_model(net, params, "neymotin_2020_model")
+
     # source of synapse is always at soma
 
     # layer2 Pyr -> layer2 Pyr
@@ -412,12 +500,20 @@ def law_2021_model(
            Perception." Cerebral Cortex, 32, 668–688 (2022).
     """
 
-    net = jones_2009_model(
+    hnn_core_root = Path(hnn_core.__file__).parent
+    if params is None:
+        params_fname = hnn_core_root / "param" / "default.json"
+        params = read_params(params_fname)
+
+    net = neymotin_2020_model(
         params,
         add_drives_from_params,
         legacy_mode,
         mesh_shape=mesh_shape,
     )
+    # Ensure model_variant and params' cell types match current model (same cell types
+    # as 'neymotin_2020_model')
+    net._model_variant = _validate_params_for_model(net, params, "law_2021_model")
 
     # Update biophysics (increase gabab duration of inhibition)
     net.cell_types["L2_pyramidal"]["cell_object"].synapses["gabab"]["tau1"] = 45.0
@@ -500,8 +596,8 @@ def calcium_model(
            Brain Topography, 35, 19–35 (2022).
     """
     hnn_core_root = Path(hnn_core.__file__).parent
-    params_fname = hnn_core_root / "param" / "default.json"
     if params is None:
+        params_fname = hnn_core_root / "param" / "default.json"
         params = read_params(params_fname)
 
     net = jones_2009_model(
@@ -511,6 +607,10 @@ def calcium_model(
         mesh_shape=mesh_shape,
     )
 
+    # Ensure model_variant and params' cell types match current model (same cell types
+    # as 'neymotin_2020_model')
+    net._model_variant = _validate_params_for_model(net, params, "calcium_model")
+
     # Replace L5 pyramidal cell template with updated calcium
     cell_name = "L5_pyramidal"
     pos = net.cell_types[cell_name]["cell_object"].pos
@@ -527,13 +627,15 @@ def duecker_ET_model(
     """ "Initiate like old calcium model and then replace with new cells"""
 
     hnn_core_root = Path(hnn_core.__file__).parent
-    params_fname = hnn_core_root / "param" / "default_duecker_ET.json"
     if params is None:
+        params_fname = hnn_core_root / "param" / "default_duecker_ET.json"
         params = read_params(params_fname)
 
     cell_types = {
         "L2_inhibitory": {
-            "cell_object": human_gen_interneuron(cell_name="L2Inh", layer=2),
+            "cell_object": human_gen_interneuron(
+                cell_name=_short_name("L2_inhibitory"), layer=2
+            ),
             "cell_metadata": {
                 "morpho_type": "interneuron",
                 "electro_type": "inhibitory",
@@ -546,7 +648,7 @@ def duecker_ET_model(
             },
         },
         "L2_pyramidal": {
-            "cell_object": pyramidal_humanL23(cell_name="L2Pyr"),
+            "cell_object": pyramidal_humanL23(cell_name=_short_name("L2_pyramidal")),
             "cell_metadata": {
                 "morpho_type": "pyramidal",
                 "electro_type": "excitatory",
@@ -559,7 +661,9 @@ def duecker_ET_model(
             },
         },
         "L5_inhibitory": {
-            "cell_object": human_gen_interneuron(cell_name="L5Inh", layer=5),
+            "cell_object": human_gen_interneuron(
+                cell_name=_short_name("L5_inhibitory"), layer=5
+            ),
             "cell_metadata": {
                 "morpho_type": "interneuron",
                 "electro_type": "inhibitory",
@@ -572,7 +676,7 @@ def duecker_ET_model(
             },
         },
         "L5_pyramidal": {
-            "cell_object": pyramidal_humanL5ET(cell_name="L5Pyr"),
+            "cell_object": pyramidal_humanL5ET(cell_name=_short_name("L5_pyramidal")),
             "cell_metadata": {
                 "morpho_type": "pyramidal",
                 "electro_type": "excitatory",
@@ -613,6 +717,16 @@ def duecker_ET_model(
         cell_types=cell_types,
     )
 
+    # check variant and cell types. Basket cells are replaced by
+    # interneurons in duecker_ET_model, so their parameters are rejected
+    net._model_variant = _validate_params_for_model(
+        net,
+        params,
+        "duecker_ET_model",
+        require_variant=True,
+        excluded_cells=("L2Basket", "L5Basket"),
+    )
+
     delay = net.delay
 
     # layer2 Pyr -> layer2 Pyr

hnn_core/parallel_backends.py

@@ -80,9 +80,9 @@ def _gather_trial_data(sim_data, net, n_trials, postproc, bsl_cor="jones"):
 
         N_pyr_x = net._N_pyr_x
         N_pyr_y = net._N_pyr_y
-        if bsl_cor == "jones":
+        if bsl_cor == "jones" or bsl_cor == "neymotin":
             if net._verbose:
-                print("Applying Jones baseline correction", flush=True)
+                print("Applying Neymotin, 2020 baseline correction", flush=True)
             dpl._baseline_renormalize(N_pyr_x, N_pyr_y)  # XXX cf. #270
 
         dpl._convert_fAm_to_nAm()  # always applied, cf. #264

hnn_core/param/default_duecker_ET.json

@@ -1,4 +1,5 @@
-{   "sim_prefix": "duecker_ET_model",
+{   
+    "model_variant": "duecker_ET_model",
     "tstop": 170,
     "dt": 0.025,
     "celsius": 37.0,

hnn_core/param/duecker_ET_gamma_50Hz_bursty.json

@@ -1,4 +1,4 @@
-{   "sim_prefix": "duecker_ET_model",
+{   "model_variant": "duecker_ET_model",
     "tstop": 170,
     "dt": 0.025,
     "celsius": 37.0,