diff --git a/GROK-IDA-INIT.md b/GROK-IDA-INIT.md new file mode 100644 index 0000000000..30c06d96cc --- /dev/null +++ b/GROK-IDA-INIT.md @@ -0,0 +1,271 @@ +# Grok handoff: IDA consistent initialization (`nrnida`) + +Use this file when starting a **new** Grok session rooted in `~/neuron/nrnida`. + +--- + +## Status (2026-08-12) + +| Phase | Status | Notes | +|-------|--------|--------| +| **0** Mode API + `IDA_Y_INIT` plumbing | **Done** | `CVode.dae_init_mode` 0–2; default still heuristic (0) | +| **R** Industry / circuit reinit survey | **Done** | `~/neuron/notes/ida_phase_R_industry_map.md` | +| **Battery IC** (mode **3**) LM + extracellular hold | **Done** | C→V hold \(\Delta v\); L hold; op-amp \(\tau\); ext \(V_m\)/layers | +| Mode 3 **\(C y'=f\)** for \(y'\) | **Done** | Continuous limit; no `dteps` bias; nano-step only on residual fallback | +| Soft sparse13 factor fail | **Done** | Mode 1/3 can fall back without abort | +| **Three-panel IC audit** | **Done** | `dae_init_audit` / `dae_init_audit_file`; audit suppresses mode-3 fallback | +| Singular / algebraic residual diagnosis | **Done** | Top residual eqs classified (algebraic / near-singular \(c\)) on mode-3 fail | +| Automated tests | **Partial** | `test_ida_init_mode.py` (mode 3, SEClamp, forcing suite A3, A5 stats) | +| Manual GUI circuits | **Local** | `external/tests/nrntest/nrniv/ida/*.ses` (gitignored); `~/models/nrndc1sim` | +| **Plan A** forcing \(t^+\) for free \(y'\) | **A0–A5 done** | continuous `Vector.play` + LM `dforce` / FD | +| **MOD density `PROCEDURE dforce`** | **Parked** | tip `210c43c56` on `hines-grok/ida-mod-dforce-park` only | +| **(b) Source-current discontinuities** | **Done (v1+polish)** | Plan: `~/neuron/notes/ida_plan_b_source_currents.md`. PWLClamp; E0–E4; end free-\(y\); electrode–xc seed. **Docs:** `dae_init_mode` electrode/xtral notes in `cvode.rst`. **E2:** true Vm is `seg.v` (held); `vext` may jump when `xc=0`. **CI:** `ctest -R hoctests::test_ida_source_current`. Smoke `nrndc1sim`: mode-3 OK at `finitialize`; per-step `re_init` after transfer often falls back (algebraic clamp/current). Default still mode 0. | + +**Working tip:** `7688c6238` — *A5: IDA IC path stats, clearer mode-3 fallback messages* +**Branch:** `hines-grok/ida-init` (tracks `origin/hines-grok/ida-init`) +**Default remains mode 0.** Mode 3 ready for broader validation; residual fail → heuristic (except audit armed). + +### Parked: density MOD `dforce` + +| Item | Value | +|------|--------| +| Branch | `hines-grok/ida-mod-dforce-park` (also on origin) | +| Commit | `210c43c56` *Call density PROCEDURE dforce at IDA IC* | +| Content | `src/nrnoc/mod_dforce.cpp`, `test/hoctests/mod_dforce.mod`, hook in `nrndaspk.cpp`, docs, test | +| On `ida-init` | **Not present** (reset to A5 after park) | +| Intent | Variable-capacitance / assigned-rate MOD hook; **do not revive** until param-in-C / charge work is scheduled | +| Restore | `git cherry-pick 210c43c56` or merge that branch | + +Motivation for park: that commit was a rush on **time-dependent capacitance** rates. Preferred next slice is **source currents**, not density `dforce`. Defer param-in-C and charge-conservation patterns (`dcmdt` / `NET_RECEIVE` charge jumps). + +### Plan A (forcing \(t^+\) info) — done + +**Forcing \(t^+\) info:** right-limit value \(u(t^+)\) and classical derivative \(u'(t^+)\) of exogenous drives after a discontinuity (or at `finitialize`). Geometric / DAE literature: **1-jet** of \(u\) at \(t^+\). + +| Slice | Status | Notes | +|-------|--------|--------| +| **A0** Spec + continuous-play \(t^+\) oracle | **Done** | `src/nrniv/vecplay_tplus.h`; `VecPlayContinuous::forcing_tplus` | +| **A1** Wire play instances at IC | **Done** | `nrn_collect_forcing_tplus` at `Daspk::init`; `Daspk::last_forcing_tplus` | +| **A2** Mode 3 free \(y'\) from \(u'\) (LM / series CR) | **Done** | null(C) via \(Z^\top G y' = Z^\top b'\); ramp \(V_1'=1,V_2'=0.5\) | +| **A3** Multi-event + finitialize suite | **Done** | istep, kink, end extrap, flat end, finitialize slope, multi-event | +| **A4** LM `dforce` thin API | **Done** | `LinearMechanism.dforce(callable, bdot)`; FD fallback | +| **A5** Polish / diagnostics | **Done** | `dae_init_stats`; clearer fallback msgs; counters | + +Continuous play \(t^+\) rules: hold for \(t < t_0\); **outgoing** slope at a knot (including \(t=t_0\)); **linear extrapolation of the last two points** past the end (flat last segment ⇒ \(u'=0\)). + +--- + +## Why this worktree + +Focused DAE/IDA initial-condition work. Branch from **master**; do not revive permanent VMX / `v12` state doubling from historical `origin/idainit`. + +| Worktree / path | Branch | Purpose | +|-----------------|--------|---------| +| `~/neuron/nrnida` | `hines-grok/ida-init` | IDA IC modes, battery project, audit, Plan A | +| park | `hines-grok/ida-mod-dforce-park` | Density MOD dforce (local/remote park) | +| notes | `~/neuron/notes/` | Design notes outside the repo | + +--- + +## North star (success criteria) + +After `finitialize` and after each discontinuity (`NET_RECEIVE`, `Vector.play`, `at_time`, clamps, …): + +1. **Residual:** \(g(y', y, t) = 0\) within integrator WRMS / ewt tolerance. +2. **Physical \(dt\to 0\):** hold continuous content; free absolute algebraics as network requires. + - **Hold:** ODE states; capacitor \(\Delta v\) / charge; inductor current; membrane \(V_m\) with extracellular; \(xc>0\) layer drops. + - **Free to jump:** absolute node voltages (e.g. series \(C\)–\(R\) + \(I\) step). + +Nano-\(dt\) fully implicit Euler is a **fallback**, not the definition of (2). Mode 3 recovers \(y'\) from \(C y' = f(y,t)\) at fixed projected \(y\). + +### Non-goals (still) + +- Permanent `USE_VMX` / doubled continuous state for integration (`origin/idainit`). +- Block on SUNDIALS 3 / PR #1960 (`slds`) for IC design (revalidate later). +- Pure speculative \(g_{\mathrm{ic}}\) stacks without industry map (Phase R already done). +- **Now deferred:** density `PROCEDURE dforce`; param-in-C; charge-conservation / `dcmdt` productization. + +--- + +## Architecture (current) + +```text +Discontinuity / finitialize + │ + ▼ +Daspk::init() [src/nrncvode/nrndaspk.cpp] + │ collect continuous Vector.play forcing t+ (A1) + │ mode 0 → heuristic nano-step (dae_init_dteps) + │ mode 1 → IDACalcIC(IDA_Y_INIT) then heuristic fallback + │ mode 2 → IDA_Y_INIT only + │ mode 3 → battery y-project → y' from C*y'=f(y) + │ + null(C) free y' from play b' / LM.dforce (A2–A4) + │ residual fail → heuristic + │ (audit armed: no fallback; diagnose top residual eqs) + ▼ +Cvode::res → G ≈ C y' − f(y) (membrane cm, xc layers, nrndae_dkres) +``` + +| Piece | Location | +|-------|----------| +| IC modes / audit / \(C y'=f\) / free \(y'\) | `src/nrncvode/nrndaspk.{h,cpp}` | +| HOC API | `src/nrncvode/cvodeobj.cpp` (`dae_init_mode`, `dae_init_dteps`, `dae_init_audit`, `dae_init_audit_file`, `dae_init_stats`) | +| Battery project (LM) | `src/nrniv/linmod.cpp` `battery_ic_project()` | +| LM `dforce` / bdot | `src/nrniv/linmod1.cpp`, model side in linmod | +| Continuous play \(t^+\) | `src/nrniv/vecplay_tplus.h`, `vrecitem` | +| NrnDAE entry + LM \(y'\) seed | `src/nrniv/nrndae.cpp` | +| Extracellular battery | `src/nrnoc/extcelln.cpp` `nrn_extracellular_battery_ic()` | +| Soft factor fail | `src/nrnoc/solve.cpp` + `nrn_sparse13_soft_fail` | +| Built-in IClamp | `src/nrnoc/stim.mod` (`ELECTRODE_CURRENT i`, `at_time(del)`, `at_time(del+dur)`) | +| Docs | `docs/progref/simctrl/cvode.rst` | +| Tests | `test/hoctests/tests/test_ida_init_mode.py` | +| Design notes | `~/neuron/notes/ida_y_init_adoption.md`, `ida_phase_R_industry_map.md` | + +### Mode 3 idea (productized) + +1. **Project \(y\):** hold continuous content (LM floating \(C\), diagonal \(L\), op-amp lag; extracellular \(V_m\) / \(xc>0\) drops); free absolute algebraics as the network requires. +2. **Recover \(y'\):** diagonal / simple-mass solve \(C y' = f(y,t)\) (membrane \(10^{-3} c_m\), mechanism identity, LM single-column/difference stamps). **No `dteps`.** +3. **Free algebraic rates:** when continuous play (or LM force) has classical \(u'\), adjust \(y'\) in null(\(C\)) via \(Z^\top G y' = Z^\top b'\). +4. **Residual check:** on failure, classify top residual eqs; fall back to nano-step unless audit is armed. + +`finitialize` still chooses \(y\) via `v_init` + `INITIAL` + defaults (`vext=0`). Mode 3 does **not** invent a consistent algebraic \(y\) when that recipe is wrong. + +### Three-panel audit + +```text +CVode().dae_init_audit(2, T) # level 0–2; arm first reinit with t >= T (one-shot) +CVode().dae_init_audit_file("ic_audit.txt") # append; empty → stdout +``` + +| Panel | Meaning | +|-------|---------| +| **A pre** | Continuous \((y,y')\) + residual at **integrator retreat** (`interpolate` after possible overshoot). Continuous play is at event \(t\). Do **not** re-eval after the jump. | +| **B post-event pre-IC** | After discontinuity, before projector | +| **C post-IC** | After mode 0/1/2/3 path (mode 3 + audit: pure mode 3 even if residual fails) | + +--- + +## Next: (b) Source-current discontinuities + +### Scope (preferred) + +Focus on **electrode / source current** discontinuities and their classical \(t^+\) rates, **not** variable \(C\) / charge-conserving MOD capacitance. + +| In scope | Out of scope (defer) | +|----------|----------------------| +| `IClamp` step (`at_time` + piecewise constant \(i\)) | Density `PROCEDURE dforce` (parked) | +| Sinusoid / continuous \(i(t)\) with kinks or known \(i'\) | Param-in-\(C\) / time-dependent `cm` productization | +| Interaction with **extracellular** (electrode current → `vext` network) | `dcmdt` / charge-jump PP patterns productization | +| How \(i\) and \(i'\) enter residual / free \(y'\) | Thread B (`vext=0` INITIAL abandonment) as primary | + +### Framing already accepted (discussion) + +1. **Built-in source pattern (NMODL)** is already the right shape for steps: + - `ELECTRODE_CURRENT i` (not transmembrane — matters for extracellular). + - `at_time(t_event)` so DASPK reinit sees breakpoints. + - Piecewise `BREAKPOINT` assignment of \(i\) (e.g. `stim.mod`). +2. **Jump in \(i\)** (step on/off): primarily a **\(y\) algebraic** problem after hold — free voltages jump; continuous content held. Residual \(G \approx C y' - f\) absorbs the new \(i\) into \(f\) once \(y\) is on the manifold; classical \(i'\) is zero almost everywhere between steps. +3. **Classical \(i'\)** matters when free algebraic \(y'\) must track a **continuous** source with nonzero derivative (ramps, sinusoids, continuous play into a force term) — same role as Plan A \(b'\) for LM: only the **null(\(C\))** free directions need \(i'\) / \(b'\). +4. **Extracellular** is “harder \(y\)”, not new keywords: electrode current couples into the extracellular network; battery IC already holds \(V_m\) / \(xc>0\) drops. Source (b) should validate that path under IClamp ± `extracellular`, not invent new NMODL. +5. **Do not** treat density MOD `dforce` as the vehicle for (b). That park was for assigned **capacitance rates**. Source work should reuse **play \(t^+\)**, **LM.dforce**, and/or **explicit electrode \(i,i'\) policy** as appropriate. + +### Suggested order of attack (when coding starts) + +Prefer **discussion → minimal spike → tests** over large productization. + +1. **Catalog** how `IClamp` / `at_time` already force reinit; what residual looks like mode 0 vs 3 on a pure istep into a compartment (with/without LM series \(C\)–\(R\)). +2. **Spike:** IClamp step at known \(t\) under mode 3 — residual + physical \(\Delta v\) hold; compare audit panels A/B/C. +3. **Spike:** continuous sinusoid (or play into amp) where \(i'\) is known — does free \(y'\) need an electrode analogue of A2? +4. **Spike:** same + `extracellular` — battery hold vs free `vext` nodes. +5. Only then: thin API if something is missing (do **not** default to density `dforce`). + +### Relation to residual vs “charge conservation” + +- Residual success: \(g(y',y,t)\approx 0\) after IC. +- Physical success: capacitor content continuous; absolute nodes free to jump when \(I\) steps. +- Charge-conserving **parameter jumps in \(C\)** are a **different** problem (parked with density dforce / dcmdt). Do not conflate with electrode \(i\) steps. + +--- + +## Other deferred / optional work + +1. Broader model validation; event-reinit matrix of models. +2. Thread **B:** abandon fixed `vext=0` INITIAL when needed (e.g. forced `e_extracellular`). +3. Policy: when to recommend mode 3 over 0; `cm→0` / ideal clamp as algebraic. +4. SUNDIALS 3 revalidation; do **not** resurrect permanent VMX. +5. Untracked local `docs/index.rst` in worktree — leave alone unless docs build work. + +--- + +## Build + +Typical local build: `~/neuron/nrnida/build/`. + +```bash +cd ~/neuron/nrnida +mkdir -p build && cd build +cmake .. -G Ninja \ + -DCMAKE_BUILD_TYPE=RelWithDebInfo \ + -DNRN_ENABLE_TESTS=ON \ + -DNRN_ENABLE_MPI=OFF +ninja -j$(nproc) +``` + +**Note:** build tree may still contain objects from the parked `mod_dforce` commit (e.g. `mod_dforce.cpp.o`). After switching to A5 tip, a clean rebuild is safer if linking oddities appear: `ninja -t clean` then `ninja`, or reconfigure. + +--- + +## Tests + +```bash +cd build +export PYTHONPATH="$PWD/lib/python${PYTHONPATH:+:$PYTHONPATH}" +export LD_LIBRARY_PATH="$PWD/lib${LD_LIBRARY_PATH:+:$LD_LIBRARY_PATH}" +# Prefer install tree or PATH that matches the built nrnivmodl if tests compile mod: +export PATH="$PWD/bin:$PATH" + +python ../test/hoctests/tests/test_ida_init_mode.py +``` + +Manual: `~/models/nrndc1sim` with `h.cvode.dae_init_mode(3)` before `GUI()` / run. + +--- + +## Starting prompt (new session) + +``` +Read ~/neuron/nrnida/GROK-IDA-INIT.md and, if present, +~/neuron/notes/ida_y_init_adoption.md and ida_phase_R_industry_map.md. + +Repo: ~/neuron/nrnida, branch hines-grok/ida-init @ 7688c6238 (A5). +Parked density dforce: hines-grok/ida-mod-dforce-park @ 210c43c56 — do not revive by default. +Build: prefer existing build/; clean rebuild if parked objects linger. + +Goal: consistent IDA IC after finitialize and discontinuities — +residual ~0 AND physical dt→0 (hold continuous content). + +Done: dae_init_mode 0–3; mode 3 = battery y-hold + y' from C*y'=f + Plan A +forcing t+ (play + LM.dforce) for free y'; audit; A5 stats. + +Next (b): source-current discontinuities (IClamp / sinusoid + at_time, +± extracellular). Discussion-first. Defer param-in-C, charge conservation, +density PROCEDURE dforce. Reuse Plan A framing: i jumps → free y; i' only +for free algebraic y' when classical derivative is nonzero. + +Do not revive permanent VMX/v12. Default stays mode 0 until mode 3 is +broadly validated. +``` + +--- + +## Historical context + +| Branch / path | Lesson | +|---------------|--------| +| `origin/idainit` | `IDA_YA_YDP_INIT` + permanent v12/VMX; LinearCircuit better; extracellular never fully worked; doubled state | +| Pure `IDA_Y_INIT` (mode 1/2) | Algebraic LM OK; folded caps → singular \(\partial g/\partial y\) at \(c_j=0\) | +| Nano-step heuristic | Makes \(\|g\|\) small by finite \(C/\mathrm{d}t\); \(O(\mathrm{dteps}\cdot\|y'\|)\) residual bias (fails tiny \(c_m\)) | +| Mode 3 \(C y'=f\) | Clears residual when \(y\) already on manifold (SEClamp + tiny \(c_m\) scm2eem) | +| Density MOD dforce (`210c43c56`) | Parked; wrong priority vs source currents; capacitance-rate hook | + +Industry map (Phase R): hold differentials / content; solve algebraics; recover \(y'\) from mass equation. diff --git a/docs/progref/modelspec/programmatic/linmod.rst b/docs/progref/modelspec/programmatic/linmod.rst index 9b71bb5bad..ca498604cd 100644 --- a/docs/progref/modelspec/programmatic/linmod.rst +++ b/docs/progref/modelspec/programmatic/linmod.rst @@ -130,7 +130,67 @@ LinearMechanism callable can change the elements of b and g (but do not introduce new elements into g) as a function of time and states. It may be useful for stability and performance to place the linearized part of b into g. - Consider the following pendulum.py with equations + Consider the following pendulum.py with equations + + .. method:: LinearMechanism.dforce + + Syntax: + ``lm.dforce(bdot)`` + + ``lm.dforce(dforce_callable, bdot)`` + + Description: + Supplies :math:`db/dt` for IDA consistent initialization when + :meth:`CVode.dae_init_mode` is 3 (forcing :math:`t^+` info). The + ``bdot`` Vector must have the same size as ``b``. + + With a callable, that function is invoked at each IDA reinit with + ``t`` set to the IC time; it should fill ``bdot``. Without a + callable, ``bdot`` is used as-is (the user may update it before + ``re_init``). + + If neither ``dforce`` nor continuous :meth:`Vector.play` into ``b`` + provides :math:`b'`, but a force callable was passed to the + constructor, a one-sided finite difference of that callable is used + as a fallback. + + Example (series :math:`C`–:math:`R` with sinusoid current into node 0): + + .. code-block:: + python + + import math + from neuron import n + + A, w = 1.0, 2 * math.pi + c = n.Matrix(2, 2) + g = n.Matrix(2, 2) + y = n.Vector(2) + b = n.Vector(2) + bdot = n.Vector(2) + c.setval(0, 0, 1); c.setval(0, 1, -1) + c.setval(1, 0, -1); c.setval(1, 1, 1) + g.setval(1, 1, 1) + + def force(): + b.x[0] = A * math.sin(w * n.t) + + def dforce(): + bdot.x[0] = A * w * math.cos(w * n.t) + + lm = n.LinearMechanism(force, c, g, y, b) + lm.dforce(dforce, bdot) + n.CVode().active(True) + n.CVode().use_daspk(True) + n.CVode().dae_init_mode(3) + n.finitialize(0) + + See also: + :meth:`CVode.dae_init_mode`, :meth:`CVode.dae_init_audit` + + .. + + Description (continued, pendulum): Example: diff --git a/docs/progref/simctrl/cvode.rst b/docs/progref/simctrl/cvode.rst index 09f7982c33..76d428227f 100644 --- a/docs/progref/simctrl/cvode.rst +++ b/docs/progref/simctrl/cvode.rst @@ -2111,6 +2111,242 @@ CVode +.. method:: CVode.dae_init_mode + + .. tab:: Python + + Syntax: + ``mode = cvode.dae_init_mode()`` + + ``mode = cvode.dae_init_mode(mode)`` + + Description: + Selects how the DAE (IDA) integrator obtains a consistent initial + condition after ``finitialize`` or a discontinuity, see + :meth:`CVode.use_daspk`. + + * ``0`` (default) — legacy short fully-implicit step heuristic + controlled by :meth:`CVode.dae_init_dteps`. + * ``1`` — try Sundials ``IDACalcIC`` with ``IDA_Y_INIT`` (solve all + of ``y`` given ``y'``); on failure fall back to the heuristic. + * ``2`` — pure ``IDA_Y_INIT`` only (no heuristic fallback). + * ``3`` — **content hold + continuous** :math:`y'` recovery for + :math:`C(t)\,y' = f(y,t)`: + + 1. **Hold continuous content** and free absolute algebraics: + LinearMechanism floating capacitors → voltage sources holding + branch :math:`\Delta v`; diagonal mass (inductor current) held; + op-amp lag (:math:`\tau v_k'`) holds the differentiated voltage; + extracellular holds :math:`V_m` and capacitive layer drops. + 2. **Recover** :math:`y'` from :math:`C\,y' = f(y,t)` at that + fixed :math:`y` (diagonal membrane/mechanism mass, simple LM + mass stamps). This is the continuous :math:`dt\to 0` limit and + does **not** use :meth:`CVode.dae_init_dteps`. + + On residual failure, prints a short equation classification + (algebraic vs near-singular :math:`c`) and falls back to the + mode-``0`` heuristic (unless an IC audit is armed, which keeps + pure mode-``3`` state for diagnosis). + + ``IDA_Y_INIT`` does not require a differential/algebraic ``id`` vector. + It works well for invertible algebraic LinearMechanism systems. + For folded capacitor equations (``C*(v1'-v2')``) the Newton matrix + can be singular under ``IDA_Y_INIT``; mode ``1`` then falls back to + the heuristic. Default remains ``0`` until mode ``3`` is fully + validated across models. + + Mode ``3`` assumes the post-event (or ``finitialize``) :math:`y` is + already on the algebraic manifold for rows with :math:`C=0`. It + does not invent a consistent :math:`y` when user/``INITIAL`` state + is algebraically inconsistent (e.g. clamp amp ≠ :math:`v`, or + ``xc=0`` with an inconsistent ``vext`` gauge). + + **Electrode / source currents** (e.g. :class:`IClamp`, continuous + :meth:`Vector.play` into a force, piecewise-linear test stimuli): + + * Use :keyword:`ELECTRODE_CURRENT` and :func:`at_time` breakpoints + for step sources so IDA reinit lands on the event. + * A **jump** in electrode current is primarily an algebraic + :math:`y` problem after hold: free absolute levels (and free + end-node voltages when a point process sits at location ``0`` or + ``1``); continuous content (:math:`V_m`, floating :math:`\Delta v`, + capacitive :math:`xc` drops) is held. Classical :math:`i'` is not + required for residual success between flat segments. + * A **kink** or ramp (nonzero classical :math:`i'`) matters for free + algebraic rates when mass leaves a null space (e.g. series + :math:`C`–:math:`R` with play into :math:`b`). Prefer continuous + :meth:`Vector.play` or :meth:`LinearMechanism.dforce` for the + 1-jet of the drive; mode ``3`` reuses that Plan-A path. + * With **extracellular**, electrode current is not transmembrane: + it couples into the ``vext`` network. Mode ``3`` holds membrane + :math:`V_m` (``seg.v``) and capacitive layer drops; absolute + ``vext`` may jump. For coupled membrane/``xc`` mass, rates + include electrode current in the :math:`C y' = f` seed. + * **Default remains mode ``0``** until you validate mode ``3`` on + your models. Use :meth:`CVode.dae_init_stats` and + :meth:`CVode.dae_init_audit` when diagnosing reinits. + + **Forcing** :math:`t^+` **info:** the right-limit value + :math:`u(t^+)` and classical derivative :math:`u'(t^+)` of + exogenous drives after a discontinuity (or at ``finitialize``). + In the geometric / DAE literature this pair is the **1-jet** of + :math:`u` at :math:`t^+`. Continuous :meth:`Vector.play` supplies + this from piecewise-linear samples (outgoing segment slope at a + knot; linear extrapolation of the last two points past the end of + the ``t`` vector — see :meth:`Vector.play`). At each IDA reinit, + those plays are sampled and listed in :meth:`CVode.dae_init_audit` + under ``forcing t+ info``. Mode ``3`` uses the sample for + LinearMechanism free :math:`y'` (common mode of floating + capacitors, etc.): after seeding :math:`C y' = b - G y`, free + directions are adjusted so differentiated algebraics hold when + :math:`b'` is known from continuous play (e.g. series + :math:`C`–:math:`R` ramp: :math:`V_R' = R I'`) or from + :meth:`LinearMechanism.dforce` (analytic :math:`b'`, or a finite- + difference fallback when only the force callable is provided). + + .. tab:: HOC + + Syntax: + ``mode = cvode.dae_init_mode()`` + + ``mode = cvode.dae_init_mode(mode)`` + + Description: + Same as the Python tab: ``0`` heuristic (default), ``1`` + ``IDA_Y_INIT`` with heuristic fallback, ``2`` pure ``IDA_Y_INIT``, + ``3`` battery content hold plus :math:`C y' = f(y)` for :math:`y'`, + including electrode/source and extracellular notes above. + +---- + + + +.. method:: CVode.dae_init_stats + + .. tab:: Python + + Syntax: + ``n = cvode.dae_init_stats()`` + + ``cvode.dae_init_stats(1)`` + + ``n = cvode.dae_init_stats(vec)`` + + Description: + IDA consistent-initialization path counters (Plan A5). + + * No argument — print a short summary (also included in + :meth:`CVode.statistics`) and return the number of IDA reinits. + * ``1`` — reset all IC counters. + * ``vec`` — fill a :class:`Vector` of length 8: + + 0. total IDA reinits + 1. mode 3 successes + 2. mode 3 residual failures that fell back to nano-step + 3. reinits that applied free :math:`y'` from continuous play + 4. reinits that applied free :math:`y'` from :meth:`LinearMechanism.dforce` + 5. reinits that applied free :math:`y'` from FD of the force callable + 6. last IC path mode (0 heuristic, 1/2 CalcIC, 3 battery+forcing) + 7. last IC forcing source flags (bit 1=play, 2=dforce, 4=fd, 8=applied) + + .. tab:: HOC + + Same as Python: ``cvode.dae_init_stats()``, ``cvode.dae_init_stats(1)``, + or ``cvode.dae_init_stats(vec)``. + +---- + + + +.. method:: CVode.dae_init_audit + + .. tab:: Python + + Syntax: + ``level = cvode.dae_init_audit()`` + + ``level = cvode.dae_init_audit(level)`` + + ``level = cvode.dae_init_audit(level, t)`` + + Description: + Diagnostic **three-panel audit** of DAE (IDA) consistent initialization + after ``finitialize`` or a discontinuity reinit. Intended for + development: after a run where something looks wrong near a known + time, re-run with the audit armed at that time. + + * ``level = 0`` — off (default). + * ``level = 1`` — summary only (WRMS / max residual A·B·C, jump sizes). + * ``level = 2`` — summary plus top residual rows per panel + (``y``, ``y'``, residual ≈ ``c*y' - f(y)``). + + With the optional second argument ``t``, the next reinit with + simulation time ``>= t`` produces one dump, then the audit disarms + (one-shot). Typical workflow: notice an issue near time ``T``, then + ``cvode.dae_init_audit(2, T)`` and re-run. Output goes to stdout, or + to a file set by :meth:`CVode.dae_init_audit_file`. + + Panels: + + * **A pre** — continuous ``(y, y')`` and residual at the integrator + **retreat** to the discontinuity time (``interpolate`` after a + step that may overshoot). Continuous play is synchronized to that + time when the residual is captured. Prefer this over the raw step + endpoint. Unavailable at the first ``finitialize``. + * **B post-event pre-IC** — state after the discontinuity and before + the IC projector; residual is usually large on affected equations. + * **C post-IC** — after heuristic / ``IDA_Y_INIT`` / mode-3 + (content hold + :math:`C y'=f`). When mode ``3`` fails residual + and an audit is armed, fallback is suppressed so panel C shows + pure mode ``3``. + + Here “event” means any discontinuity reinit (``NET_RECEIVE``, + ``Vector.play``, ``at_time``, clamps, etc.), not only network events. + + .. tab:: HOC + + Syntax: + ``level = cvode.dae_init_audit()`` + + ``level = cvode.dae_init_audit(level)`` + + ``level = cvode.dae_init_audit(level, t)`` + + Description: + Same as the Python tab. + +---- + + + +.. method:: CVode.dae_init_audit_file + + .. tab:: Python + + Syntax: + ``cvode.dae_init_audit_file()`` + + ``cvode.dae_init_audit_file(path)`` + + Description: + Destination for :meth:`CVode.dae_init_audit` text. With no argument + or an empty string, write to stdout. With a path, append each audit + dump to that file. Returns ``1`` if a file path is set, else ``0``. + + .. tab:: HOC + + Syntax: + ``cvode.dae_init_audit_file()`` + + ``cvode.dae_init_audit_file("path")`` + + Description: + Same as the Python tab. + +---- + + + .. method:: CVode.dae_init_dteps .. tab:: Python @@ -2128,7 +2364,8 @@ CVode The size of the "infinitesimal" fixed fully implicit step used for initialization of the DAE solver, see :func:`use_daspk` , in order to meet the the initial condition requirement of f(y',y,t)=0. The default - is 1e-9 ms. + is 1e-9 ms. Used when :meth:`CVode.dae_init_mode` is ``0``, or as the + fallback when mode is ``1``. The default heuristic for meeting the initial condition requirement based on the pre-initialization value of all the states and an initialization time @@ -2194,8 +2431,9 @@ CVode The size of the "infinitesimal" fixed fully implicit step used for initialization of the DAE solver, see :func:`use_daspk` , in order to meet the the initial condition requirement of f(y',y,t)=0. The default - is 1e-9 ms. + is 1e-9 ms. Also used as fallback when :meth:`CVode.dae_init_mode` is 1. + The default heuristic for meeting the initial condition requirement based on the pre-initialization value of all the states and an initialization time diff --git a/src/nrncvode/cvodeobj.cpp b/src/nrncvode/cvodeobj.cpp index 228b7ac2a5..133b068088 100644 --- a/src/nrncvode/cvodeobj.cpp +++ b/src/nrncvode/cvodeobj.cpp @@ -20,6 +20,7 @@ extern void (*nrn_multisplit_setup_)(); #include "nrndaspk.h" #include "nrniv_mf.h" #include "nrnpy.h" +#include "ivocvect.h" #include "tqueue.hpp" #include "mymath.h" #include @@ -323,6 +324,69 @@ static double dae_init_dteps(void* v) { return Daspk::dteps_; } +// 0 = heuristic only (default); 1 = IDA_Y_INIT then heuristic fallback; +// 2 = IDA_Y_INIT only; 3 = battery content hold + y' from C*y'=f(y). +static double dae_init_mode(void* v) { + hoc_return_type_code = HocReturnType::integer; + if (ifarg(1)) { + Daspk::init_mode_ = (int) chkarg(1, 0, 3); + } + return (double) Daspk::init_mode_; +} + +// Three-panel IDA IC audit: 0 off, 1 summary, 2 residual rows. +// Optional second arg t arms first reinit with time >= t (one-shot). +static double dae_init_audit(void* v) { + hoc_return_type_code = HocReturnType::integer; + if (ifarg(1)) { + int level = (int) chkarg(1, 0, 2); + Daspk::audit_set_level(level); + if (ifarg(2)) { + Daspk::audit_arm_at(*getarg(2)); + } + } + return (double) Daspk::audit_level(); +} + +// Append audit text to path; no arg or empty string → stdout. +static double dae_init_audit_file(void* v) { + if (ifarg(1)) { + Daspk::audit_set_file(gargstr(1)); + } else { + Daspk::audit_set_file(nullptr); + } + return Daspk::audit_path_.empty() ? 0. : 1.; +} + +// A5: IDA IC path statistics. +// dae_init_stats() — print summary (same as part of statistics()) +// dae_init_stats(1) — reset counters +// dae_init_stats(vec) — fill Vector: [n_init, mode3_ok, mode3_fb, +// play, dforce, fd, last_path, last_flags] +static double dae_init_stats(void* v) { + if (ifarg(1)) { + if (hoc_is_object_arg(1) && is_vector_arg(1)) { + IvocVect* vec = vector_arg(1); + vec->resize(8); + vec->elem(0) = (double) Daspk::ic_init_count_; + vec->elem(1) = (double) Daspk::ic_mode3_ok_count(); + vec->elem(2) = (double) Daspk::ic_mode3_fallback_count(); + vec->elem(3) = (double) Daspk::ic_forcing_play_inits_; + vec->elem(4) = (double) Daspk::ic_forcing_dforce_inits_; + vec->elem(5) = (double) Daspk::ic_forcing_fd_inits_; + vec->elem(6) = (double) Daspk::last_ic_path_mode(); + vec->elem(7) = (double) Daspk::last_ic_forcing_flags(); + return 8.; + } + if (chkarg(1, 0, 1) == 1.) { + Daspk::reset_ic_stats(); + return 0.; + } + } + Daspk::print_ic_stats(); + return (double) Daspk::ic_init_count_; +} + static double use_mxb(void* v) { hoc_return_type_code = HocReturnType::boolean; if (ifarg(1)) { @@ -604,6 +668,10 @@ static Member_func members[] = {{"solve", solve}, {"store_events", store_events}, {"condition_order", condition_order}, {"dae_init_dteps", dae_init_dteps}, + {"dae_init_mode", dae_init_mode}, + {"dae_init_audit", dae_init_audit}, + {"dae_init_audit_file", dae_init_audit_file}, + {"dae_init_stats", dae_init_stats}, {"simgraph_remove", simgraph_remove}, {"state_magnitudes", state_magnitudes}, {"ncs_netcons", ncs_netcons}, @@ -653,7 +721,13 @@ static void destruct(void* v) { void Cvode_reg() { class2oc("CVode", cons, destruct, members, omembers, nullptr); net_cvode_instance = new NetCvode(1); - Daspk::dteps_ = 1e-9; // change with cvode.dae_init_dteps(newval) + Daspk::dteps_ = 1e-9; // change with cvode.dae_init_dteps(newval) + Daspk::init_mode_ = 0; // heuristic IC; use dae_init_mode(1|2|3) for other paths + Daspk::audit_level_ = 0; + Daspk::audit_armed_ = 0; + Daspk::audit_serial_ = 0; + Daspk::audit_path_.clear(); + Daspk::reset_ic_stats(); } /* Functions Called by the CVODE Solver */ diff --git a/src/nrncvode/nrndaspk.cpp b/src/nrncvode/nrndaspk.cpp index 48fa792ec7..40cc7688c5 100644 --- a/src/nrncvode/nrndaspk.cpp +++ b/src/nrncvode/nrndaspk.cpp @@ -6,12 +6,17 @@ #include #include +#include +#include +#include #include "spmatrix.h" #include "nrnoc2iv.h" #include "cvodeobj.h" #include "nrndaspk.h" #include "netcvode.h" #include "nrn_ansi.h" +#include "vecplay_tplus.h" +#include "nrndae.h" #include "ida/ida.h" #include "ida/ida_impl.h" #include "mymath.h" @@ -30,8 +35,18 @@ double Daspk::dteps_; extern void nrndae_dkres(double*, double*, double*); extern void nrndae_dkpsol(double); +extern int nrndae_battery_ic_project(); +extern void nrndae_seed_yp_from_f(double* f, double* yp); +extern int nrndae_complete_yp_from_forcing(double* yp, const std::vector& forcing); +extern void nrndae_append_dforce_to_forcing_list(double tt, std::vector& out); extern void nrn_solve(NrnThread*); +extern int nrn_sparse13_soft_fail; +extern int nrn_sparse13_factor_error(); void nrn_daspk_init_step(double, double, int); +void nrn_cable_battery_ic(); +#if EXTRACELLULAR +void nrn_extracellular_battery_ic(); +#endif // this is private in ida.cpp but we want to check if our initialization // is good. Unfortunately ewt is set on the first call to solve which // is too late for us. @@ -129,7 +144,12 @@ static int msolve(IDAMem mem, N_Vector b, N_Vector w, N_Vector ycur, N_Vector, N nvec_y = ycur; nvec_yp = b; thread_cj = mem->ida_cj; + thread_ier = 0; nrn_multithread_job(msolve_thread); + if (nrn_sparse13_factor_error()) { + // Recoverable linear solve failure (e.g. singular J during IDA_Y_INIT). + return 1; + } return thread_ier; } @@ -146,6 +166,12 @@ Daspk::Daspk(Cvode* cv, int neq) { use_parasite_ = false; spmat_ = nullptr; mem_ = nullptr; + audit_pre_t_ = 0.; + audit_pre_max_abs_ = 0.; + audit_pre_wrms_ = -1.; + audit_pre_neq_ = 0; + audit_pre_valid_ = false; + audit_pre_res_valid_ = false; } Daspk::~Daspk() { @@ -190,6 +216,80 @@ void Daspk::info() {} int Daspk::init_failure_style_; int Daspk::init_try_again_; int Daspk::first_try_init_failures_; +int Daspk::init_mode_; +int Daspk::calcic_fallback_count_; +int Daspk::ic_init_count_ = 0; +int Daspk::ic_mode3_ok_count_ = 0; +int Daspk::ic_mode3_fallback_count_ = 0; +int Daspk::ic_forcing_play_inits_ = 0; +int Daspk::ic_forcing_dforce_inits_ = 0; +int Daspk::ic_forcing_fd_inits_ = 0; +int Daspk::last_ic_path_mode_ = -1; +int Daspk::last_ic_forcing_flags_ = 0; +int Daspk::audit_level_ = 0; +double Daspk::audit_t_select_ = 0.; +int Daspk::audit_armed_ = 0; +int Daspk::audit_serial_ = 0; +std::string Daspk::audit_path_; +std::vector Daspk::last_forcing_tplus_; +double Daspk::last_forcing_t_ = 0.; + +const std::vector& Daspk::last_forcing_tplus() { + return last_forcing_tplus_; +} + +double Daspk::last_forcing_t() { + return last_forcing_t_; +} + +void Daspk::reset_ic_stats() { + ic_init_count_ = 0; + ic_mode3_ok_count_ = 0; + ic_mode3_fallback_count_ = 0; + ic_forcing_play_inits_ = 0; + ic_forcing_dforce_inits_ = 0; + ic_forcing_fd_inits_ = 0; + calcic_fallback_count_ = 0; + last_ic_path_mode_ = -1; + last_ic_forcing_flags_ = 0; +} + +int Daspk::last_ic_path_mode() { + return last_ic_path_mode_; +} +int Daspk::last_ic_forcing_flags() { + return last_ic_forcing_flags_; +} +int Daspk::ic_mode3_ok_count() { + return ic_mode3_ok_count_; +} +int Daspk::ic_mode3_fallback_count() { + return ic_mode3_fallback_count_; +} + +void Daspk::print_ic_stats() { + if (!ic_init_count_ && !calcic_fallback_count_ && !first_try_init_failures_) { + return; + } + Printf(" IDA IC: %d reinit(s)", ic_init_count_); + if (ic_mode3_ok_count_ || ic_mode3_fallback_count_) { + Printf("; mode3 ok=%d fallback=%d", ic_mode3_ok_count_, ic_mode3_fallback_count_); + } + if (ic_forcing_play_inits_ || ic_forcing_dforce_inits_ || ic_forcing_fd_inits_) { + Printf("; free y' from play=%d dforce=%d fd=%d", + ic_forcing_play_inits_, + ic_forcing_dforce_inits_, + ic_forcing_fd_inits_); + } + Printf("\n"); + if (calcic_fallback_count_) { + Printf(" %d IDA IC mode 1/3 residual failure(s) fell back to heuristic\n", + calcic_fallback_count_); + } + if (first_try_init_failures_) { + Printf(" %d First try Initialization failures\n", first_try_init_failures_); + } +} static void do_ode_thread(neuron::model_sorted_token const& sorted_token, NrnThread& ntr) { auto* const nt = &ntr; @@ -221,76 +321,20 @@ N_VGetArrayPointer(ida->delta_)[i]); } #endif -int Daspk::init() { +int Daspk::check_init_residual() { extern double t; -#if 0 -printf("Daspk_init t_=%20.12g t-t_=%g t0_-t_=%g\n", -cv_->t_, t-cv_->t_, cv_->t0_-cv_->t_); -#endif - N_VConst(0., yp_); - - // the new initial condition is based on a dteps_ step backward euler - // linear solution with respect to the old state in order to - // start the following initial condition calculation with a "valid" - // (in a linear system sense) initial state. - - double tt = cv_->t_; - double dtinv = 1. / dteps_; - if (init_failure_style_ & 010) { - cv_->play_continuous(tt); - nrn_daspk_init_step(tt, dteps_, 1); - nrn_daspk_init_step(tt, dteps_, 1); - cv_->daspk_gather_y(yp_); - cv_->play_continuous(tt); - nrn_daspk_init_step(tt, dteps_, 1); - cv_->daspk_gather_y(cv_->y_); - N_VLinearSum(dtinv, cv_->y_, -dtinv, yp_, yp_); - } else { -#if 0 - cv_->play_continuous(tt); - nrn_daspk_init_step(tt, dteps_, 1); - cv_->daspk_gather_y(cv_->y_); - tt = cv_->t_ + dteps_; - cv_->play_continuous(tt); - nrn_daspk_init_step(tt, dteps_, 1); - cv_->daspk_gather_y(yp_); - N_VLinearSum(dtinv, yp_, -dtinv, cv_->y_, yp_); - cv_->daspk_scatter_y(cv_->y_); -#else - cv_->play_continuous(tt); - nrn_daspk_init_step(tt, dteps_, 1); // first approx to y (and maybe good enough) - nrn_daspk_init_step(tt, dteps_, 1); // 2nd approx to y (trouble with 2sramp.hoc) - - cv_->daspk_gather_y(cv_->y_); - tt = cv_->t_ + dteps_; - cv_->play_continuous(tt); - nrn_daspk_init_step(tt, dteps_, 0); // rhs contains delta y (for v, vext, linmod - cv_->gather_ydot(yp_); - N_VScale(dtinv, yp_, yp_); -#endif - } - thread_cv = cv_; - nvec_yp = yp_; - nrn_multithread_job(nrn_ensure_model_data_are_sorted(), do_ode_thread); - ida_init(); t = cv_->t_; -#if 1 - // test - // printf("test\n"); if (!IDAEwtSet((IDAMem) mem_, cv_->y_)) { hoc_execerror("Bad Ida error weight vector", 0); } use_parasite_ = false; - // check(cv_->t_, this); res_gvardt(cv_->t_, cv_->y_, yp_, parasite_, cv_); double norm = N_VWrmsNorm(parasite_, ((IDAMem) mem_)->ida_ewt); - // printf("norm=%g at t=%g\n", norm, t); if (norm > 1.) { switch (init_failure_style_ & 03) { case 0: Printf("IDA initialization failure, weighted norm of residual=%g\n", norm); return IDA_ERR_FAIL; - break; case 1: Printf("IDA initialization warning, weighted norm of residual=%g\n", norm); break; @@ -301,11 +345,6 @@ cv_->t_, t-cv_->t_, cv_->t0_-cv_->t_); Printf(" subtracting (for next 1e-6 ms): f(y', y, %g)*exp(-1e7*(t-%g))\n", nt_t, nt_t); break; } -#if 0 -for (int i=0; i < cv_->neq_; ++i) { - printf(" %d %g %g %g %g\n", i, nt_t, N_VGetArrayPointer(cv_->y_)[i], N_VGetArrayPointer(yp_)[i], N_VGetArrayPointer(delta_)[i]); -} -#endif if (init_try_again_ < 0) { ++first_try_init_failures_; init_try_again_ += 1; @@ -313,11 +352,787 @@ for (int i=0; i < cv_->neq_; ++i) { init_try_again_ = 0; return err; } + // style 1 or 2: accept with warning return 0; } + return 0; +} + +// Fill y_ and yp_ using the legacy nano-step heuristic (no residual check). +static void seed_y_yp_heuristic(Daspk* d) { + Cvode* cv = d->cv_; + N_Vector yp = d->yp_; + N_VConst(0., yp); + + double tt = cv->t_; + double dtinv = 1. / Daspk::dteps_; + if (Daspk::init_failure_style_ & 010) { + cv->play_continuous(tt); + nrn_daspk_init_step(tt, Daspk::dteps_, 1); + nrn_daspk_init_step(tt, Daspk::dteps_, 1); + cv->daspk_gather_y(yp); + cv->play_continuous(tt); + nrn_daspk_init_step(tt, Daspk::dteps_, 1); + cv->daspk_gather_y(cv->y_); + N_VLinearSum(dtinv, cv->y_, -dtinv, yp, yp); + } else { + cv->play_continuous(tt); + nrn_daspk_init_step(tt, Daspk::dteps_, 1); + nrn_daspk_init_step(tt, Daspk::dteps_, 1); + + cv->daspk_gather_y(cv->y_); + tt = cv->t_ + Daspk::dteps_; + cv->play_continuous(tt); + nrn_daspk_init_step(tt, Daspk::dteps_, 0); + cv->gather_ydot(yp); + N_VScale(dtinv, yp, yp); + } + thread_cv = cv; + nvec_yp = yp; + nrn_multithread_job(nrn_ensure_model_data_are_sorted(), do_ode_thread); +} + +int Daspk::init_heuristic() { + extern double t; + seed_y_yp_heuristic(this); + ida_init(); + t = cv_->t_; + return check_init_residual(); +} + +int Daspk::init_ida_y_init() { + extern double t; + // Seed: voltage/DAE yp = 0; mechanism ODE yp = f(y) via do_ode. + // Suitable when residual uniquely determines y at that yp (e.g. pure + // algebraic LinearMechanism with invertible g). For folded capacitors + // C*(v1'-v2'), dF/dy is singular when cj=0 (IDA_Y_INIT), so CalcIC may + // soft-fail and mode 1 falls back to the heuristic. + N_VConst(0., yp_); + double tt = cv_->t_; + cv_->play_continuous(tt); + cv_->daspk_gather_y(cv_->y_); + thread_cv = cv_; + nvec_yp = yp_; + nrn_multithread_job(nrn_ensure_model_data_are_sorted(), do_ode_thread); + ida_init(); + t = cv_->t_; + + // tout1 only sets integration direction / rough t scale for IDACalcIC. + realtype tout1 = tt + 1.0; + if (tout1 == tt) { + tout1 = tt + 1e-3; + } + nrn_sparse13_soft_fail = 1; + int ier = IDACalcIC(mem_, IDA_Y_INIT, tout1); + nrn_sparse13_soft_fail = 0; + if (ier != IDA_SUCCESS) { + Printf("IDACalcIC(IDA_Y_INIT) failed, err=%d\n", ier); + return ier; + } + // Corrected y is already in cv_->y_ (IDA y0); scatter into NEURON structures. + cv_->daspk_scatter_y(cv_->y_); + t = cv_->t_; + nt_t = cv_->t_; + return check_init_residual(); +} + +// Mode 3 step after y-project: set yp from C*yp = f(y) at fixed y (continuous +// limit). Residual G = C*yp - f (same assembly as Cvode::res). Unlike the +// dteps companion nano-step, there is no O(dteps) bias in y'. +// +// Handles: diagonal membrane cm, identity mechanism ODEs, simple LM mass +// (via nrndae_seed_yp_from_f), and single-layer xc when present. Algebraic +// rows (c=0) leave yp=0 — residual can only clear if f already vanishes. +static void seed_yp_from_Cy_eq_f(Daspk* d) { + Cvode* cv = d->cv_; + NrnThread* nt = nrn_threads; + CvodeThreadData& z = cv->ctd_[0]; + double tt = cv->t_; + nt->_t = tt; + cv->play_continuous(tt); + auto const sorted = nrn_ensure_model_data_are_sorted(); + // f(y) for G = C*yp - f (same first half as Cvode::res) + nrn_rhs(sorted, *nt); + cv->do_ode(sorted, *nt); + cv->gather_ydot(d->delta_); + double* F = cv->n_vector_data(d->delta_, 0); + double* yp = cv->n_vector_data(d->yp_, 0); + + N_VConst(0., d->yp_); + +#if EXTRACELLULAR + // Extracellular + membrane mass are coupled (see Cvode::res): + // delta[vi] = F[vi] - c*(yp_vi - yp_vx0) + // delta[vx0] = F[vx0] + c*(yp_vi - yp_vx0) - (xc chain) + // Electrode current sits only in F[vi]. Solving C yp = F then requires + // c*(yp_vi - yp_vx0) = F[vi] + // and for 1-layer xc: cx*yp_vx0 = F[vi] + F[vx0] + // so yp_vx0 = (F[vi]+F[vx0])/cx, not F[vx0]/cx alone (the old seed left + // residual ±I_electrode on the vext row after a source step). + if (z.cmlext_) { + assert(z.cmlext_->ml.size() == 1); + Memb_list* mlx = &z.cmlext_->ml[0]; + Memb_list* mlc = (z.cmlcap_ && z.cmlcap_->ml.size() == 1) ? &z.cmlcap_->ml[0] : nullptr; + int n = mlx->nodecount; + for (int i = 0; i < n; ++i) { + Node* nd = mlx->nodelist[i]; + int jx = nd->eqn_index_; // vext[0] + int jv = nd->eqn_index_ - 1; // vi + double c = 0.; + if (mlc) { + // same node order as cmlcap list when both present + for (int k = 0; k < mlc->nodecount; ++k) { + if (mlc->nodelist[k] == nd) { + c = 1e-3 * mlc->data(k, 0); + break; + } + } + } + if (nrn_nlayer_extracellular == 1) { + double cx = 1e-3 * mlx->data(i, neuron::extracellular::xc_index, 0); + if (cx > 0.) { + if (c > 0.) { + // Coupled cm+xc: electrode in F[vi] must charge both + // (yp_vx = (F_vi+F_vx)/cx, yp_vi = yp_vx + F_vi/c) + yp[jx] = (F[jv] + F[jx]) / cx; + yp[jv] = yp[jx] + F[jv] / c; + } else { + // No membrane mass (e.g. zero-area node): no cm coupling + // in res; use diagonal xc only. Algebraic vi needs F[vi]=0. + yp[jx] = F[jx] / cx; + } + } else if (c > 0.) { + // xc algebraic, membrane capacitive: only relative rate from F[vi] + yp[jv] = F[jv] / c; + } + } else { + // Multi-layer (default nlayer is often 2). res couples membrane + // only to vext[0], and layer caps in a chain to ground: + // c*(yp_vi-yp0)=F[vi] + // cx_k*(yp[k]-yp[k+1]) carries F[vi]+F[0]+...+F[k] toward ground + // so the outermost rate is + // yp[last]=(F[vi]+sum_k F[vext[k]])/cx_last + // then walk inward. (Old seed used F[layer]/cx alone and left + // electrode residual on outer rows.) + int nlay = nrn_nlayer_extracellular; + // cumulative F from membrane + all layers + double f_cum = F[jv]; + for (int k = 0; k < nlay; ++k) { + f_cum += F[jx + k]; + } + // outermost + int k = nlay - 1; + int jj = jx + k; + double cx = 1e-3 * mlx->data(i, neuron::extracellular::xc_index, k); + if (cx > 0.) { + yp[jj] = f_cum / cx; + } + // remove this layer's F and step inward + for (k = nlay - 2; k >= 0; --k) { + f_cum -= F[jx + k + 1]; + jj = jx + k; + cx = 1e-3 * mlx->data(i, neuron::extracellular::xc_index, k); + if (cx > 0.) { + yp[jj] = yp[jj + 1] + f_cum / cx; + } + } + if (c > 0.) { + yp[jv] = yp[jx] + F[jv] / c; + } + } + } + } else #endif + // Capacitive membrane without extracellular: c*vm' = f + if (z.cmlcap_) { + assert(z.cmlcap_->ml.size() == 1); + Memb_list* ml = &z.cmlcap_->ml[0]; + int n = ml->nodecount; + for (int i = 0; i < n; ++i) { + Node* nd = ml->nodelist[i]; + int j = nd->eqn_index_ - 1; + double c = 1e-3 * ml->data(i, 0); + if (c == 0.) { + continue; // algebraic membrane node + } + yp[j] = F[j] / c; + } + } - return 0; + // Mechanism ODEs: identity mass y' = f_ode + for (int i = z.neq_v_; i < z.nvsize_; ++i) { + yp[i] = F[i]; + } + + // LinearMechanism mass (diagonal / lag / simple floating difference) + nrndae_seed_yp_from_f(F, yp); + + // A2/A4: free y' in null(C) from play / dforce / FD (db/dt). + // Uses Daspk::last_forcing_tplus_ collected at the start of Daspk::init. + const int fflags = nrndae_complete_yp_from_forcing(yp, Daspk::last_forcing_tplus()); + Daspk::last_ic_forcing_flags_ = fflags; + if (fflags & NRN_IC_FORCING_PLAY) { + ++Daspk::ic_forcing_play_inits_; + } + if (fflags & NRN_IC_FORCING_DFORCE) { + ++Daspk::ic_forcing_dforce_inits_; + } + if (fflags & NRN_IC_FORCING_FD) { + ++Daspk::ic_forcing_fd_inits_; + } +} + +// After residual failure: classify largest residual equations (algebraic vs +// tiny-cm differential) to aid singular / inconsistent IC diagnosis. +static void diagnose_battery_ic_residual(Daspk* d) { + Cvode* cv = d->cv_; + if (cv->neq_ <= 0) { + return; + } + // Residual left in parasite_ by check_init_residual + double* r = cv->n_vector_data(d->parasite_, 0); + double* yp = cv->n_vector_data(d->yp_, 0); + CvodeThreadData& z = cv->ctd_[0]; + + // Per-eq diagonal c estimate for voltage equations (0 = algebraic) + std::vector cdiag(cv->neq_, 0.); + for (int i = z.neq_v_; i < z.nvsize_ && i < cv->neq_; ++i) { + cdiag[i] = 1.; // mechanism ODE identity + } + if (z.cmlcap_) { + Memb_list* ml = &z.cmlcap_->ml[0]; + for (int i = 0; i < ml->nodecount; ++i) { + Node* nd = ml->nodelist[i]; + int j = nd->eqn_index_ - 1; + if (j >= 0 && j < cv->neq_) { + cdiag[j] = 1e-3 * ml->data(i, 0); + } + } + } + + // Top few residual equations by |r| + constexpr int ntop = 5; + int idx[ntop]; + double ar[ntop]; + for (int k = 0; k < ntop; ++k) { + idx[k] = -1; + ar[k] = -1.; + } + for (int i = 0; i < cv->neq_; ++i) { + double a = std::fabs(r[i]); + for (int k = 0; k < ntop; ++k) { + if (a > ar[k]) { + for (int m = ntop - 1; m > k; --m) { + ar[m] = ar[m - 1]; + idx[m] = idx[m - 1]; + } + ar[k] = a; + idx[k] = i; + break; + } + } + } + + constexpr double tiny_c = 1e-12; // 1e-3 * cm with cm ~ 1e-9 uF/cm2 scale + Printf(" battery IC residual diagnosis (top residual eqs):\n"); + for (int k = 0; k < ntop && idx[k] >= 0; ++k) { + if (ar[k] <= 0.) { + break; + } + int i = idx[k]; + double c = (i < (int) cdiag.size()) ? cdiag[i] : 0.; + const char* kind = "unknown"; + if (i >= z.neq_v_) { + kind = "mechanism-ODE"; + } else if (c == 0.) { + kind = "algebraic (c=0): y may need to change; y' cannot clear residual"; + } else if (c < tiny_c) { + kind = "near-singular c: huge |y'| or treat as algebraic / ideal clamp"; + } else { + kind = "capacitive"; + } + Printf(" eq %d |res|=%.6g c~%.3g y'=%.6g — %s\n", i, ar[k], c, yp[i], kind); + } +} + +int Daspk::init_battery() { + extern double t; + // Mode 3: (1) hold continuous content (LM + extracellular), (2) y' from + // C*y' = f(y) at fixed y. Nano-step heuristic is only a fallback (caller). + double tt = cv_->t_; + cv_->play_continuous(tt); + // Sync Node <-> IDA y (vi/vext transform) before projectors read Node.v. + // Also evaluate residual once so POINT_PROCESS BREAKPOINT (IClamp/PWL) + // and play side effects match the post-event t+ world before project. + cv_->daspk_gather_y(cv_->y_); + N_VConst(0., yp_); + res_gvardt(tt, cv_->y_, yp_, delta_, cv_); + cv_->daspk_scatter_y(cv_->y_); + int berr = nrndae_battery_ic_project(); + if (berr != 0) { + Printf("nrndae_battery_ic_project failed, err=%d\n", berr); + return berr; + } + // Free zero-area cable nodes (electrode at 0/1); hold CAP-list voltages. + nrn_cable_battery_ic(); +#if EXTRACELLULAR + // Hold Vm on CAP nodes and xc content; free zero-area Vm. + nrn_extracellular_battery_ic(); +#endif + // Projected states → IDA N_Vector (vi,vext transform in gather). + cv_->daspk_gather_y(cv_->y_); + cv_->daspk_scatter_y(cv_->y_); + seed_yp_from_Cy_eq_f(this); + + ida_init(); + t = cv_->t_; + nt_t = cv_->t_; + cv_->daspk_gather_y(cv_->y_); + cv_->daspk_scatter_y(cv_->y_); + int err = check_init_residual(); + if (err != 0) { + diagnose_battery_ic_residual(this); + } + return err; +} + +void Daspk::audit_set_level(int level) { + audit_level_ = level; + if (level <= 0) { + audit_armed_ = 0; + } +} + +int Daspk::audit_level() { + return audit_level_; +} + +void Daspk::audit_arm_at(double t) { + audit_t_select_ = t; + audit_armed_ = 1; +} + +double Daspk::audit_t_select() { + return audit_t_select_; +} + +void Daspk::audit_set_file(const char* path) { + if (!path || !path[0]) { + audit_path_.clear(); + } else { + audit_path_ = path; + } +} + +const char* Daspk::audit_file() { + return audit_path_.c_str(); +} + +FILE* Daspk::audit_open_out() { + if (audit_path_.empty()) { + return stdout; + } + FILE* f = fopen(audit_path_.c_str(), "a"); + if (!f) { + Printf("dae_init_audit: cannot open '%s' for append; using stdout\n", audit_path_.c_str()); + return stdout; + } + return f; +} + +static void audit_copy_nv_out(Cvode* cv, N_Vector nv, std::vector& out) { + out.resize(cv->neq_); + int k = 0; + for (int tid = 0; tid < cv->nctd_; ++tid) { + double* p = cv->n_vector_data(nv, tid); + int n = cv->ctd_[tid].nvsize_; + for (int i = 0; i < n; ++i) { + out[k++] = p[i]; + } + } +} + +static void audit_copy_nv_in(Cvode* cv, const std::vector& in, N_Vector nv) { + int k = 0; + for (int tid = 0; tid < cv->nctd_; ++tid) { + double* p = cv->n_vector_data(nv, tid); + int n = cv->ctd_[tid].nvsize_; + for (int i = 0; i < n; ++i) { + p[i] = in[k++]; + } + } +} + +bool Daspk::audit_should_fire() const { + if (audit_level_ <= 0 || !audit_armed_) { + return false; + } + // First reinit with t >= select (user knows the time of interest). + return cv_->t_ >= audit_t_select_ - NetCvode::eps(audit_t_select_); +} + +// Summarize residual already in delta_ (no res call). +static void audit_residual_stats(Cvode* cv, + N_Vector delta, + N_Vector y_for_ewt, + void* mem, + double* max_abs, + double* wrms) { + *max_abs = 0.; + for (int tid = 0; tid < cv->nctd_; ++tid) { + double* d = cv->n_vector_data(delta, tid); + int n = cv->ctd_[tid].nvsize_; + for (int i = 0; i < n; ++i) { + double a = std::fabs(d[i]); + if (a > *max_abs) { + *max_abs = a; + } + } + } + *wrms = -1.; + if (mem) { + if (IDAEwtSet((IDAMem) mem, y_for_ewt)) { + *wrms = N_VWrmsNorm(delta, ((IDAMem) mem)->ida_ewt); + } + } +} + +void Daspk::audit_save_pre_from_delta() { + if (audit_level_ <= 0) { + return; + } + // Caller has just evaluated res into delta_ at cv_->t_ with continuous play + // synchronized to that t (advance end, or interpolate/retreat to event). + audit_copy_nv_out(cv_, cv_->y_, audit_y_pre_); + audit_copy_nv_out(cv_, yp_, audit_yp_pre_); + audit_copy_nv_out(cv_, delta_, audit_r_pre_); + audit_pre_t_ = cv_->t_; + audit_pre_neq_ = cv_->neq_; + audit_residual_stats(cv_, delta_, cv_->y_, mem_, &audit_pre_max_abs_, &audit_pre_wrms_); + audit_pre_valid_ = (audit_pre_neq_ > 0 && (int) audit_y_pre_.size() == audit_pre_neq_); + audit_pre_res_valid_ = audit_pre_valid_ && ((int) audit_r_pre_.size() == audit_pre_neq_); +} + +void Daspk::audit_eval_residual(N_Vector y, N_Vector yp, double* max_abs, double* wrms) { + res_gvardt(cv_->t_, y, yp, delta_, cv_); + audit_residual_stats(cv_, delta_, y, mem_, max_abs, wrms); +} + +void Daspk::audit_dump_panel(FILE* f, + const char* title, + N_Vector y, + N_Vector yp, + N_Vector delta, + double max_abs, + double wrms, + int max_rows) { + fprintf(f, "--- %s ---\n", title); + if (wrms >= 0.) { + fprintf(f, " WRMS(residual)=%.6g max|residual|=%.6g neq=%d\n", wrms, max_abs, cv_->neq_); + } else { + fprintf(f, " WRMS(residual)=n/a max|residual|=%.6g neq=%d\n", max_abs, cv_->neq_); + } + if (audit_level_ < 2 || cv_->neq_ <= 0) { + return; + } + // Rank equations by |residual|; print top max_rows (or all if neq small). + struct Row { + int i; + double y, yp, r, ar; + }; + std::vector rows; + rows.reserve(cv_->neq_); + int k = 0; + for (int tid = 0; tid < cv_->nctd_; ++tid) { + double* py = cv_->n_vector_data(y, tid); + double* pyp = cv_->n_vector_data(yp, tid); + double* pr = cv_->n_vector_data(delta, tid); + int n = cv_->ctd_[tid].nvsize_; + for (int i = 0; i < n; ++i, ++k) { + Row row; + row.i = k; + row.y = py[i]; + row.yp = pyp[i]; + row.r = pr[i]; + row.ar = std::fabs(pr[i]); + rows.push_back(row); + } + } + int nprint = cv_->neq_; + if (max_rows > 0 && nprint > max_rows) { + std::partial_sort(rows.begin(), + rows.begin() + max_rows, + rows.end(), + [](const Row& a, const Row& b) { return a.ar > b.ar; }); + nprint = max_rows; + fprintf(f, + " (top %d of %d eqs by |residual|; form: residual ~ c*y' - f(y))\n", + nprint, + cv_->neq_); + } else { + fprintf(f, " (all eqs; residual ~ c*y' - f(y))\n"); + } + fprintf(f, " %6s %16s %16s %16s\n", "eq", "y", "y'", "residual"); + for (int i = 0; i < nprint; ++i) { + const Row& r = rows[i]; + fprintf(f, " %6d %16.8g %16.8g %16.8g\n", r.i, r.y, r.yp, r.r); + } +} + +int Daspk::init() { +#if 0 +printf("Daspk_init t_=%20.12g t-t_=%g t0_-t_=%g mode=%d\n", +cv_->t_, t-cv_->t_, cv_->t0_-cv_->t_, init_mode_); +#endif + const bool do_audit = audit_should_fire(); + + // A1/A4: continuous Vector.play forcing t+ plus optional LM.dforce / FD b'. + // Play state (ubound_index_) is already post-event when reinit runs. + last_forcing_t_ = cv_->t_; + nrn_collect_forcing_tplus(cv_->t_, last_forcing_tplus_); + nrndae_append_dforce_to_forcing_list(cv_->t_, last_forcing_tplus_); + // Stdout when audit level >= 1 and this reinit is not already writing a + // three-panel dump (that dump includes the same block). + if (audit_level_ >= 1 && !do_audit && !last_forcing_tplus_.empty()) { + nrn_dump_forcing_tplus(stdout, last_forcing_t_, last_forcing_tplus_); + } + + // Capture panel B *before* the projector changes y (and possibly yp). + std::vector yB, ypB, rB; + double maxB = 0., wrmsB = -1.; + bool have_B = false; + // Panel A: continuous pre-reinit snapshot (prefer post-interpolate residual). + std::vector yA, ypA, rA; + double tA = 0., maxA_saved = 0., wrmsA_saved = -1.; + bool have_A = false; + bool have_A_res = false; + if (do_audit) { + if (audit_pre_valid_ && audit_pre_neq_ == cv_->neq_) { + yA = audit_y_pre_; + ypA = audit_yp_pre_; + tA = audit_pre_t_; + have_A = true; + if (audit_pre_res_valid_ && (int) audit_r_pre_.size() == audit_pre_neq_) { + rA = audit_r_pre_; + maxA_saved = audit_pre_max_abs_; + wrmsA_saved = audit_pre_wrms_; + have_A_res = true; + } + } + cv_->play_continuous(cv_->t_); + cv_->daspk_gather_y(cv_->y_); + if (have_A) { + audit_copy_nv_in(cv_, ypA, yp_); + } else { + N_VConst(0., yp_); + } + audit_eval_residual(cv_->y_, yp_, &maxB, &wrmsB); + audit_copy_nv_out(cv_, cv_->y_, yB); + audit_copy_nv_out(cv_, yp_, ypB); + audit_copy_nv_out(cv_, delta_, rB); + have_B = true; + // Scatter post-event y back so projector sees the same model state. + cv_->daspk_scatter_y(cv_->y_); + } + + int path_mode = init_mode_; + bool fell_back = false; + int err = 0; + last_ic_forcing_flags_ = 0; + ++ic_init_count_; + if (init_mode_ == 0) { + err = init_heuristic(); + path_mode = 0; + } else if (init_mode_ == 3) { + // Hold continuous content (LM / extracellular), then y' from C*y'=f(y) + // plus free y' from forcing t+ (play / dforce / FD). + // On residual failure, fall back to nano-step heuristic (unless audit + // is armed — keep pure mode-3 state for panel C). + err = init_battery(); + path_mode = 3; + if (err != 0) { + if (do_audit) { + Printf( + "mode 3 IC residual failed (err=%d, forcing_flags=0x%x); audit armed — not " + "falling back to heuristic\n", + err, + last_ic_forcing_flags_); + } else { + Printf( + "mode 3 IC residual failed (err=%d, forcing_flags=0x%x); falling back to " + "nano-step heuristic IC\n", + err, + last_ic_forcing_flags_); + ++calcic_fallback_count_; + ++ic_mode3_fallback_count_; + fell_back = true; + err = init_heuristic(); + path_mode = 0; + } + } else { + ++ic_mode3_ok_count_; + } + } else { + err = init_ida_y_init(); + path_mode = init_mode_; + if (err != 0 && init_mode_ == 1) { + Printf( + "IDACalcIC residual/project failed (err=%d); falling back to nano-step " + "heuristic IC\n", + err); + ++calcic_fallback_count_; + fell_back = true; + err = init_heuristic(); + path_mode = 0; + } + } + last_ic_path_mode_ = path_mode; + + if (do_audit && have_B) { + // Preserve post-IC (y, yp) before temporarily loading A/B for printing. + std::vector yC, ypC; + audit_copy_nv_out(cv_, cv_->y_, yC); + audit_copy_nv_out(cv_, yp_, ypC); + + FILE* f = audit_open_out(); + const bool close_f = (f != stdout); + ++audit_serial_; + fprintf(f, + "\n=== IDA IC three-panel audit serial=%d t=%.15g requested_mode=%d " + "path_mode=%d fallback=%d err=%d ===\n", + audit_serial_, + cv_->t_, + init_mode_, + path_mode, + fell_back ? 1 : 0, + err); + if (have_A) { + fprintf(f, + " note: reinit after discontinuity (NET_RECEIVE / Vector.play / at_time / " + "...)\n"); + } else { + fprintf(f, " note: no continuous pre-snapshot (typical of finitialize)\n"); + } + if (init_mode_ == 3 && path_mode == 3 && do_audit) { + fprintf(f, + " note: mode 3 panel C: C*y'=f(y) seed + free y' from forcing t+ " + "(null(C) / Z'G y'=Z'b') flags=0x%x" + "%s\n", + last_ic_forcing_flags_, + err != 0 ? "; residual failed; fallback suppressed (audit armed)" : ""); + if (last_ic_forcing_flags_ & NRN_IC_FORCING_APPLIED) { + fprintf(f, + " note: free y' sources:%s%s%s\n", + (last_ic_forcing_flags_ & NRN_IC_FORCING_PLAY) ? " play" : "", + (last_ic_forcing_flags_ & NRN_IC_FORCING_DFORCE) ? " dforce" : "", + (last_ic_forcing_flags_ & NRN_IC_FORCING_FD) ? " fd" : ""); + } else if (err == 0) { + fprintf(f, + " note: no free-y' forcing applied (C*y'=f seed only; ok if u'=0 or no " + "singular C)\n"); + } + } + // A1: always show forcing t+ in the three-panel dump when present + nrn_dump_forcing_tplus(f, last_forcing_t_, last_forcing_tplus_); + + double maxA = 0., wrmsA = -1.; + if (have_A) { + // Use residual captured with continuous play at tA (typically the + // interpolate/retreat residual). Do not re-call res after the jump. + audit_copy_nv_in(cv_, yA, cv_->y_); + audit_copy_nv_in(cv_, ypA, yp_); + if (have_A_res) { + audit_copy_nv_in(cv_, rA, delta_); + maxA = maxA_saved; + wrmsA = wrmsA_saved; + } else { + // Fallback only if residual was not stored (should be rare). + double t_save = cv_->t_; + cv_->t_ = tA; + audit_eval_residual(cv_->y_, yp_, &maxA, &wrmsA); + cv_->t_ = t_save; + } + char title[192]; + snprintf(title, + sizeof title, + "A pre (continuous at t=%.15g; residual from retreat/step, not re-eval after " + "jump)", + tA); + audit_dump_panel(f, title, cv_->y_, yp_, delta_, maxA, wrmsA, 40); + } else { + fprintf(f, + "--- A pre --- (unavailable: no prior continuous snapshot; e.g. " + "finitialize)\n"); + } + + // Panel B: vectors captured before the projector (no re-eval). + audit_copy_nv_in(cv_, yB, cv_->y_); + audit_copy_nv_in(cv_, ypB, yp_); + audit_copy_nv_in(cv_, rB, delta_); + audit_dump_panel(f, + "B post-event pre-IC (y after discontinuity; y' from pre or 0)", + cv_->y_, + yp_, + delta_, + maxB, + wrmsB, + 40); + + // Panel C: post-IC state (saved above) + audit_copy_nv_in(cv_, yC, cv_->y_); + audit_copy_nv_in(cv_, ypC, yp_); + double maxC = 0., wrmsC = -1.; + audit_eval_residual(cv_->y_, yp_, &maxC, &wrmsC); + audit_dump_panel(f, "C post-IC", cv_->y_, yp_, delta_, maxC, wrmsC, 40); + + double max_dy = 0.; + int i_max = -1; + for (size_t i = 0; i < yB.size() && i < yC.size(); ++i) { + double d = std::fabs(yC[i] - yB[i]); + if (d > max_dy) { + max_dy = d; + i_max = (int) i; + } + } + fprintf(f, "--- summary ---\n"); + fprintf(f, " max|res| A/B/C = %.6g / %.6g / %.6g\n", maxA, maxB, maxC); + if (wrmsA >= 0. || wrmsB >= 0. || wrmsC >= 0.) { + fprintf(f, " WRMS A/B/C = %.6g / %.6g / %.6g\n", wrmsA, wrmsB, wrmsC); + } + fprintf(f, " max|y_C - y_B| = %.6g (eq %d)\n", max_dy, i_max); + if (have_A && (int) yC.size() == (int) yA.size()) { + double max_from_pre = 0.; + int j_max = -1; + for (size_t i = 0; i < yA.size(); ++i) { + double d = std::fabs(yC[i] - yA[i]); + if (d > max_from_pre) { + max_from_pre = d; + j_max = (int) i; + } + } + fprintf(f, " max|y_C - y_A| = %.6g (eq %d)\n", max_from_pre, j_max); + } + fprintf(f, "=== end IDA IC audit ===\n\n"); + if (close_f) { + fclose(f); + } + audit_armed_ = 0; + // Restore post-IC into model and IDA vectors. + audit_copy_nv_in(cv_, yC, cv_->y_); + audit_copy_nv_in(cv_, ypC, yp_); + cv_->daspk_scatter_y(cv_->y_); + } + + // After a successful IC: continuous residual at this IC time for a later A + // if the next reinit has no intervening interpolate (unusual). + if (err == 0 && audit_level_ > 0) { + // check_init_residual left residual in parasite_ / via res; re-eval once + // so delta_ matches current y,yp at cv_->t_. + double max_abs = 0., wrms = -1.; + audit_eval_residual(cv_->y_, yp_, &max_abs, &wrms); + audit_save_pre_from_delta(); + } + return err; } int Daspk::advance_tn(double tstop) { @@ -342,6 +1157,11 @@ int Daspk::advance_tn(double tstop) { #endif cv_->t0_ = tn; cv_->tn_ = cv_->t_; + // Continuous residual at step end. If an event retreats via interpolate(), + // that call overwrites this snapshot with the true pre-event point. + if (ier >= 0) { + audit_save_pre_from_delta(); + } // printf("Daspk::advance_tn complete.\n"); return ier; } @@ -356,7 +1176,12 @@ int Daspk::interpolate(double tt) { } cv_->t_ = tt; // interpolation does not call res. So we have to. + // Continuous play is synchronized to tt inside res — this is the natural + // pre-event residual for panel A (play / NetCon / at_time retreat). res_gvardt(cv_->t_, cv_->y_, yp_, delta_, cv_); + if (ier >= 0) { + audit_save_pre_from_delta(); + } // if(MyMath::eq(t, cv_->t_, NetCvode::eps(cv_->t_))) { // printf("t=%.15g t_=%.15g\n", t, cv_->t_); //} @@ -379,9 +1204,7 @@ void Daspk::statistics() { printf("nonlinear conv. failures = %d\n", iwork_[15-1]); printf("linear conv. failures = %d\n", iwork_[16-1]); #endif - if (first_try_init_failures_) { - Printf(" %d First try Initialization failures\n", first_try_init_failures_); - } + print_ic_stats(); } static void* daspk_scatter_thread(NrnThread* nt) { @@ -615,8 +1438,15 @@ printf(" %g", b[i]); printf("\n"); #endif - _nt->cj = cj; - _nt->_dt = 1. / cj; + // IDACalcIC(IDA_Y_INIT) sets cj = 0 (solve all y given yp; J = dG/dy only). + // Avoid 1/cj and skip the mechanism-ODE 1/cj scaling in that case. + if (cj == 0.0) { + _nt->cj = 0.0; + _nt->_dt = 1.0; + } else { + _nt->cj = cj; + _nt->_dt = 1. / cj; + } _nt->_vcv = this; daspk_scatter_y(y, _nt->id); // I'm not sure this is necessary. @@ -639,6 +1469,11 @@ for (i=0; i < z.neq_v_; ++i) { } #endif daspk_nrn_solve(_nt); // not the cvode one + if (nrn_sparse13_factor_error()) { + // Force next psol to rebuild via nrn_lhs (spClear resets sparse Error). + solve_state_ = INVALID; + return 1; // recoverable; IDACalcIC may fail soft when J singular + } #if 0 //printf("after nrn_solve matrix\n"); //spPrint(sp13mat_, 1,1,1); @@ -653,8 +1488,12 @@ for (i=0; i < neq_v_; ++i) { // the ode's of the form m' = (minf - m)/mtau in model descriptions compute // b = b/(1 + dt*mtau) since cvode required J = 1 - gam*df/dy // so we need to scale those results by 1/cj. - for (i = z.neq_v_; i < z.nvsize_; ++i) { - b[i] *= _nt->_dt; + // When cj==0 (IDA_Y_INIT), leave solvemem results unscaled; mechanism-ODE + // quality in that limit may still need refinement (see ida_y_init notes). + if (cj != 0.0) { + for (i = z.neq_v_; i < z.nvsize_; ++i) { + b[i] *= _nt->_dt; + } } #if 0 for (i=0; i < z.nvsize_; ++i) { diff --git a/src/nrncvode/nrndaspk.h b/src/nrncvode/nrndaspk.h index 6826a25b69..f2604582d1 100644 --- a/src/nrncvode/nrndaspk.h +++ b/src/nrncvode/nrndaspk.h @@ -4,6 +4,10 @@ #include "nvector_nrnthread.h" #include "nvector_nrnthread_ld.h" #include "nvector_nrnserial_ld.h" +#include "vecplay_tplus.h" +#include +#include +#include class Cvode; @@ -18,9 +22,54 @@ class Daspk { N_Vector ewtvec(); N_Vector acorvec(); + // Three-panel IC audit (stdout or file). See dae_init_audit docs. + static void audit_set_level(int level); + static int audit_level(); + static void audit_arm_at(double t); + static double audit_t_select(); + static void audit_set_file(const char* path); // empty/null → stdout + static const char* audit_file(); + + // Last IC-time continuous Vector.play forcing t+ collection (A1; for A2 y'). + static const std::vector& last_forcing_tplus(); + static double last_forcing_t(); + + // A5: cumulative IC path statistics (reset on CVode construction / clear). + static void reset_ic_stats(); + static void print_ic_stats(); // used by statistics() + // last IC: path_mode and forcing source bitmask (NRN_IC_FORCING_*) + static int last_ic_path_mode(); + static int last_ic_forcing_flags(); + static int ic_mode3_ok_count(); + static int ic_mode3_fallback_count(); + private: void ida_init(); void info(); + // Heuristic nano-step IC (legacy). Returns 0 on success / accepted warn path. + int init_heuristic(); + // IDACalcIC(IDA_Y_INIT) with yp=0 + ODE f(y). Returns 0 on residual success. + int init_ida_y_init(); + // Mode 3: battery content hold (LM / extracellular) + y' from C*y'=f(y). + int init_battery(); + // Shared residual WRMS check and parasite / style handling after y,yp ready. + int check_init_residual(); + + // IC audit helpers + bool audit_should_fire() const; + // Save continuous (y, yp) and residual already in delta_ (caller just evaluated res). + // Prefer calling after interpolate/retreat so panel A is at t_event, not step end tn. + void audit_save_pre_from_delta(); + void audit_eval_residual(N_Vector y, N_Vector yp, double* max_abs, double* wrms); + void audit_dump_panel(FILE* f, + const char* title, + N_Vector y, + N_Vector yp, + N_Vector delta, + double max_abs, + double wrms, + int max_rows); + FILE* audit_open_out(); public: void* mem_; @@ -35,4 +84,44 @@ class Daspk { static double dteps_; static int init_try_again_; static int first_try_init_failures_; + // 0 = heuristic only (default); 1 = IDA_Y_INIT then heuristic fallback; + // 2 = IDA_Y_INIT only; 3 = battery content hold + y' from C*y'=f(y), + // with heuristic fallback on residual failure. + static int init_mode_; + // Count of mode 1/3 residual failures that fell back to the heuristic. + static int calcic_fallback_count_; + // A5: mode-3 specific + forcing-source tallies + static int ic_init_count_; + static int ic_mode3_ok_count_; + static int ic_mode3_fallback_count_; + static int ic_forcing_play_inits_; // inits that used play b' + static int ic_forcing_dforce_inits_; // inits that used dforce/bdot + static int ic_forcing_fd_inits_; // inits that used FD of force callable + static int last_ic_path_mode_; + static int last_ic_forcing_flags_; + + // Audit control (process-wide; one IDA path typically). + static int audit_level_; // 0 off, 1 summary, 2 three-panel (top residual rows) + static double audit_t_select_; // first reinit with t >= this (when armed) + static int audit_armed_; // 1 = waiting for t match + static int audit_serial_; // reinit count (all reinits) + static std::string audit_path_; // empty → stdout + + // Forcing t+ from continuous Vector.play at last Daspk::init (A1). + static std::vector last_forcing_tplus_; + static double last_forcing_t_; + + private: + // Continuous pre-reinit snapshot for panel A (play/NetCon/at_time retreat). + // Residual is captured when res was evaluated with continuous play at that t + // (do not re-eval after the discontinuity has been applied). + std::vector audit_y_pre_; + std::vector audit_yp_pre_; + std::vector audit_r_pre_; + double audit_pre_t_; + double audit_pre_max_abs_; + double audit_pre_wrms_; + int audit_pre_neq_; + bool audit_pre_valid_; + bool audit_pre_res_valid_; // residual vector matches y/yp/t above }; diff --git a/src/nrncvode/vrecitem.h b/src/nrncvode/vrecitem.h index adaafb674e..1644acab1b 100644 --- a/src/nrncvode/vrecitem.h +++ b/src/nrncvode/vrecitem.h @@ -292,6 +292,8 @@ class VecPlayContinuous: public PlayRecord { void continuous(double tt); double interpolate(double tt); + /** Forcing t⁺ info: u(tt) and classical u'(tt) matching continuous play. */ + void forcing_tplus(double tt, double* value, double* deriv) const; double interp(double th, double x0, double x1) { return x0 + (x1 - x0) * th; } diff --git a/src/nrniv/linmod.cpp b/src/nrniv/linmod.cpp index 03e7e0bd0f..c7fe1fe2f3 100644 --- a/src/nrniv/linmod.cpp +++ b/src/nrniv/linmod.cpp @@ -23,8 +23,12 @@ // daspk equation indices. #include +#include +#include #include "linmod.h" #include "nrnpy.h" +#include "ocmatrix.h" +#include "nrn_ansi.h" // for extern t if needed LinearModelAddition::LinearModelAddition(Matrix* cmat, Matrix* gmat, @@ -37,16 +41,109 @@ LinearModelAddition::LinearModelAddition(Matrix* cmat, Object* f_callable) : NrnDAE(cmat, yvec, y0, nnode, nodes, elayer) , b_(*bvec) - , f_callable_(f_callable) { + , f_callable_(f_callable) + , dforce_callable_(nullptr) + , bdot_(nullptr) { // printf("LinearModelAddition %p\n", this); g_ = new MatrixMap(gmat); } LinearModelAddition::~LinearModelAddition() { // printf("~LinearModelAddition %p\n", this); + if (dforce_callable_) { + hoc_obj_unref(dforce_callable_); + dforce_callable_ = nullptr; + } delete g_; } +void LinearModelAddition::set_dforce(Object* dforce_callable, Vect* bdot) { + if (dforce_callable && !bdot) { + hoc_execerror("LinearMechanism.dforce: callable requires a bdot Vector", 0); + } + if (dforce_callable_) { + hoc_obj_unref(dforce_callable_); + dforce_callable_ = nullptr; + } + if (dforce_callable) { + dforce_callable_ = dforce_callable; + hoc_obj_ref(dforce_callable_); + } + bdot_ = bdot; + // size_ may still be 0 before first matrix alloc; compare to b_ + if (bdot_ && bdot_->size() != b_.size()) { + hoc_execerror("LinearMechanism.dforce: bdot size must match b", 0); + } +} + +bool LinearModelAddition::fill_bdot_for_ic(double tt, double* out, double fd_h) { + if (!out || size_ <= 0) { + return false; + } + for (int i = 0; i < size_; ++i) { + out[i] = 0.; + } + // Analytic / user vector path + if (dforce_callable_) { + extern double t; + const double t_sav = t; + t = tt; + if (nrn_threads) { + nrn_threads->_t = tt; + } + if (!neuron::python::methods.hoccommand_exec || + !neuron::python::methods.hoccommand_exec(dforce_callable_)) { + t = t_sav; + if (nrn_threads) { + nrn_threads->_t = t_sav; + } + hoc_execerror("LinearMechanism.dforce callable failed", 0); + } + t = t_sav; + if (nrn_threads) { + nrn_threads->_t = t_sav; + } + } + if (bdot_) { + for (int i = 0; i < size_; ++i) { + out[i] = bdot_->elem(i); + } + return true; + } + // FD fallback: only when f_callable updates b and no bdot vector was given + if (!f_callable_ || fd_h <= 0.) { + return false; + } + extern double t; + const double t_sav = t; + std::vector b0(size_); + // Ensure b is current at tt + t = tt; + if (nrn_threads) { + nrn_threads->_t = tt; + } + f_(y_, yptmp_, size_); + for (int i = 0; i < size_; ++i) { + b0[i] = b_[i]; + } + t = tt + fd_h; + if (nrn_threads) { + nrn_threads->_t = t; + } + f_(y_, yptmp_, size_); + for (int i = 0; i < size_; ++i) { + out[i] = (b_[i] - b0[i]) / fd_h; + b_[i] = b0[i]; // restore + } + t = t_sav; + if (nrn_threads) { + nrn_threads->_t = t_sav; + } + // re-sync b via f_ at original t + f_(y_, yptmp_, size_); + return true; +} + void LinearModelAddition::alloc_(int size, int start, int nnode, Node** nodes, int* elayer) { // printf("LinearModelAddition::alloc_ %p\n", this); assert(b_.size() == size); @@ -82,3 +179,269 @@ double LinearModelAddition::jacobian_multiplier_() { MatrixMap* LinearModelAddition::jacobian_(Vect& y) { return g_; } + +int LinearModelAddition::battery_ic_project() { + // Consistent IC: hold continuous content, solve algebraics. + // + // 1) Floating capacitors (both Cij and Cji nonzero): replace with voltage + // sources yi - yj = hold, free battery current (classic C→battery). + // 2) Diagonal C_ii only (inductor current state): hold yi, drop eqn i + // (L→current source of value yi). + // 3) One-sided C_row,col (op-amp lag tau*vk' in row o): hold y_col (vk), + // drop dynamic equation row (not a floating two-terminal cap). + // + // Solves augmented [G stamps + constraints] for y. + if (assumed_identity_) { + return 0; + } + const int n = size_; + if (n <= 0) { + return -1; + } + + constexpr double ctol = 1e-18; + + struct FloatingCap { + int i; + int j; + double hold; + }; + struct StateHold { + int var; // y index held continuous + int drop_eqn; // dynamic residual row replaced by hold + double hold; + }; + + std::vector caps; + std::vector row_used(n, 0); + std::vector col_in_floating(n, 0); + + // --- floating capacitors: mutual off-diagonal pair --- + for (int i = 0; i < n; ++i) { + for (int j = i + 1; j < n; ++j) { + const double cij = (*c_)(i, j); + const double cji = (*c_)(j, i); + if (std::fabs(cij) > ctol && std::fabs(cji) > ctol) { + FloatingCap cap; + cap.i = i; + cap.j = j; + cap.hold = y_[i] - y_[j]; + caps.push_back(cap); + row_used[i] = 1; + row_used[j] = 1; + col_in_floating[i] = 1; + col_in_floating[j] = 1; + } + } + } + + // --- remaining C rows: inductor (diag) or op-amp lag (one-sided) --- + std::vector holds; + for (int r = 0; r < n; ++r) { + if (row_used[r]) { + continue; + } + int nz_cols[8]; + int nnz = 0; + for (int j = 0; j < n && nnz < 8; ++j) { + if (std::fabs((*c_)(r, j)) > ctol) { + nz_cols[nnz++] = j; + } + } + if (nnz == 0) { + continue; + } + StateHold h; + h.drop_eqn = r; + if (nnz == 1 && nz_cols[0] == r) { + // Diagonal mass: inductor current (or similar) — hold y_r + h.var = r; + h.hold = y_[r]; + } else { + // Prefer off-diagonal column (variable being differentiated), e.g. + // OpAmp: C[o][k]=tau ⇒ hold output voltage y_k, drop eqn o. + int jhold = nz_cols[0]; + for (int k = 0; k < nnz; ++k) { + if (nz_cols[k] != r) { + jhold = nz_cols[k]; + break; + } + } + h.var = jhold; + h.hold = y_[jhold]; + } + holds.push_back(h); + row_used[r] = 1; + } + + const int nc = static_cast(caps.size()); + const int nh = static_cast(holds.size()); + const int m = n + nc; // unknowns: y[0..n) and is[0..nc) + + // Equations: (n - nh) KCL/alg rows + nc floating constraints + nh holds = n+nc + Matrix* A = Matrix::instance(m, m); + Vect rhs(m); + Vect sol(m); + for (int i = 0; i < m; ++i) { + rhs[i] = 0.; + for (int j = 0; j < m; ++j) { + (*A)(i, j) = 0.; + } + } + + std::vector drop(n, 0); + for (const auto& h: holds) { + drop[h.drop_eqn] = 1; + } + + int eq = 0; + // Keep algebraic/KCL rows (not dropped dynamic equations) + for (int r = 0; r < n; ++r) { + if (drop[r]) { + continue; + } + for (int j = 0; j < n; ++j) { + (*A)(eq, j) = (*g_)(r, j); + } + // battery currents on floating caps + for (int k = 0; k < nc; ++k) { + const int i = caps[k].i; + const int j = caps[k].j; + if (r == i) { + (*A)(eq, n + k) = 1.; + } else if (r == j) { + (*A)(eq, n + k) = -1.; + } + } + rhs[eq] = b_[r]; + ++eq; + } + // Floating capacitor constraints + for (int k = 0; k < nc; ++k) { + (*A)(eq, caps[k].i) = 1.; + (*A)(eq, caps[k].j) = -1.; + rhs[eq] = caps[k].hold; + ++eq; + } + // State holds (inductor current, op-amp output voltage, …) + for (int k = 0; k < nh; ++k) { + (*A)(eq, holds[k].var) = 1.; + rhs[eq] = holds[k].hold; + ++eq; + } + if (eq != m) { + // Structure mismatch — refuse rather than solve a bad system + delete A; + return -1; + } + + A->solv(&rhs, &sol, true); + for (int i = 0; i < n; ++i) { + y_[i] = sol[i]; + } + delete A; + return 0; +} + +bool LinearModelAddition::b_element_is(int i, double* play_target) const { + if (!play_target || i < 0 || i >= size_) { + return false; + } + // Vect stores contiguous doubles; play into b._ref_x[i] targets &b[i] + return play_target == &b_[i]; +} + +void LinearModelAddition::complete_yp_from_bdot(const double* bdot, double* yp_global) { + // C y' = b - G y is underdetermined when C is singular. Existing seed sets a + // particular solution; free components live in null(C). Differentiated + // left-null constraints Z^T G y' = Z^T b' fix those when G has a component + // on the free mode (e.g. R to ground on a floating C–R series pair). + if (!bdot || !yp_global || size_ <= 0) { + return; + } + if (assumed_identity_) { + return; // C = I, no free y' + } + if (!c_ || !g_) { + return; + } + const int n = size_; + constexpr double ctol = 1e-12; + constexpr double gtol = 1e-18; + + std::vector yp(n); + for (int i = 0; i < n; ++i) { + yp[i] = yp_global[bmap_[i] - 1]; + } + + // Floating capacitor pairs: right null n = e_i + e_j, left null Z = e_i + e_j + // when C is a pure mutual difference stamp on (i,j). + std::vector used(n, 0); + for (int i = 0; i < n; ++i) { + for (int j = i + 1; j < n; ++j) { + const double cij = (*c_)(i, j); + const double cji = (*c_)(j, i); + if (std::fabs(cij) <= ctol || std::fabs(cji) <= ctol) { + continue; + } + // Confirm Z = e_i+e_j is ~left null of C + bool z_null = true; + for (int k = 0; k < n; ++k) { + const double ztc = (*c_)(i, k) + (*c_)(j, k); + if (std::fabs(ztc) > + ctol * (1. + std::fabs((*c_)(i, k)) + std::fabs((*c_)(j, k)))) { + z_null = false; + break; + } + } + if (!z_null) { + continue; + } + // Z^T G yp and Z^T G n with n = e_i + e_j + double ztg_yp = 0.; + for (int k = 0; k < n; ++k) { + ztg_yp += ((*g_)(i, k) + (*g_)(j, k)) * yp[k]; + } + const double ztg_n = (*g_)(i, i) + (*g_)(i, j) + (*g_)(j, i) + (*g_)(j, j); + const double ztbdot = bdot[i] + bdot[j]; + if (std::fabs(ztg_n) < gtol) { + continue; // free mode invisible to G — leave yp + } + const double alpha = (ztbdot - ztg_yp) / ztg_n; + yp[i] += alpha; + yp[j] += alpha; + used[i] = 1; + used[j] = 1; + } + } + + // Pure algebraic rows (C row ~ 0): enforce (G y')_r = bdot[r] by adjusting + // free algebraic variables where possible (simple diagonal G_rr). + for (int r = 0; r < n; ++r) { + bool crow0 = true; + for (int k = 0; k < n; ++k) { + if (std::fabs((*c_)(r, k)) > ctol) { + crow0 = false; + break; + } + } + if (!crow0) { + continue; + } + const double grr = (*g_)(r, r); + if (std::fabs(grr) < gtol) { + continue; + } + // residual of differentiated row: bdot[r] - sum_k G_rk yp[k] + double gyp = 0.; + for (int k = 0; k < n; ++k) { + gyp += (*g_)(r, k) * yp[k]; + } + // Adjust yp[r] only (local diagonal completion) + yp[r] += (bdot[r] - gyp) / grr; + } + + for (int i = 0; i < n; ++i) { + yp_global[bmap_[i] - 1] = yp[i]; + } +} diff --git a/src/nrniv/linmod.h b/src/nrniv/linmod.h index 479ceb6d51..e9dcb8919e 100644 --- a/src/nrniv/linmod.h +++ b/src/nrniv/linmod.h @@ -20,6 +20,51 @@ class LinearModelAddition: public NrnDAE { Object* f_callable = NULL); virtual ~LinearModelAddition(); + /** + * Spike: consistent IC by replacing capacitors with voltage sources that + * hold branch Δv = y_i - y_j (battery replacement). Solves the augmented + * algebraic system [G, B; B^T, 0] [y; i_s] = [b; hold] and writes y back. + * @return 0 on success, nonzero on failure / nothing to do. + */ + int battery_ic_project(); + + /** + * After y is fixed and a particular C*yp ≈ b-Gy is seeded, adjust yp in + * null(C) so differentiated algebraic constraints hold: Z^T G yp = Z^T bdot + * for left-null vectors Z of C (e.g. common mode of a floating capacitor). + * bdot is db/dt from forcing t+ (same layout as b_). yp_global is the full + * IDA y' vector (bmap_ applied inside). + */ + void complete_yp_from_bdot(const double* bdot, double* yp_global); + + /** True if play_target points at b_[i]; used to map Vector.play → b'. */ + bool b_element_is(int i, double* play_target) const; + + int size() const { + return size_; + } + + /** + * A4: optional db/dt for IC free y'. + * @param dforce_callable if non-null, executed at IC (e.g. Python) before reading bdot + * @param bdot vector same length as b; holds db/dt (not owned) + * Either or both may be set. Clear with set_dforce(nullptr, nullptr). + */ + void set_dforce(Object* dforce_callable, Vect* bdot); + + /** Call dforce (if any) and fill out[0..size) with bdot; FD fallback if needed. */ + bool fill_bdot_for_ic(double tt, double* out, double fd_h = 1e-8); + + Object* f_callable() const { + return f_callable_; + } + Object* dforce_callable() const { + return dforce_callable_; + } + Vect* bdot_vec() const { + return bdot_; + } + private: void f_(Vect& y, Vect& yprime, int size); MatrixMap* jacobian_(Vect& y); @@ -29,4 +74,6 @@ class LinearModelAddition: public NrnDAE { MatrixMap* g_; Vect& b_; Object* f_callable_; + Object* dforce_callable_; // A4: optional analytic db/dt provider + Vect* bdot_; // A4: user vector for db/dt (not owned) }; diff --git a/src/nrniv/linmod1.cpp b/src/nrniv/linmod1.cpp index 053ff068c3..95f0eae63c 100644 --- a/src/nrniv/linmod1.cpp +++ b/src/nrniv/linmod1.cpp @@ -47,7 +47,42 @@ static double valid(void* v) { return double(((LinearMechanism*) v)->valid()); } -static Member_func members[] = {{"valid", valid}, {nullptr, nullptr}}; +// dforce(bdot_vector) +// dforce(python_callable, bdot_vector) +// Optional analytic db/dt for IDA mode-3 free y' (Plan A4). Callable is +// invoked at IC with t set to the IC time; it should fill bdot_vector. +static double dforce(void* v) { + LinearMechanism* m = (LinearMechanism*) v; + if (!m->model_) { + return 0.; + } + Object* callable = nullptr; + Vect* bdot = nullptr; + if (!ifarg(1)) { + m->model_->set_dforce(nullptr, nullptr); + return 0.; + } + if (hoc_is_object_arg(1)) { + Object* o = *hoc_objgetarg(1); + if (o && o->ctemplate && strcmp(o->ctemplate->sym->name, "PythonObject") == 0) { + callable = o; + if (!ifarg(2) || !is_vector_arg(2)) { + hoc_execerror("LinearMechanism.dforce: need dforce(callable, bdot_Vector)", 0); + } + bdot = vector_arg(2); + } else if (is_vector_arg(1)) { + bdot = vector_arg(1); + } else { + hoc_execerror("LinearMechanism.dforce: arg must be Vector or (callable, Vector)", 0); + } + } else { + hoc_execerror("LinearMechanism.dforce: arg must be Vector or (callable, Vector)", 0); + } + m->model_->set_dforce(callable, bdot); + return 1.; +} + +static Member_func members[] = {{"valid", valid}, {"dforce", dforce}, {nullptr, nullptr}}; static void* cons(Object*) { LinearMechanism* m = new LinearMechanism(); diff --git a/src/nrniv/nrndae.cpp b/src/nrniv/nrndae.cpp index cf1d27619d..034ec6db29 100644 --- a/src/nrniv/nrndae.cpp +++ b/src/nrniv/nrndae.cpp @@ -1,10 +1,19 @@ #include <../../nrnconf.h> +#include #include +#include +#include #include "nrndae.h" #include "nrndae_c.h" #include "nrnoc2iv.h" #include "treeset.h" #include "utils/enumerate.h" +#include "linmod.h" +#include "netcvode.h" +#include "vrecitem.h" +#include "vecplay_tplus.h" + +extern NetCvode* net_cvode_instance; extern int secondorder; @@ -103,6 +112,21 @@ void nrndae_dkres(double* y, double* yprime, double* delta) { } } +int nrndae_battery_ic_project() { + int err = 0; + for (NrnDAE* item: nrndae_list) { + auto* lm = dynamic_cast(item); + if (!lm) { + continue; + } + const int e = lm->battery_ic_project(); + if (e != 0) { + err = e; + } + } + return err; +} + inline void NrnDAE::alloc_(int size, int start, int nnode, Node** nodes, int* elayer) {} void NrnDAE::alloc(int start_index) { @@ -163,14 +187,13 @@ NrnDAE::NrnDAE(Matrix* cmat, cmat = assumed_identity_; } c_ = new MatrixMap(cmat); - Vect& elay = *elayer; nnode_ = nnode; nodes_ = nodes; if (nnode_ > 0) { elayer_ = new int[nnode_]; if (elayer) { for (int i = 0; i < nnode_; ++i) { - elayer_[i] = int(elay[i]); + elayer_[i] = int((*elayer)[i]); } } else { for (int i = 0; i < nnode_; ++i) { @@ -236,20 +259,17 @@ void NrnDAE::update() { void NrnDAE::init() { // printf("NrnDAE::init %lx\n", (long)this); // printf("init size_=%d %d %d %d\n", size_, y_->size(), y0_->size(), b_->size()); - Vect& y0 = *y0_; v2y(); if (f_init_) { f_init_(data_); + } else if (y0_) { + for (int i = nnode_; i < size_; ++i) { + y_[i] = (*y0_)[i]; + } } else { - if (y0_) { - for (int i = nnode_; i < size_; ++i) { - y_[i] = y0[i]; - } - } else { - for (int i = nnode_; i < size_; ++i) { - y_[i] = 0.; - } + for (int i = nnode_; i < size_; ++i) { + y_[i] = 0.; } } // for (i=0; i < nnode_; ++i) printf(" i=%d y[i]=%g\n", i, y[i]); @@ -301,6 +321,171 @@ void NrnDAE::dkres(double* y, double* yprime, double* delta) { } } +void NrnDAE::seed_yp_from_f(double* f, double* yp) { + // Set yp on mapped equations so C*yp ≈ f for simple mass structure. + // Residual uses delta -= C*yp with the same map (see dkres). + constexpr double ctol = 1e-18; + if (assumed_identity_) { + for (int i = 0; i < size_; ++i) { + yp[bmap_[i] - 1] = f[bmap_[i] - 1]; + } + return; + } + if (!c_ || size_ <= 0) { + return; + } + for (int r = 0; r < size_; ++r) { + int jnz[8]; + double cval[8]; + int nnz = 0; + for (int j = 0; j < size_ && nnz < 8; ++j) { + const double cij = (*c_)(r, j); + if (std::fabs(cij) > ctol) { + jnz[nnz] = j; + cval[nnz] = cij; + ++nnz; + } + } + if (nnz == 0) { + continue; // algebraic row: yp free / leave existing + } + if (nnz == 1) { + // Diagonal mass or one-sided lag: c_rj * yp_j = f_r + yp[bmap_[jnz[0]] - 1] = f[bmap_[r] - 1] / cval[0]; + continue; + } + if (nnz == 2 && std::fabs(cval[0] + cval[1]) < ctol * (1. + std::fabs(cval[0]))) { + // Pure difference stamp C*(yp_a - yp_b) = f_r (floating capacitor row). + // Gauge: leave the more negative column's yp unchanged (often 0), set the other. + const int ja = jnz[0]; + const int jb = jnz[1]; + // c_ra * yp_a + c_rb * yp_b = f_r with c_rb ≈ -c_ra + // => yp_a - yp_b = f_r / c_ra + const double scale = cval[0]; + const double dyp = f[bmap_[r] - 1] / scale; + // Keep yp[jb] as is (0 unless set by another row), set yp[ja] + yp[bmap_[ja] - 1] = yp[bmap_[jb] - 1] + dyp; + } + // denser rows: leave yp; residual check / heuristic fallback may apply + } +} + +void nrndae_seed_yp_from_f(double* f, double* yp) { + for (NrnDAE* item: nrndae_list) { + item->seed_yp_from_f(f, yp); + } +} + +int nrndae_complete_yp_from_forcing(double* yp, const std::vector& forcing) { + int flags = 0; + if (!yp) { + return 0; + } + std::vector* prl = net_cvode_instance ? net_cvode_instance->playrec_list() + : nullptr; + extern double t; + + for (NrnDAE* item: nrndae_list) { + auto* lm = dynamic_cast(item); + if (!lm) { + continue; + } + const int n = lm->size(); + if (n <= 0) { + continue; + } + std::vector bdot(n, 0.); + bool have_bdot = false; + int src = 0; + + // A1/A2: continuous Vector.play → components of b + if (prl && !forcing.empty()) { + for (const auto& e: forcing) { + if (e.playrec_index < 0 || e.playrec_index >= (int) prl->size()) { + continue; + } + PlayRecord* pr = (*prl)[e.playrec_index]; + if (!pr || pr->type() != VecPlayContinuousType) { + continue; + } + auto* vpc = static_cast(pr); + double* target = nullptr; + if (vpc->pd_) { + target = static_cast(vpc->pd_); + } + if (!target) { + continue; + } + for (int i = 0; i < n; ++i) { + if (lm->b_element_is(i, target)) { + bdot[i] = e.deriv; + have_bdot = true; + src |= NRN_IC_FORCING_PLAY; + } + } + } + } + + // A4: dforce / bdot vector (overrides play); else FD if f_callable and no play + const bool have_dforce = lm->bdot_vec() || lm->dforce_callable(); + if (have_dforce || (lm->f_callable() && !have_bdot)) { + std::vector bdot_df(n, 0.); + if (lm->fill_bdot_for_ic(t, bdot_df.data())) { + for (int i = 0; i < n; ++i) { + bdot[i] = bdot_df[i]; + } + have_bdot = true; + if (have_dforce) { + src |= NRN_IC_FORCING_DFORCE; + } else { + src |= NRN_IC_FORCING_FD; + } + } + } + + if (have_bdot) { + lm->complete_yp_from_bdot(bdot.data(), yp); + flags |= src | NRN_IC_FORCING_APPLIED; + } + } + return flags; +} + +void nrndae_append_dforce_to_forcing_list(double tt, std::vector& out) { + int lm_i = 0; + for (NrnDAE* item: nrndae_list) { + auto* lm = dynamic_cast(item); + if (!lm || lm->size() <= 0) { + ++lm_i; + continue; + } + if (!lm->bdot_vec() && !lm->dforce_callable() && !lm->f_callable()) { + ++lm_i; + continue; + } + std::vector bdot(lm->size(), 0.); + if (!lm->fill_bdot_for_ic(tt, bdot.data())) { + ++lm_i; + continue; + } + for (int i = 0; i < lm->size(); ++i) { + NrnForcingTPlus e{}; + e.deriv = bdot[i]; + e.playrec_index = -1 - lm_i; + e.ubound_index = i; + if (lm->dforce_callable()) { + std::snprintf(e.label, sizeof e.label, "LM[%d].dforce b'[%d]", lm_i, i); + } else if (lm->bdot_vec()) { + std::snprintf(e.label, sizeof e.label, "LM[%d].bdot[%d]", lm_i, i); + } else { + std::snprintf(e.label, sizeof e.label, "LM[%d].bdot_fd[%d]", lm_i, i); + } + out.push_back(e); + } + ++lm_i; + } +} + void NrnDAE::rhs() { NrnThread* _nt = nrn_threads; diff --git a/src/nrniv/nrndae.h b/src/nrniv/nrndae.h index 6b3bbe765c..8bf75515dd 100644 --- a/src/nrniv/nrndae.h +++ b/src/nrniv/nrndae.h @@ -12,6 +12,7 @@ #pragma once #include "ivocvect.h" #include "matrixmap.h" +#include "vecplay_tplus.h" #include "neuron/container/data_handle.hpp" @@ -59,6 +60,16 @@ class NrnDAE { */ void dkres(double* y, double* yprime, double* delta); + /** + * Seed local contributions to yp from C*yp = f for simple C structure + * (identity, pure diagonal, or single-column lag). Floating mutual C + * uses a zero common-mode gauge when the row is a pure difference stamp. + * + * @param f full-system rhs f(y) (same layout as IDA residual f-part) + * @param yp full-system y' to fill (may already hold membrane rates) + */ + void seed_yp_from_f(double* f, double* yp); + /** * Initialize the dynamics. * @@ -158,6 +169,7 @@ class NrnDAE { */ virtual void alloc_(int size, int start, int nnode, Node** nodes, int* elayer); + protected: /// the matrix \f[$C$ in $C y' = f(y)$\f] MatrixMap* c_; @@ -218,5 +230,37 @@ void nrndae_register(NrnDAE* n); */ void nrndae_deregister(NrnDAE* n); +/** + * Battery-style IC: for each LinearModelAddition, project states + * (capacitors → held Δv voltage sources, etc.). Returns 0 if all projections OK. + */ +int nrndae_battery_ic_project(); + +/** + * After y is fixed, seed yp from diagonal / single-column C rows so that + * C*yp ≈ f on LinearMechanism equations (used by dae_init_mode 3). + * f and yp are full IDA state vectors (neq). + */ +void nrndae_seed_yp_from_f(double* f, double* yp); + +/** + * A2: using continuous Vector.play forcing t+ (u'), map db/dt into each + * LinearMechanism and complete free yp components (null space of C). + * forcing may be empty (no-op). yp is the full IDA y' vector. + */ +// Returns bitmask of sources used for free y': 1=play, 2=dforce/bdot, 4=FD, 8=applied. +int nrndae_complete_yp_from_forcing(double* yp, const std::vector& forcing); + +/** Append LinearMechanism.dforce / FD b' entries for IC audit (A4). */ +void nrndae_append_dforce_to_forcing_list(double tt, std::vector& out); + +// Source bits for nrndae_complete_yp_from_forcing +enum { + NRN_IC_FORCING_PLAY = 1, + NRN_IC_FORCING_DFORCE = 2, + NRN_IC_FORCING_FD = 4, + NRN_IC_FORCING_APPLIED = 8 +}; + typedef std::list NrnDAEPtrList; typedef NrnDAEPtrList::const_iterator NrnDAEPtrListIterator; diff --git a/src/nrniv/vecplay_tplus.h b/src/nrniv/vecplay_tplus.h new file mode 100644 index 0000000000..c0f4d672bf --- /dev/null +++ b/src/nrniv/vecplay_tplus.h @@ -0,0 +1,138 @@ +#pragma once + +/** + * @file vecplay_tplus.h + * @brief Continuous Vector.play right-limit value and classical derivative. + * + * **Forcing \(t^+\) info:** the right-limit value \(u(t^+)\) and classical + * derivative \(u'(t^+)\) of an exogenous drive after a discontinuity (or at + * finitialize). In the geometric / DAE literature this pair is the **1-jet** + * of \(u\) at \(t^+\). + * + * Geometry matches `VecPlayContinuous::interpolate` for a given active upper + * knot index `ubound_index` (use `n - 1` when the full `t` vector is active): + * - \(t < t_0\): hold \(y_0\), derivative 0 + * - interior / at a knot (including \(t = t_0\)): linear segment; **outgoing** slope + * - \(t \ge t_{\mathrm{ubound}}\): linear extrapolation of the last two points + * on that bound (constant last value only if that segment is flat, or if + * `ubound_index == 0`) + * - coincident knot times \(t_i = t_{i+1}\): value average, derivative 0 + */ + +#include +#include +#include +#include +#include + +/** One continuous Vector.play forcing \(t^+\) sample (value + classical deriv). */ +struct NrnForcingTPlus { + double value{}; + double deriv{}; + int playrec_index{-1}; + int ubound_index{-1}; + char label[128]{}; // Vector object name if available +}; + +/** + * Collect forcing \(t^+\) info from every VecPlayContinuous in the playrec list. + * Call at IDA IC time with post-event continuous play state (ubound_index_ current). + * Clears and fills `out`. Returns (int)out.size(). + */ +int nrn_collect_forcing_tplus(double tt, std::vector& out); + +/** Print collected forcing \(t^+\) lines to f (audit / debug). */ +void nrn_dump_forcing_tplus(FILE* f, double tt, const std::vector& entries); + +/** + * Continuous play value and classical derivative at time `tt`. + * + * @param n number of samples (length of y and t) + * @param y sample values + * @param t sample times (nondecreasing) + * @param tt query time + * @param ubound_index active upper knot (same role as VecPlayContinuous::ubound_index_) + * @param value output u(tt); may be null + * @param deriv output u'(tt) classical; may be null + * @return 0 on success, -1 if n < 1 or y/t null + */ +inline int nrn_vecplay_continuous_tplus(int n, + const double* y, + const double* t, + double tt, + int ubound_index, + double* value, + double* deriv) { + if (n < 1 || !y || !t) { + return -1; + } + if (ubound_index < 0) { + ubound_index = 0; + } + if (ubound_index >= n) { + ubound_index = n - 1; + } + + auto set = [&](double v, double d) { + if (value) { + *value = v; + } + if (deriv) { + *deriv = d; + } + }; + + // Match VecPlayContinuous::interpolate when tt >= t[ubound] + if (tt >= t[ubound_index]) { + if (ubound_index == 0) { + set(y[0], 0.); + return 0; + } + const double t0 = t[ubound_index - 1]; + const double t1 = t[ubound_index]; + const double x0 = y[ubound_index - 1]; + const double x1 = y[ubound_index]; + if (t0 == t1) { + set(0.5 * (x0 + x1), 0.); + return 0; + } + const double slope = (x1 - x0) / (t1 - t0); + // linear extrapolation when tt > t1 + set(x0 + slope * (tt - t0), slope); + return 0; + } + + // Strictly before the first knot: hold (right-limit at t0 uses outgoing segment). + if (tt < t[0]) { + set(y[0], 0.); + return 0; + } + + // Find last like VecPlayContinuous::search: first index with t[last] > tt + // (unique times ⇒ outgoing segment at a knot, including tt == t[0]). + int last = 1; + while (last < ubound_index && tt >= t[last]) { + ++last; + } + while (last > 0 && tt < t[last - 1]) { + --last; + } + if (last < 1) { + last = 1; + } + if (last > ubound_index) { + last = ubound_index; + } + + const double t0 = t[last - 1]; + const double t1 = t[last]; + const double x0 = y[last - 1]; + const double x1 = y[last]; + if (t0 == t1) { + set(0.5 * (x0 + x1), 0.); + return 0; + } + const double slope = (x1 - x0) / (t1 - t0); + set(x0 + slope * (tt - t0), slope); + return 0; +} diff --git a/src/nrniv/vrecord.cpp b/src/nrniv/vrecord.cpp index 131595070c..a44117d08c 100644 --- a/src/nrniv/vrecord.cpp +++ b/src/nrniv/vrecord.cpp @@ -11,6 +11,7 @@ #include "ocpointer.h" #include "vrecitem.h" +#include "vecplay_tplus.h" #include "netcvode.h" #include "cvodeobj.h" @@ -353,36 +354,107 @@ void VecPlayContinuous::continuous(double tt) { } } +void VecPlayContinuous::forcing_tplus(double tt, double* value, double* deriv) const { + // Full sample arrays; active upper knot is ubound_index_ (same as interpolate). + const int n = t_ ? static_cast(t_->size()) : 0; + // IvocVect stores doubles contiguously via data() / &elem(0) + const double* y = (y_ && n > 0) ? &y_->elem(0) : nullptr; + const double* tv = (t_ && n > 0) ? &t_->elem(0) : nullptr; + if (nrn_vecplay_continuous_tplus(n, y, tv, tt, ubound_index_, value, deriv) != 0) { + if (value) { + *value = 0.; + } + if (deriv) { + *deriv = 0.; + } + } +} + +int nrn_collect_forcing_tplus(double tt, std::vector& out) { + out.clear(); + if (!net_cvode_instance) { + return 0; + } + std::vector* prl = net_cvode_instance->playrec_list(); + if (!prl) { + return 0; + } + for (std::size_t i = 0; i < prl->size(); ++i) { + PlayRecord* pr = (*prl)[i]; + if (!pr || pr->type() != VecPlayContinuousType) { + continue; + } + auto* vpc = static_cast(pr); + NrnForcingTPlus e; + e.playrec_index = static_cast(i); + e.ubound_index = vpc->ubound_index_; + vpc->forcing_tplus(tt, &e.value, &e.deriv); + e.label[0] = '\0'; + if (vpc->y_ && vpc->y_->obj_) { + const char* nm = hoc_object_name(vpc->y_->obj_); + if (nm) { + std::snprintf(e.label, sizeof e.label, "%s", nm); + } + } + if (e.label[0] == '\0') { + std::snprintf(e.label, sizeof e.label, "VecPlayContinuous[%d]", e.playrec_index); + } + out.push_back(e); + } + return static_cast(out.size()); +} + +void nrn_dump_forcing_tplus(FILE* f, double tt, const std::vector& entries) { + if (!f) { + return; + } + fprintf(f, + "--- forcing t+ info (t=%.15g) n_play_continuous=%d ---\n", + tt, + (int) entries.size()); + fprintf(f, " (right-limit u and classical u'; 1-jet of exogenous continuous Vector.play)\n"); + if (entries.empty()) { + fprintf(f, " (none)\n"); + return; + } + fprintf(f, " %4s %10s %16s %16s %s\n", "idx", "ubound", "u(t+)", "u'(t+)", "label"); + for (const auto& e: entries) { + fprintf(f, + " %4d %10d %16.8g %16.8g %s\n", + e.playrec_index, + e.ubound_index, + e.value, + e.deriv, + e.label); + } +} + double VecPlayContinuous::interpolate(double tt) { + // Keep last_index_ cache in sync with historical search side effects, then + // evaluate with the shared t⁺ geometry (value only). Value at/after t0 + // matches prior play (linear segments); only the t⁺ *derivative* treatment + // of t==t0 differs for forcing (outgoing slope). if (tt >= t_->elem(ubound_index_)) { last_index_ = ubound_index_; if (last_index_ == 0) { - // printf("return last tt=%g ubound=%g y=%g\n", tt, t_->elem(ubound_index_), - // y_->elem(last_index_)); return y_->elem(last_index_); } - } else if (tt <= t_->elem(0)) { + } else if (tt < t_->elem(0)) { last_index_ = 0; - // printf("return elem(0) tt=%g t0=%g y=%g\n", tt, t_->elem(0), y_->elem(0)); return y_->elem(0); + } else if (tt == t_->elem(0)) { + last_index_ = (ubound_index_ > 0) ? 1 : 0; + if (last_index_ == 0) { + return y_->elem(0); + } + // fall through to shared evaluator via tplus (value = y0 on first segment) } else { search(tt); } - double x0 = y_->elem(last_index_ - 1); - double x1 = y_->elem(last_index_); - double t0 = t_->elem(last_index_ - 1); - double t1 = t_->elem(last_index_); - // printf("IvocVectRecorder::continuous tt=%g t0=%g t1=%g theta=%g x0=%g x1=%g\n", tt, t0, t1, - // (tt - t0)/(t1 - t0), x0, x1); - if (t0 == t1) { - return (x0 + x1) / 2.; - } - return interp((tt - t0) / (t1 - t0), x0, x1); -#if 0 - // dt - double theta = tt/dt_ - last_index_; - interp(theta, x0, x1); -#endif + double value = 0.; + const int n = static_cast(t_->size()); + nrn_vecplay_continuous_tplus(n, &y_->elem(0), &t_->elem(0), tt, ubound_index_, &value, nullptr); + return value; } void VecPlayContinuous::search(double tt) { diff --git a/src/nrnoc/extcelln.cpp b/src/nrnoc/extcelln.cpp index 3ad4db7fa7..ee27515f31 100644 --- a/src/nrnoc/extcelln.cpp +++ b/src/nrnoc/extcelln.cpp @@ -3,10 +3,13 @@ #include #include +#include #include "section.h" #include "nrniv_mf.h" #include "hocassrt.h" #include "parse.hpp" +#include "nrn_ansi.h" +#include "multicore.h" extern int nrn_use_daspk_; @@ -740,5 +743,125 @@ ndesave[i].rhs[j] -= ndesave[i].v[k]*ndesave[i].m[j][(nlayer)+k-j]; #endif #endif +/* + * Battery-style consistent IC for extracellular + membrane capacitance. + * Hold continuous content (Vm and xc layer drops / grounded xc voltages) via + * large conductances on the cj=0 (resistive) matrix, solve, update, then + * exactly restore holds while preserving algebraic common-mode shift when + * the outermost layer has no capacitance to ground. + * + * Called from Daspk battery IC (dae_init_mode 3) after LinearMechanism project. + */ +void nrn_extracellular_battery_ic() { + if (!use_sparse13) { + return; + } + auto const sorted_token = nrn_ensure_model_data_are_sorted(); + for (int it = 0; it < nrn_nthread; ++it) { + NrnThread* nt = nrn_threads + it; + Memb_list* ml = nt->_ecell_memb_list; + if (!ml || ml->nodecount == 0) { + continue; + } + const int cnt = ml->nodecount; + const int nl = nrn_nlayer_extracellular; + + // Save continuous holds from current state + std::vector vm_hold(cnt); + std::vector vext_hold(cnt * nl); + for (int i = 0; i < cnt; ++i) { + Node* nd = ml->nodelist[i]; + Extnode* nde = nd->extnode; + vm_hold[i] = nd->v(); + for (int j = 0; j < nl; ++j) { + vext_hold[i * nl + j] = nde->v[j]; + } + } + + double const cj_sav = nt->cj; + double const dt_sav = nt->_dt; + nt->cj = 0.0; + nt->_dt = 1e9; + + setup_tree_matrix(sorted_token, *nt); + + // Stiff springs on actual_d / sparse entries. nrn_solve copies + // actual_d → sparse13 diagonals, so node holds must update actual_d. + // Skip Vm hold on zero-area nodes (algebraic; electrode at loc 0/1). + constexpr double ghold = 1e9; + for (int i = 0; i < cnt; ++i) { + Node* nd = ml->nodelist[i]; + Extnode* nde = nd->extnode; + const int vi = nd->v_node_index; + const bool hold_vm = (NODEAREA(nd) > 0.); + if (hold_vm) { + // Membrane: continuous Vm ⇒ dvi - dvx0 = 0 + nt->actual_d(vi) += ghold; + if (nd->_d_matelm) { + *nd->_d_matelm += ghold; + } + *nde->_d[0] += ghold; + *nde->_x12[0] -= ghold; + *nde->_x21[0] -= ghold; + } + + for (int j = 0; j < nl; ++j) { + double const xc = *nde->param[xc_index_ext(j)]; + if (xc <= 0.0) { + continue; + } + if (j == nl - 1) { + // Grounded layer capacitance: continuous vext[last] + *nde->_d[j] += ghold; + } else { + // Inter-layer drop continuous + *nde->_d[j] += ghold; + *nde->_d[j + 1] += ghold; + *nde->_x12[j + 1] -= ghold; + *nde->_x21[j + 1] -= ghold; + } + } + } + + nrn_solve(nt); + nrn_update_voltage(sorted_token, *nt); + + // Exact hold restore; allow common-mode shift only if outer xc == 0 + for (int i = 0; i < cnt; ++i) { + Node* nd = ml->nodelist[i]; + Extnode* nde = nd->extnode; + const bool hold_vm = (NODEAREA(nd) > 0.); + double const xc_last = *nde->param[xc_index_ext(nl - 1)]; + double alpha = 0.0; + if (xc_last <= 0.0) { + // Algebraic outer layer: keep solved absolute level of vext[0] + alpha = nde->v[0] - vext_hold[i * nl + 0]; + } + // Rebuild vext from held drops (or held absolutes) + alpha + // Start from outermost continuous or algebraic anchor + if (xc_last > 0.0) { + nde->v[nl - 1] = vext_hold[i * nl + (nl - 1)]; + } else { + nde->v[nl - 1] = vext_hold[i * nl + (nl - 1)] + alpha; + } + for (int j = nl - 2; j >= 0; --j) { + double const xc = *nde->param[xc_index_ext(j)]; + if (xc > 0.0) { + double const drop = vext_hold[i * nl + j] - vext_hold[i * nl + j + 1]; + nde->v[j] = nde->v[j + 1] + drop; + } else { + nde->v[j] = vext_hold[i * nl + j] + alpha; + } + } + if (hold_vm) { + nd->v() = vm_hold[i]; + } + // else: leave solved Vm (algebraic zero-area end free) + } + + nt->cj = cj_sav; + nt->_dt = dt_sav; + } +} #endif /*EXTRACELLULAR*/ diff --git a/src/nrnoc/nrn_ansi.h b/src/nrnoc/nrn_ansi.h index 0ac243ec6d..5266546274 100644 --- a/src/nrnoc/nrn_ansi.h +++ b/src/nrnoc/nrn_ansi.h @@ -92,6 +92,12 @@ extern void hoc_level_pushsec(Section*); void nrn_ba(neuron::model_sorted_token const&, NrnThread&, int); extern void nrn_rhs_ext(NrnThread*); extern void nrn_setup_ext(NrnThread*); +/** Battery-style IC: hold positive-area node V; free zero-area algebraics (end PP). */ +void nrn_cable_battery_ic(); +#if EXTRACELLULAR +/** Battery-style IC: hold Vm and xc continuous content when extracellular present. */ +void nrn_extracellular_battery_ic(); +#endif void nrn_cap_jacob(neuron::model_sorted_token const&, NrnThread*, Memb_list*); void nrn_div_capacity(neuron::model_sorted_token const&, NrnThread*, Memb_list*); void nrn_mul_capacity(neuron::model_sorted_token const&, NrnThread*, Memb_list*); @@ -99,6 +105,9 @@ extern void clear_point_process_struct(Prop* p); extern void ext_con_coef(void); extern void nrn_multisplit_ptr_update(void); extern void nrn_use_daspk(int); +// When non-zero, sparse13 factor failures do not hoc_execerror; see nrn_sparse13_factor_error(). +extern int nrn_sparse13_soft_fail; +extern int nrn_sparse13_factor_error(); extern void nrn_update_ps2nt(void); neuron::model_sorted_token nrn_ensure_model_data_are_sorted(); extern void activstim_rhs(void); diff --git a/src/nrnoc/solve.cpp b/src/nrnoc/solve.cpp index a42fbe960a..f181731128 100644 --- a/src/nrnoc/solve.cpp +++ b/src/nrnoc/solve.cpp @@ -330,6 +330,13 @@ static void dashes(Section* sec, int offset, int first) { } /* solve the matrix equation */ +int nrn_sparse13_soft_fail = 0; +static int nrn_sparse13_factor_error_ = 0; + +int nrn_sparse13_factor_error() { + return nrn_sparse13_factor_error_; +} + void nrn_solve(NrnThread* _nt) { #if 0 printf("\nnrn_solve enter %lx\n", (long)_nt); @@ -359,6 +366,10 @@ void nrn_solve(NrnThread* _nt) { update_sp13_mat_based_on_actual_d(_nt); e = spFactor(_nt->_sp13mat); if (e != spOKAY) { + if (nrn_sparse13_soft_fail) { + nrn_sparse13_factor_error_ = e; + return; + } switch (e) { case spZERO_DIAG: hoc_execerror("spFactor error:", "Zero Diagonal"); @@ -368,6 +379,7 @@ void nrn_solve(NrnThread* _nt) { hoc_execerror("spFactor error:", "Singular"); } } + nrn_sparse13_factor_error_ = 0; update_sp13_rhs_based_on_actual_rhs(_nt); spSolve(_nt->_sp13mat, _nt->_sp13_rhs, _nt->_sp13_rhs); update_actual_d_based_on_sp13_mat(_nt); diff --git a/src/nrnoc/treeset.cpp b/src/nrnoc/treeset.cpp index e24a45b196..1216bbfd15 100644 --- a/src/nrnoc/treeset.cpp +++ b/src/nrnoc/treeset.cpp @@ -29,6 +29,7 @@ #include #include +#include #include @@ -604,6 +605,108 @@ void setup_tree_matrix(neuron::model_sorted_token const& cache_token, NrnThread& nrn_nonvint_block_conductance(nt.end, nt.node_d_storage(), nt.id); } +/* + * Battery-style free-y IC for cable algebraics (mode 3). + * + * Hold voltages on nodes with positive area (capacitive membrane content); + * leave zero-area nodes free so electrode current at loc 0/1 can set absolute + * levels through the axial network. Uses cj=0 resistive matrix + stiff + * springs (delta-v ≈ 0) on held nodes, then exact restore of those voltages. + * + * Called from Daspk battery IC before extracellular battery projection. + * No-op when not using sparse13 (IDA/DASPK path). + */ +void nrn_cable_battery_ic() { + if (!use_sparse13) { + return; + } + auto const sorted_token = nrn_ensure_model_data_are_sorted(); + for (int it = 0; it < nrn_nthread; ++it) { + NrnThread* nt = nrn_threads + it; + // Extracellular uses a separate battery path (Vm / layer holds). + // Applying absolute-v holds here fights vi/vext coordinates. + if (nt->_ecell_memb_list && nt->_ecell_memb_list->nodecount > 0) { + continue; + } + const int n = nt->end; + if (n <= 0) { + continue; + } + std::vector v_hold(n); + std::vector held(n, 0); + // Hold nodes that carry membrane capacitance in the CAP mechanism list. + // Section ends (loc 0/1) are omitted from CAP ⇒ free algebraics for + // electrode current. NODEAREA is not a reliable zero-end test here + // (can be non-zero on connection nodes while CAP still omits them). + for (NrnThreadMembList* tml = nt->tml; tml; tml = tml->next) { + if (tml->index != CAP) { + continue; + } + Memb_list* ml = tml->ml; + for (int j = 0; j < ml->nodecount; ++j) { + Node* nd = ml->nodelist[j]; + int i = nd->v_node_index; + if (i >= 0 && i < n) { + held[i] = 1; + v_hold[i] = nd->v(); + } + } + break; + } + int nheld = 0, nfree = 0; + for (int i = 0; i < n; ++i) { + if (held[i]) { + ++nheld; + } else { + ++nfree; + } + } + // Need CAP content to hold and free algebraics to move. Pure LM + // models (no CAP) must not run this path — cj=0 solve can singularize. + if (nheld == 0 || nfree == 0) { + continue; + } + + double const cj_sav = nt->cj; + double const dt_sav = nt->_dt; + nt->cj = 0.0; + nt->_dt = 1e9; + + setup_tree_matrix(sorted_token, *nt); + + // Add hold springs on actual_d (nrn_solve copies actual_d → sparse13). + constexpr double ghold = 1e9; + for (int i = 0; i < n; ++i) { + if (!held[i]) { + continue; + } + nt->actual_d(i) += ghold; + } + + int soft_sav = nrn_sparse13_soft_fail; + nrn_sparse13_soft_fail = 1; + nrn_solve(nt); + int ferr = nrn_sparse13_factor_error(); + nrn_sparse13_soft_fail = soft_sav; + if (ferr) { + // Leave voltages unchanged on factor failure (e.g. LM+cable mix). + nt->cj = cj_sav; + nt->_dt = dt_sav; + continue; + } + nrn_update_voltage(sorted_token, *nt); + + for (int i = 0; i < n; ++i) { + if (held[i]) { + nt->_v_node[i]->v() = v_hold[i]; + } + } + + nt->cj = cj_sav; + nt->_dt = dt_sav; + } +} + /* membrane mechanisms needed by other mechanisms (such as Eion by HH) may require that the needed mechanism appear before it in the list. (because ina must be initialized before it is incremented by HH) diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index a87f557a01..b56839b5e7 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -26,7 +26,9 @@ add_executable( unit_tests/utils/enumerate.cpp unit_tests/utils/Sprintf.cpp unit_tests/oc/hoc_interpreter.cpp + unit_tests/vecplay_tplus.cpp cover/unit_tests/cover.cpp) +target_include_directories(testneuron PRIVATE ${PROJECT_SOURCE_DIR}/src/nrniv) set(catch2_targets testneuron) if(NRN_ENABLE_THREADS) add_executable(nrn-benchmarks common/catch2_main.cpp benchmarks/threads/test_multicore.cpp) diff --git a/test/hoctests/pwlclamp.mod b/test/hoctests/pwlclamp.mod new file mode 100644 index 0000000000..fcbb748bf7 --- /dev/null +++ b/test/hoctests/pwlclamp.mod @@ -0,0 +1,122 @@ +COMMENT +Test-only piecewise-linear (PWL) electrode current stimulus for IDA IC work. + +Mirrors continuous Vector.play tables used in test_ida_init_mode.py +(A1 ramp, A3 istep/kink/...) so Section models can receive the same jump and +kink waveforms as LinearMechanism force play. + +IClamp-like: ELECTRODE_CURRENT i, at_time at each knot, BREAKPOINT assigns i(t). + +Right-continuous at jumps (last knot index with tk[j] <= t). Past the last +knot: linear extrapolation of the last non-zero-length segment when possible; +otherwise hold the last value. (Matches continuous-play end rules used in A3.) + +NOT a product mechanism. Does NOT register classical di/dt into Plan A +forcing t+ (play / LM.dforce remain the free-y' 1-jet sources). + +Max NKT knots; set nkt and tk[i], ik[i] from Python/HOC. +ENDCOMMENT + +NEURON { + POINT_PROCESS PWLClamp + RANGE i, nkt + ELECTRODE_CURRENT i +} + +UNITS { + (nA) = (nanoamp) +} + +PARAMETER { + nkt = 0 + tk[16] (ms) + ik[16] (nA) +} + +ASSIGNED { + i (nA) +} + +INITIAL { + i = ival(t) +} + +BREAKPOINT { + LOCAL j + if (nkt > 0) { + FROM j = 0 TO nkt-1 { + at_time(tk[j]) + } + } + i = ival(t) +} + +PROCEDURE set_knot(j, tval (ms), ival_ (nA)) { + : Python/HOC cannot assign PARAMETER arrays directly on this PP; + : use set_knot(j, t, i) for j = 0 .. nkt-1 after setting nkt. + if (j >= 0 && j < 16) { + tk[j] = tval + ik[j] = ival_ + } +} + +FUNCTION get_tk(j) (ms) { + if (j >= 0 && j < 16) { + get_tk = tk[j] + } else { + get_tk = 0 + } +} + +FUNCTION get_ik(j) (nA) { + if (j >= 0 && j < 16) { + get_ik = ik[j] + } else { + get_ik = 0 + } +} + +FUNCTION ival(tt (ms)) (nA) { + LOCAL j, jmax, t0, t1, i0, i1 + if (nkt < 1) { + ival = 0 + } else if (tt < tk[0]) { + : hold first knot value before the first time (play-like) + ival = ik[0] + } else if (tt >= tk[nkt-1]) { + : at/after last knot: extrap last non-degenerate segment if any + if (nkt >= 2) { + j = nkt - 2 + WHILE (j > 0 && tk[j+1] <= tk[j]) { + j = j - 1 + } + if (tk[j+1] > tk[j]) { + ival = ik[j+1] + (ik[j+1] - ik[j]) * (tt - tk[j+1]) / (tk[j+1] - tk[j]) + } else { + ival = ik[nkt-1] + } + } else { + ival = ik[0] + } + } else { + : last index j with tk[j] <= tt (right-continuous at coincident knots) + jmax = 0 + FROM j = 0 TO nkt-1 { + if (tk[j] <= tt) { + jmax = j + } + } + if (jmax >= nkt-1) { + ival = ik[nkt-1] + } else if (tk[jmax+1] > tk[jmax]) { + t0 = tk[jmax] + t1 = tk[jmax+1] + i0 = ik[jmax] + i1 = ik[jmax+1] + ival = i0 + (i1 - i0) * (tt - t0) / (t1 - t0) + } else { + : degenerate (jump) segment: sit on right value at this index + ival = ik[jmax] + } + } +} diff --git a/test/hoctests/tests/test_ida_init_mode.py b/test/hoctests/tests/test_ida_init_mode.py new file mode 100644 index 0000000000..1e789f2b34 --- /dev/null +++ b/test/hoctests/tests/test_ida_init_mode.py @@ -0,0 +1,786 @@ +"""Phase 0: CVode.dae_init_mode and IDA_Y_INIT with LinearMechanism. + +Each model case runs in a subprocess so LinearMechanism teardown does not +interact across tests (NrnDAE / sparse13 matrix reallocation). +""" + +from neuron import h +import math +import os +import subprocess +import sys + +h.load_file("stdrun.hoc") +cvode = h.CVode() + + +def test_dae_init_mode_api(): + assert cvode.dae_init_mode() == 0 + assert cvode.dae_init_mode(1) == 1 + assert cvode.dae_init_mode(2) == 2 + assert cvode.dae_init_mode(3) == 3 + assert cvode.dae_init_mode(0) == 0 + + +def test_dae_init_audit_api(): + assert cvode.dae_init_audit() == 0 + assert cvode.dae_init_audit(1) == 1 + assert cvode.dae_init_audit(2, 5.0) == 2 + assert cvode.dae_init_audit(0) == 0 + assert cvode.dae_init_audit_file() == 0 + # empty path → stdout; non-empty sets file mode flag + assert cvode.dae_init_audit_file("") == 0 + + +def test_dae_init_stats_api(): + """A5: dae_init_stats reset / vector fill.""" + code = r""" +from neuron import h +h.load_file('stdrun.hoc') +cv = h.CVode() +cv.dae_init_stats(1) # reset +v = h.Vector() +n = cv.dae_init_stats(v) +assert n == 8 and v.size() == 8, (n, v.size()) +assert v[0] == 0 # no reinits yet +# series CR + play ramp once +c = h.Matrix(2, 2) +g = h.Matrix(2, 2) +y = h.Vector(2) +b = h.Vector(2) +c.setval(0, 0, 1); c.setval(0, 1, -1) +c.setval(1, 0, -1); c.setval(1, 1, 1) +g.setval(1, 1, 1) +lm = h.LinearMechanism(c, g, y, b) +tvec = h.Vector([0, 1, 1, 2]) +ivec = h.Vector([0, 0, 0.5, 1.0]) +ivec.play(b._ref_x[0], tvec, True) +h.cvode_active(True) +cv.use_daspk(1) +cv.dae_init_mode(3) +h.finitialize(0) +h.continuerun(1.0) +cv.re_init() +cv.dae_init_stats(v) +assert v[0] >= 2, list(v) # finitialize + re_init at least +assert v[1] >= 1, list(v) # mode3 ok +assert v[3] >= 1, list(v) # play free y' +assert int(v[6]) == 3, list(v) # last path mode 3 +assert int(v[7]) & 8, list(v) # APPLIED bit +print('ok') +""" + assert "ok" in _run_isolated(code) + + +def _sanitizer_child_env(env=None): + """Re-apply sanitizer preload for macOS SIP (see NeuronTestHelper.cmake).""" + env = os.environ.copy() if env is None else env + try: + env[os.environ["NRN_SANITIZER_PRELOAD_VAR"]] = os.environ[ + "NRN_SANITIZER_PRELOAD_VAL" + ] + except KeyError: + pass + return env + + +def _run_isolated(code: str): + env = _sanitizer_child_env() + exe = os.environ.get("NRN_PYTHON_EXECUTABLE", sys.executable) + r = subprocess.run( + [exe, "-c", code], + env=env, + capture_output=True, + text=True, + ) + if r.returncode != 0: + raise AssertionError(f"stdout:\n{r.stdout}\nstderr:\n{r.stderr}") + return r.stdout + + +def test_pure_resistive_all_modes_isolated(): + code = r""" +from neuron import h +import math +h.load_file('stdrun.hoc') +cvode = h.CVode() +c = h.Matrix(2, 2) +g = h.Matrix(2, 2) +y = h.Vector(2) +y0 = h.Vector(2) +b = h.Vector([1.0, 2.0]) +g.setval(0, 0, 1.0) +g.setval(1, 1, 0.5) +lm = h.LinearMechanism(c, g, y, y0, b) +h.cvode_active(True) +for mode in (0, 1, 2, 3): + cvode.dae_init_mode(mode) + h.finitialize(0.0) + assert math.isclose(y[0], 1.0, rel_tol=1e-6, abs_tol=1e-6), mode + assert math.isclose(y[1], 4.0, rel_tol=1e-6, abs_tol=1e-6), mode +print('ok') +""" + assert "ok" in _run_isolated(code) + + +def test_series_cr_mode0_mode1_isolated(): + code = r""" +from neuron import h +import math +h.load_file('stdrun.hoc') +cvode = h.CVode() +I, R, C = 1.0, 2.0, 1.0 +c = h.Matrix(2, 2) +g = h.Matrix(2, 2) +y = h.Vector(2) +y0 = h.Vector(2) +b = h.Vector([I, 0.0]) +c.setval(0, 0, C) +c.setval(0, 1, -C) +c.setval(1, 0, -C) +c.setval(1, 1, C) +g.setval(1, 1, 1.0 / R) +lm = h.LinearMechanism(c, g, y, y0, b) +h.cvode_active(True) +IR = I * R +for mode in (0, 1): + cvode.dae_init_mode(mode) + h.finitialize(0.0) + assert math.isclose(y[0] - y[1], 0.0, abs_tol=1e-6), mode + assert math.isclose(y[0], IR, rel_tol=1e-4, abs_tol=1e-4), (mode, y[0]) + assert math.isclose(y[1], IR, rel_tol=1e-4, abs_tol=1e-4), (mode, y[1]) +cvode.dae_init_mode(1) +h.finitialize(0.0) +h.continuerun(0.5) +assert math.isclose(y[1], IR, rel_tol=1e-3, abs_tol=1e-3) +assert y[0] > y[1] +print('ok') +""" + assert "ok" in _run_isolated(code) + + +def test_series_cr_battery_mode3_isolated(): + """Battery IC: C→V hold; absolute V jumps to I*R, Δv held at 0.""" + code = r""" +from neuron import h +import math +h.load_file('stdrun.hoc') +cvode = h.CVode() +I, R, Cval = 1.0, 2.0, 1.0 +c = h.Matrix(2, 2) +g = h.Matrix(2, 2) +y = h.Vector(2) +y0 = h.Vector(2) +b = h.Vector([I, 0.0]) +c.setval(0, 0, Cval) +c.setval(0, 1, -Cval) +c.setval(1, 0, -Cval) +c.setval(1, 1, Cval) +g.setval(1, 1, 1.0 / R) +lm = h.LinearMechanism(c, g, y, y0, b) +h.cvode_active(True) +cvode.dae_init_mode(3) +h.finitialize(0.0) +IR = I * R +assert math.isclose(y[0] - y[1], 0.0, abs_tol=1e-9), (y[0], y[1]) +assert math.isclose(y[0], IR, rel_tol=1e-6, abs_tol=1e-6), y[0] +assert math.isclose(y[1], IR, rel_tol=1e-6, abs_tol=1e-6), y[1] +h.continuerun(0.5) +assert math.isclose(y[1], IR, rel_tol=1e-3, abs_tol=1e-3) +assert y[0] > y[1] - 1e-9 +print('ok') +""" + assert "ok" in _run_isolated(code) + + +def test_opamp_tau_battery_holds_output_voltage(): + """OpAmp lag (tau>0): C[o][k]=tau. Hold v_out, not bogus v_out-I. + + Before the inductor/opamp extension, mode 3 treated (v,I) as a floating + capacitor and drove v→~1 (algebraic follower). Continuous IC from rest + must keep v_out≈0 when + is stepped via b. + """ + code = r""" +from neuron import h +import math +h.load_file('stdrun.hoc') +cvode = h.CVode() +# LinearCircuit OpAmp lag stamp (tau>0): y0=v_out, y1=I_out +# tau*v' + (gain+1)*v = gain*Vplus +# I + v/R = 0 +gain, tau, R, Vplus = 1e6, 1.0, 1.0, 1.0 +c = h.Matrix(2, 2) +g = h.Matrix(2, 2) +y = h.Vector(2) +y0 = h.Vector(2) +b = h.Vector([0.0, gain * Vplus]) +c.setval(1, 0, tau) +g.setval(0, 0, 1.0 / R) +g.setval(0, 1, 1.0) +g.setval(1, 0, gain + 1.0) +lm = h.LinearMechanism(c, g, y, y0, b) +h.cvode_active(True) +cvode.dae_init_dteps(1e-9, 1) # warn on residual, still check y +cvode.dae_init_mode(3) +h.finitialize(0.0) +# Continuous from rest: v_out held at 0, I=0 (not algebraic steady v~1) +assert abs(y[0]) < 1e-6, f'v_out should be held near 0, got {y[0]}' +assert abs(y[1]) < 1e-6, f'I should be ~0, got {y[1]}' +print('ok') +""" + assert "ok" in _run_isolated(code) + + +def test_extracellular_battery_mode3_holds_vm(): + """Cable + extracellular: mode 3 holds Vm continuous after finitialize.""" + code = r""" +from neuron import h +import math +h.load_file('stdrun.hoc') +cvode = h.CVode() +s = h.Section(name='s') +s.L = s.diam = 10 +s.insert('pas') +s.insert('extracellular') +# modest xg so outer layer is not pure capacitive float +for seg in s: + seg.xg[0] = 1e-3 + seg.xc[0] = 1.0 +ic = h.IClamp(s(0.5)) +ic.delay = 0 +ic.dur = 1e9 +ic.amp = 0.01 +h.cvode_active(True) +cvode.use_daspk(1) +cvode.dae_init_dteps(1e-9, 1) +cvode.dae_init_mode(3) +h.finitialize(-65) +# At t=0 from finitialize, Vm should remain near v_init +assert abs(s(0.5).v - (-65)) < 1e-3, s(0.5).v +h.continuerun(0.5) +# Should integrate without falling back hard-fail +assert math.isfinite(s(0.5).v) +print('ok') +""" + assert "ok" in _run_isolated(code) + + +def test_inductor_battery_holds_current(): + """Inductor: diagonal C[k][k]=L. Hold current; voltage may jump.""" + code = r""" +from neuron import h +import math +h.load_file('stdrun.hoc') +cvode = h.CVode() +# Node voltage y0, inductor current y1; R to ground; inject Is into node. +# LinearCircuit Inductor stamp (j=ground eliminated): +# KCL: -I + v/R = Is => G[0][0]=1/R, G[0][1]=-1, b[0]=Is +# L*I' - v = 0 => C[1][1]=L, G[1][0]=-1 +L, R, Is = 1.0, 2.0, 1.0 +c = h.Matrix(2, 2) +g = h.Matrix(2, 2) +y = h.Vector(2) +y0 = h.Vector(2) +b = h.Vector([Is, 0.0]) +c.setval(1, 1, L) +g.setval(0, 0, 1.0 / R) +g.setval(0, 1, -1.0) +g.setval(1, 0, -1.0) +lm = h.LinearMechanism(c, g, y, y0, b) +h.cvode_active(True) +cvode.dae_init_dteps(1e-9, 1) +cvode.dae_init_mode(3) +h.finitialize(0.0) +# Hold I=0 (rest); KCL => v = Is*R +assert abs(y[1]) < 1e-6, f'I_L should be held ~0, got {y[1]}' +assert math.isclose(y[0], Is * R, rel_tol=1e-5, abs_tol=1e-5), y[0] +print('ok') +""" + assert "ok" in _run_isolated(code) + + +def test_ida_ic_three_panel_audit_isolated(): + """Smoke: three-panel audit on finitialize and after a later reinit.""" + code = r""" +from neuron import h +import tempfile, os +h.load_file('stdrun.hoc') +cvode = h.CVode() +c = h.Matrix(2, 2) +g = h.Matrix(2, 2) +y = h.Vector(2) +y0 = h.Vector(2) +b = h.Vector([1.0, 0.0]) +c.setval(0, 0, 1.0) +c.setval(0, 1, -1.0) +c.setval(1, 0, -1.0) +c.setval(1, 1, 1.0) +g.setval(1, 1, 0.5) +lm = h.LinearMechanism(c, g, y, y0, b) +h.cvode_active(True) +cvode.dae_init_mode(3) +path = tempfile.mktemp(prefix='ida_ic_audit_', suffix='.txt') +cvode.dae_init_audit_file(path) +# finitialize audit +cvode.dae_init_audit(2, 0.0) +h.finitialize(0.0) +# second audit: after short run, force reinit at t via second finitialize-like path +# Arm for t>=0.1 then advance with at_time-style stop via continuerun + re-init +cvode.dae_init_audit(2, 0.0) # next reinit (finitialize again) +h.finitialize(0.0) +cvode.dae_init_audit(0) +cvode.dae_init_audit_file('') +text = open(path).read() +os.remove(path) +assert 'IDA IC three-panel audit' in text, text +assert 'B post-event pre-IC' in text, text +assert 'C post-IC' in text, text +assert 'summary' in text, text +assert 'C*y' in text or "C*y'=f" in text or 'diagonal' in text or 'content' in text or 'mode 3' in text, text +print('ok') +""" + assert "ok" in _run_isolated(code) + + +def test_forcing_tplus_play_ramp_at_reinit(): + """A1/A2: forcing t+ in audit; free y' from u' on series C–R ramp. + + After jump at t=1 onto slope I'=0.5 (I=0.5, R=C=1): + V2' = R*I' = 0.5, V1' = V2' + I/C = 1.0 + """ + code = r""" +from neuron import h +import tempfile, os +h.load_file('stdrun.hoc') +cvode = h.CVode() +# Series floating C + R to ground; I into node 0 via b[0] +c = h.Matrix(2, 2) +g = h.Matrix(2, 2) +y = h.Vector(2) +y0 = h.Vector(2) +b = h.Vector([0.0, 0.0]) +c.setval(0, 0, 1.0) +c.setval(0, 1, -1.0) +c.setval(1, 0, -1.0) +c.setval(1, 1, 1.0) +g.setval(1, 1, 1.0) +lm = h.LinearMechanism(c, g, y, y0, b) +tvec = h.Vector([0, 1, 1, 2, 2, 5]) +ivec = h.Vector([0, 0, 0.5, 1.0, 0, 0]) +ivec.play(b._ref_x[0], tvec, True) +h.cvode_active(True) +cvode.use_daspk(1) +cvode.dae_init_mode(3) +h.finitialize(0.0) +h.continuerun(1.0) +path = tempfile.mktemp(prefix='ida_forcing_tplus_', suffix='.txt') +cvode.dae_init_audit_file(path) +cvode.dae_init_audit(2, 0.0) +cvode.re_init() +text = open(path).read() +os.remove(path) +cvode.dae_init_audit(0) +cvode.dae_init_audit_file('') +assert 'forcing t+ info' in text, text +found_u = False +for line in text.splitlines(): + parts = line.split() + if len(parts) < 4: + continue + try: + u = float(parts[2]) + up = float(parts[3]) + except ValueError: + continue + if abs(u - 0.5) < 1e-6 and abs(up - 0.5) < 1e-6: + found_u = True + break +assert found_u, 'expected u≈0.5 and u′≈0.5 in forcing t+ dump:\n' + text +# A2: panel C y' should be V1'≈1, V2'≈0.5 (eq order = LM y order) +# Parse "C post-IC" table rows: eq y y' residual +in_c = False +yp = {} +for line in text.splitlines(): + if 'C post-IC' in line: + in_c = True + continue + if in_c and line.startswith('---'): + break + if not in_c: + continue + parts = line.split() + if len(parts) >= 4 and parts[0].isdigit(): + eq = int(parts[0]) + yp[eq] = float(parts[2]) +assert 0 in yp and 1 in yp, text +assert abs(yp[0] - 1.0) < 1e-4, (yp, text) +assert abs(yp[1] - 0.5) < 1e-4, (yp, text) +print('ok') +""" + assert "ok" in _run_isolated(code) + + +def test_seclamp_tiny_cm_mode3_no_init_failure(): + """SEClamp + tiny cm: mode 3 C*y'=f clears residual (scm2eem-style). + + Legacy nano-step leaves O(dteps*|y'|) residual and fails WRMS under + default style 0. Mode 3 should report err=0 and WRMS C ~ 0. + """ + code = r""" +from neuron import h +import math +import tempfile, os +h.load_file('stdrun.hoc') +cvode = h.CVode() +s = h.Section(name='eem') +s.L = 10 +s.diam = 100 / s.L / h.PI +s.insert('hh') +s.cm = 0.001 +vc = h.SEClamp(s(0.5)) +vc.rs = 0.1 +vc.dur1 = 1e9 +vc.amp1 = -65.0 +h.cvode_active(True) +cvode.use_daspk(1) +cvode.atol(1e-6) +cvode.dae_init_mode(3) +path = tempfile.mktemp(prefix='ida_seclamp_', suffix='.txt') +cvode.dae_init_audit_file(path) +cvode.dae_init_audit(2, 0.0) +h.finitialize(-65) +cvode.dae_init_audit(0) +cvode.dae_init_audit_file('') +text = open(path).read() +os.remove(path) +assert 'requested_mode=3' in text, text +assert 'path_mode=3' in text, text +assert 'fallback=0' in text, text +assert 'err=0' in text, text +# summary WRMS C is the third number — expect exact 0 after C*y'=f +assert 'WRMS A/B/C' in text, text +# last token on that line should be 0 (or ~0) +for line in text.splitlines(): + if 'WRMS A/B/C' in line: + parts = line.split('=')[-1].split('/') + wrms_c = float(parts[-1].strip()) + assert wrms_c < 1e-6, (wrms_c, line) + break +assert math.isclose(s(0.5).v, -65.0, abs_tol=1e-6) +vc.amp1 = -60.0 +cvode.re_init() +h.continuerun(0.1) +assert math.isfinite(s(0.5).v) +print('ok') +""" + assert "ok" in _run_isolated(code) + + +def test_mode3_forcing_tplus_suite_a3(): + """A3: istep / kink / end extrapolate / finitialize / multi-event (series C–R). + + Algebraic oracles (R=C=1): V2' = I', V1' = I' + I (from V2'=R*I', V1'=V2'+I/C). + Distills external nrntest/nrniv/ida iramp/istep-style cases into CI tests. + """ + code = r""" +from neuron import h +import tempfile, os + +def parse_forcing_and_yp(text): + u = up = None + yp = {} + in_c = False + wrms_c = None + in_forcing = False + for line in text.splitlines(): + if 'forcing t+ info' in line: + in_forcing = True + continue + if in_forcing and line.startswith('---'): + in_forcing = False + if in_forcing: + parts = line.split() + if len(parts) >= 4 and parts[0].isdigit(): + try: + int(parts[1]) + u, up = float(parts[2]), float(parts[3]) + except ValueError: + pass + if 'WRMS' in line and 'A/B/C' in line: + wrms_c = float(line.split('=')[-1].split('/')[-1].strip()) + if 'C post-IC' in line: + in_c = True + continue + if in_c and line.startswith('---'): + in_c = False + continue + if in_c: + parts = line.split() + if len(parts) >= 4 and parts[0].isdigit(): + yp[int(parts[0])] = float(parts[2]) + return u, up, yp, wrms_c + +def check_case(name, u, up, yp, wrms_c, eu, eup, eyp0, eyp1, tol=1e-4): + assert wrms_c is not None and abs(wrms_c) < 1e-9, (name, 'wrms', wrms_c) + assert u is not None and abs(u - eu) < tol, (name, 'u', u, eu) + assert up is not None and abs(up - eup) < tol, (name, 'up', up, eup) + assert abs(yp.get(0, 1e9) - eyp0) < tol, (name, 'yp0', yp, eyp0) + assert abs(yp.get(1, 1e9) - eyp1) < tol, (name, 'yp1', yp, eyp1) + +h.load_file('stdrun.hoc') +cvode = h.CVode() +# Single series C–R LM for the whole suite (avoid stacking mechanisms) +c = h.Matrix(2, 2) +g = h.Matrix(2, 2) +y = h.Vector(2) +y0 = h.Vector(2) +b = h.Vector([0.0, 0.0]) +c.setval(0, 0, 1.0) +c.setval(0, 1, -1.0) +c.setval(1, 0, -1.0) +c.setval(1, 1, 1.0) +g.setval(1, 1, 1.0) +lm = h.LinearMechanism(c, g, y, y0, b) +h.cvode_active(True) +cvode.use_daspk(1) +cvode.dae_init_mode(3) + +play_iv = None +play_tv = None + +def set_play(tvec, ivec): + global play_iv, play_tv + if play_iv is not None: + play_iv.play_remove() + play_tv = h.Vector(tvec) + play_iv = h.Vector(ivec) + play_iv.play(b._ref_x[0], play_tv, True) + +def audit_reinit(): + path = tempfile.mktemp(prefix='ida_a3_', suffix='.txt') + cvode.dae_init_audit_file(path) + cvode.dae_init_audit(2, 0.0) + cvode.re_init() + text = open(path).read() + os.remove(path) + cvode.dae_init_audit(0) + cvode.dae_init_audit_file('') + return text + +# --- T_istep: jump to I=0.5 then flat (I'=0) --- +set_play([0, 1, 1, 5], [0, 0, 0.5, 0.5]) +h.finitialize(0.0) +h.continuerun(1.0) +u, up, yp, wrms = parse_forcing_and_yp(audit_reinit()) +check_case('istep', u, up, yp, wrms, 0.5, 0.0, 0.5, 0.0) + +# --- T_kink: continuous I, slope becomes 1 at t=1 (I=0) --- +set_play([0, 1, 2], [0, 0, 1]) +h.finitialize(0.0) +h.continuerun(1.0) +u, up, yp, wrms = parse_forcing_and_yp(audit_reinit()) +check_case('kink', u, up, yp, wrms, 0.0, 1.0, 1.0, 1.0) + +# --- T_end_extrap: past last knot, last segment slope 2 --- +set_play([0, 1, 2], [0, 1, 3]) +h.finitialize(0.0) +h.continuerun(2.5) +u, up, yp, wrms = parse_forcing_and_yp(audit_reinit()) +# I(2.5)=4, I'=2 → V2'=2, V1'=6 +check_case('end_extrap', u, up, yp, wrms, 4.0, 2.0, 6.0, 2.0, tol=1e-3) + +# --- T_flat_end: last two y equal → I'=0 after end --- +set_play([0, 1, 2], [0, 1, 1]) +h.finitialize(0.0) +h.continuerun(3.0) +u, up, yp, wrms = parse_forcing_and_yp(audit_reinit()) +check_case('flat_end', u, up, yp, wrms, 1.0, 0.0, 1.0, 0.0) + +# --- T_finitialize: first segment slope 0.5 --- +set_play([0, 2], [0, 1]) +path = tempfile.mktemp(prefix='ida_a3_fini_', suffix='.txt') +cvode.dae_init_audit_file(path) +cvode.dae_init_audit(2, 0.0) +h.finitialize(0.0) +text = open(path).read() +os.remove(path) +cvode.dae_init_audit(0) +cvode.dae_init_audit_file('') +u, up, yp, wrms = parse_forcing_and_yp(text) +check_case('finitialize_slope', u, up, yp, wrms, 0.0, 0.5, 0.5, 0.5) + +# --- T_multi_event: iramp-like jump+ramp then off --- +set_play([0, 1, 1, 2, 2, 5], [0, 0, 0.5, 1.0, 0, 0]) +h.finitialize(0.0) +h.continuerun(1.0) +u, up, yp, wrms = parse_forcing_and_yp(audit_reinit()) +check_case('multi_t1_ramp', u, up, yp, wrms, 0.5, 0.5, 1.0, 0.5) +h.continuerun(2.0) +u, up, yp, wrms = parse_forcing_and_yp(audit_reinit()) +check_case('multi_t2_off', u, up, yp, wrms, 0.0, 0.0, 0.0, 0.0) + +if play_iv is not None: + play_iv.play_remove() +print('ok') +""" + assert "ok" in _run_isolated(code) + + +def test_mode3_dforce_sinusoid_a4(): + """A4: LinearMechanism.dforce supplies analytic b' for free y' (sinusoid I).""" + code = r""" +from neuron import h +import math +import tempfile, os +h.load_file('stdrun.hoc') +cvode = h.CVode() +A, w = 1.0, 2.0 * math.pi # I = A sin(w t), I' = A w cos(w t) +c = h.Matrix(2, 2) +g = h.Matrix(2, 2) +y = h.Vector(2) +y0 = h.Vector(2) +b = h.Vector([0.0, 0.0]) +bdot = h.Vector([0.0, 0.0]) +c.setval(0, 0, 1.0) +c.setval(0, 1, -1.0) +c.setval(1, 0, -1.0) +c.setval(1, 1, 1.0) +g.setval(1, 1, 1.0) + +def force(): + b.x[0] = A * math.sin(w * h.t) + +def dforce(): + bdot.x[0] = A * w * math.cos(w * h.t) + +lm = h.LinearMechanism(force, c, g, y, y0, b) +lm.dforce(dforce, bdot) +h.cvode_active(True) +cvode.use_daspk(1) +cvode.dae_init_mode(3) +# IC at t=0: I=0, I'=A*w → V2'=A*w, V1'=A*w + I = A*w +path = tempfile.mktemp(prefix='ida_dforce_', suffix='.txt') +cvode.dae_init_audit_file(path) +cvode.dae_init_audit(2, 0.0) +h.finitialize(0.0) +text = open(path).read() +os.remove(path) +assert 'dforce' in text or 'bdot' in text, text +assert 'err=0' in text, text +in_c = False +yp = {} +for line in text.splitlines(): + if 'C post-IC' in line: + in_c = True + continue + if in_c and line.startswith('---'): + break + if in_c: + parts = line.split() + if len(parts) >= 4 and parts[0].isdigit(): + yp[int(parts[0])] = float(parts[2]) +expect = A * w +assert abs(yp.get(1, 1e9) - expect) < 1e-4, (yp, expect, text) +assert abs(yp.get(0, 1e9) - expect) < 1e-4, (yp, expect, text) # I=0 so V1'=V2' +# mid-ramp time: re_init at t=0.25 without play +h.t = 0.25 +force() +dforce() +cvode.dae_init_audit_file(path) +cvode.dae_init_audit(2, 0.0) +cvode.re_init() +text = open(path).read() +os.remove(path) +I = A * math.sin(w * 0.25) +Ip = A * w * math.cos(w * 0.25) +yp = {} +in_c = False +for line in text.splitlines(): + if 'C post-IC' in line: + in_c = True + continue + if in_c and line.startswith('---'): + break + if in_c: + parts = line.split() + if len(parts) >= 4 and parts[0].isdigit(): + yp[int(parts[0])] = float(parts[2]) +assert abs(yp.get(1, 1e9) - Ip) < 1e-3, (yp, Ip, text) +assert abs(yp.get(0, 1e9) - (Ip + I)) < 1e-3, (yp, Ip, I, text) +print('ok') +""" + assert "ok" in _run_isolated(code) + + +def test_mode3_dforce_fd_fallback_a4(): + """A4: without dforce, FD of f_callable estimates b' for free y'.""" + code = r""" +from neuron import h +import math +import tempfile, os +h.load_file('stdrun.hoc') +cvode = h.CVode() +A, w = 1.0, 2.0 # slower so FD with h=1e-8 is accurate +c = h.Matrix(2, 2) +g = h.Matrix(2, 2) +y = h.Vector(2) +y0 = h.Vector(2) +b = h.Vector([0.0, 0.0]) +c.setval(0, 0, 1.0) +c.setval(0, 1, -1.0) +c.setval(1, 0, -1.0) +c.setval(1, 1, 1.0) +g.setval(1, 1, 1.0) + +def force(): + b.x[0] = A * math.sin(w * h.t) + +lm = h.LinearMechanism(force, c, g, y, y0, b) +# no lm.dforce — FD fallback +h.cvode_active(True) +cvode.use_daspk(1) +cvode.dae_init_mode(3) +path = tempfile.mktemp(prefix='ida_fd_', suffix='.txt') +cvode.dae_init_audit_file(path) +cvode.dae_init_audit(2, 0.0) +h.finitialize(0.0) +text = open(path).read() +os.remove(path) +assert 'bdot_fd' in text or 'err=0' in text, text +yp = {} +in_c = False +for line in text.splitlines(): + if 'C post-IC' in line: + in_c = True + continue + if in_c and line.startswith('---'): + break + if in_c: + parts = line.split() + if len(parts) >= 4 and parts[0].isdigit(): + yp[int(parts[0])] = float(parts[2]) +expect = A * w # cos(0)=1 +assert abs(yp.get(1, 1e9) - expect) < 1e-2, (yp, expect, text) +assert abs(yp.get(0, 1e9) - expect) < 1e-2, (yp, expect, text) +print('ok') +""" + assert "ok" in _run_isolated(code) + + +if __name__ == "__main__": + test_dae_init_mode_api() + test_dae_init_audit_api() + test_dae_init_stats_api() + test_pure_resistive_all_modes_isolated() + test_series_cr_mode0_mode1_isolated() + test_series_cr_battery_mode3_isolated() + test_opamp_tau_battery_holds_output_voltage() + test_inductor_battery_holds_current() + test_extracellular_battery_mode3_holds_vm() + test_ida_ic_three_panel_audit_isolated() + test_forcing_tplus_play_ramp_at_reinit() + test_seclamp_tiny_cm_mode3_no_init_failure() + test_mode3_forcing_tplus_suite_a3() + test_mode3_dforce_sinusoid_a4() + test_mode3_dforce_fd_fallback_a4() + print("ok") diff --git a/test/hoctests/tests/test_ida_source_current.py b/test/hoctests/tests/test_ida_source_current.py new file mode 100644 index 0000000000..d782bbebba --- /dev/null +++ b/test/hoctests/tests/test_ida_source_current.py @@ -0,0 +1,956 @@ +"""Plan (b): source-current discontinuities — PWL MOD + Section↔LM parity. + +WP0b–WP3: PWLClamp electrode stimulus; mode-3 IClamp/PWL jump/kink; +Section↔LM parity (E0/E1); end-ri LM twins (Z0/Z1). + +Run via ctest (preferred for CI — loads PWLClamp from hoctests nrnivmodl hash):: + + ctest -R hoctests::test_ida_source_current --output-on-failure + +Or standalone (auto nrnivmodl of test/hoctests/pwlclamp.mod if needed):: + + export PATH=build/bin:$PATH LD_LIBRARY_PATH=build/lib:$LD_LIBRARY_PATH + export PYTHONPATH=build/lib/python:$PYTHONPATH + python test/hoctests/tests/test_ida_source_current.py + +Each model case runs in a subprocess (LinearMechanism teardown hygiene). +""" + +from __future__ import annotations + +import os +import subprocess +import sys +import tempfile + +# Canonical waveform tables (same as test_ida_init_mode A1/A3 play) +W_ISTEP = ([0.0, 1.0, 1.0, 5.0], [0.0, 0.0, 0.5, 0.5]) +W_KINK = ([0.0, 1.0, 2.0], [0.0, 0.0, 1.0]) +W_JUMP_RAMP = ([0.0, 1.0, 1.0, 2.0, 2.0, 5.0], [0.0, 0.0, 0.5, 1.0, 0.0, 0.0]) + + +def _repo_hoctests(): + return os.path.abspath(os.path.join(os.path.dirname(__file__), "..")) + + +def _find_nrnivmodl() -> str: + env = os.environ.get("NRNIVMODL") + if env and os.path.isfile(env): + return env + # Prefer build tree next to a PYTHONPATH neuron install + for d in os.environ.get("PATH", "").split(os.pathsep): + cand = os.path.join(d, "nrnivmodl") + if os.path.isfile(cand) and os.access(cand, os.X_OK): + # Skip broken pyenv shims that point at missing .data/bin + try: + r = subprocess.run( + [cand, "-help"], + capture_output=True, + text=True, + timeout=10, + ) + # nrnivmodl -help may exit nonzero; existence of binary is enough if + # it is the cmake-installed script under build/bin + if "build/bin" in cand or r.returncode in (0, 1, 2): + if os.path.isfile(cand): + return cand + except Exception: + continue + return "nrnivmodl" + + +_PWL_MECH_DIR = None + + +def _build_pwlclamp_mech_dir() -> str: + """Once-per-process compile of pwlclamp.mod; return dir with libnrnmech.so.""" + global _PWL_MECH_DIR + if _PWL_MECH_DIR and os.path.isfile(os.path.join(_PWL_MECH_DIR, "libnrnmech.so")): + return _PWL_MECH_DIR + # Reuse prior local build if present + cached = "/tmp/pwlclamp_build/x86_64" + if os.path.isfile(os.path.join(cached, "libnrnmech.so")): + _PWL_MECH_DIR = cached + return cached + mod = os.path.join(_repo_hoctests(), "pwlclamp.mod") + build = tempfile.mkdtemp(prefix="pwlclamp_") + nrnivmodl = _find_nrnivmodl() + r = subprocess.run( + [nrnivmodl, mod], + cwd=build, + capture_output=True, + text=True, + env=os.environ.copy(), + ) + if r.returncode != 0: + raise RuntimeError(f"nrnivmodl failed ({nrnivmodl}):\n{r.stdout}\n{r.stderr}") + arch = None + for name in os.listdir(build): + p = os.path.join(build, name) + if os.path.isdir(p) and os.path.isfile(os.path.join(p, "libnrnmech.so")): + arch = p + break + if not arch: + raise RuntimeError( + f"libnrnmech.so not found under {build}: {os.listdir(build)}" + ) + _PWL_MECH_DIR = arch + return arch + + +def _sanitizer_child_env(env=None): + """Re-apply sanitizer preload for macOS SIP (see NeuronTestHelper.cmake).""" + env = os.environ.copy() if env is None else env + try: + env[os.environ["NRN_SANITIZER_PRELOAD_VAR"]] = os.environ[ + "NRN_SANITIZER_PRELOAD_VAL" + ] + except KeyError: + pass + return env + + +def _run_isolated(code: str, timeout: float = 120.0) -> str: + """Run code in a subprocess with the same interpreter/env; return stdout.""" + env = _sanitizer_child_env() + # Point children at a ready PWLClamp mechanism directory + try: + env["NRN_PWLCLAMP_MECH"] = _build_pwlclamp_mech_dir() + except Exception as e: + # IClamp-only tests do not need the mech + env.pop("NRN_PWLCLAMP_MECH", None) + env["NRN_PWLCLAMP_BUILD_ERR"] = str(e) + exe = os.environ.get("NRN_PYTHON_EXECUTABLE", sys.executable) + proc = subprocess.run( + [exe, "-c", code], + capture_output=True, + text=True, + timeout=timeout, + env=env, + ) + out = (proc.stdout or "") + (proc.stderr or "") + if proc.returncode != 0: + raise AssertionError( + f"subprocess failed rc={proc.returncode}\n--- stdout+stderr ---\n{out}" + ) + return out + + +_LOADER = r""" +import os +from neuron import h + +def ensure_pwlclamp(): + if hasattr(h, 'PWLClamp'): + return + candidates = [] + mech = os.environ.get('NRN_PWLCLAMP_MECH') + if mech: + candidates.append(mech) + for d in os.environ.get('PATH', '').split(os.pathsep): + if d and os.path.isdir(d): + candidates.append(d) + for d in candidates: + so = os.path.join(d, 'libnrnmech.so') + if os.path.isfile(so): + h.nrn_load_dll(so) + if hasattr(h, 'PWLClamp'): + return + err = os.environ.get('NRN_PWLCLAMP_BUILD_ERR', '') + raise RuntimeError( + 'PWLClamp not available; set NRN_PWLCLAMP_MECH to dir with libnrnmech.so. ' + + err + ) + +ensure_pwlclamp() +""" + + +_SET_PWL = r""" +def set_pwl(stim, tvec, ivec): + assert len(tvec) == len(ivec) and len(tvec) <= 16 + stim.nkt = len(tvec) + for j, (tt, ii) in enumerate(zip(tvec, ivec)): + stim.set_knot(j, tt, ii) +""" + + +def test_pwlclamp_ival_right_continuous(): + """WP0b: PWLClamp.ival matches jump/kink tables at knot t+.""" + code = ( + _LOADER + + _SET_PWL + + r""" +from neuron import h +s = h.Section(name='s') +stim = h.PWLClamp(s(0.5)) +# istep at t=1+: I=0.5 +set_pwl(stim, [0, 1, 1, 5], [0, 0, 0.5, 0.5]) +assert abs(stim.ival(1.0) - 0.5) < 1e-12, stim.ival(1.0) +assert abs(stim.ival(0.5) - 0.0) < 1e-12, stim.ival(0.5) +assert abs(stim.ival(3.0) - 0.5) < 1e-12, stim.ival(3.0) +# kink at t=1: I=0, then ramp +set_pwl(stim, [0, 1, 2], [0, 0, 1]) +assert abs(stim.ival(1.0) - 0.0) < 1e-12, stim.ival(1.0) +assert abs(stim.ival(1.5) - 0.5) < 1e-12, stim.ival(1.5) +# jump-ramp A1 at t=1+: I=0.5 +set_pwl(stim, [0, 1, 1, 2, 2, 5], [0, 0, 0.5, 1.0, 0, 0]) +assert abs(stim.ival(1.0) - 0.5) < 1e-12, stim.ival(1.0) +print('ok') +""" + ) + assert "ok" in _run_isolated(code) + + +def test_iclamp_step_mode3_section(): + """WP1: pure capacitive Section + IClamp step; mode 3 path, Vm hold.""" + code = r""" +from neuron import h +import math +h.load_file('stdrun.hoc') +cvode = h.CVode() +s = h.Section(name='s') +s.L = s.diam = 10 +s.insert('pas') +s.g_pas = 1e-4 +s.e_pas = -65 +ic = h.IClamp(s(0.5)) +ic.delay = 1.0 +ic.dur = 1e9 +ic.amp = 0.1 +h.cvode_active(True) +cvode.use_daspk(1) +cvode.dae_init_mode(3) +cvode.dae_init_stats(1) +h.finitialize(-65) +h.continuerun(0.999) +v_before = s(0.5).v +h.continuerun(1.0) +# at/after step: path_mode 3, residual ok via stats +v = h.Vector() +n = cvode.dae_init_stats(v) +assert n == 8 and v.size() == 8, list(v) +# After continuerun past del, at least one reinit should have used mode 3 +# v[6] last path mode, v[1] mode3 ok count (see A5 layout) +assert int(v[6]) == 3 or v[1] >= 1, list(v) +# Vm continuous across the step (charge hold); allow tiny float drift +v_after = s(0.5).v +assert abs(v_after - v_before) < 1e-3, (v_before, v_after) +assert math.isfinite(v_after) +# integrate a bit further +h.continuerun(1.5) +assert math.isfinite(s(0.5).v) +print('ok') +""" + assert "ok" in _run_isolated(code) + + +def test_pwl_istep_mode3_section(): + """WP1: W_istep via PWLClamp on capacitive Section; mode 3.""" + code = ( + _LOADER + + _SET_PWL + + r""" +from neuron import h +import math +h.load_file('stdrun.hoc') +cvode = h.CVode() +s = h.Section(name='s') +s.L = s.diam = 10 +s.insert('pas') +s.g_pas = 1e-4 +s.e_pas = -65 +stim = h.PWLClamp(s(0.5)) +set_pwl(stim, [0, 1, 1, 5], [0, 0, 0.5, 0.5]) +h.cvode_active(True) +cvode.use_daspk(1) +cvode.dae_init_mode(3) +cvode.dae_init_stats(1) +h.finitialize(-65) +h.continuerun(0.999) +v_before = s(0.5).v +h.continuerun(1.0) +v = h.Vector() +cvode.dae_init_stats(v) +assert int(v[6]) == 3 or v[1] >= 1, list(v) +assert abs(s(0.5).v - v_before) < 1e-3, (v_before, s(0.5).v) +assert abs(stim.ival(1.0) - 0.5) < 1e-12 +h.continuerun(1.2) +assert math.isfinite(s(0.5).v) +print('ok') +""" + ) + assert "ok" in _run_isolated(code) + + +def test_pwl_kink_mode3_section(): + """WP1b: W_kink via PWLClamp; continuous I, slope change; mode 3.""" + code = ( + _LOADER + + _SET_PWL + + r""" +from neuron import h +import math +h.load_file('stdrun.hoc') +cvode = h.CVode() +s = h.Section(name='s') +s.L = s.diam = 10 +s.insert('pas') +s.g_pas = 1e-4 +s.e_pas = -65 +stim = h.PWLClamp(s(0.5)) +set_pwl(stim, [0, 1, 2], [0, 0, 1]) +h.cvode_active(True) +cvode.use_daspk(1) +cvode.dae_init_mode(3) +cvode.dae_init_stats(1) +h.finitialize(-65) +h.continuerun(0.999) +v_before = s(0.5).v +h.continuerun(1.0) +v = h.Vector() +cvode.dae_init_stats(v) +assert int(v[6]) == 3 or v[1] >= 1, list(v) +assert abs(s(0.5).v - v_before) < 1e-3, (v_before, s(0.5).v) +assert abs(stim.ival(1.0) - 0.0) < 1e-12 +h.continuerun(1.5) +assert math.isfinite(s(0.5).v) +print('ok') +""" + ) + assert "ok" in _run_isolated(code) + + +def test_e0_section_lm_parity_istep_kink(): + """WP1b / E0: Section+PWL vs LM+play for W_istep and W_kink (capacitive). + + LM: single C to ground + leak R + I into the node (play). + Section: pas + PWLClamp mid-compartment. + Compare relative Vm change across the event and mode-3 success. + Absolute IR drop scaling differs (area units); compare continuity and path. + """ + code = ( + _LOADER + + _SET_PWL + + r""" +from neuron import h +import math +h.load_file('stdrun.hoc') + +def run_section(tvec, ivec, t_event): + cvode = h.CVode() + s = h.Section(name='s') + s.L = s.diam = 10 + s.insert('pas') + s.g_pas = 0.001 + s.e_pas = 0.0 + s.cm = 1.0 + stim = h.PWLClamp(s(0.5)) + set_pwl(stim, tvec, ivec) + h.cvode_active(True) + cvode.use_daspk(1) + cvode.dae_init_mode(3) + cvode.dae_init_stats(1) + h.finitialize(0.0) + h.continuerun(t_event - 1e-4) + v0 = s(0.5).v + h.continuerun(t_event) + v1 = s(0.5).v + stats = h.Vector() + cvode.dae_init_stats(stats) + return v0, v1, list(stats), stim.ival(t_event) + +def run_lm(tvec, ivec, t_event): + cvode = h.CVode() + # C=1, R=1e3 (weak leak), I into node 0; single voltage unknown + # Use 1-state: C v' + v/R = I + c = h.Matrix(1, 1) + g = h.Matrix(1, 1) + y = h.Vector(1) + y0 = h.Vector(1) + b = h.Vector([0.0]) + c.setval(0, 0, 1.0) + g.setval(0, 0, 1.0 / 1e3) + lm = h.LinearMechanism(c, g, y, y0, b) + tv = h.Vector(tvec) + iv = h.Vector(ivec) + iv.play(b._ref_x[0], tv, True) + h.cvode_active(True) + cvode.use_daspk(1) + cvode.dae_init_mode(3) + cvode.dae_init_stats(1) + h.finitialize(0.0) + h.continuerun(t_event - 1e-4) + v0 = y[0] + h.continuerun(t_event) + v1 = y[0] + stats = h.Vector() + cvode.dae_init_stats(stats) + return v0, v1, list(stats), b[0] + +def check_case(name, tvec, ivec, t_event): + sv0, sv1, sstats, si = run_section(tvec, ivec, t_event) + lv0, lv1, lstats, li = run_lm(tvec, ivec, t_event) + # mode 3 used on both + assert int(sstats[6]) == 3 or sstats[1] >= 1, (name, 'sec stats', sstats) + assert int(lstats[6]) == 3 or lstats[1] >= 1, (name, 'lm stats', lstats) + # continuous content: voltage continuous across event + assert abs(sv1 - sv0) < 1e-3, (name, 'sec dV', sv0, sv1) + assert abs(lv1 - lv0) < 1e-3, (name, 'lm dV', lv0, lv1) + # drive value at t+ agrees + assert abs(si - li) < 1e-9, (name, 'I', si, li) + print(name, 'ok', 'I=', si, 'sec dV=', sv1-sv0, 'lm dV=', lv1-lv0) + +check_case('istep', [0, 1, 1, 5], [0, 0, 0.5, 0.5], 1.0) +check_case('kink', [0, 1, 2], [0, 0, 1], 1.0) +print('ok') +""" + ) + assert "ok" in _run_isolated(code) + + +def test_pwl_jump_ramp_section_mode3(): + """WP1b: W_jump_ramp on Section via PWLClamp (A1-style table).""" + code = ( + _LOADER + + _SET_PWL + + r""" +from neuron import h +import math +h.load_file('stdrun.hoc') +cvode = h.CVode() +s = h.Section(name='s') +s.L = s.diam = 10 +s.insert('pas') +s.g_pas = 1e-4 +s.e_pas = -65 +stim = h.PWLClamp(s(0.5)) +set_pwl(stim, [0, 1, 1, 2, 2, 5], [0, 0, 0.5, 1.0, 0, 0]) +h.cvode_active(True) +cvode.use_daspk(1) +cvode.dae_init_mode(3) +cvode.dae_init_stats(1) +h.finitialize(-65) +h.continuerun(0.999) +v_before = s(0.5).v +h.continuerun(1.0) +assert abs(stim.ival(1.0) - 0.5) < 1e-12 +v = h.Vector() +cvode.dae_init_stats(v) +assert int(v[6]) == 3 or v[1] >= 1, list(v) +assert abs(s(0.5).v - v_before) < 1e-3, (v_before, s(0.5).v) +h.continuerun(1.5) +assert math.isfinite(s(0.5).v) +print('ok') +""" + ) + assert "ok" in _run_isolated(code) + + +def test_e1_extracellular_pwl_istep_mode3(): + """WP3 E1: 1-layer extracellular + PWL istep; mode 3 holds Vm, path_mode=3. + + Uses nlayer_extracellular(1). Sets xg/xc on all segments (including ends). + Requires coupled cm+xc seed (electrode current in F[vi] charges both). + """ + code = ( + _LOADER + + _SET_PWL + + r""" +from neuron import h +import math +h.load_file('stdrun.hoc') +h.nlayer_extracellular(1) +cvode = h.CVode() +s = h.Section(name='s') +s.L = s.diam = 10 +s.nseg = 1 +s.insert('pas') +s.g_pas = 1e-5 +s.e_pas = 0.0 +s.insert('extracellular') +for seg in s: + seg.xg[0] = 1e-3 + seg.xc[0] = 1.0 +stim = h.PWLClamp(s(0.5)) +set_pwl(stim, [0, 1, 1, 5], [0, 0, 0.1, 0.1]) +h.cvode_active(True) +cvode.use_daspk(1) +cvode.dae_init_mode(3) +cvode.dae_init_stats(1) +h.finitialize(0.0) +h.continuerun(0.999) +vm0 = s(0.5).v # transmembrane Vm with extracellular +h.continuerun(1.0) +vm1 = s(0.5).v # transmembrane Vm with extracellular +st = h.Vector() +cvode.dae_init_stats(st) +assert int(st[6]) == 3, list(st) +assert st[2] == 0, ('fallback', list(st)) +assert abs(vm1 - vm0) < 1e-6, (vm0, vm1) +assert math.isfinite(s(0.5).v) +print('ok') +""" + ) + assert "ok" in _run_isolated(code) + + +def test_e1_section_lm_parity_istep_kink(): + """WP3 E1 parity: Section+xtral+PWL vs LM floating Cm + xg + play force. + + LM twin (2-node): y0=v_int, y1=vext; floating C between them; xg on vext; + I into y0 via play. Compare Vm continuity and mode-3 path on both. + """ + code = ( + _LOADER + + _SET_PWL + + r""" +from neuron import h +import math +h.load_file('stdrun.hoc') + +def run_section(tvec, ivec, t_event): + h.nlayer_extracellular(1) + cvode = h.CVode() + s = h.Section(name='s') + s.L = s.diam = 10 + s.nseg = 1 + s.insert('pas') + s.g_pas = 1e-5 + s.e_pas = 0.0 + s.insert('extracellular') + for seg in s: + seg.xg[0] = 1e-3 + seg.xc[0] = 1.0 + stim = h.PWLClamp(s(0.5)) + set_pwl(stim, tvec, ivec) + h.cvode_active(True) + cvode.use_daspk(1) + cvode.dae_init_mode(3) + cvode.dae_init_stats(1) + h.finitialize(0.0) + h.continuerun(t_event - 1e-4) + vm0 = s(0.5).v # transmembrane Vm with extracellular + h.continuerun(t_event) + vm1 = s(0.5).v # transmembrane Vm with extracellular + st = h.Vector() + cvode.dae_init_stats(st) + return vm0, vm1, list(st), stim.ival(t_event) + +def run_lm(tvec, ivec, t_event): + cvode = h.CVode() + Cm, xg = 1.0, 1e-3 + c = h.Matrix(2, 2) + g = h.Matrix(2, 2) + y = h.Vector(2) + y0 = h.Vector(2) + b = h.Vector([0.0, 0.0]) + c.setval(0, 0, Cm); c.setval(0, 1, -Cm) + c.setval(1, 0, -Cm); c.setval(1, 1, Cm) + g.setval(1, 1, xg) + lm = h.LinearMechanism(c, g, y, y0, b) + tv = h.Vector(tvec) + iv = h.Vector(ivec) + iv.play(b._ref_x[0], tv, True) + h.cvode_active(True) + cvode.use_daspk(1) + cvode.dae_init_mode(3) + cvode.dae_init_stats(1) + h.finitialize(0.0) + h.continuerun(t_event - 1e-4) + vm0 = y[0] - y[1] + h.continuerun(t_event) + vm1 = y[0] - y[1] + st = h.Vector() + cvode.dae_init_stats(st) + return vm0, vm1, list(st), b[0] + +def check(name, tvec, ivec, te): + svm0, svm1, sst, si = run_section(tvec, ivec, te) + lvm0, lvm1, lst, li = run_lm(tvec, ivec, te) + assert int(sst[6]) == 3 and sst[2] == 0, (name, 'sec', sst) + assert int(lst[6]) == 3 and lst[2] == 0, (name, 'lm', lst) + assert abs(svm1 - svm0) < 1e-6, (name, 'sec dVm', svm0, svm1) + assert abs(lvm1 - lvm0) < 1e-6, (name, 'lm dVm', lvm0, lvm1) + assert abs(si - li) < 1e-9, (name, 'I', si, li) + print(name, 'ok') + +check('istep', [0, 1, 1, 5], [0, 0, 0.1, 0.1], 1.0) +check('kink', [0, 1, 2], [0, 0, 1], 1.0) +print('ok') +""" + ) + assert "ok" in _run_isolated(code) + + +def test_end_electrode_mode3_path(): + """A: IClamp at loc 0 and 1; mode 3 path_mode==3 (cable free-y).""" + code = r""" +from neuron import h +h.load_file('stdrun.hoc') +for loc in (0.0, 1.0): + cv = h.CVode() + s = h.Section() + s.L = 100 + s.diam = 10 + s.nseg = 5 + s.Ra = 100 + s.insert('pas') + s.g_pas = 1e-5 + s.e_pas = 0.0 + ic = h.IClamp(s(loc)) + ic.delay = 1.0 + ic.dur = 1e9 + ic.amp = 0.05 + h.cvode_active(True) + cv.use_daspk(1) + cv.dae_init_mode(3) + cv.dae_init_stats(1) + h.finitialize(0.0) + h.continuerun(0.999) + vmid0 = s(0.5).v + h.continuerun(1.0) + st = h.Vector() + cv.dae_init_stats(st) + assert int(st[6]) == 3 and st[2] == 0, (loc, list(st)) + assert abs(s(0.5).v - vmid0) < 1e-6, (loc, vmid0, s(0.5).v) +print('ok') +""" + assert "ok" in _run_isolated(code) + + +def test_z0_end_ri_lm_parity_istep(): + """WP3 Z0: end loc 0; LM twin uses sec(0.0001).ri(); Section path_mode=3.""" + code = ( + _LOADER + + _SET_PWL + + r""" +from neuron import h +import math +h.load_file('stdrun.hoc') + +s = h.Section(name='s') +s.L = 100 +s.diam = 10 +s.nseg = 5 +s.Ra = 100 +s.insert('pas') +s.g_pas = 1e-5 +s.e_pas = 0.0 +h.finitialize(0.0) +R_end = s(0.0001).ri() +assert R_end > 0, R_end +Iamp = 0.05 +# Drop cable before LM so sparse13 is pure linear mechanism +for sec in list(h.allsec()): + h.delete_section(sec=sec) + +# --- LM twin --- +cv2 = h.CVode() +c = h.Matrix(2, 2) +g = h.Matrix(2, 2) +y = h.Vector(2) +y0 = h.Vector(2) +b = h.Vector([0.0, 0.0]) +c.setval(1, 1, 1.0) +g.setval(0, 0, 1.0 / R_end) +g.setval(0, 1, -1.0 / R_end) +g.setval(1, 0, -1.0 / R_end) +g.setval(1, 1, 1.0 / R_end) +lm = h.LinearMechanism(c, g, y, y0, b) +tv = h.Vector([0, 1, 1, 5]) +iv = h.Vector([0, 0, Iamp, Iamp]) +iv.play(b._ref_x[0], tv, True) +h.cvode_active(True) +cv2.use_daspk(1) +cv2.dae_init_mode(3) +cv2.dae_init_stats(1) +h.finitialize(0.0) +h.continuerun(1.0) +cv2.re_init() +st2 = h.Vector() +cv2.dae_init_stats(st2) +assert int(st2[6]) == 3 and st2[2] == 0, list(st2) +assert abs(y[1]) < 1e-6, y[1] +assert math.isclose(y[0], Iamp * R_end, rel_tol=1e-4, abs_tol=1e-4), (y[0], Iamp * R_end) +print('ok', 'R_end', R_end, 'I*R', Iamp * R_end) +""" + ) + assert "ok" in _run_isolated(code) + + # Section end path in its own process (no LM) + code_sec = ( + _LOADER + + _SET_PWL + + r""" +from neuron import h +h.load_file('stdrun.hoc') +s = h.Section(name='s') +s.L = 100 +s.diam = 10 +s.nseg = 5 +s.Ra = 100 +s.insert('pas') +s.g_pas = 1e-5 +s.e_pas = 0.0 +stim = h.PWLClamp(s(0.0)) +set_pwl(stim, [0, 1, 1, 5], [0, 0, 0.05, 0.05]) +cvode = h.CVode() +h.cvode_active(True) +cvode.use_daspk(1) +cvode.dae_init_mode(3) +cvode.dae_init_stats(1) +h.finitialize(0.0) +seg = s(0.1) +h.continuerun(0.999) +vm0 = seg.v +h.continuerun(1.0) +vm1 = seg.v +st = h.Vector() +cvode.dae_init_stats(st) +assert int(st[6]) == 3 and st[2] == 0, list(st) +assert abs(vm1 - vm0) < 1e-3, (vm0, vm1) +print('ok') +""" + ) + assert "ok" in _run_isolated(code_sec) + + +def test_z1_end_ri_lm_parity_istep(): + """WP3 Z1: loc 1; LM series sec(1.0).ri(); Section path_mode=3.""" + code = ( + _LOADER + + _SET_PWL + + r""" +from neuron import h +import math +h.load_file('stdrun.hoc') +s = h.Section(name='s') +s.L = 100 +s.diam = 10 +s.nseg = 5 +s.Ra = 100 +s.insert('pas') +s.g_pas = 1e-5 +s.e_pas = 0.0 +h.finitialize(0.0) +R_end = s(1.0).ri() +Iamp = 0.05 +for sec in list(h.allsec()): + h.delete_section(sec=sec) + +cv2 = h.CVode() +c = h.Matrix(2, 2) +g = h.Matrix(2, 2) +y = h.Vector(2) +y0 = h.Vector(2) +b = h.Vector([0.0, 0.0]) +c.setval(1, 1, 1.0) +g.setval(0, 0, 1.0 / R_end) +g.setval(0, 1, -1.0 / R_end) +g.setval(1, 0, -1.0 / R_end) +g.setval(1, 1, 1.0 / R_end) +lm = h.LinearMechanism(c, g, y, y0, b) +tv = h.Vector([0, 1, 1, 5]) +iv = h.Vector([0, 0, Iamp, Iamp]) +iv.play(b._ref_x[0], tv, True) +h.cvode_active(True) +cv2.use_daspk(1) +cv2.dae_init_mode(3) +cv2.dae_init_stats(1) +h.finitialize(0.0) +h.continuerun(1.0) +cv2.re_init() +st2 = h.Vector() +cv2.dae_init_stats(st2) +assert int(st2[6]) == 3 and st2[2] == 0, list(st2) +assert abs(y[1]) < 1e-6, y[1] +assert math.isclose(y[0], Iamp * R_end, rel_tol=1e-4, abs_tol=1e-4), (y[0], Iamp * R_end) +print('ok') +""" + ) + assert "ok" in _run_isolated(code) + + code_sec = ( + _LOADER + + _SET_PWL + + r""" +from neuron import h +h.load_file('stdrun.hoc') +s = h.Section(name='s') +s.L = 100 +s.diam = 10 +s.nseg = 5 +s.Ra = 100 +s.insert('pas') +s.g_pas = 1e-5 +s.e_pas = 0.0 +stim = h.PWLClamp(s(1.0)) +set_pwl(stim, [0, 1, 1, 5], [0, 0, 0.05, 0.05]) +cvode = h.CVode() +h.cvode_active(True) +cvode.use_daspk(1) +cvode.dae_init_mode(3) +cvode.dae_init_stats(1) +h.finitialize(0.0) +seg = s(0.9) +h.continuerun(0.999) +vm0 = seg.v +h.continuerun(1.0) +vm1 = seg.v +st = h.Vector() +cvode.dae_init_stats(st) +assert int(st[6]) == 3 and st[2] == 0, list(st) +assert abs(vm1 - vm0) < 1e-3, (vm0, vm1) +print('ok') +""" + ) + assert "ok" in _run_isolated(code_sec) + + +def test_e2_xc0_extracellular_istep(): + """B E2: xc=0 (resistive xtral) + IClamp step; mode 3 holds Vm (=seg.v). + + With extracellular, seg.v is transmembrane Vm; vext may jump when xc=0. + Do not use v-vext as Vm. + """ + code = r""" +from neuron import h +h.load_file('stdrun.hoc') +h.nlayer_extracellular(1) +cv = h.CVode() +s = h.Section() +s.L = s.diam = 10 +s.nseg = 1 +s.insert('pas') +s.g_pas = 1e-5 +s.e_pas = 0.0 +s.insert('extracellular') +for seg in s: + seg.xg[0] = 1e-3 + seg.xc[0] = 0.0 +ic = h.IClamp(s(0.5)) +ic.delay = 1.0 +ic.dur = 1e9 +ic.amp = 0.1 +h.cvode_active(True) +cv.use_daspk(1) +cv.dae_init_mode(3) +cv.dae_init_stats(1) +h.finitialize(0.0) +h.continuerun(0.999) +vm0 = s(0.5).v # transmembrane +vext0 = s(0.5).vext[0] +h.continuerun(1.0) +vm1 = s(0.5).v +vext1 = s(0.5).vext[0] +st = h.Vector() +cv.dae_init_stats(st) +assert int(st[6]) == 3 and st[2] == 0, list(st) +assert abs(vm1 - vm0) < 1e-6, (vm0, vm1, vext0, vext1) +# algebraic outer layer: vext free to jump under electrode step +assert abs(vext1 - vext0) > 1e-6 or abs(ic.amp) < 1e-12 +print('ok', 'dVm', vm1 - vm0, 'dvext', vext1 - vext0) +""" + assert "ok" in _run_isolated(code) + + +def test_e3_multilayer_xc_istep(): + """B E3: default nlayer=2 with xc>0 on both layers; mode 3 istep.""" + code = r""" +from neuron import h +h.load_file('stdrun.hoc') +# default nlayer is 2 +cv = h.CVode() +s = h.Section() +s.L = s.diam = 10 +s.nseg = 1 +s.insert('pas') +s.g_pas = 1e-5 +s.e_pas = 0.0 +s.insert('extracellular') +nl = int(h.nlayer_extracellular()) +assert nl >= 2, nl +for seg in s: + for j in range(nl): + seg.xg[j] = 1e-3 + seg.xc[j] = 1.0 +ic = h.IClamp(s(0.5)) +ic.delay = 1.0 +ic.dur = 1e9 +ic.amp = 0.1 +h.cvode_active(True) +cv.use_daspk(1) +cv.dae_init_mode(3) +cv.dae_init_stats(1) +h.finitialize(0.0) +h.continuerun(0.999) +vm0 = s(0.5).v # transmembrane Vm with extracellular +h.continuerun(1.0) +vm1 = s(0.5).v # transmembrane Vm with extracellular +st = h.Vector() +cv.dae_init_stats(st) +assert int(st[6]) == 3 and st[2] == 0, list(st) +assert abs(vm1 - vm0) < 1e-6, (vm0, vm1) +print('ok', 'nlayer', nl) +""" + assert "ok" in _run_isolated(code) + + +def test_e4_xg_scale_istep(): + """B E4: weak and strong xg; mode 3 istep holds Vm.""" + code = r""" +from neuron import h +h.load_file('stdrun.hoc') +h.nlayer_extracellular(1) +for xg in (1e-6, 1e-3, 1.0): + cv = h.CVode() + s = h.Section() + s.L = s.diam = 10 + s.nseg = 1 + s.insert('pas') + s.g_pas = 1e-5 + s.e_pas = 0.0 + s.insert('extracellular') + for seg in s: + seg.xg[0] = xg + seg.xc[0] = 1.0 + ic = h.IClamp(s(0.5)) + ic.delay = 1.0 + ic.dur = 1e9 + ic.amp = 0.1 + h.cvode_active(True) + cv.use_daspk(1) + cv.dae_init_mode(3) + cv.dae_init_stats(1) + h.finitialize(0.0) + h.continuerun(0.999) + vm0 = s(0.5).v # transmembrane Vm with extracellular + h.continuerun(1.0) + vm1 = s(0.5).v # transmembrane Vm with extracellular + st = h.Vector() + cv.dae_init_stats(st) + assert int(st[6]) == 3 and st[2] == 0, (xg, list(st)) + assert abs(vm1 - vm0) < 1e-6, (xg, vm0, vm1) + print('xg', xg, 'ok') +print('ok') +""" + assert "ok" in _run_isolated(code) + + +if __name__ == "__main__": + test_pwlclamp_ival_right_continuous() + test_iclamp_step_mode3_section() + test_pwl_istep_mode3_section() + test_pwl_kink_mode3_section() + test_pwl_jump_ramp_section_mode3() + test_e0_section_lm_parity_istep_kink() + test_e1_extracellular_pwl_istep_mode3() + test_e1_section_lm_parity_istep_kink() + test_end_electrode_mode3_path() + test_z0_end_ri_lm_parity_istep() + test_z1_end_ri_lm_parity_istep() + test_e2_xc0_extracellular_istep() + test_e3_multilayer_xc_istep() + test_e4_xg_scale_istep() + print("all ok") diff --git a/test/unit_tests/vecplay_tplus.cpp b/test/unit_tests/vecplay_tplus.cpp new file mode 100644 index 0000000000..3febfd7b90 --- /dev/null +++ b/test/unit_tests/vecplay_tplus.cpp @@ -0,0 +1,144 @@ +#include +#include + +#include +#include + +#include "vecplay_tplus.h" + +using Catch::Matchers::WithinAbs; +using Catch::Matchers::WithinRel; + +static double fd_deriv(const std::vector& y, + const std::vector& t, + double tt, + int ubound, + double h = 1e-9) { + double v0 = 0., v1 = 0.; + REQUIRE(nrn_vecplay_continuous_tplus( + (int) y.size(), y.data(), t.data(), tt, ubound, &v0, nullptr) == 0); + REQUIRE(nrn_vecplay_continuous_tplus( + (int) y.size(), y.data(), t.data(), tt + h, ubound, &v1, nullptr) == 0); + return (v1 - v0) / h; +} + +TEST_CASE("vecplay continuous t+ before first knot", "[vecplay][tplus]") { + std::vector y{1., 2., 3.}; + std::vector t{0., 1., 2.}; + const int ub = 2; + double v = 0., d = 0.; + REQUIRE(nrn_vecplay_continuous_tplus(3, y.data(), t.data(), -1.0, ub, &v, &d) == 0); + REQUIRE_THAT(v, WithinAbs(1.0, 1e-15)); + REQUIRE_THAT(d, WithinAbs(0.0, 1e-15)); + // At t0: value y0, classical right derivative = outgoing slope (1 here) + REQUIRE(nrn_vecplay_continuous_tplus(3, y.data(), t.data(), 0.0, ub, &v, &d) == 0); + REQUIRE_THAT(v, WithinAbs(1.0, 1e-15)); + REQUIRE_THAT(d, WithinAbs(1.0, 1e-15)); +} + +TEST_CASE("vecplay continuous t+ interior and kink uses outgoing slope", "[vecplay][tplus]") { + // Piecewise linear: slope 1 on [0,1], slope 2 on [1,2] + std::vector y{0., 1., 3.}; + std::vector t{0., 1., 2.}; + const int ub = 2; + double v = 0., d = 0.; + + REQUIRE(nrn_vecplay_continuous_tplus(3, y.data(), t.data(), 0.5, ub, &v, &d) == 0); + REQUIRE_THAT(v, WithinAbs(0.5, 1e-15)); + REQUIRE_THAT(d, WithinAbs(1.0, 1e-15)); + + // At knot t=1: outgoing segment slope 2, value 1 + REQUIRE(nrn_vecplay_continuous_tplus(3, y.data(), t.data(), 1.0, ub, &v, &d) == 0); + REQUIRE_THAT(v, WithinAbs(1.0, 1e-15)); + REQUIRE_THAT(d, WithinAbs(2.0, 1e-15)); + + REQUIRE(nrn_vecplay_continuous_tplus(3, y.data(), t.data(), 1.5, ub, &v, &d) == 0); + REQUIRE_THAT(v, WithinAbs(2.0, 1e-15)); + REQUIRE_THAT(d, WithinAbs(2.0, 1e-15)); +} + +TEST_CASE("vecplay continuous t+ extrapolates last two points past end", "[vecplay][tplus]") { + std::vector y{0., 1., 3.}; // last slope (3-1)/(2-1) = 2 + std::vector t{0., 1., 2.}; + const int ub = 2; + double v = 0., d = 0.; + REQUIRE(nrn_vecplay_continuous_tplus(3, y.data(), t.data(), 2.0, ub, &v, &d) == 0); + REQUIRE_THAT(v, WithinAbs(3.0, 1e-15)); + REQUIRE_THAT(d, WithinAbs(2.0, 1e-15)); + REQUIRE(nrn_vecplay_continuous_tplus(3, y.data(), t.data(), 3.0, ub, &v, &d) == 0); + REQUIRE_THAT(v, WithinAbs(5.0, 1e-15)); // 3 + 2*(3-2) + REQUIRE_THAT(d, WithinAbs(2.0, 1e-15)); +} + +TEST_CASE("vecplay continuous t+ flat last segment hold after end", "[vecplay][tplus]") { + std::vector y{0., 1., 1.}; + std::vector t{0., 1., 2.}; + const int ub = 2; + double v = 0., d = 0.; + REQUIRE(nrn_vecplay_continuous_tplus(3, y.data(), t.data(), 5.0, ub, &v, &d) == 0); + REQUIRE_THAT(v, WithinAbs(1.0, 1e-15)); + REQUIRE_THAT(d, WithinAbs(0.0, 1e-15)); +} + +TEST_CASE("vecplay continuous t+ single sample", "[vecplay][tplus]") { + std::vector y{7.}; + std::vector t{0.}; + double v = 0., d = 0.; + REQUIRE(nrn_vecplay_continuous_tplus(1, y.data(), t.data(), 10.0, 0, &v, &d) == 0); + REQUIRE_THAT(v, WithinAbs(7.0, 1e-15)); + REQUIRE_THAT(d, WithinAbs(0.0, 1e-15)); +} + +TEST_CASE("vecplay continuous t+ coincident times average value zero deriv", "[vecplay][tplus]") { + std::vector y{0., 2., 4.}; + std::vector t{0., 1., 1.}; // jump disc at t=1 + const int ub = 2; + double v = 0., d = 0.; + // On the coincident segment at/after t=1 with ubound=2: average, deriv 0 + REQUIRE(nrn_vecplay_continuous_tplus(3, y.data(), t.data(), 1.0, ub, &v, &d) == 0); + REQUIRE_THAT(v, WithinAbs(3.0, 1e-15)); // (2+4)/2 + REQUIRE_THAT(d, WithinAbs(0.0, 1e-15)); +} + +TEST_CASE("vecplay continuous t+ matches finite difference of value", "[vecplay][tplus]") { + // Ramp like iramp: 0 until 1, then 0.5 at 1, 1.0 at 2 + std::vector y{0., 0., 0.5, 1.0, 0., 0.}; + std::vector t{0., 1., 1., 2., 2., 5.}; + // After disc at index of second t=1, active ubound for ramp segment is 3 (t=2,y=1) + // Full vector ubound = 5 + const int ub = 5; + for (double tt: {0.5, 1.0, 1.5, 2.0, 3.0, 4.0}) { + double v = 0., d = 0.; + REQUIRE(nrn_vecplay_continuous_tplus((int) y.size(), y.data(), t.data(), tt, ub, &v, &d) == + 0); + // Right FD of value should match classical deriv where smooth; at kinks + // FD is one-sided outgoing if h>0 + const double d_fd = fd_deriv(y, t, tt, ub); + REQUIRE_THAT(d, WithinAbs(d_fd, 1e-5)); + } +} + +TEST_CASE("vecplay continuous t+ iramp segment after jump at t=1", "[vecplay][tplus]") { + // Values: hold 0, jump to 0.5 at t=1, ramp to 1 at t=2, jump to 0 + // t: 0, 1, 1, 2, 2, 5 y: 0, 0, 0.5, 1, 0, 0 + std::vector y{0., 0., 0.5, 1.0, 0., 0.}; + std::vector t{0., 1., 1., 2., 2., 5.}; + double v = 0., d = 0.; + + // During ramp: ubound typically at least the end of ramp knot (index 3) + // With full ubound=5, at tt=1.5 we are on segment between t[2]=1,y=0.5 and t[3]=2,y=1 + REQUIRE(nrn_vecplay_continuous_tplus(6, y.data(), t.data(), 1.5, 5, &v, &d) == 0); + REQUIRE_THAT(v, WithinAbs(0.75, 1e-15)); + REQUIRE_THAT(d, WithinAbs(0.5, 1e-15)); + + // At t=1+ with full vector: after coincident pair, outgoing is ramp slope 0.5 + // search at tt=1: first t[j] > 1 → j=3 (t[3]=2), segment (2,3): 0.5→1, slope 0.5 + REQUIRE(nrn_vecplay_continuous_tplus(6, y.data(), t.data(), 1.0, 5, &v, &d) == 0); + REQUIRE_THAT(v, WithinAbs(0.5, 1e-15)); + REQUIRE_THAT(d, WithinAbs(0.5, 1e-15)); +} + +TEST_CASE("vecplay continuous t+ bad args", "[vecplay][tplus]") { + double v = 0., d = 0.; + REQUIRE(nrn_vecplay_continuous_tplus(0, nullptr, nullptr, 0., 0, &v, &d) == -1); +}