Take a gradient at a collar psi you supply

[dfalster]

Stacked on #89 — review that one first; this branch targets it.

Both implementations called find_root_collar_psi() and read opt_root_psi_
back, so a model that TRACKS the optimum rather than finding it — TF24f carries
the collar as an ODE state — got a confident answer about the re-solved optimum
instead, with nothing saying so. psi evaluates there; dpsi_dtheta (default
zero) is the indirect term, since an imposed collar's response is the caller's to
know. M, H and dY_dpsi come back as the coefficients of their sensitivity
ODE.

stationarity no longer routes but still decides profit's column, so
psi = psi* with dpsi_dtheta = -M/H reproduces the solving path bit-for-bit.
A clamped psi reports "clamped" and no gradient: the collar moves with the
bound, not with dpsi_dtheta. Closes #88.

[dfalster]

The shape of it

psi and dpsi_dtheta on leaf_gradient() and leaf_gradient_batch();
gradient::Prescribed plus psi/dpsi_dtheta arguments to gradient::batch
and an optional const Prescribed* on gradient::at. Status gains
Prescribed and Clamped; Result gains M, dY_dpsi and psi.

One composite serves both paths. gradient_ift / .gradient_ift() take
dpsi_dtheta: null means "derive it by the IFT", which is right for a collar
that was solved for; non-null means the caller imposed the collar and knows how
it moves. That is deliberate rather than economical — it is what makes the
bit-for-bit equivalence below a real assertion instead of two implementations
that happen to line up.

The stationarity test: refined from what #88 proposed

The issue said the prescribed path should "skip the stationarity test", and the
routing half of that is right and is what the code does — method is refused,
there is no ift-vs-fd choice, and differencing the solve would answer about the
optimum instead of about the point.

But the number keeps its meaning and gains a better one: how far the collar
you supplied sits from the optimum, in MPa. It is reported, and it makes exactly
one decision — profit's. The rule, one sentence for both paths:

use the best available value of dprofit/dpsi.

At a stationary point that is the analytic zero (#89's envelope). Away from one
it is dprofit_droot_collar_psi, which is exact in ψ and already computed —
strictly better than differencing the same quantity, which is what the other four
columns must do.

⚠️ The test is stationarity, not status. Routing on status — which
carries "prescribed" on one path and "interior" on the other — would look
equivalent and silently break the equality below.

Verification

• psi = <the solved psi*> with dpsi_dtheta = -M/H reproduces the solving
  path bit-for-bit. identical() on gradient, M, dY_dpsi and psi, in
  both leaf_gradient() and leaf_gradient_batch(), and the batch's prescribed
  rows are identical() to leaf_gradient() called alone. Also asserted the
  other way: with dpsi_dtheta left at its default the answers differ by more
  than 1e-12, so the equality is not an artefact of the default.
• Off-optimum, against a central difference of A and profit taken the long
  way with the collar held fixed. And dY_dpsi[["profit"]] is identical() to
  dprofit_droot_collar_psi at that ψ, so "exact, not differenced" is asserted
  rather than described.
• Clamp: outside the interval, exactly at the bound, and one step inside it.
• Full R suite 0 failures, gradient_golden.tsv: worst relative difference 0.
  C++ 463 checks / 0 failures, 288 golden points bit-identical, make bench
  builds. The solving path is untouched — no golden cell moved.

⚠️ Two bugs the tests found, both real

1. Null dpsi_dtheta meant two different things and C++ used the wrong one.
gradient_ift's null is "derive -M/H". On the prescribed path a null means the
opposite — the collar does not move with θ, so the answer is the partial. The
first version passed the null straight through, so a batch asking for a
prescribed gradient with no dpsi_dtheta silently got the solving path's
indirect term added on top. Caught by the batch-vs-leaf_gradient() equality
test, which is exactly the kind of thing it exists for: the R side had the
semantics right, so only a cross-implementation comparison could see it. Fixed in
at() rather than in batch() so every C++ consumer gets it — plant included.

2. dpsi_dtheta was validated inside the composite, which a clamped point
never reaches. So a wrong-length vector errored at some ψ and was silently
accepted at others. Moved up front, next to the pars check. Argument checking
cannot be conditional on what the model does with the arguments.

A third, smaller one, found while fixing the second: the length-1 recycling
branch overwrote the caller's names with pars before the name check, so
dpsi_dtheta = c(stem_b = 1) against pars = "vcmax_25" was applied to
vcmax_25 — a named argument quietly meaning a different parameter, which is the
failure naming is supposed to prevent. Names are checked first now.

The clamp policy, which is the design decision

evaluate_root_collar_psi clamps, and outputs_at() already rejected clamped
evaluations on exact equality. On the prescribed path clamping is routine —
TF24f depends on it to pull an out-of-range tracked state back inside — so
"reject" needed a policy rather than an error.

The chosen one: status = "clamped", gradient/M/dY_dpsi all NA, and
value/H/stationarity/psi kept, because those describe the point rather
than the derivative. Not thrown, because a fit visits these routinely and a
per-row status is how a batch says which rows to distrust. Not the direct term
either: the collar used is min(max(psi, a(θ)), b(θ)), so it moves with the
bound, and the direct term alone is the same plausible-and-wrong answer the
active-set guard exists to refuse. It fires for a ψ within one step of an end
too, where dY/dpsi cannot be centred.

⚠️ NA payload differs between the two implementations, deliberately. R uses
NA_real_ and C++ util::na_value (a quiet NaN), because gradient.hpp is
R-free and cannot reach NA_REAL. This matches the existing per-row error path.
The contract is is.na(), not the payload, and the tests assert it that way.

Not done here

No plain-C++ test of gradient::at with a Prescribed — the namespace has never
had one, and the R↔C++ bit-for-bit comparison covers the numerics. Worth adding
when traitecoevo/plant#614 becomes the first C++ consumer, since that is the
point at which the API shape needs a guard of its own.

[aornugent]

The design is right and the load-bearing test is the right one. Keeping a single
composite and asserting that psi = psi* with dpsi_dtheta = -M/H reproduces the
solving path with identical() is what makes this a second entry point rather than a
second implementation, and refusing a gradient at a clamped psi — instead of handing
back the direct term, which would be plausible and wrong — is the correct call for the
same reason the pinned case gets the fallback. Returning M, H and dY_dpsi so a
caller can integrate ds/dt = k(M + H s) themselves is a better boundary than trying
to guess the caller's dynamics.

I read the R and C++ at() against each other line by line and they agree
structurally: resid is in scope where the prescribed branch uses it, dY_dpsi is
filled on every path that later reads it, the zeros materialisation in at() is in
the right place, and the dpsi buffer in batch() outlives the Prescribed pointing
into it. Four things below, the first of which is a defect.

1. An integer psi is silently reported as clamped

leaf_gradient() calls

.gradient_check_psi(psi, dpsi_dtheta, method)

and discards the return; the helper ends in invisible(NULL) and never coerces. The
clamp test is then

clamped <- !identical(l$opt_root_psi_, psi)

and identical(3, 3L) is FALSE. So leaf_gradient(..., psi = 3L) — or the entirely
ordinary for (p in 2:5) leaf_gradient(..., psi = p), since 2:5 yields integers —
returns status = "clamped", an all-NA gradient, and a psi in the result equal to
what was asked for, so the "where it was pulled to" diagnostic shows no movement and
there is nothing to notice.

Three things make this worth fixing rather than documenting:

• the batch counterpart does not have it — .gradient_check_psi_batch() returns
  as.numeric(psi), so the two entry points disagree on the same input;
• C++ does not have it either — util::identical(double, double) is a == b, so the
  R reference and the implementation it is required to match bit-for-bit differ in
  which inputs they call clamped;
• the sibling this cites as precedent uses the type-tolerant form:
  .gradient_outputs_at() tests l$opt_root_psi_ != psi.

One line either way: assign the checker's return and have it as.numeric(), or use
!=. Worth a test with a non-double psi — nothing in the new tests passes anything
but a double literal or a value read back out of a previous result, which is why this
survives the suite.

(Found by reading plus checking identical(3, 3L) at the R prompt, not from a built
run — so worth confirming against a build before fixing.)

2. The prescribed path has no no-gradient, and adopts the shut-down sentinel

On the solving path usable == FALSE is a hard stop: dprofit_droot_collar_psi
returns a sentinel zero at a shut-down point rather than a derivative, and the code
says so. On the prescribed path use_ift = !clamped regardless of usable, and then

} else if (prescribed) {
  dY_dpsi[["profit"]] <- resid
}

adopts that sentinel as if it were dprofit/dpsi. With the default dpsi_dtheta = 0
it costs nothing. With a non-zero one — which is the case this feature exists for — the
profit column silently loses its indirect term at exactly the operating points where a
tracking model spends its dry-patch time. Either the prescribed path needs its own
no-gradient, or there needs to be a test at a shut-down driver with a non-clamped
psi showing the state is unreachable. Right now the Prescribed/Clamped pair is
the whole vocabulary, and usable is computed and then not consulted.

Relatedly, .gradient_dpsi_dtheta checks is.numeric(dpsi_dtheta) && !anyNA(...) but
not is.finite. Inf therefore reaches the composite, where #89's profit line is
direct + rounded(0 * d_psi) — giving ±Inf in four columns and NaN in profit alone.
Cheapest fix is on the #89 side (assign profit's column rather than multiply); a
finiteness check here would also do it.

3. Status reordering, and the default: arm

Error moves from 3 to 5. Source-compatible, not value-compatible, so any consumer
holding the integer breaks quietly. Separately, status_name's default: return "error" means the next member added anywhere before Error gets labelled "error"
rather than failing to compile — an explicit case Status::Error: with no default:
makes the switch exhaustive and turns that into a compiler diagnostic.

4. Minor: method is per call, and its comment says per row

In src/gradient.cpp,

method_used[i] = psi_p != nullptr ? "prescribed" : (r.used_ift ? "ift" : "fd");

psi_p is call-level, so the ternary does not vary across i — but the comment above
it justifies the placement as "per row rather than per call because a batch can carry a
psi for every row and still have one throw". A row that threw reports
status = "error" with method = "prescribed". Defensible, but the comment is
describing behaviour the expression does not have.

On the stack

This targets #89's branch, so it inherits that PR's rebase problem: 26574f4 moved
every cell of gradient_golden.tsv and landed a competing expect_golden() tolerance
design. Worth settling underneath before this one moves.

[dfalster]

All four addressed, pushed as 1461b48 (rebased onto the updated #89, which is
now on top of master).

1. Integer psi — confirmed against a build, and fixed

You were right to flag it as unconfirmed; it reproduces. identical(3, 3L) is
FALSE, the checker ended in invisible(NULL), and the return was discarded —
so leaf_gradient(..., psi = 3L) came back status = "clamped" with an all-NA
gradient and a psi in the result showing no movement.

Fixed the way you suggested and by the route the two siblings already used:
.gradient_check_psi() now ends in as.numeric(psi) and the caller assigns it.
I chose coercion over != so that the value stored in the result is the same
double the comparison used — with != an integer psi would pass the clamp test
but still be reported back as an integer.

There is a test, and it has teeth in both entry points: it asserts
psi = 3 and psi = 3L give identical() gradients and both report
"prescribed", and that leaf_gradient_batch() agrees on the same input.

2. The shut-down sentinel — this was the real one, and it needed a new binding

You are right that usable is computed and then not consulted, and right that
use_ift = !clamped lets the sentinel through. What I could not do was gate on
usable: it demands H < 0, which is a maximum test, and a prescribed collar
away from the optimum may legitimately sit where profit is convex. Nothing on
that path divides by H, so that is fine and must stay fine.

What actually disables the point is infeasibility — and the flag for it has
existed since #79 while being unreachable from R, because RcppR6 has no form for
a bool* and the generated binding silently drops it. So every R-side composite
was the "composite that ignores it inherits the bug" case the C++ vignette
warns about; #90 just made it reachable.

Added Leaf$dprofit_droot_collar_psi_checked(psi) returning
{dprofit, feasible}, used by R; C++ reads the out-parameter directly. An
infeasible prescribed point is now "no-gradient" with an all-NA gradient.

The test is the case clamped cannot catch, which is the one you asked for: it
hands back exactly the collar the shut-down state seated, verifies it lands
unclamped (identical(g$psi, solved$psi)), and asserts "no-gradient".
Confirmed reachable and previously wrong:

prescribed psi lands unclamped: TRUE
dprofit there: 0   feasible: FALSE     <- previously adopted as dprofit/dpsi

On the Inf half — fixed on the #89 side as you suggested, by assigning profit's
column rather than multiplying. I did not add an is.finite check to
.gradient_dpsi_dtheta: a non-finite dpsi_dtheta is now handled rather than
rejected, and an infinite sensitivity is a legitimate thing for a caller's ODE to
be carrying into a report.

3. Status — Error stays at 3, and the switch is exhaustive

Prescribed and Clamped are appended after Error, so no existing integer
value moves. Out of narrative order, with a comment saying why so the next person
does not tidy it.

default: is gone; every member has a case. There is a trailing return "error" after the switch for -Wreturn-type on an out-of-range cast, which
keeps the diagnostic you want on a new member while staying warning-clean.

4. method comment — corrected

You are right that the expression is call-level. The comment now says so, and
says what falls out of it: a row that threw reports status = "error" with
method = "prescribed", because method names the route that was asked for and
status says what happened to the row.

---

⚠️ One thing worth flagging about the rebase

The first rebase attempt committed conflict markers into six files — I ran
git add -A && git rebase --continue and only read the tail of git's output, so
I resolved the one conflict I saw and staged five I did not. Caught by grepping
for markers before building, fixed in the amended commit, and the suites are the
proof rather than my say-so. Not a resolution I would ask anyone to re-derive
from the diff, so: every conflicted hunk kept both sides — #89's envelope
assignment and this PR's caller-supplied d_psi — and the composite signature is
now (..., dpsi_dtheta, envelope, M_out, out).

Local: R suite 0 failures, gradient_golden.tsv: worst relative difference 0
(so R and C++ still agree bit-for-bit on all five columns), C++ 486 checks / 0
failures, 576 golden points bit-identical, make bench builds.