From e924db2b4ed2736afb27da6b7d1c5807d98aeb2b Mon Sep 17 00:00:00 2001 From: Johannes Mahl Date: Mon, 6 Jul 2026 16:03:40 -0700 Subject: [PATCH 1/4] update benchmark skill include variable initial guess fit recovery here rather than pytest because static initial guesses are meaningless and variable/random ones could break test pipeline in a non-reproducable way --- .claude/skills/benchmark/SKILL.md | 6 +- .claude/skills/benchmark/benchmark_gir.py | 171 +++++++++++++++++++++- TODO.md | 33 +++-- docs/ai/benchmark.md | 31 ++++ 4 files changed, 221 insertions(+), 20 deletions(-) diff --git a/.claude/skills/benchmark/SKILL.md b/.claude/skills/benchmark/SKILL.md index 981aabe..1d4ff4c 100644 --- a/.claude/skills/benchmark/SKILL.md +++ b/.claude/skills/benchmark/SKILL.md @@ -1,7 +1,7 @@ --- name: benchmark -description: 'Benchmark GIR compiled evaluator vs interpreter. Args: example number (default: 2), `--fit` for full fit, `-n N` for repetitions, `--nfev` for residual-evaluation counts, `--plan-time` for planning vs fit wall time.' -argument-hint: '[example_num] [--fit] [-n N] [--nfev] [--plan-time]' +description: 'Benchmark GIR compiled evaluator vs interpreter. Args: example number (default: 2), `--fit` for full fit, `-n N` for repetitions, `--nfev` for residual-evaluation counts, `--plan-time` for planning vs fit wall time, `--par-variability` for fit robustness vs perturbed initial guesses (`--starts N`).' +argument-hint: '[example_num] [--fit] [-n N] [--nfev] [--plan-time] [--par-variability] [--starts N]' disable-model-invocation: true allowed-tools: Bash(.venv/bin/python .claude/skills/benchmark/benchmark_gir.py *) --- @@ -11,6 +11,6 @@ source of truth for this skill. Pass through the user's optional benchmark arguments unchanged: -- `[example_num] [--fit] [-n N] [--nfev] [--plan-time]` +- `[example_num] [--fit] [-n N] [--nfev] [--plan-time] [--par-variability] [--starts N]` When this wrapper and the shared doc differ, follow the shared doc. diff --git a/.claude/skills/benchmark/benchmark_gir.py b/.claude/skills/benchmark/benchmark_gir.py index a23f127..9239341 100644 --- a/.claude/skills/benchmark/benchmark_gir.py +++ b/.claude/skills/benchmark/benchmark_gir.py @@ -427,6 +427,125 @@ def timed_sched(*args, **kwargs): return plan_total, fit_time +# ------------------------------------------------------------------ +# Parameter variability (fit robustness vs. perturbed initial guesses) +# ------------------------------------------------------------------ + + +# +def capture_par_variability(example_num, *, n_starts=4): + """Run the standard pipeline from perturbed initial guesses. + + Runs the baseline + fit_2d pipeline once from the example's own + initial values (reference) and ``n_starts`` more times with every + free parameter's init scaled by a fixed factor ladder, clipped + into bounds. Deterministic by construction -- no RNG -- so + repeated invocations give identical numbers. + + This is a diagnostic (robustness report), not a pass/fail test: + a parameter whose fitted value depends on its start indicates a + shallow or multi-modal objective, or init-dependent machinery + (e.g. conv-kernel support sized from the initial value). + """ + + ladder = [0.6, 0.8, 1.25, 1.6] + fitted_runs = [] + labels = [] + redchis = [] + + for run in range(n_starts + 1): + file, dynamics_calls = load_example(example_num, add_dynamics=False) + file.define_baseline( + time_start=0, time_stop=10, time_type="ind", show_plot=False + ) + file.fit_baseline(model_name="2D", stages=2, try_ci=0) + for call in dynamics_calls: + file.add_time_dependence(**call) + + model = file.model_active + assert model is not None + free_names = [ + name for name, par in model.lmfit_pars.items() if par.vary and not par.expr + ] + + label = "reference" if run == 0 else f"start {run}" + if run > 0: + for i, name in enumerate(free_names): + par = model.lmfit_pars[name] + factor = ladder[(i + run - 1) % len(ladder)] + value = par.value * factor + if abs(par.value) < 1e-12: + # multiplicative perturbation is a no-op at 0; nudge + # by a fraction of the bound range when available + if np.isfinite(par.min) and np.isfinite(par.max): + value = par.value + 0.1 * (par.max - par.min) + else: + continue + if np.isfinite(par.min) and np.isfinite(par.max): + span = par.max - par.min + value = min( + max(value, par.min + 0.01 * span), + par.max - 0.01 * span, + ) + elif value <= par.min or value >= par.max: + continue # would leave bounds; keep original init + par.value = value + + t0 = time.perf_counter() + file.fit_2d(model_name="2D", stages=2, try_ci=0) + wall = time.perf_counter() - t0 + + df = file.get_fit_results(fit_type="2d") + fitted = dict(zip(df["name"], df["value"], strict=True)) + fitted_runs.append({name: fitted[name] for name in free_names}) + + try: + redchi = model.result[1].redchi + except (AttributeError, IndexError, TypeError): + redchi = float("nan") + labels.append(label) + redchis.append(redchi) + print(f" {label:12s} redchi={redchi:10.4g} wall={wall:6.2f} s") + + # Separate the two failure signals: a run that converged to a worse + # optimum (secondary minimum) vs. spread among runs that reached the + # best optimum (identifiability / init-dependent machinery). + redchi_arr = np.array(redchis) + best_redchi = np.nanmin(redchi_arr) + ok_idx = [ + i + for i in range(len(redchis)) + if np.isfinite(redchi_arr[i]) and redchi_arr[i] <= best_redchi * 1.01 + ] + off_idx = [i for i in range(len(redchis)) if i not in ok_idx] + + print() + if off_idx: + off_str = ", ".join(f"{labels[i]} (redchi {redchis[i]:.4g})" for i in off_idx) + print(f" off-optimum runs (excluded from spread): {off_str}") + print(f" {'parameter':32s}{'reference':>12s}{'spread':>12s}{'rel':>8s}") + worst_rel = 0.0 + worst_name = "-" + for name in fitted_runs[0]: + values = np.array([fitted_runs[i][name] for i in ok_idx]) + spread = float(values.max() - values.min()) + scale = max(abs(float(values.mean())), 1e-300) + rel = spread / scale + if rel > worst_rel: + worst_rel, worst_name = rel, name + flag = " <-- start-sensitive" if rel > 0.01 else "" + print( + f" {name:32s}{fitted_runs[0][name]:12.6g}{spread:12.3g}" + f"{rel * 100:7.2f}%{flag}" + ) + print() + print( + f" {len(ok_idx)}/{len(redchis)} runs at best optimum; " + f"worst spread among them: {worst_name} ({worst_rel * 100:.2f}%)" + ) + return worst_rel, worst_name, len(off_idx) + + # ------------------------------------------------------------------ # Full-fit benchmark # ------------------------------------------------------------------ @@ -536,11 +655,29 @@ def bench_fit(example_num, dynamics_calls, *, n_reps=3): "wall time. If --example is 0, runs all examples." ), ) + parser.add_argument( + "--par-variability", + action="store_true", + help=( + "Fit from deterministically perturbed initial guesses and " + "report the spread of fitted values across starts (no RNG, " + "reproducible). --starts controls the number of perturbed " + "starts. If --example is 0, runs all examples." + ), + ) + parser.add_argument( + "--starts", + type=int, + default=4, + help="Perturbed starts for --par-variability (default: 4)", + ) args = parser.parse_args() - # --nfev / --plan-time with example 0 means "all examples"; skip the + # Multi-example modes with example 0 mean "all examples"; skip the # single-example preamble so the capture fn can load each one itself. - skip_preamble = (args.nfev or args.plan_time) and args.example == 0 + skip_preamble = ( + args.nfev or args.plan_time or args.par_variability + ) and args.example == 0 if not skip_preamble: folder = _find_example_folder(args.example) print(f"Example: {folder.name}") @@ -615,6 +752,36 @@ def bench_fit(example_num, dynamics_calls, *, n_reps=3): print(f"PLAN-TIME CAPTURE -- example {args.example:02d}") print("=" * 60) capture_plan_time(args.example) + elif args.par_variability: + if args.example == 0: + summaries: dict[int, str] = {} + for n in range(1, 5): + print() + print("=" * 60) + print(f"PARAMETER VARIABILITY -- example {n:02d}") + print("=" * 60) + try: + worst_rel, worst_name, n_off = capture_par_variability( + n, n_starts=args.starts + ) + summaries[n] = ( + f"worst {worst_rel * 100:6.2f}% ({worst_name})" + f"{f' [{n_off} off-optimum]' if n_off else ''}" + ) + except FileNotFoundError as e: + print(f" skipped: {e}") + summaries[n] = "skipped" + print() + print("=" * 60) + print("PARAMETER VARIABILITY SUMMARY") + for n, summary_str in summaries.items(): + print(f" example {n:02d}: {summary_str}") + print("=" * 60) + else: + print("=" * 60) + print(f"PARAMETER VARIABILITY -- example {args.example:02d}") + print("=" * 60) + capture_par_variability(args.example, n_starts=args.starts) else: bench_per_call(file, n_calls=args.calls) diff --git a/TODO.md b/TODO.md index 806b3c4..e6c3c22 100644 --- a/TODO.md +++ b/TODO.md @@ -2,15 +2,15 @@ ## Fitting -- [ ] **Mismatched initial guesses**: round-trip tests — one each for basic, profile, profile+dynamics. -- [ ] **Conv kernel support is sized from the initial parameter value**: `Component.create_t_kernel` builds the kernel time axis once at model construction (`t_range = par_init * kernel_width`, e.g. ±4·SD for `gaussCONV`) and never rebuilds it. When the fitted width grows past its init, the kernel is silently truncated, biasing the recovered width (observed while building `04_parameter_profiles`, 2026-06-11: truth SD=10, init 5 → fitted SD≈10.8 at any SNR; 03's SD 0.148 vs truth 0.15 with init 0.1 is likely the same effect, mild). The GIR path snapshots the same static axis (`graph_ir.py`, `kernel_time`). Fix candidates: rebuild the kernel axis when the kernel parameter value changes; size the support from the parameter's max bound; or at minimum warn when `fitted_value * kernel_width` exceeds the kernel range. Workaround used in the examples: initialize conv widths generously above the expected value (commented in the YAMLs). +- [ ] **Enable 1D time-trace fitting (post-SbS kinetics)**: wire standalone `TIME_1D` dynamics models (already evaluable on the mcp path per `docs/design/supported_models.md`, but connected to no fit method) to a fit entry point, so parameter-vs-time traces extracted from SbS results can be fit to `functions/time.py` dynamics (`expFun` sums, IRF convolution) inside the package instead of in external scripts. Frame it as the diagnostic/initialization rung before `fit_2d`, not a statistical equivalent (per-slice correlations and unpropagated SbS uncertainties make two-step inferior). Design direction: promote SbS results into a time-axis `File` — the promoted object inherits limits/CI/MCMC/save/export machinery, and trace fits land in a `SavedFitSlot` like every other fit type (no parallel fitting pipeline). Requires `File` to support time as the primary axis, which today it assumes is energy. Weighting traces by per-slice stderr ties into the `sigma_type` expansion item below. +- [ ] **Conv kernel support is sized from the initial parameter value**: `Component.create_t_kernel` builds the kernel time axis once at model construction (`t_range = par_init * kernel_width`, e.g. ±4·SD for `gaussCONV`) and never rebuilds it. When the fitted width grows past its init, the kernel is silently truncated, biasing the recovered width (observed while building `04_parameter_profiles`, 2026-06-11: truth SD=10, init 5 → fitted SD≈10.8 at any SNR; 03's SD 0.148 vs truth 0.15 with init 0.1 is likely the same effect, mild). The GIR path snapshots the same static axis (`graph_ir.py`, `kernel_time`). Fix candidates: rebuild the kernel axis when the kernel parameter value changes; size the support from the parameter's max bound; or at minimum warn when `fitted_value * kernel_width` exceeds the kernel range. Workaround used in the examples: initialize conv widths generously above the expected value (commented in the YAMLs). When fixing, add a deterministic regression test that inits the conv width well below truth (the failure mode above); broad robustness-vs-start checking deliberately lives in the benchmark skill (`--par-variability`, added 2026-07-06), not the test suite, to avoid weak-or-flaky convergence asserts. ## Noise and simulation - [ ] **Simulator noise-language cleanup**: align simulator docs/metadata with the fit-results noise schema. Keep simulator `noise_type` meaning "noise distribution / random generator" (`gaussian`, `poisson`, `none`), not sigma shape. Fix the stale `Simulator.set_noise_type()` docstring that mentions `uniform`; clarify `detection` vs. `noise_type` vs. `noise_level`; and, for analog Gaussian simulations, consider saving the derived `sigma_data = noise_level * max(abs(clean_data))` alongside existing metadata. For parameter sweeps, store derived `sigma_data` per configuration when it depends on each clean dataset. - [ ] **Warn on subcycle-boundary time samples**: in `Dynamics.normalize_time`, emit a warning when a time sample lands within epsilon of a subcycle boundary (`|t * f_eff - round(t * f_eff)| < eps`). Subcycle masks switch discretely there, so the sample's assignment — and hence the model prediction for that whole row — is sensitive to floating-point representation of the time axis. This silently biased a multi-cycle fit when synthetic data was generated on `np.arange` axes but fit against their `%.6e`-rounded CSV reload (see `03_multi_cycle_dynamics/data/generate_data.ipynb`, 2026-06-11). A warning turns that silent bias into a one-line diagnosis; it also flags the physically ambiguous case of measuring exactly at the switching instant. -- [ ] **Future `sigma_type` expansion in FitResults**: after the first constant, user-supplied sigma schema lands, extend uncertainty handling beyond scalar `sigma_data`. Keep `noise_type` for the statistical assumption/distribution and use `sigma_type` for sigma shape: initially `constant`, later `per_spectrum` and `per_point`. Add HDF5 storage, validation, baseline/SBS/2D alignment, `compare_models()` behavior, and tests for vector/matrix sigma. Defer automatic Poisson-derived sigma until residual-space variance propagation is explicit. +- [ ] **Future `sigma_type` expansion in FitResults**: the constant, user-supplied sigma schema has landed (`SIGMA_TYPE_CONSTANT` in [fit_io.py](src/trspecfit/utils/fit_io.py) ~L54; `validate_noise_metadata` hard-locks `sigma_type` to `"constant"`), so this is now unblocked: extend uncertainty handling beyond scalar `sigma_data`. Keep `noise_type` for the statistical assumption/distribution and use `sigma_type` for sigma shape: initially `constant`, later `per_spectrum` and `per_point`. Add HDF5 storage, validation, baseline/SBS/2D alignment, `compare_models()` behavior, and tests for vector/matrix sigma. Defer automatic Poisson-derived sigma until residual-space variance propagation is explicit. ## Performance & architecture @@ -21,28 +21,31 @@ - **Persist `correl` and `acceptance_fraction` into the slots**: 2026-06 added live-only accessors (`get_correlations`, `get_conf_intervals`, `get_mcmc`) reading `model.result` as a stopgap for notebook 12, so these are NOT yet saved. Add per-parameter correlations to the slot `params` payload and `acceptance_fraction` to the slot `mcmc` payload, with `.fit.h5` read/write support and save/load round-trip tests, so they survive persistence like the rest of the slot. - **Relocate the live accessors to `FitResults`**: 2026-06 added `File.get_correlations`, `File.get_conf_intervals`, `File.get_mcmc` (and the private `File._result_model` resolver) reading `model.result[...]` directly. These conceptually belong on `FitResults` (like `compare_models`, which already lives there with `File.compare_models` as thin sugar). The existing `File.get_fit_results` is in the same boat. Decide whether all of these should move into `FitResults` (with thin `File.*` sugar that delegates), and whether they read live `model.result` or persisted slots — then move them and update callers (notebook 12 reads them). - The raw list-index access (`result[1..4]`) and the deeper unified-results-object question are deferred to this item. -- [ ] **Disentangle plotting from saving/conversion in the fit pipeline**: figure rendering is currently entangled with data conversion and file IO. `fitlib.results_to_df` (results → DataFrame) and `File._save_2d_fit_legacy` / `_save_sbs_fit_legacy` (CSV writers) also render figures, and `fit_baseline` / `fit_slice_by_slice` / `fit_2d` reach plotting only by calling that save-legacy path. Split into (a) compute/convert, (b) explicit save/export, (c) an explicit plotting API — cf. the `_save_img_flag` helper and `FitResults.plot_residuals`. Until then, the fit methods gate display on `show_output` and saving on `auto_export` via a `save_files` flag through the legacy methods. Scope: beyond the examples-upgrade branch. +- [ ] **Disentangle plotting from saving/conversion in the fit pipeline**: figure rendering is currently entangled with data conversion and file IO. `fitlib.results_to_df` (results → DataFrame) and `File._save_2d_fit_legacy` / `_save_sbs_fit_legacy` (CSV writers) also render figures, and `fit_slice_by_slice` / `fit_2d` reach plotting only by calling that save-legacy path (`fit_baseline` is already disentangled — it calls `fitlib.plt_fit_res_1d` directly, gated by `_save_img_flag`; use it as the template). Blocks the legacy-shim removal item in the v1.0.0 checklist, since the `_save_*_legacy` impls are the live SbS/2D plotting path. Split into (a) compute/convert, (b) explicit save/export, (c) an explicit plotting API — cf. the `_save_img_flag` helper and `FitResults.plot_residuals`. Until then, the fit methods gate display on `show_output` and saving on `auto_export` via a `save_files` flag through the legacy methods. Scope: beyond the examples-upgrade branch. - [ ] **Decouple directory creation from path computation (`create_model_path`)**: `Project.create_model_path` ([trspecfit.py](src/trspecfit/trspecfit.py) ~L2167) `mkdir`s the `{path_results}/{file}/{fit_type}/{model}/` tree as a side effect of *computing* the path. Every fit method calls it to build the `save_path` handed to `fit_wrapper` (`fit_baseline` ~L2572, `fit_2d` ~L3983, spectrum ~L2804, sbs ~L3120), and the project plot-grid display builds image paths through it too (`fit_baselines` ~L1070, `fit_2d` ~L1432) — so empty `*_fits/` dir trees appear even when `auto_export: False` and nothing is written (observed in 21, 2026-06-27; every fitting example already does this). `fit_wrapper` only writes when `save_output=1 if auto_export else 0`. Fix: make `create_model_path` return the path *without* `mkdir`, and `mkdir`-on-write at the actual write sites (`fit_wrapper`'s save branch + the explicit `save_*`/`export_*` functions); the plot-grid readers then just build a path and `.exists()`-check. Do **not** gate on `auto_export` directly — explicit `save_fit()`/`export_fit()` legitimately need the dir and aren't `auto_export`-driven. Cosmetic only (the dirs are gitignored, no functional bug); blast radius spans baseline/spectrum/sbs/2d + savers, needs no-regression tests. Related to the plotting/saving disentanglement item above; scope: beyond the examples-upgrade branch. ## Testing -- [ ] **Pylance/Pyright `| None` noise in tests**: ~280 Pyright errors across tests, all from accessing `File`/`Model` attributes typed as `ndarray | None`. Current `# type guard` asserts are inconsistent and don't propagate through helper methods. Find a cleaner pattern (e.g. `TypeGuard`, narrowing wrapper, or Pyright config) and apply consistently. +- [ ] **Pylance/Pyright `| None` noise in tests**: several hundred Pyright errors across tests, mostly from accessing `File`/`Model` attributes typed as `ndarray | None` (`[tool.pyright]` only includes `src`, so these surface in editors, not CI). Current `# type guard` asserts (~110) are inconsistent and don't propagate through helper methods. Find a cleaner pattern (e.g. `TypeGuard`, narrowing wrapper, or Pyright config) and apply consistently. ## User and AI ergonomics -- [ ] **Curate the public API surface before v1.0.0**: audit user-facing classes (`File`, `Project`, `Simulator`, `Model`, etc.) and decide which methods/attributes should be discoverable in notebooks and docs. Add curated `__dir__()` output for autocomplete, keep `__all__`/API docs aligned, and gradually rename or deprecate internal helper methods that should not look like primary user workflows. This should improve both human notebook ergonomics and AI/LLM efficiency by making the intended workflow surface smaller, clearer, and easier to infer. -- [ ] **Document API tiers**: add a short guide that separates stable user API (`Project`, `File`, `Simulator`, `PlotConfig`), advanced public API (`mcp.Model`, `Component`, `Par`, `ParameterSweep`, `MC`), and internal implementation modules (`graph_ir`, `eval_1d`, `eval_2d`, low-level parsing/HDF5 helpers). Use this as the source of truth for docs, tests, examples, and AI-agent guidance. - [ ] **Add tool-neutral agent orientation**: add `AGENTS.md` or `docs/ai/agent-orientation.md` pointing agents to `CLAUDE.md`, `TODO.md`, `PLAN.md`, `docs/design/repo_architecture.md`, supported-model docs, common commands, and API-change guardrails. Keep it concise so any LLM can quickly find the intended workflow and repo boundaries. -- [ ] **Add more AI-friendly task recipes**: extend `docs/ai/` with checklists for common repo changes, such as adding YAML syntax, adding plotting options, changing fitting workflows, modifying GIR/evaluator behavior, extending save/load fields, and preparing a release. -- [ ] **`trspec.fit` project website**: stand up the site as the home for fully rendered example notebooks (plots, fit tables) so users can browse without installing — preferred over Read the Docs, which has build timeouts/execution limits for notebooks. Room to grow into interactive tooling: a model-builder UI (GUI that outputs the YAML model files) and other future work. Keep git notebook sources stripped (`nbstripout`); render/execute at publish time so outputs never drift from the code. +- [ ] **Add more AI-friendly task recipes**: `docs/ai/` already covers add-function, check-example, check-docs, code-review, changelog, and benchmark. Extend with checklists for the remaining common changes: adding YAML syntax, adding plotting options, changing fitting workflows, modifying GIR/evaluator behavior, extending save/load fields, and preparing a release. - [ ] **Revisit default `Project.name` for saving examples**: `Project.name` defaults to `"test"` ([trspecfit.py](src/trspecfit/trspecfit.py) ~L202), so a bare `save_fit()` / `save_fits()` with no path writes `./fit_results/test.fit.h5`. The example `project.yaml` files set no `name:`. When working through the saving notebooks (`11_save_load_export`, `20_multi_file_independent_fit`), decide on a convention — set a meaningful `name:` per example so the default archive path is self-describing, and/or keep using explicit content-named paths (`comparison.fit.h5`, `batch.fit.h5`). Pick one and apply consistently. - [ ] **Add minimal runnable workflow examples**: supplement notebooks with small script-like examples or docs snippets for the canonical public workflows: load data, load a model, set limits, fit baseline, fit 2D, inspect results, simulate data, and run a parameter sweep. -- [ ] **Improve public validation errors**: make user-facing errors state what failed, where it failed (file/model/component/parameter when applicable), and what the user or agent should change next. Prioritize YAML parsing, model loading, fit setup, and unsupported-model fallback paths. -- [ ] **Tighten public type hints and aliases**: reduce ambiguous `Any` on public APIs, document key aliases such as `ModelRef`, and keep return types crisp for IDEs, Pyright, generated docs, and LLM code navigation. +- [ ] **Improve public validation errors**: YAML parsing and model loading already meet the bar — `ModelValidationError` messages carry model/component/parameter context plus hints (verified 2026-07-06). Remaining scope: audit fit setup and unsupported-model fallback paths and bring them to the same standard (state what failed, where, and what to change next). - [ ] **Document fast verification slices**: add focused pytest commands for common edits (public workflow/API, YAML parser, functions, GIR/evaluator, plotting) so contributors and agents can validate changes quickly before running the full suite. -## Build & release +## Road to v1.0.0 & first adopters -- [ ] **Automate tagging and pushing**: automate `git tag v1.2.3` + `git push v1.2.3` as part of the release workflow. -- [ ] **Remove legacy/backwards-compat code**: before v1.0.0 release, audit codebase for legacy fallbacks and backwards compatibility shims and consider removing. Known shims slated for removal: - - `File.save_sbs_fit` / `File.save_2d_fit` (deprecated wrappers — replace internal callers and drop the public methods plus their `_save_*_fit_legacy` impls; users should migrate to `File.export_fit(fit_type=...)`). +Rationale (2026-07-06): lock-in comes from schemas and public names, not from users — so stabilize the contract surfaces (1–5) before the website (6) and outreach (tracked outside the repo: `~/Desktop/trspecfit_outreach.md`). Early TR-XPS/TR-XAS adopters are then a feature, not a risk: they supply the usage signal that feature decisions currently lack, and their needs are a strict subset of the general "parameterized lineshapes vs. control axis" problem, so serving them does not distort the architecture. + +1. [ ] **Version the model YAML schema**: model YAML files are the artifact user groups accumulate, and currently the only unversioned contract — the fit archive already carries `SCHEMA_VERSION` ([fit_io.py](src/trspecfit/utils/fit_io.py) ~L46) and refuses mismatched reads/appends. Add an optional top-level format-version key in `utils/parsing.py` (absent = current version) so future syntax changes can warn or migrate instead of silently misparsing old model files. +2. [ ] **Curate the public API surface**: audit user-facing classes (`File`, `Project`, `Simulator`, `Model`, etc.) and decide which methods/attributes should be discoverable in notebooks and docs. Add curated `__dir__()` output for autocomplete, keep `__all__`/API docs aligned, and gradually rename or deprecate internal helper methods that should not look like primary user workflows. This should improve both human notebook ergonomics and AI/LLM efficiency by making the intended workflow surface smaller, clearer, and easier to infer. Prerequisite for outreach: once external groups have notebooks, every public-looking name is frozen in practice. +3. [ ] **Document API tiers**: add a short guide that separates stable user API (`Project`, `File`, `Simulator`, `PlotConfig`), advanced public API (`mcp.Model`, `Component`, `Par`, `ParameterSweep`, `MC`), and internal implementation modules (`graph_ir`, `eval_1d`, `eval_2d`, low-level parsing/HDF5 helpers). Use this as the source of truth for docs, tests, examples, and AI-agent guidance. +4. [ ] **Remove legacy/backwards-compat code**: audit codebase for legacy fallbacks and backwards-compatibility shims and consider removing before v1.0.0. Full removal depends on the plotting/saving disentanglement item (Performance & architecture) — the `_save_*_fit_legacy` impls are still the live SbS/2D plotting path. Known shims slated for removal: + - `File.save_sbs_fit` / `File.save_2d_fit` (deprecated wrappers — replace internal callers and drop the public methods plus their `_save_*_fit_legacy` impls; users should migrate to `File.export_fit(fit_type=...)`). + - Legacy pre-rename format branch in [sweep.py](src/trspecfit/utils/sweep.py) ~L624 and the legacy `save_img` int-mapping helper in [plot.py](src/trspecfit/utils/plot.py) ~L813. +5. [ ] **State a pre-1.0 stability & deprecation policy**: short docs note establishing the social contract for early adopters — pre-1.0 public API may change, but only with a `DeprecationWarning` and a changelog entry; versioned schemas (`.fit.h5` archive, model YAML) refuse or migrate rather than misread. Cheap to write, and it is what makes early adoption safe for both sides. +6. [ ] **`trspec.fit` project website**: stand up the site as the home for fully rendered example notebooks (plots, fit tables) so users can browse without installing — preferred over Read the Docs, which has build timeouts/execution limits for notebooks. Room to grow into interactive tooling: a model-builder UI (GUI that outputs the YAML model files) and other future work. Keep git notebook sources stripped (`nbstripout`); render/execute at publish time so outputs never drift from the code. diff --git a/docs/ai/benchmark.md b/docs/ai/benchmark.md index 1378bef..1d2163c 100644 --- a/docs/ai/benchmark.md +++ b/docs/ai/benchmark.md @@ -41,6 +41,8 @@ Parse the arguments: - First positional integer -> `--example N` (default: `2`) - `--fit` -> include full-fit benchmark - `-n N` -> fit repetitions (default: `3`) +- `--par-variability` -> fit-robustness report (see below); `--starts N` sets + the number of perturbed starts (default: `4`) Run: @@ -72,6 +74,35 @@ table at the end. .venv/bin/python .claude/skills/benchmark/benchmark_gir.py --example --plan-time ``` +## Parameter variability (fit robustness vs. initial guesses) + +`--par-variability` runs the baseline + `fit_2d` pipeline once from the +example's own initial values (reference) and `--starts` more times (default 4) +with every free parameter's init scaled by a fixed factor ladder +(x0.6 / x0.8 / x1.25 / x1.6, clipped into bounds). No RNG — repeated +invocations give identical numbers. + +The report separates the two failure signals: + +- **off-optimum runs** — a start that converged to a worse redchi (secondary + minimum). Excluded from the spread statistics and listed separately. +- **spread among best-optimum runs** — fitted-value spread across starts that + reached the same optimum. Nonzero spread here indicates a flat objective + direction (the parameter is not identifiable from the data) or + init-dependent machinery: state derived once from initial parameter values + and never rebuilt during the fit, such as a convolution-kernel support. + Parameters above 1% relative spread are flagged `start-sensitive`. + +This is a diagnostic, not a pass/fail test — it deliberately lives here rather +than in the pytest suite, where perturbed-start convergence asserts would be +either weak (tiny perturbations) or flaky (aggressive ones). Accepts +`--example 0` for an all-examples summary. + +```bash +.venv/bin/python .claude/skills/benchmark/benchmark_gir.py --example --par-variability +.venv/bin/python .claude/skills/benchmark/benchmark_gir.py --example 0 --par-variability --starts 6 +``` + ## Profiling (GIR path only) For flamegraphs of the GIR hot path, use `--profile` to run a GIR-only loop From 17c966bb3e4fc3c43143d9dd171537b8ee62b13d Mon Sep 17 00:00:00 2001 From: Johannes Mahl Date: Wed, 8 Jul 2026 11:59:09 -0700 Subject: [PATCH 2/4] fix Voigt amplitude dependence on energy grid --- src/trspecfit/functions/energy.py | 7 ++++--- tests/test_functions_energy.py | 17 +++++++++++++++++ 2 files changed, 21 insertions(+), 3 deletions(-) diff --git a/src/trspecfit/functions/energy.py b/src/trspecfit/functions/energy.py index 94f7554..5e9b330 100644 --- a/src/trspecfit/functions/energy.py +++ b/src/trspecfit/functions/energy.py @@ -255,7 +255,7 @@ def Voigt(x: np.ndarray, A: float, x0: float, SD: float, W: float) -> np.ndarray x : ndarray Energy axis A : float - Peak amplitude (maximum value, approximately, for narrow peaks) + Peak amplitude (maximum value, at x = x0) x0 : float Peak center position SD : float @@ -270,8 +270,9 @@ def Voigt(x: np.ndarray, A: float, x0: float, SD: float, W: float) -> np.ndarray """ voigt = np.real(wofz(((x - x0) + 1j * (W / 2)) / SD / np.sqrt(2))) - max_voigt = np.max(voigt, axis=-1, keepdims=True) - return np.asarray(A * voigt / max_voigt) + # analytic peak value (profile at x = x0), keeps A grid-independent + peak_voigt = np.real(wofz(1j * (W / 2) / SD / np.sqrt(2))) + return np.asarray(A * voigt / peak_voigt) # diff --git a/tests/test_functions_energy.py b/tests/test_functions_energy.py index 2e03a15..9147871 100644 --- a/tests/test_functions_energy.py +++ b/tests/test_functions_energy.py @@ -368,6 +368,23 @@ def test_2d_broadcast_normalizes_per_slice(self): assert result.shape == (2, x.shape[-1]) np.testing.assert_allclose(np.max(result, axis=-1), A[:, 0], rtol=1e-3) + # + def test_grid_independent_amplitude(self): + """Value at a fixed x must not depend on the sampled window/grid. + + The analytic peak normalization (wofz at dx=0) replaced the old + max-over-grid normalization, which rescaled the whole profile when + x0 fell outside the energy window. + """ + + full_x = np.linspace(-10, 10, 2001) # includes the peak at x0=0 + idx_probe = np.argmin(np.abs(full_x - 2.0)) + x_probe = full_x[idx_probe] + tail_x = np.linspace(x_probe, 10, 18) # excludes the peak entirely + full = Voigt(full_x, A=4.0, x0=0.0, SD=1.0, W=1.0) + tail = Voigt(tail_x, A=4.0, x0=0.0, SD=1.0, W=1.0) + assert tail[0] == pytest.approx(full[idx_probe], rel=1e-9) + # # From f234568e2657214f8e64080ae43e8f615738278d Mon Sep 17 00:00:00 2001 From: Johannes Mahl Date: Wed, 8 Jul 2026 12:48:32 -0700 Subject: [PATCH 3/4] remove y0 from dynamics functions; add stepFun offset primitive Breaking change. The per-function y0 offset broke the causality convention f(t < t0) = 0 inconsistently across dynamics functions (sqrtFun/erfFun leaked y0 before t0; the rest clamped to 0) and duplicated what an additive component already expresses. Drop y0 from linFun, expFun, sinFun, sinDivX, erfFun, and sqrtFun, and add stepFun(t, A, t0) as the dedicated causal offset primitive (erfFun is its Gaussian-broadened counterpart). Offsets and plateaus are now modeled by adding a stepFun component, e.g. expFun + stepFun sharing t0 for a decay to a nonzero plateau. Wire stepFun through both evaluation paths (mcp and the lowered GIR/eval_2d dispatch, new DynFuncKind.STEPFUN) and migrate every model YAML, notebook, and doc that referenced dynamics y0 (all values were 0, so pure deletions). Stale y0 entries now fail model validation with a parameter-count error. The energy-domain Offset(x, y0) background is unchanged. Also from the functions audit: - correct docstrings: expRiseCONV (anti-causal, not causal), boxCONV (hard edges, no smoothing), GLS/GLP (valid range m in [0, 1]; GLP gives NaN/inf for m < 0), DS (asymmetric tail extends toward x < x0). - add end-to-end lowered-2D parity tests for every dynamics driver (stepFun, sinFun, linFun, sinDivX, erfFun, sqrtFun), verifying each DynFuncKind enum/dispatch entry and param count against the mcp interpreter -- previously only expFun was exercised this way, and this change decremented all those param counts. --- CHANGELOG.md | 18 ++ PLAN.md | 76 ++++++- docs/ai/add-function.md | 2 +- docs/design/lowered_evaluator.md | 31 ++- docs/design/repo_architecture.md | 6 +- .../01_basic_fitting/data/generate_data.ipynb | 2 +- .../data/models_time_truth.yaml | 1 - .../01_basic_fitting/models_time.yaml | 1 - .../02_dependent_parameters/models_time.yaml | 1 - .../data/models_time_truth.yaml | 2 - .../03_multi_cycle_dynamics/models_time.yaml | 4 +- .../data/models_time_truth.yaml | 1 - .../04_parameter_profiles/models_time.yaml | 1 - .../10_model_comparison/models_time.yaml | 2 - .../models_time.yaml | 1 - .../data/generate_data.ipynb | 3 +- .../data/models_time_truth.yaml | 1 - .../21_multi_file_shared_fit/models_time.yaml | 1 - .../01_simulator/models_time.yaml | 7 +- .../02_ml_training_data/models_time.yaml | 7 +- pyproject.toml | 2 +- src/trspecfit/eval_2d.py | 13 +- src/trspecfit/functions/energy.py | 6 + src/trspecfit/functions/time.py | 115 ++++++---- src/trspecfit/graph_ir.py | 3 + tests/models/file_time.yaml | 60 +++-- tests/models/project_time.yaml | 3 - tests/test_config_functions.py | 7 +- tests/test_evaluate_2d.py | 43 ++++ tests/test_functions_time.py | 207 ++++++++++-------- tests/test_graph_ir.py | 38 +--- tests/test_mcp_library.py | 3 - tests/test_model_parser.py | 4 - tests/test_project_fit.py | 4 - 34 files changed, 420 insertions(+), 256 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 48a85c5..47792b7 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,24 @@ and this project adheres to [Semantic Versioning](https://semver.org/). This file is maintained using the shared changelog workflow in [`docs/ai/changelog.md`](docs/ai/changelog.md). +## [0.10.0] - 2026-07-08 + +### Added + +- **`stepFun(t, A, t0)` dynamics primitive**: a causal step (0 before `t0`, `A` after). Since dynamics components combine by addition, `stepFun` is the dedicated way to model baselines and plateaus (e.g. `expFun` + `stepFun` sharing `t0` gives a decay to a nonzero plateau); `erfFun` is its Gaussian-broadened counterpart. Registered on both the mcp and compiled (GIR) evaluation paths. + +### Changed + +- **Breaking: `y0` removed from all dynamics functions** (`linFun`, `expFun`, `sinFun`, `sinDivX`, `erfFun`, `sqrtFun`). The per-function offset broke the causality convention `f(t < t0) = 0` inconsistently (`sqrtFun`/`erfFun` leaked `y0` before `t0`; the others clamped to 0) and duplicated what an additive component expresses directly. Migration for model YAMLs: delete `y0:` lines that were `0` (the common case); replace a nonzero `y0` with a `stepFun` component sharing the function's `t0`. Stale `y0:` entries now fail model validation with a parameter-count error. The energy-domain `Offset(x, y0)` background is unaffected. +- **Breaking: `Voigt` amplitude is now grid-independent.** The profile was normalized by its maximum over the sampled energy window, so the value at a fixed energy depended on the window/grid, and a peak center pushed outside the fit window rescaled the in-window tail up to amplitude `A`. Normalization now uses the analytic peak value `wofz(i·W/(2·SD·√2))` at `dx = 0`; `A` is the true peak height regardless of the sampled grid. Fitted `A` values may shift slightly relative to earlier releases (the two normalizations agree only when the sampled grid contains the exact peak). `voigtCONV` is unchanged (its kernel axis is symmetric around the peak, where max-over-grid is exact). + +### Fixed + +- **`expRiseCONV` docstring claimed the kernel is causal.** It is deliberately anti-causal — the mirror of `expDecayCONV`, nonzero only for `x ≤ 0`, so the convolved response rises before the excitation and saturates at `t0`. The docstring now says so. +- **`boxCONV` docstring claimed "smooth edges"** — the kernel is a hard `|x| ≤ width/2` threshold; docstring corrected. +- **`GLS`/`GLP` docstrings** now state the valid mixing range `m ∈ [0, 1]` (for `GLP`, `m < 0` makes the denominator `1 + 4·m·u²` cross zero, producing NaN/inf). +- **`DS` docstring** now documents that the asymmetric tail extends toward `x < x0` (correct for a kinetic-energy axis; mirrored on a binding-energy axis). + ## [0.9.3] - 2026-06-29 ### Added diff --git a/PLAN.md b/PLAN.md index 0d24ddf..c3371b6 100644 --- a/PLAN.md +++ b/PLAN.md @@ -1,11 +1,71 @@ -# Active Plan +# PLAN: time-functions fixes + y0 overhaul -No active multi-step feature in progress. +Branch: `fix-time-functions`. Source: audit report (`~/Desktop/functions bugs.txt`), +all findings verified. Expanded scope (agreed 2026-07-06): remove `y0` from all +dynamics primitives in favor of a dedicated step function, plus the Voigt +normalization fix and the remaining docstring corrections. -- Backlog lives in [`TODO.md`](TODO.md). -- Shipped-feature design notes live in [`docs/design/`](docs/design/) (and - archived deep-dives in [`docs/design/archive/`](docs/design/archive/)). +## Design decisions -Populate this file when starting the next multi-step feature; clear it on -completion per `CLAUDE.md` (archive the design into `docs/design/` if the -decisions are durable, otherwise let the changelog stand as the record). +- **Remove `y0`** from all six dynamics functions: `linFun`, `expFun`, `sinFun`, + `sinDivX`, `erfFun`, `sqrtFun`. Rationale: `y0` breaks the causality convention + (`f(t 24 PARAM_INPUT(pos=0) # A -> dynamics - 21 -> 24 PARAM_INPUT(pos=1) # tau -> dynamics - 22 -> 24 PARAM_INPUT(pos=2) # t0 -> dynamics - 23 -> 24 PARAM_INPUT(pos=3) # y0 -> dynamics - 2 -> 25 BASE_INPUT # GLP_01_A base value - 24 -> 25 TRACE_INPUT # dynamics trace - 25 -> 10 PARAM_INPUT(pos=0) # resolved A -> GLP_01 (replaces edge 2->10) + 20 -> 23 PARAM_INPUT(pos=0) # A -> dynamics + 21 -> 23 PARAM_INPUT(pos=1) # tau -> dynamics + 22 -> 23 PARAM_INPUT(pos=2) # t0 -> dynamics + 2 -> 24 BASE_INPUT # GLP_01_A base value + 23 -> 24 TRACE_INPUT # dynamics trace + 24 -> 10 PARAM_INPUT(pos=0) # resolved A -> GLP_01 (replaces edge 2->10) ``` #### 1.6 Expression model as a graph @@ -870,9 +868,9 @@ from `functions/time.py`: ```python DYNAMICS_DISPATCH = { - 0: fcts_time.expFun, # (t, A, tau, t0, y0) -> (n_time,) - 1: fcts_time.sinFun, # (t, A, f, phi, t0, y0) -> (n_time,) - 2: fcts_time.linFun, # (t, m, t0, y0) -> (n_time,) + 0: fcts_time.expFun, # (t, A, tau, t0) -> (n_time,) + 1: fcts_time.sinFun, # (t, A, f, phi, t0) -> (n_time,) + 2: fcts_time.linFun, # (t, m, t0) -> (n_time,) ... } ``` @@ -1030,8 +1028,9 @@ This replaces `2 * n_free_params + 1` evaluator calls per iteration with 1 evaluator call + 1 Jacobian call. For 4 free params, that's 9 -> 2 calls, ~4.5x fewer evaluations per iteration. -**Variable projection (VARPRO):** Linear parameters (amplitudes `A`, -offset `y0`, slope `m`) can be solved in closed form given the nonlinear +**Variable projection (VARPRO):** Linear parameters (amplitudes `A`, the +`Offset` background `y0`, slope `m`) can be solved in closed form given +the nonlinear parameters. Reduces optimizer dimensionality. The graph makes identifying linear params straightforward: any param that appears as a linear factor in its component's function. diff --git a/docs/design/repo_architecture.md b/docs/design/repo_architecture.md index 8a12509..b7e34d1 100644 --- a/docs/design/repo_architecture.md +++ b/docs/design/repo_architecture.md @@ -212,8 +212,10 @@ current peak sum (e.g. Shirley). Add new peak or background shapes here. ### `functions/time.py` Dynamics and convolution kernels. Dynamics functions (e.g. `expFun`, -`sinFun`, `linFun`, `erfFun`, `sqrtFun`) share signature -`func(t, par1, ..., t0, y0)` with the invariant `f(t < t0) = 0`. +`sinFun`, `linFun`, `erfFun`, `sqrtFun`, `stepFun`) share signature +`func(t, par1, ..., t0)` with the invariant `f(t < t0) = 0`; constant +offsets are their own additive component (`stepFun`, or `erfFun` for a +broadened onset) rather than a `y0` parameter on every function. Convolution kernels are named `funcCONV` (e.g. `gaussCONV`) with a companion `funcCONV_kernel_width(...)` returning the kernel-width multiplier. Add new time-domain behavior or IRF kernels here. diff --git a/examples/fitting_workflows/01_basic_fitting/data/generate_data.ipynb b/examples/fitting_workflows/01_basic_fitting/data/generate_data.ipynb index 401e81c..5491338 100644 --- a/examples/fitting_workflows/01_basic_fitting/data/generate_data.ipynb +++ b/examples/fitting_workflows/01_basic_fitting/data/generate_data.ipynb @@ -169,7 +169,7 @@ "| `gaussCONV SD` | 3 | IRF width |\n", "| `expFun A` | 5 | Peak-shift amplitude |\n", "| `expFun tau` | 50 | Shift decay time constant |\n", - "| `expFun t0` / `y0` | 0 / 0 | Time zero / offset |" + "| `expFun t0` | 0 | Time zero |" ] } ], diff --git a/examples/fitting_workflows/01_basic_fitting/data/models_time_truth.yaml b/examples/fitting_workflows/01_basic_fitting/data/models_time_truth.yaml index d3dcaa1..360be80 100644 --- a/examples/fitting_workflows/01_basic_fitting/data/models_time_truth.yaml +++ b/examples/fitting_workflows/01_basic_fitting/data/models_time_truth.yaml @@ -8,4 +8,3 @@ MonoExpIRF: A: [5, True] tau: [50, True] t0: [0, True] - y0: [0, True] diff --git a/examples/fitting_workflows/01_basic_fitting/models_time.yaml b/examples/fitting_workflows/01_basic_fitting/models_time.yaml index 8fa610d..640f464 100644 --- a/examples/fitting_workflows/01_basic_fitting/models_time.yaml +++ b/examples/fitting_workflows/01_basic_fitting/models_time.yaml @@ -8,4 +8,3 @@ MonoExpIRF: A: [5, True, 0, 15] tau: [10, True, 1, 100] t0: [0, False] - y0: [0, False] diff --git a/examples/fitting_workflows/02_dependent_parameters/models_time.yaml b/examples/fitting_workflows/02_dependent_parameters/models_time.yaml index 9475358..580cfb7 100644 --- a/examples/fitting_workflows/02_dependent_parameters/models_time.yaml +++ b/examples/fitting_workflows/02_dependent_parameters/models_time.yaml @@ -7,4 +7,3 @@ sine_wave: f: [0.25, True, 1E-3, 10] phi: [0, True, -0.75, 0.75] t0: [0, False] - y0: [0, False] diff --git a/examples/fitting_workflows/03_multi_cycle_dynamics/data/models_time_truth.yaml b/examples/fitting_workflows/03_multi_cycle_dynamics/data/models_time_truth.yaml index eca07ad..c6a84cf 100644 --- a/examples/fitting_workflows/03_multi_cycle_dynamics/data/models_time_truth.yaml +++ b/examples/fitting_workflows/03_multi_cycle_dynamics/data/models_time_truth.yaml @@ -13,7 +13,6 @@ MonoExpNeg: A: [-2, True] tau: [0.4, True] t0: [0, True] - y0: [0, True] # subcycle 2: amplitude mirrors subcycle 1 via expression, tau independent MonoExpPos: @@ -21,4 +20,3 @@ MonoExpPos: A: ["-expFun_01_A"] tau: [0.8, True] t0: [0, True] - y0: [0, True] diff --git a/examples/fitting_workflows/03_multi_cycle_dynamics/models_time.yaml b/examples/fitting_workflows/03_multi_cycle_dynamics/models_time.yaml index 0fccaea..5076e97 100644 --- a/examples/fitting_workflows/03_multi_cycle_dynamics/models_time.yaml +++ b/examples/fitting_workflows/03_multi_cycle_dynamics/models_time.yaml @@ -12,13 +12,12 @@ IRF: SD: [0.1, True, 0, 1] # Subcycle 1: negative exponential shift, relaxing back to 0. -# t0/y0 are on the subcycle's local clock, which resets at each subcycle start. +# t0 is on the subcycle's local clock, which resets at each subcycle start. MonoExpNeg: expFun: A: [-1.5, True, -5, -1E-3] tau: [1, True, 0.05, 5] t0: [0, False] - y0: [0, False] # Subcycle 2: amplitude mirrors subcycle 1 via an expression (no extra free # parameter), while tau is fit independently — expressions work across @@ -28,4 +27,3 @@ MonoExpPos: A: ["-expFun_01_A"] tau: [1, True, 0.05, 5] t0: [0, False] - y0: [0, False] diff --git a/examples/fitting_workflows/04_parameter_profiles/data/models_time_truth.yaml b/examples/fitting_workflows/04_parameter_profiles/data/models_time_truth.yaml index 0af788e..6435757 100644 --- a/examples/fitting_workflows/04_parameter_profiles/data/models_time_truth.yaml +++ b/examples/fitting_workflows/04_parameter_profiles/data/models_time_truth.yaml @@ -14,4 +14,3 @@ BandBendingRecovery: A: [0.5, True] tau: [100, True] t0: [0, True] - y0: [0, True] diff --git a/examples/fitting_workflows/04_parameter_profiles/models_time.yaml b/examples/fitting_workflows/04_parameter_profiles/models_time.yaml index 7129d22..70c9b11 100644 --- a/examples/fitting_workflows/04_parameter_profiles/models_time.yaml +++ b/examples/fitting_workflows/04_parameter_profiles/models_time.yaml @@ -26,4 +26,3 @@ BandBendingRecovery: A: [0.3, True, 0, 1] # collapse amplitude (eV/nm) tau: [50, True, 5, 500] # recovery time constant (ps) t0: [0, False] - y0: [0, False] diff --git a/examples/fitting_workflows/10_model_comparison/models_time.yaml b/examples/fitting_workflows/10_model_comparison/models_time.yaml index b30d354..48469f4 100644 --- a/examples/fitting_workflows/10_model_comparison/models_time.yaml +++ b/examples/fitting_workflows/10_model_comparison/models_time.yaml @@ -10,7 +10,6 @@ MonoExpPosIRF: A: [2, True, 0, 10] tau: [30, True, 1, 100] t0: [0, False] - y0: [0, False] # MonoExpPos — same exponential rise, no IRF convolution. Sharp turn-on # at t = t0; will undershoot the smoothed rise around t = 0. @@ -19,4 +18,3 @@ MonoExpPos: A: [2, True, 0, 10] tau: [30, True, 1, 100] t0: [0, False] - y0: [0, False] diff --git a/examples/fitting_workflows/20_multi_file_independent_fit/models_time.yaml b/examples/fitting_workflows/20_multi_file_independent_fit/models_time.yaml index ceac58f..1bed166 100644 --- a/examples/fitting_workflows/20_multi_file_independent_fit/models_time.yaml +++ b/examples/fitting_workflows/20_multi_file_independent_fit/models_time.yaml @@ -8,4 +8,3 @@ MonoExpPosIRF: A: [5, file, 0, 15] tau: [10, project, 1, 100] t0: [0, False] - y0: [0, False] diff --git a/examples/fitting_workflows/21_multi_file_shared_fit/data/generate_data.ipynb b/examples/fitting_workflows/21_multi_file_shared_fit/data/generate_data.ipynb index ac8f331..ec2134a 100644 --- a/examples/fitting_workflows/21_multi_file_shared_fit/data/generate_data.ipynb +++ b/examples/fitting_workflows/21_multi_file_shared_fit/data/generate_data.ipynb @@ -177,8 +177,7 @@ "| `gaussCONV SD` | 3 | IRF width |\n", "| `expFun A` | varies | Shift amplitude |\n", "| `expFun tau` | 50 | Decay time constant |\n", - "| `expFun t0` | 0 | Time zero |\n", - "| `expFun y0` | 0 | Offset |" + "| `expFun t0` | 0 | Time zero |" ] } ], diff --git a/examples/fitting_workflows/21_multi_file_shared_fit/data/models_time_truth.yaml b/examples/fitting_workflows/21_multi_file_shared_fit/data/models_time_truth.yaml index 6c30ffe..9a48413 100644 --- a/examples/fitting_workflows/21_multi_file_shared_fit/data/models_time_truth.yaml +++ b/examples/fitting_workflows/21_multi_file_shared_fit/data/models_time_truth.yaml @@ -7,4 +7,3 @@ MonoExpPosIRF: A: [5, True, 0, 15] tau: [50, True, 1, 100] t0: [0, False, 0, 1] - y0: [0, False, 0, 1] diff --git a/examples/fitting_workflows/21_multi_file_shared_fit/models_time.yaml b/examples/fitting_workflows/21_multi_file_shared_fit/models_time.yaml index ceac58f..1bed166 100644 --- a/examples/fitting_workflows/21_multi_file_shared_fit/models_time.yaml +++ b/examples/fitting_workflows/21_multi_file_shared_fit/models_time.yaml @@ -8,4 +8,3 @@ MonoExpPosIRF: A: [5, file, 0, 15] tau: [10, project, 1, 100] t0: [0, False] - y0: [0, False] diff --git a/examples/synthetic_data/01_simulator/models_time.yaml b/examples/synthetic_data/01_simulator/models_time.yaml index 2efeb06..fbad753 100644 --- a/examples/synthetic_data/01_simulator/models_time.yaml +++ b/examples/synthetic_data/01_simulator/models_time.yaml @@ -8,20 +8,18 @@ IRF: SD: [5.0E-2, True, 0, 1] MonoExpPos: - # expFun(t, A, tau, t0, y0) + # expFun(t, A, tau, t0) expFun: A: [1, True, 0, 5] tau: [2.5, True, 1, 10] t0: [0, False, 0, 1] - y0: [0, False, 0, 1] MonoExpNeg: - # expFun(t, A, tau, t0, y0) + # expFun(t, A, tau, t0) expFun: A: [-1, True, -5, 0] tau: [2.5, True, 1, 10] t0: [0, False, 0, 1] - y0: [0, False, 0, 1] MonoExpPosIRF: gaussCONV: @@ -30,4 +28,3 @@ MonoExpPosIRF: A: [5, True, 0, 15] tau: [50, True, 1, 100] t0: [0, False, 0, 1] - y0: [0, False, 0, 1] diff --git a/examples/synthetic_data/02_ml_training_data/models_time.yaml b/examples/synthetic_data/02_ml_training_data/models_time.yaml index 4fd2839..9af37f8 100644 --- a/examples/synthetic_data/02_ml_training_data/models_time.yaml +++ b/examples/synthetic_data/02_ml_training_data/models_time.yaml @@ -8,20 +8,18 @@ IRF: SD: [5.0E-2, True, 0, 1] MonoExpPos: - # expFun(t, A, tau, t0, y0) + # expFun(t, A, tau, t0) expFun: A: [1, True, 0, 5] tau: [2.5, True, 1, 10] t0: [0, False, 0, 1] - y0: [0, False, 0, 1] MonoExpNeg: - # expFun(t, A, tau, t0, y0) + # expFun(t, A, tau, t0) expFun: A: [-1, True, -5, 0] tau: [2.5, True, 1, 10] t0: [0, False, 0, 1] - y0: [0, False, 0, 1] MonoExpPosIRF: gaussCONV: @@ -30,4 +28,3 @@ MonoExpPosIRF: A: [5, True, 0, 15] tau: [50, True, 1, 100] t0: [0, False, 0, 1] - y0: [0, False, 0, 1] diff --git a/pyproject.toml b/pyproject.toml index c27ad97..df02c51 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta" [project] name = "trspecfit" -version = "0.9.3" +version = "0.10.0" authors = [ {name = "Johannes Mahl", email = "johannes.a.mahl@gmail.com"}, ] diff --git a/src/trspecfit/eval_2d.py b/src/trspecfit/eval_2d.py index b92f496..88310bf 100644 --- a/src/trspecfit/eval_2d.py +++ b/src/trspecfit/eval_2d.py @@ -29,12 +29,13 @@ # --------------------------------------------------------------------------- DYNAMICS_DISPATCH: dict[int, tuple] = { - DynFuncKind.EXPFUN: (fcts_time.expFun, 4), - DynFuncKind.SINFUN: (fcts_time.sinFun, 5), - DynFuncKind.LINFUN: (fcts_time.linFun, 3), - DynFuncKind.SINDIVX: (fcts_time.sinDivX, 4), - DynFuncKind.ERFFUN: (fcts_time.erfFun, 4), - DynFuncKind.SQRTFUN: (fcts_time.sqrtFun, 3), + DynFuncKind.EXPFUN: (fcts_time.expFun, 3), + DynFuncKind.SINFUN: (fcts_time.sinFun, 4), + DynFuncKind.LINFUN: (fcts_time.linFun, 2), + DynFuncKind.SINDIVX: (fcts_time.sinDivX, 3), + DynFuncKind.ERFFUN: (fcts_time.erfFun, 3), + DynFuncKind.SQRTFUN: (fcts_time.sqrtFun, 2), + DynFuncKind.STEPFUN: (fcts_time.stepFun, 2), } # Convolution kernel dispatch: kernel function evaluated on the frozen diff --git a/src/trspecfit/functions/energy.py b/src/trspecfit/functions/energy.py index 5e9b330..e8646c1 100644 --- a/src/trspecfit/functions/energy.py +++ b/src/trspecfit/functions/energy.py @@ -296,6 +296,8 @@ def GLS(x: np.ndarray, A: float, x0: float, F: float, m: float) -> np.ndarray: - m = 1: Pure Lorentzian - 0 < m < 1: Weighted mixture Typical value: m ≈ 0.3 + Keep m in [0, 1] (set parameter bounds accordingly); values + outside give negative mixture weights and unphysical shapes. Returns ------- @@ -328,6 +330,8 @@ def GLP(x: np.ndarray, A: float, x0: float, F: float, m: float) -> np.ndarray: - m = 1: Pure Lorentzian - 0 < m < 1: Hybrid shape Typical value: m ≈ 0.3 + Keep m in [0, 1] (set parameter bounds accordingly); m < 0 makes + the denominator 1 + 4*m*u² cross zero, producing NaN/inf. Returns ------- @@ -360,6 +364,8 @@ def DS(x: np.ndarray, A: float, x0: float, F: float, alpha: float) -> np.ndarray - 0 < alpha < 0.3: Typical for metals (e.g., Al: 0.10-0.15) - Larger alpha: Stronger asymmetry, more pronounced tail - Range: typically 0-0.5 for physical systems + The asymmetric tail extends toward x < x0, correct for a + kinetic-energy axis; on a binding-energy axis it appears mirrored. Returns ------- diff --git a/src/trspecfit/functions/time.py b/src/trspecfit/functions/time.py index 4b36a22..001bfc8 100644 --- a/src/trspecfit/functions/time.py +++ b/src/trspecfit/functions/time.py @@ -6,11 +6,10 @@ Use CamelCase naming (UpperCamelCase or lowerCamelCase) for function names. **Dynamics Functions:** -Signature: func(t, par1, par2, ..., t0, y0) +Signature: func(t, par1, par2, ..., t0) - t: Time axis (numpy array) - par1, par2, ...: Function-specific parameters - t0: Time zero (function starts at this time) -- y0: Offset value (baseline) - Returns: f(t) = 0 for t < t0, dynamics for t >= t0 **Convolution Kernels:** @@ -23,12 +22,15 @@ **Time Zero Convention:** All dynamics functions are zero before t0 and activate at t >= t0. This reflects physical causality: response begins after excitation. +Exception: erfFun is a smoothed step, so it is nonzero (approaching 0) +shortly before t0; it crosses A/2 at t0. -**Offset Convention:** -Parameter y0 sets the asymptotic value or baseline. -- Decays: approach y0 as t → ∞ -- Rises: start from 0, reach y0 + A -- Oscillations: oscillate around y0 +**Offsets and Baselines:** +Dynamics components combine by addition, so constant offsets are their +own component rather than a parameter of every function: +- stepFun: sharp onset to a constant value A at t0 +- erfFun: Gaussian-broadened onset, rises from 0 to A around t0 +For example, a decay to a nonzero plateau is expFun + stepFun sharing t0. **Time Resolution:** Functions inherit time axis from Dynamics model. Consider: @@ -42,7 +44,6 @@ - A: Amplitude (change in signal) - tau: Time constant (decay/rise time, 1/e point) - t0: Time zero (start of dynamics) -- y0: Offset/baseline value - f: Frequency (for oscillations) - phi: Phase (for oscillations) - SD: Standard deviation (for Gaussian kernels) @@ -93,7 +94,35 @@ def none(t: np.ndarray) -> np.ndarray: # -def linFun(t: np.ndarray, m: float, t0: float, y0: float) -> np.ndarray: +def stepFun(t: np.ndarray, A: float, t0: float) -> np.ndarray: + """ + Step function (constant offset switching on at t0). + + The causal offset primitive: dynamics components combine by addition, + so add stepFun to model baselines or plateaus (e.g. expFun + stepFun + sharing t0 gives a decay to a nonzero plateau). For a Gaussian-broadened + onset use erfFun instead. + + Parameters + ---------- + t : ndarray + Time axis + A : float + Step height (constant value for t >= t0) + t0 : float + Time zero (onset of the step) + + Returns + ------- + ndarray + Step function: 0 for t=t0 + """ + + return np.where(t < t0, 0.0, A) + + +# +def linFun(t: np.ndarray, m: float, t0: float) -> np.ndarray: """ Linear dynamics (constant rate of change). @@ -107,20 +136,18 @@ def linFun(t: np.ndarray, m: float, t0: float, y0: float) -> np.ndarray: - m < 0: Linear decrease t0 : float Time zero (start of linear change) - y0 : float - Offset value at t0 (initial value) Returns ------- ndarray - Linear function: 0 for t=t0 + Linear function: 0 for t=t0 """ - return np.where(t < t0, 0.0, m * (t - t0) + y0) + return np.where(t < t0, 0.0, m * (t - t0)) # -def expFun(t: np.ndarray, A: float, tau: float, t0: float, y0: float) -> np.ndarray: +def expFun(t: np.ndarray, A: float, tau: float, t0: float) -> np.ndarray: """ Exponential decay or rise dynamics. @@ -130,29 +157,25 @@ def expFun(t: np.ndarray, A: float, tau: float, t0: float, y0: float) -> np.ndar Time axis A : float Amplitude (initial change at t0). - - A > 0: Decay from y0+A to y0 - - A < 0: Rise from y0 to y0+|A| + - A > 0: Jumps to A at t0, decays toward 0 + - A < 0: Jumps to -|A| at t0, rises toward 0 tau : float Time constant (1/e time). Units: [time units] At t = t0 + tau, signal changes by factor of e (≈2.718) t0 : float Time zero (start of exponential) - y0 : float - Asymptotic value (baseline as t → ∞) Returns ------- ndarray - Exponential: 0 for t=t0 + Exponential: 0 for t=t0 """ - return np.where(t < t0, 0.0, A * np.exp(-1 / tau * (t - t0)) + y0) + return np.where(t < t0, 0.0, A * np.exp(-1 / tau * (t - t0))) # -def sinFun( - t: np.ndarray, A: float, f: float, phi: float, t0: float, y0: float -) -> np.ndarray: +def sinFun(t: np.ndarray, A: float, f: float, phi: float, t0: float) -> np.ndarray: """ Sinusoidal oscillations (coherent dynamics). @@ -172,20 +195,19 @@ def sinFun( - phi = π: Starts at zero (negative slope) t0 : float Time zero (start of oscillation) - y0 : float - Offset (center line of oscillation) Returns ------- ndarray - Sinusoid: 0 for t=t0 + Sinusoid: 0 for t=t0 + Oscillates around 0; add stepFun to shift the center line. """ - return np.where(t < t0, 0.0, A * np.sin(2 * np.pi * f * (t - t0) + phi) + y0) + return np.where(t < t0, 0.0, A * np.sin(2 * np.pi * f * (t - t0) + phi)) # -def sinDivX(t: np.ndarray, A: float, f: float, t0: float, y0: float) -> np.ndarray: +def sinDivX(t: np.ndarray, A: float, f: float, t0: float) -> np.ndarray: """ Damped sinc function: sin(x)/x oscillation. @@ -199,50 +221,49 @@ def sinDivX(t: np.ndarray, A: float, f: float, t0: float, y0: float) -> np.ndarr Frequency in [1/time units] t0 : float Time zero (start of oscillation) - y0 : float - Offset value Returns ------- ndarray - Sinc oscillation: 0 for t=t0 + Sinc oscillation: 0 for t=t0 """ # np.sinc(u) = sin(pi*u)/(pi*u), so u=2*f*(t-t0) gives sin(2*pi*f*dt)/(2*pi*f*dt) - return np.where(t < t0, 0.0, A * np.sinc(2 * f * (t - t0)) + y0) + return np.where(t < t0, 0.0, A * np.sinc(2 * f * (t - t0))) # -def erfFun(t: np.ndarray, A: float, SD: float, t0: float, y0: float) -> np.ndarray: +def erfFun(t: np.ndarray, A: float, SD: float, t0: float) -> np.ndarray: """ Error function rise (step with Gaussian broadening). - erfFun ≈ step ⊗ Gaussian(SD) + erfFun ≈ stepFun ⊗ Gaussian(SD) + + As a smoothed step this is the one dynamics function that is nonzero + (approaching 0) shortly before t0; it crosses A/2 at t0 and rises to A. Parameters ---------- t : ndarray Time axis A : float - Amplitude (total change from initial to final value) + Amplitude (final value, asymptote as t → ∞) SD : float Standard deviation of Gaussian broadening (rise time ~2.355*SD) Smaller SD → sharper rise t0 : float Center of rise (50% point) - y0 : float - Final value (asymptote as t → ∞) Returns ------- ndarray - Error function: A/2 * (1 + erf((t-t0)/(SD*√2))) + y0 + Error function: A/2 * (1 + erf((t-t0)/(SD*√2))) """ - return np.asarray(A / 2 * (1 + erf((t - t0) / (SD * np.sqrt(2)))) + y0) + return np.asarray(A / 2 * (1 + erf((t - t0) / (SD * np.sqrt(2))))) # -def sqrtFun(t: np.ndarray, A: float, t0: float, y0: float) -> np.ndarray: +def sqrtFun(t: np.ndarray, A: float, t0: float) -> np.ndarray: """ Square root rise (diffusion dynamics). @@ -254,17 +275,15 @@ def sqrtFun(t: np.ndarray, A: float, t0: float, y0: float) -> np.ndarray: Amplitude scaling factor t0 : float Time zero (start of diffusion) - y0 : float - Offset value Returns ------- ndarray - Square root rise: 0 for t=t0 + Square root rise: 0 for t=t0 """ # numpy array .clip sets all t int: # def expRiseCONV(x: np.ndarray, tau: float) -> np.ndarray: """ - Causal exponential rise kernel. + Anti-causal exponential rise kernel (mirror of expDecayCONV). + + The kernel is nonzero only for x <= 0, so the convolved response at + time t draws from the signal at later times: it rises before the + excitation and saturates at t0. Parameters ---------- @@ -474,7 +497,7 @@ def boxCONV(x: np.ndarray, width: float) -> np.ndarray: Returns ------- ndarray - Rectangular function: 1 inside width, 0 outside (with smooth edges) + Rectangular function: 1 inside width, 0 outside (hard edges) """ return np.where(np.abs(x) <= width / 2, 1.0, 0.0) diff --git a/src/trspecfit/graph_ir.py b/src/trspecfit/graph_ir.py index 804f590..64f5866 100644 --- a/src/trspecfit/graph_ir.py +++ b/src/trspecfit/graph_ir.py @@ -128,6 +128,7 @@ class DynFuncKind(IntEnum): SINDIVX = 3 ERFFUN = 4 SQRTFUN = 5 + STEPFUN = 6 # @@ -157,6 +158,7 @@ class ParamSourceKind(IntEnum): "sinDivX": DynFuncKind.SINDIVX, "erfFun": DynFuncKind.ERFFUN, "sqrtFun": DynFuncKind.SQRTFUN, + "stepFun": DynFuncKind.STEPFUN, } _FUNCTION_NAME_TO_PROFILE_FUNC: dict[str, ProfileFuncKind] = { @@ -2894,6 +2896,7 @@ def schedule_2d(graph: GraphIR) -> ScheduledPlan2D: int(DynFuncKind.SINDIVX): fcts_time.sinDivX, int(DynFuncKind.ERFFUN): fcts_time.erfFun, int(DynFuncKind.SQRTFUN): fcts_time.sqrtFun, + int(DynFuncKind.STEPFUN): fcts_time.stepFun, } _CONV_KERNEL_DISPATCH: dict[int, Callable[..., Any]] = { diff --git a/tests/models/file_time.yaml b/tests/models/file_time.yaml index 74e1cdc..034defe 100644 --- a/tests/models/file_time.yaml +++ b/tests/models/file_time.yaml @@ -19,12 +19,11 @@ IRF: SD: [5.0E-2, True, 0, 1] MonoExpPos: - # expFun(t, A, tau, t0, y0) + # expFun(t, A, tau, t0) expFun: A: [1, True, 0, 5] tau: [2.5, True, 1, 10] t0: [0, False, 0, 1] - y0: [0, False, 0, 1] # stronger amplitude for profile round-trip tests (~20% modulation) MonoExpPosStrong: @@ -32,15 +31,55 @@ MonoExpPosStrong: A: [2, True, 0, 10] tau: [2.5, True, 1, 10] t0: [0, False, 0, 1] - y0: [0, False, 0, 1] MonoExpNeg: - # expFun(t, A, tau, t0, y0) + # expFun(t, A, tau, t0) expFun: A: [-1, True, -5, 0] tau: [2.5, True, 1, 10] t0: [0, False, 0, 1] - y0: [0, False, 0, 1] + +MonoStep: + # stepFun(t, A, t0) + stepFun: + A: [1.5, True, 0, 5] + t0: [0, False, 0, 1] + +MonoSin: + # sinFun(t, A, f, phi, t0) + # phi initialized nonzero (with t0=0) so a phi/t0 ordering bug in the + # lowered path diverges from the interpreter even at initial theta. + sinFun: + A: [1, True, 0, 5] + f: [0.05, True, 0.01, 0.2] + phi: [0.3, True, -3.15, 3.15] + t0: [0, False, 0, 1] + +MonoLin: + # linFun(t, m, t0) + linFun: + m: [0.02, True, -1, 1] + t0: [0, False, 0, 1] + +MonoSinDivX: + # sinDivX(t, A, f, t0) + sinDivX: + A: [1, True, 0, 5] + f: [0.05, True, 0.01, 0.2] + t0: [0, False, 0, 1] + +MonoErf: + # erfFun(t, A, SD, t0) + erfFun: + A: [1, True, 0, 5] + SD: [5, True, 1, 20] + t0: [20, False, 0, 50] + +MonoSqrt: + # sqrtFun(t, A, t0) + sqrtFun: + A: [0.2, True, 0, 5] + t0: [0, False, 0, 1] MonoExpPosIRF: gaussCONV: @@ -49,7 +88,6 @@ MonoExpPosIRF: A: [1, True, 0, 5] tau: [2.5, True, 1, 10] t0: [0, False, 0, 1] - y0: [0, False, 0, 1] MonoExpPosLorentzIRF: lorentzCONV: @@ -58,7 +96,6 @@ MonoExpPosLorentzIRF: A: [1, True, 0, 5] tau: [2.5, True, 1, 10] t0: [0, False, 0, 1] - y0: [0, False, 0, 1] MonoExpPosVoigtIRF: voigtCONV: @@ -68,7 +105,6 @@ MonoExpPosVoigtIRF: A: [1, True, 0, 5] tau: [2.5, True, 1, 10] t0: [0, False, 0, 1] - y0: [0, False, 0, 1] MonoExpPosExpSymIRF: expSymCONV: @@ -77,7 +113,6 @@ MonoExpPosExpSymIRF: A: [1, True, 0, 5] tau: [2.5, True, 1, 10] t0: [0, False, 0, 1] - y0: [0, False, 0, 1] MonoExpPosExpDecayIRF: expDecayCONV: @@ -86,7 +121,6 @@ MonoExpPosExpDecayIRF: A: [1, True, 0, 5] tau: [2.5, True, 1, 10] t0: [0, False, 0, 1] - y0: [0, False, 0, 1] MonoExpPosExpRiseIRF: expRiseCONV: @@ -95,7 +129,6 @@ MonoExpPosExpRiseIRF: A: [1, True, 0, 5] tau: [2.5, True, 1, 10] t0: [0, False, 0, 1] - y0: [0, False, 0, 1] MonoExpPosBoxIRF: boxCONV: @@ -104,7 +137,6 @@ MonoExpPosBoxIRF: A: [1, True, 0, 5] tau: [2.5, True, 1, 10] t0: [0, False, 0, 1] - y0: [0, False, 0, 1] # Empty subcycle placeholder (applies to all times in multi-cycle) ModelNone: @@ -118,7 +150,6 @@ MonoExpPosExpr: A: ["-expFun_01_A"] tau: ["expFun_01_tau"] t0: [0, False, 0, 1] - y0: [0, False, 0, 1] # Bi-exponential: two expFun, second shares t0 with first via expression. # Single-cycle, no convolution -- fully lowerable. @@ -127,12 +158,10 @@ BiExpSharedT0: A: [3, True, 0, 10] tau: [5, True, 1, 20] t0: [0, False, 0, 1] - y0: [0, False, 0, 1] expFun: A: [-1, True, -5, 0] tau: [20, True, 5, 50] t0: ["expFun_01_t0"] - y0: [0, False, 0, 1] # convolution as last component (should fail ordering validation) conv_last: @@ -140,6 +169,5 @@ conv_last: A: [1, True, 0, 5] tau: [2.5, True, 1, 10] t0: [0, False, 0, 1] - y0: [0, False, 0, 1] gaussCONV: SD: [5.0E-2, True, 0, 1] \ No newline at end of file diff --git a/tests/models/project_time.yaml b/tests/models/project_time.yaml index e0326b8..084d814 100644 --- a/tests/models/project_time.yaml +++ b/tests/models/project_time.yaml @@ -6,7 +6,6 @@ MonoExpProject: A: [5, "file", 0, 15] tau: [5, "project", 1, 100] t0: [0, "static", 0, 1] - y0: [0, "static", 0, 1] # Bi-exponential decay with Gaussian IRF — t0 of second component is # an expression referencing t0 of first (project-vary). Tests that @@ -20,9 +19,7 @@ BiExpProject: A: [2, "file", -20, 20] tau: [2, "project", 0.1, 50] t0: [3, "project", -5, 10] - y0: [0, "static", 0, 1] expFun: A: [1, "file", -20, 20] tau: [20, "project", 1, 200] t0: ["expFun_01_t0"] - y0: [0, "static", 0, 1] diff --git a/tests/test_config_functions.py b/tests/test_config_functions.py index 2117fef..864b8c6 100644 --- a/tests/test_config_functions.py +++ b/tests/test_config_functions.py @@ -51,10 +51,13 @@ def test_LinBack_strips_spectrum(self): # time: dynamics functions — first arg (t) stripped # def test_expFun(self): - assert get_function_parameters("expFun") == ["A", "tau", "t0", "y0"] + assert get_function_parameters("expFun") == ["A", "tau", "t0"] def test_sinFun(self): - assert get_function_parameters("sinFun") == ["A", "f", "phi", "t0", "y0"] + assert get_function_parameters("sinFun") == ["A", "f", "phi", "t0"] + + def test_stepFun(self): + assert get_function_parameters("stepFun") == ["A", "t0"] def test_none_has_no_params(self): """none() takes only t, so parameters should be empty.""" diff --git a/tests/test_evaluate_2d.py b/tests/test_evaluate_2d.py index 6e82737..6aefea4 100644 --- a/tests/test_evaluate_2d.py +++ b/tests/test_evaluate_2d.py @@ -249,6 +249,49 @@ def test_glp_time_dep_amplitude(self): F_idx = plan.opt_param_names.index("GLP_01_F") _perturb_theta(plan, model, theta, [A_idx, F_idx], [5.0, 0.3]) + # Per-function parity for every non-expFun dynamics driver: each must + # match MCP at initial theta AND under perturbation of its params. + # expFun is already exercised pervasively via MonoExpPos; this list + # covers the remaining DynFuncKind enum/dispatch entries and their + # param counts (all decremented when y0 was removed). Perturbing more + # than one param where available also catches param-order bugs. + _DYN_DRIVER_PARITY_CASES = [ + ("MonoStep", [("stepFun_01_A", 1.0)]), + ( + "MonoSin", + [("sinFun_01_A", 0.5), ("sinFun_01_f", 0.02), ("sinFun_01_phi", 0.1)], + ), + ("MonoLin", [("linFun_01_m", 0.01)]), + ("MonoSinDivX", [("sinDivX_01_A", 0.5), ("sinDivX_01_f", 0.02)]), + ("MonoErf", [("erfFun_01_A", 0.5), ("erfFun_01_SD", 1.0)]), + ("MonoSqrt", [("sqrtFun_01_A", 0.1)]), + ] + + # + @pytest.mark.parametrize( + ("dyn_model", "param_perturbations"), + _DYN_DRIVER_PARITY_CASES, + ids=[c[0] for c in _DYN_DRIVER_PARITY_CASES], + ) + def test_dynamics_driver_parity(self, dyn_model, param_perturbations): + """evaluate_2d == interpreter for every lowerable dynamics driver.""" + + _file, model = _make_2d_model( + ["glp_only"], + [("GLP_01_A", [dyn_model])], + ) + graph = build_graph(model) + assert can_lower_2d(graph) + plan = schedule_2d(graph) + + theta = _compare_evaluator_vs_interpreter(model, plan) + perturb_idx = [ + plan.opt_param_names.index(f"GLP_01_A_{suffix}") + for suffix, _ in param_perturbations + ] + perturb_deltas = [delta for _, delta in param_perturbations] + _perturb_theta(plan, model, theta, perturb_idx, perturb_deltas) + def test_gls_time_dep_mixing(self): """GLS with mixing parameter m time-dependent.""" diff --git a/tests/test_functions_time.py b/tests/test_functions_time.py index 07f2003..fe203e4 100644 --- a/tests/test_functions_time.py +++ b/tests/test_functions_time.py @@ -17,6 +17,7 @@ sinDivX, sinFun, sqrtFun, + stepFun, ) @@ -38,47 +39,88 @@ def test_returns_zeros(self): assert result.shape == t.shape +# +# +class TestStepFun: + # + def test_zero_before_t0(self): + t = make_time_axis() + result = stepFun(t, A=2.0, t0=0.0) + np.testing.assert_allclose(result[t < 0], 0.0) + + # + def test_constant_after_t0(self): + t = make_time_axis() + A = 2.5 + result = stepFun(t, A=A, t0=0.0) + np.testing.assert_allclose(result[t >= 0], A) + + # + def test_shifted_t0(self): + t = make_time_axis() + result = stepFun(t, A=1.0, t0=5.0) + np.testing.assert_allclose(result[t < 5.0], 0.0) + np.testing.assert_allclose(result[t >= 5.0], 1.0) + + # + def test_zero_amplitude(self): + """A=0 gives all zeros.""" + + t = make_time_axis() + result = stepFun(t, A=0.0, t0=0.0) + np.testing.assert_allclose(result, 0.0, atol=1e-15) + + # + def test_offset_composition(self): + """expFun + stepFun sharing t0 reproduces a decay to a plateau.""" + + t = make_time_axis() + combined = expFun(t, A=5.0, tau=2.0, t0=0.0) + stepFun(t, A=1.0, t0=0.0) + np.testing.assert_allclose(combined[t < 0], 0.0) + assert combined[-1] == pytest.approx(1.0, abs=1e-6) + + # # class TestLinFun: # def test_zero_before_t0(self): t = make_time_axis() - result = linFun(t, m=2.0, t0=0.0, y0=1.0) + result = linFun(t, m=2.0, t0=0.0) np.testing.assert_allclose(result[t < 0], 0.0) # def test_value_at_t0(self): t = make_time_axis() - result = linFun(t, m=2.0, t0=0.0, y0=1.0) + result = linFun(t, m=2.0, t0=0.0) idx = np.argmin(np.abs(t - 0.0)) - assert result[idx] == pytest.approx(1.0, abs=0.01) + assert result[idx] == pytest.approx(0.0, abs=0.01) # def test_slope(self): t = make_time_axis() m = 3.0 - result = linFun(t, m=m, t0=0.0, y0=0.0) + result = linFun(t, m=m, t0=0.0) # At t=10, value should be m*10 = 30 idx = np.argmin(np.abs(t - 10.0)) assert result[idx] == pytest.approx(30.0, abs=0.1) # - def test_offset_t0(self): + def test_shifted_t0(self): t = make_time_axis() - result = linFun(t, m=1.0, t0=5.0, y0=2.0) + result = linFun(t, m=1.0, t0=5.0) # Zero before t0=5 np.testing.assert_allclose(result[t < 5.0], 0.0) - # At t=10: m*(10-5) + y0 = 7 + # At t=10: m*(10-5) = 5 idx = np.argmin(np.abs(t - 10.0)) - assert result[idx] == pytest.approx(7.0, abs=0.1) + assert result[idx] == pytest.approx(5.0, abs=0.1) # def test_zero_slope(self): - """m=0 with y0=0 gives all zeros.""" + """m=0 gives all zeros.""" t = make_time_axis() - result = linFun(t, m=0.0, t0=0.0, y0=0.0) + result = linFun(t, m=0.0, t0=0.0) np.testing.assert_allclose(result, 0.0, atol=1e-15) # @@ -86,7 +128,7 @@ def test_monotonic_positive_slope(self): """Positive m: monotonically increasing after t0.""" t = make_time_axis() - result = linFun(t, m=2.0, t0=0.0, y0=0.0) + result = linFun(t, m=2.0, t0=0.0) active = result[t >= 0] assert np.all(np.diff(active) >= -1e-12) @@ -97,50 +139,50 @@ class TestExpFun: # def test_zero_before_t0(self): t = make_time_axis() - result = expFun(t, A=1.0, tau=5.0, t0=0.0, y0=0.0) + result = expFun(t, A=1.0, tau=5.0, t0=0.0) np.testing.assert_allclose(result[t < 0], 0.0) # def test_value_at_t0(self): t = make_time_axis() - result = expFun(t, A=3.0, tau=5.0, t0=0.0, y0=1.0) + result = expFun(t, A=3.0, tau=5.0, t0=0.0) idx = np.argmin(np.abs(t - 0.0)) - assert result[idx] == pytest.approx(4.0, abs=0.01) # A + y0 + assert result[idx] == pytest.approx(3.0, abs=0.01) # A # - def test_decay_to_y0(self): - """At t >> tau, value approaches y0.""" + def test_decay_to_zero(self): + """At t >> tau, value approaches 0.""" t = make_time_axis() - result = expFun(t, A=5.0, tau=2.0, t0=0.0, y0=1.0) - assert result[-1] == pytest.approx(1.0, abs=1e-6) + result = expFun(t, A=5.0, tau=2.0, t0=0.0) + assert result[-1] == pytest.approx(0.0, abs=1e-6) # def test_value_at_one_tau(self): - """At t = t0 + tau, value = A*exp(-1) + y0.""" + """At t = t0 + tau, value = A*exp(-1).""" t = make_time_axis() tau = 5.0 - result = expFun(t, A=1.0, tau=tau, t0=0.0, y0=0.0) + result = expFun(t, A=1.0, tau=tau, t0=0.0) idx = np.argmin(np.abs(t - tau)) assert result[idx] == pytest.approx(np.exp(-1), abs=1e-3) # def test_negative_amplitude_rise(self): - """A < 0 gives a rise from y0+A toward y0.""" + """A < 0 jumps to -|A| at t0 and rises toward 0.""" t = make_time_axis() - result = expFun(t, A=-2.0, tau=5.0, t0=0.0, y0=0.0) + result = expFun(t, A=-2.0, tau=5.0, t0=0.0) idx_t0 = np.argmin(np.abs(t - 0.0)) assert result[idx_t0] == pytest.approx(-2.0, abs=0.01) assert result[-1] == pytest.approx(0.0, abs=1e-3) # def test_zero_amplitude(self): - """A=0 with y0=0 gives all zeros.""" + """A=0 gives all zeros.""" t = make_time_axis() - result = expFun(t, A=0.0, tau=5.0, t0=0.0, y0=0.0) + result = expFun(t, A=0.0, tau=5.0, t0=0.0) np.testing.assert_allclose(result, 0.0, atol=1e-15) # @@ -148,7 +190,7 @@ def test_monotonic_decay(self): """Positive A: monotonically decreasing after t0.""" t = make_time_axis() - result = expFun(t, A=3.0, tau=5.0, t0=0.0, y0=0.0) + result = expFun(t, A=3.0, tau=5.0, t0=0.0) active = result[t >= 0] assert np.all(np.diff(active) <= 1e-12) @@ -159,7 +201,7 @@ class TestSinFun: # def test_zero_before_t0(self): t = make_time_axis() - result = sinFun(t, A=1.0, f=0.5, phi=0.0, t0=0.0, y0=0.0) + result = sinFun(t, A=1.0, f=0.5, phi=0.0, t0=0.0) np.testing.assert_allclose(result[t < 0], 0.0) # @@ -168,7 +210,7 @@ def test_frequency(self): t = make_time_axis() f = 0.5 - result = sinFun(t, A=1.0, f=f, phi=0.0, t0=0.0, y0=0.0) + result = sinFun(t, A=1.0, f=f, phi=0.0, t0=0.0) # At t = 1/(4f), should be at maximum (A) idx = np.argmin(np.abs(t - 1 / (4 * f))) assert result[idx] == pytest.approx(1.0, abs=0.02) @@ -178,29 +220,26 @@ def test_phase_shift(self): """Phi = pi/2 turns sin into cos (starts at maximum).""" t = make_time_axis() - result = sinFun(t, A=1.0, f=0.5, phi=np.pi / 2, t0=0.0, y0=0.0) + result = sinFun(t, A=1.0, f=0.5, phi=np.pi / 2, t0=0.0) idx_t0 = np.argmin(np.abs(t - 0.0)) assert result[idx_t0] == pytest.approx(1.0, abs=0.02) # def test_zero_amplitude(self): - """A=0 should return y0 for t>=t0, 0 for t= 0], 2.0) + result = sinFun(t, A=0.0, f=0.5, phi=0.0, t0=0.0) + np.testing.assert_allclose(result, 0.0, atol=1e-15) # - def test_offset(self): - """y0 shifts the oscillation center.""" + def test_zero_at_t0(self): + """Phi = 0: oscillation starts at zero and centers on zero.""" t = make_time_axis() - result = sinFun(t, A=1.0, f=0.5, phi=0.0, t0=0.0, y0=3.0) - # Mean of oscillation should be ~y0 over full cycles - # At t0, sin=0 so value = y0 + result = sinFun(t, A=1.0, f=0.5, phi=0.0, t0=0.0) idx_t0 = np.argmin(np.abs(t - 0.0)) - assert result[idx_t0] == pytest.approx(3.0, abs=0.02) + assert result[idx_t0] == pytest.approx(0.0, abs=0.02) # @@ -209,7 +248,7 @@ class TestSinDivX: # def test_zero_before_t0(self): t = make_time_axis() - result = sinDivX(t, A=1.0, f=0.5, t0=0.0, y0=0.0) + result = sinDivX(t, A=1.0, f=0.5, t0=0.0) np.testing.assert_allclose(result[t < 0], 0.0) # @@ -217,7 +256,7 @@ def test_decaying_envelope(self): """Amplitude should decrease over time (sinc envelope).""" t = make_time_axis() - result = sinDivX(t, A=1.0, f=0.5, t0=0.0, y0=0.0) + result = sinDivX(t, A=1.0, f=0.5, t0=0.0) active = result[t > 0.5] # skip near t0 where sinc diverges peaks = np.abs(active[1:-1])[ (active[1:-1] > active[:-2]) & (active[1:-1] > active[2:]) @@ -227,14 +266,13 @@ def test_decaying_envelope(self): # def test_value_at_t0(self): - """At t=t0, sinc(0)=1 so value should be A + y0.""" + """At t=t0, sinc(0)=1 so value should be A.""" t = make_time_axis() A = 1.5 - y0 = 0.2 - result = sinDivX(t, A=A, f=0.5, t0=0.0, y0=y0) + result = sinDivX(t, A=A, f=0.5, t0=0.0) idx_t0 = np.argmin(np.abs(t - 0.0)) - assert result[idx_t0] == pytest.approx(A + y0, abs=0.01) + assert result[idx_t0] == pytest.approx(A, abs=0.01) # def test_first_zero_location(self): @@ -243,27 +281,26 @@ def test_first_zero_location(self): t = make_time_axis() f = 0.5 t0 = 0.0 - y0 = 0.3 - result = sinDivX(t, A=1.0, f=f, t0=t0, y0=y0) + result = sinDivX(t, A=1.0, f=f, t0=t0) t_zero = t0 + 1.0 / (2.0 * f) idx_zero = np.argmin(np.abs(t - t_zero)) - assert result[idx_zero] == pytest.approx(y0, abs=0.01) + assert result[idx_zero] == pytest.approx(0.0, abs=0.01) # def test_zero_amplitude(self): - """A=0 with y0=0 gives all zeros.""" + """A=0 gives all zeros.""" t = make_time_axis() - result = sinDivX(t, A=0.0, f=0.5, t0=0.0, y0=0.0) + result = sinDivX(t, A=0.0, f=0.5, t0=0.0) np.testing.assert_allclose(result, 0.0, atol=1e-15) # def test_asymptotic_value(self): - """For t >> t0, sinc → 0 so value → y0.""" + """For t >> t0, sinc → 0 so value → 0.""" t = np.linspace(-10, 1000, 10000) - result = sinDivX(t, A=2.0, f=0.5, t0=0.0, y0=3.0) - assert result[-1] == pytest.approx(3.0, abs=0.01) + result = sinDivX(t, A=2.0, f=0.5, t0=0.0) + assert result[-1] == pytest.approx(0.0, abs=0.01) # @@ -271,52 +308,52 @@ def test_asymptotic_value(self): class TestErfFun: # def test_midpoint_value(self): - """At t = t0, erf(0) = 0, so value = A/2 + y0.""" + """At t = t0, erf(0) = 0, so value = A/2.""" t = make_time_axis() - result = erfFun(t, A=4.0, SD=1.0, t0=10.0, y0=1.0) + result = erfFun(t, A=4.0, SD=1.0, t0=10.0) idx = np.argmin(np.abs(t - 10.0)) - assert result[idx] == pytest.approx(3.0, abs=0.01) # 4/2 + 1 + assert result[idx] == pytest.approx(2.0, abs=0.01) # A/2 # def test_asymptotic_low(self): - """For t << t0, erf → -1, so value → y0.""" + """For t << t0, erf → -1, so value → 0.""" t = make_time_axis() - result = erfFun(t, A=4.0, SD=1.0, t0=10.0, y0=1.0) - assert result[0] == pytest.approx(1.0, abs=1e-6) + result = erfFun(t, A=4.0, SD=1.0, t0=10.0) + assert result[0] == pytest.approx(0.0, abs=1e-6) # def test_asymptotic_high(self): - """For t >> t0, erf → 1, so value → A + y0.""" + """For t >> t0, erf → 1, so value → A.""" t = make_time_axis() - result = erfFun(t, A=4.0, SD=1.0, t0=10.0, y0=1.0) - assert result[-1] == pytest.approx(5.0, abs=1e-3) + result = erfFun(t, A=4.0, SD=1.0, t0=10.0) + assert result[-1] == pytest.approx(4.0, abs=1e-3) # def test_monotonic_increase(self): """Error function rise should be monotonically increasing.""" t = make_time_axis() - result = erfFun(t, A=2.0, SD=1.0, t0=10.0, y0=0.0) + result = erfFun(t, A=2.0, SD=1.0, t0=10.0) assert np.all(np.diff(result) >= -1e-12) # def test_zero_amplitude(self): - """A=0 with y0=0 gives all zeros.""" + """A=0 gives all zeros.""" t = make_time_axis() - result = erfFun(t, A=0.0, SD=1.0, t0=10.0, y0=0.0) + result = erfFun(t, A=0.0, SD=1.0, t0=10.0) np.testing.assert_allclose(result, 0.0, atol=1e-15) # def test_note_no_hard_t0_cutoff(self): """erfFun does NOT have a hard t0 cutoff — it's a smooth sigmoid. - Value at t << t0 approaches y0 but is never exactly zero.""" + Value at t << t0 approaches zero but is never exactly zero.""" t = make_time_axis() - result = erfFun(t, A=4.0, SD=1.0, t0=10.0, y0=0.0) + result = erfFun(t, A=4.0, SD=1.0, t0=10.0) assert result[0] == pytest.approx(0.0, abs=1e-6) # Slight non-zero values near t0 are expected idx_before = np.argmin(np.abs(t - 8.0)) # 2 SD before @@ -328,55 +365,47 @@ def test_note_no_hard_t0_cutoff(self): class TestSqrtFun: # def test_zero_before_t0(self): - """ - sqrtFun uses .clip(0) so t= 0] assert np.all(np.diff(active) >= -1e-12) # def test_zero_amplitude(self): - """A=0 with y0=0 gives all zeros.""" + """A=0 gives all zeros.""" t = make_time_axis() - result = sqrtFun(t, A=0.0, t0=0.0, y0=0.0) + result = sqrtFun(t, A=0.0, t0=0.0) np.testing.assert_allclose(result, 0.0, atol=1e-15) - # - def test_offset_y0(self): - """y0 shifts everything (including before t0 via clip behavior).""" - - t = make_time_axis() - result = sqrtFun(t, A=1.0, t0=0.0, y0=5.0) - # Before t0: A*sqrt(0) + y0 = y0 - np.testing.assert_allclose(result[t <= 0], 5.0, atol=0.01) - if __name__ == "__main__": pytest.main([__file__]) diff --git a/tests/test_graph_ir.py b/tests/test_graph_ir.py index eeff438..060fdf8 100644 --- a/tests/test_graph_ir.py +++ b/tests/test_graph_ir.py @@ -452,15 +452,15 @@ def test_dynamics_param_nodes(self): ) graph = build_graph(model) - # expFun has params: A, tau, t0, y0 - # A and tau are vary=True, t0 and y0 are vary=False + # expFun has params: A, tau, t0 + # A and tau are vary=True, t0 is vary=False dyn_A = _node_by_name(graph, "GLP_01_A_expFun_01_A") assert dyn_A is not None assert dyn_A.kind == NodeKind.OPT_PARAM - dyn_y0 = _node_by_name(graph, "GLP_01_A_expFun_01_y0") - assert dyn_y0 is not None - assert dyn_y0.kind == NodeKind.STATIC_PARAM + dyn_t0 = _node_by_name(graph, "GLP_01_A_expFun_01_t0") + assert dyn_t0 is not None + assert dyn_t0.kind == NodeKind.STATIC_PARAM # def test_dynamics_edges(self): @@ -477,7 +477,7 @@ def test_dynamics_edges(self): assert len(trace_nodes) >= 1 trace_nid = trace_nodes[0].id param_edges = _edges_to(graph, trace_nid, EdgeKind.PARAM_INPUT) - assert len(param_edges) == 4 # A, tau, t0, y0 + assert len(param_edges) == 3 # A, tau, t0 # PARAM_PLUS_TRACE has BASE_INPUT + TRACE_INPUT resolved = _node_by_name(graph, "GLP_01_A_resolved") @@ -1307,16 +1307,16 @@ def test_profile_par_dynamics_params_in_graph(self): _file, model = _make_time_dep_profile_model() graph = build_graph(model) - # MonoExpPos has params: A, tau, t0, y0 + # MonoExpPos has params: A, tau, t0 # The profile par is GLP_01_A_pLinear_01_m, so dynamics params # are prefixed: GLP_01_A_pLinear_01_m_expFun_01_A, etc. dyn_A = _node_by_name(graph, "GLP_01_A_pLinear_01_m_expFun_01_A") assert dyn_A is not None assert dyn_A.kind == NodeKind.OPT_PARAM - dyn_y0 = _node_by_name(graph, "GLP_01_A_pLinear_01_m_expFun_01_y0") - assert dyn_y0 is not None - assert dyn_y0.kind == NodeKind.STATIC_PARAM + dyn_t0 = _node_by_name(graph, "GLP_01_A_pLinear_01_m_expFun_01_t0") + assert dyn_t0 is not None + assert dyn_t0.kind == NodeKind.STATIC_PARAM # def test_resolved_profile_par_in_sample_edges(self): @@ -1992,7 +1992,7 @@ def test_dynamics_compiled(self): # Single substep in the group assert plan.dyn_group_indptr[1] - plan.dyn_group_indptr[0] == 1 assert plan.dyn_sub_func_id[0] == int(DynFuncKind.EXPFUN) - assert plan.dyn_sub_n_params[0] == 4 # A, tau, t0, y0 + assert plan.dyn_sub_n_params[0] == 3 # A, tau, t0 # def test_dynamics_param_rows_valid(self): @@ -2610,13 +2610,6 @@ def test_synthetic_graph_with_prefixed_expr_ref(self): source_order=3, value=0.0, ), - GraphNode( - id=4, - kind=NodeKind.STATIC_PARAM, - name="dyn_y0", - source_order=4, - value=0.0, - ), GraphNode( id=5, kind=NodeKind.DYNAMICS_TRACE, @@ -2673,7 +2666,6 @@ def test_synthetic_graph_with_prefixed_expr_ref(self): GraphEdge(source=1, target=5, kind=EdgeKind.PARAM_INPUT, position=0), GraphEdge(source=2, target=5, kind=EdgeKind.PARAM_INPUT, position=1), GraphEdge(source=3, target=5, kind=EdgeKind.PARAM_INPUT, position=2), - GraphEdge(source=4, target=5, kind=EdgeKind.PARAM_INPUT, position=3), GraphEdge(source=0, target=6, kind=EdgeKind.BASE_INPUT), GraphEdge(source=5, target=6, kind=EdgeKind.TRACE_INPUT), # EXPR_REF: B_expr references A_base_resolved @@ -2861,13 +2853,6 @@ def test_sparse_node_ids(self): source_order=3, value=0.0, ), - GraphNode( - id=104, - kind=NodeKind.STATIC_PARAM, - name="dyn_y0", - source_order=4, - value=0.0, - ), GraphNode( id=105, kind=NodeKind.DYNAMICS_TRACE, @@ -2915,7 +2900,6 @@ def test_sparse_node_ids(self): GraphEdge(source=101, target=105, kind=EdgeKind.PARAM_INPUT, position=0), GraphEdge(source=102, target=105, kind=EdgeKind.PARAM_INPUT, position=1), GraphEdge(source=103, target=105, kind=EdgeKind.PARAM_INPUT, position=2), - GraphEdge(source=104, target=105, kind=EdgeKind.PARAM_INPUT, position=3), GraphEdge(source=100, target=106, kind=EdgeKind.BASE_INPUT), GraphEdge(source=105, target=106, kind=EdgeKind.TRACE_INPUT), GraphEdge(source=106, target=109, kind=EdgeKind.PARAM_INPUT, position=0), diff --git a/tests/test_mcp_library.py b/tests/test_mcp_library.py index caeeab5..b652557 100644 --- a/tests/test_mcp_library.py +++ b/tests/test_mcp_library.py @@ -279,7 +279,6 @@ def test_dynamics_model_creation(self): "A": [2, True, 1, 1e2], "tau": [5000, True, 1e3, 1e4], "t0": [0, False, 0, 1], - "y0": [0, False, 0, 1], } ) @@ -289,7 +288,6 @@ def test_dynamics_model_creation(self): "A": [5, True, 1, 1e2], "tau": [1250, True, 1e2, 1e3], "t0": [0, False, 0, 1], - "y0": [0, False, 0, 1], } ) @@ -318,7 +316,6 @@ def test_dynamics_parameter_handling(self): "A": [1, True, 0, 5], "tau": [2.5, True, 1, 10], "t0": [0, False, 0, 1], - "y0": [0, False, 0, 1], } ) diff --git a/tests/test_model_parser.py b/tests/test_model_parser.py index a363111..b9c5475 100644 --- a/tests/test_model_parser.py +++ b/tests/test_model_parser.py @@ -215,7 +215,6 @@ def test_simple_time_model(self): assert model.components[0].par_dict["A"] == [1, True, 0, 5] assert model.components[0].par_dict["tau"] == [2.5, True, 1, 10] assert model.components[0].par_dict["t0"] == [0, False, 0, 1] - assert model.components[0].par_dict["y0"] == [0, False, 0, 1] # def test_IRF_model(self): @@ -239,7 +238,6 @@ def test_IRF_model(self): assert model.components[1].par_dict["A"] == [1, True, 0, 5] assert model.components[1].par_dict["tau"] == [2.5, True, 1, 10] assert model.components[1].par_dict["t0"] == [0, False, 0, 1] - assert model.components[1].par_dict["y0"] == [0, False, 0, 1] # def test_multi_cycle_expression_model(self): @@ -271,7 +269,6 @@ def test_multi_cycle_expression_model(self): assert model.components[2].par_dict["A"] == ["-expFun_01_A"] assert model.components[2].par_dict["tau"] == ["expFun_01_tau"] assert model.components[2].par_dict["t0"] == [0, False, 0, 1] - assert model.components[2].par_dict["y0"] == [0, False, 0, 1] # Check lmfit parameters exist with prefixed names assert "parTEST_expFun_01_A" in model.lmfit_pars @@ -351,7 +348,6 @@ def test_simple_2D_model(self): assert td_par_model.components[1].par_dict["A"] == [1, True, 0, 5] assert td_par_model.components[1].par_dict["tau"] == [2.5, True, 1, 10] assert td_par_model.components[1].par_dict["t0"] == [0, False, 0, 1] - assert td_par_model.components[1].par_dict["y0"] == [0, False, 0, 1] # end of time-dependent parameter model assert model.components[2].par_dict["F"] == [1.0, True, 0.75, 2.5] assert model.components[2].par_dict["m"] == [0.3, True, 0, 1] diff --git a/tests/test_project_fit.py b/tests/test_project_fit.py index bca00a0..e04e0b9 100644 --- a/tests/test_project_fit.py +++ b/tests/test_project_fit.py @@ -53,7 +53,6 @@ def _make_truth_file(*, amplitude=20.0, x0_shift=3.0, tau=5.0): model.lmfit_pars["GLP_01_x0_expFun_01_A"].value = x0_shift model.lmfit_pars["GLP_01_x0_expFun_01_tau"].value = tau model.lmfit_pars["GLP_01_x0_expFun_01_t0"].value = 0.0 - model.lmfit_pars["GLP_01_x0_expFun_01_y0"].value = 0.0 return file @@ -144,11 +143,9 @@ def test_biexp_expr_t0_roundtrip(self): m.lmfit_pars["GLP_01_x0_expFun_01_A"].value = dx0_1 m.lmfit_pars["GLP_01_x0_expFun_01_tau"].value = TRUE_TAU1 m.lmfit_pars["GLP_01_x0_expFun_01_t0"].value = TRUE_T0 - m.lmfit_pars["GLP_01_x0_expFun_01_y0"].value = 0.0 m.lmfit_pars["GLP_01_x0_expFun_02_A"].value = dx0_2 m.lmfit_pars["GLP_01_x0_expFun_02_tau"].value = TRUE_TAU2 m.lmfit_pars["GLP_01_x0_expFun_02_t0"].value = TRUE_T0 - m.lmfit_pars["GLP_01_x0_expFun_02_y0"].value = 0.0 truth_files.append((tf, seed)) # --- simulate and build fit files --- @@ -276,7 +273,6 @@ def test_vary_levels_map(self): assert levels["GLP_01_x0_expFun_01_A"] == "file" assert levels["GLP_01_x0_expFun_01_tau"] == "project" assert levels["GLP_01_x0_expFun_01_t0"] == "static" - assert levels["GLP_01_x0_expFun_01_y0"] == "static" # def test_vary_levels_profile_with_dynamics(self): From b78d201409de8d4ae66d0dfdfdd4987e12782d39 Mon Sep 17 00:00:00 2001 From: Johannes Mahl Date: Wed, 8 Jul 2026 13:34:50 -0700 Subject: [PATCH 4/4] improve add-function skill and housekeeping: author e-mail reset PLAN md --- CLAUDE.md | 1 + PLAN.md | 76 +++++------------------------------------ docs/ai/add-function.md | 22 ++++++++++++ pyproject.toml | 2 +- 4 files changed, 32 insertions(+), 69 deletions(-) diff --git a/CLAUDE.md b/CLAUDE.md index eac8197..2aa7aa6 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -5,6 +5,7 @@ - **Token Efficiency:** Be concise. Reference file paths and line numbers rather than quoting large code blocks. - **Subagent Protocol:** Use subagents for repo-wide scans, parallel research, or scanning large directories. Instruct them to return only concise summaries to keep the main context window lean. - **Guardrails:** Use this file for coding guardrails; heavier review checklists live in `docs/ai/code-review.md`. +- **Task recipes & skills:** Common repo tasks have step-by-step recipes in `docs/ai/*.md` (add-function, changelog, check-docs, check-example, benchmark, bump-versions, code-review), each mirrored by a thin `/`-invokable wrapper in `.claude/skills/` whose source of truth is the `docs/ai` file. Most set `disable-model-invocation: true`, so they will NOT appear in an agent's auto-invoke list—before hand-rolling one of these tasks, consult the matching `docs/ai` recipe directly. - **Renaming / API changes:** grep the entire repo—notebooks, YAML, tests, and docs all reference the public API. diff --git a/PLAN.md b/PLAN.md index c3371b6..0d24ddf 100644 --- a/PLAN.md +++ b/PLAN.md @@ -1,71 +1,11 @@ -# PLAN: time-functions fixes + y0 overhaul +# Active Plan -Branch: `fix-time-functions`. Source: audit report (`~/Desktop/functions bugs.txt`), -all findings verified. Expanded scope (agreed 2026-07-06): remove `y0` from all -dynamics primitives in favor of a dedicated step function, plus the Voigt -normalization fix and the remaining docstring corrections. +No active multi-step feature in progress. -## Design decisions +- Backlog lives in [`TODO.md`](TODO.md). +- Shipped-feature design notes live in [`docs/design/`](docs/design/) (and + archived deep-dives in [`docs/design/archive/`](docs/design/archive/)). -- **Remove `y0`** from all six dynamics functions: `linFun`, `expFun`, `sinFun`, - `sinDivX`, `erfFun`, `sqrtFun`. Rationale: `y0` breaks the causality convention - (`f(t