--- title: "Using phylloptim from C++" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Using phylloptim from C++} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} knitr::opts_chunk$set(collapse = TRUE, comment = "#>") ``` The model is a set of self-contained C++ headers. **They use no R and no Rcpp**, depend only on Boost and the header-only parts of `odelia`, and compile and run with no R installed at all. That is a guarantee rather than an accident, and `.github/workflows/cpp-tests.yml` enforces it by building the whole C++ suite on runners with no R on them. This vignette contains no runnable R. It documents the other half of the package, which R users never see and `plant` uses exclusively. ## The two umbrella headers, and why the extension matters ``` #include // the model. Plain C++, no R, no Rcpp. #include // the R binding layer. Pulls in Rcpp. ``` If you are embedding the model, you want `.hpp`. The `.h` file exists because `RcppR6` hardwires both that name and ``, so the R-facing umbrella had to take it and the model's had to move. A file extension carries the whole one-way dependency rule, which is more weight than an extension should bear, hence the banner at the top of each. The rule: `src/` and `R/` include *downward* into `inst/include/`, never the reverse. Nothing reachable from `` includes Rcpp. The CMake install rules encode the same split, shipping `*.hpp` and **excluding** `phylloptim.h` and the generated `RcppR6_*.hpp`, so a non-R consumer is never handed a header they cannot compile. ## Building Two dependencies, both header-only: Boost (root finders and the incomplete gamma function) and `odelia` (the supplied-derivative nodes and vendored XAD). `odelia` is distributed as an R package but the headers used here are plain C++, so a git checkout is enough. ```sh git clone https://github.com/traitecoevo/odelia cmake -B build -DPHYLLOPTIM_ODELIA_INCLUDE_DIR=$PWD/odelia/inst/include cmake --build build ctest --test-dir build # the C++ suite, including the golden file cmake --install build --prefix /usr/local ``` ```cmake find_package(phylloptim REQUIRED) target_link_libraries(my_program PRIVATE phylloptim::phylloptim) ``` An INTERFACE target: headers, an include path and `cxx_std_20`, nothing to link. `add_subdirectory()` works too if you would rather vendor it. **Build with optimisation on if you intend to compare against the golden file.** Measured on macOS/arm64, `-O1`, `-O2` and `-O3` all reproduce `tests/cpp/golden/operating_points.tsv` bit-for-bit and `-O0` does not, missing by about 13 ULP because it declines to contract `a*b + c` into an FMA. A debug build failing `test_golden` by ~1e-15 has found nothing. ## A solve ```cpp #include #include int main() { phylloptim::Leaf l; // default Eucalyptus saligna traits // A Leaf defaults to TF24_floor on the leaf route over one soil potential; // a root network and the collar solve, as plant runs them, are seated. l.set_supply_multilayer(); l.set_model(phylloptim::Leaf::CostCurve::TF24, /*collar_route=*/true); std::vector psi_soil{2.0}; // positive magnitudes, MPa std::vector soil_depth{1.0}; // m // The per-layer root RESISTANCES, per unit leaf area. If you have root carbon // instead, this is the architecture model that maps one to the other -- but it // is a helper you call, not something the solve does (#33). Supply the // resistances directly if that is what you have. const phylloptim::RootNetwork roots = phylloptim::root_network_from_carbon( /*root carbon per m2 LEAF*/ {20.0}, phylloptim::layer_thicknesses(soil_depth), /*beta_R_H*/ 3.4e2, /*beta_R_V*/ 9.4e3); l.set_physiology(roots, /*PPFD*/ 900, psi_soil, soil_depth, /*kmax*/ 3.14e-5, /*atm_vpd*/ 2.0, /*ca*/ 40.0, /*leaf_temp*/ 25.0, /*atm_o2_kpa*/ 21.0, /*atm_kpa*/ 101.3); l.find_root_collar_psi(); std::printf("A = %.6f, gc = %.6f, psi_leaf = %.6f\n", l.assim_colimited_, l.stom_cond_CO2_, l.opt_psi_leaf_); return 0; } ``` Every water potential is a **positive magnitude in MPa**, and that is asserted rather than documented: a negative `psi_soil` is an error, not a sign convention the model quietly accepts. ## Trait gradients from C++ ``, reached from the umbrella header, gives the trait gradient at a solved point with no R anywhere, as one forward pass of the derivative surface carrying every parameter, and batches it over observations. ```cpp #include namespace g = phylloptim::gradient; g::Drivers d; // one observation's drivers d.root_network = roots; d.PPFD = 900; d.psi_soil = {2.0}; d.soil_depth = {1.0}; d.atm_vpd = 2.0; d.ca = 40.0; d.leaf_temp = 25.0; d.atm_o2_kpa = 21.0; d.atm_kpa = 101.3; double theta[g::n_theta] = {/* 16 traits, then kmax, resistance, the three model-owned slots, the damage inputs */}; const int pars[1] = {0}; // g::par_names()[0] == "vcmax_25" phylloptim::Leaf l; l.set_supply_multilayer(); g::Settings s; // the Leaf's defaults, unless set: s.model = {phylloptim::Leaf::CostCurve::TF24, g::Route::Collar}; const std::vector out = g::batch(l, theta, 1, {d, d}, /*single=*/false, pars, 1, s); // out[i].grad[k * g::n_outputs + j] is d(output j)/d(pars[k]) // out[i].status is Interior / Pinned / NoGradient / Error ``` To evaluate at a collar potential you are carrying rather than one to be solved for, which is the case a model that *tracks* the optimum has, pass `psi`, one per observation, and optionally an `n × npars` column-major `dpsi_dtheta`: ```cpp const double psi[2] = {2.7, 2.8}; // positive magnitudes, MPa const std::vector out = g::batch(l, theta, 1, {d, d}, /*single=*/false, pars, 1, s, psi, /*dpsi_dtheta=*/nullptr); // null here means ZERO... // out[i].status is now Prescribed or Refused // out[i].M / out[i].H are the coefficients of ds/dt = k(M + H s) ``` A null `dpsi_dtheta` means zero: the collar was imposed, so it moves only as you say it does. For one observation, use `g::at()` with a `g::Prescribed`. `theta` is `g::n_theta` = **35** doubles: the sixteen traits in `set_traits`' argument order, then `leaf_specific_conductance_max`, the single-potential path's `resistance`, the three model-owned slots `CF77_lambda_`, `TF24_floor_lambda_o` and `LeastCost_beta`, and the fourteen inputs of the thermal-damage layer that reach the solve. `g::par_names()` is that order and every position has a name (`g::par_vcmax_25`, `g::par_R_d_25`, `g::par_kmax`, ...); `pars` holds positions into it. R indexes these positions, so reordering the enumeration silently differentiates the wrong parameter. `phylloptim::n_leaf_pars` is **37**, two longer: the pack the kernels read also holds `par_PPFD` and `par_NIR`, seated per observation from the drivers and never fitted. Size `theta` by `n_theta`; read it out to `n_leaf_pars` and you run past its end. **Count the initialiser against `g::n_theta` rather than against this paragraph.** A short aggregate initialiser is legal C++ and zero-fills the rest, so a `theta` written from a stale count shifts every parameter after the gap and drops the last off the end, with no diagnostic at all. That has happened: a 16-entry list under what was then an 18-element array put `kmax` in a trait's slot, and the only symptom was `dA/dvcmax = nan`. The package's own CI consumer carries a `static_assert(g::n_theta == 35)` beside its list for exactly this reason; copy it. Three things the free function gives you that a hand-rolled loop would not: * **the classification.** The solve records what kind of point it found, interior or held by a bound, and the pass places the decision variable accordingly: by the implicit function theorem at an interior optimum, on the bound at a pinned one. Where the optimum sits on a bound, `-M/H` is wrong by orders of magnitude and fails silently; `status` says which it was. * **one pass for every parameter.** The pass runs at `g::n_theta + 1` lanes, one per slot and one for the decision variable, so its cost does not depend on `pars`. * **per-row status rather than an exception**, so a proposal that reaches an unsolvable point costs that row and not the batch. **It is every trait gradient the package gives**: `leaf_gradient()` and `leaf_gradient_batch()` both call it. `Settings::model` names the curve and the route, the `Leaf`'s default (TF24_floor on the leaf route) unless set. `tests/testthat/test-gradient-referee.R` checks it against a difference of the solve for every curve on both routes, and `tests/cpp/test_leaf_surface.cpp` checks the surface itself at every kind of point. The surface is public, if you want to compose something else; plant does, to read the operating point's rows at its own scalar type: | piece | what it is | |---|---| | `supply_draw_at(collar, supply)` | the supply's draw at a collar, on any scalar | | `operating_collar(draw, pars, H, &held)` | the decision variable placed from the kind the solve recorded | | `physical_outputs_at(x, draw, pars, &held)` | A, gc, psi_leaf and profit at that placement | | `profit_curvature_at_collar()` | `H`, the curvature the implicit function theorem divides by | | `dprofit_droot_collar_psi(psi, &feasible)` | dprofit/dψ at a given collar potential, at double | **Pass the `feasible` flag.** `dprofit_droot_collar_psi` returns `0.0` on its shut-down and reversed-gradient exits, and a bare zero is indistinguishable from a stationary point, so a root-find looking for `dprofit == 0` would converge onto a shut-down state and report it as the optimum. The optional `bool*` second argument is how the two are told apart. `vignette("derivative-surface")` documents the surface: what each kind places, and what each scalar type is for. ## What is fast, and what is not `tests/cpp/bench_gradient.cpp` times the gradient beside a central difference of the whole solve, both in C++, at one interior point (TF24 on the collar route, one layer), less the seat and solve every arm pays: | arm | cost | |---|---| | the gradient, any number of parameters | about 9 µs | | a difference of the solve | about 5 µs a parameter, and tens where it rebuilds a vulnerability table (`stem_c`, `root_c`) | So the pass costs about two parameters' worth of differencing, and then nothing more: from three parameters up it is the cheaper of the two, and it is exact. That is also what an R caller gets from `leaf_gradient_batch()`: the pass runs in C++ over every observation, so the boundary is crossed once per call rather than once per observation. **Constructing a `Leaf` from R is the largest single cost in the R layer**, and most of it is R-side object construction rather than the vulnerability tables. From C++ that cost is the tables alone. Reuse the object either way. The measured figures for all of these are in the README's [Performance from R](https://github.com/traitecoevo/phylloptim#performance-from-r), the one place they are kept current. ## The golden file is the contract `tests/cpp/golden/operating_points.tsv` records 672 operating points at full precision and is compared **bit-exactly**. It is what makes a large refactor of this code checkable rather than hopeful, and any change that moves it is expected to say what moved and by how much. It is bit-exact on the platform that generated it (macOS/arm64) and compared with per-field tolerances elsewhere, because libm's `exp`/`pow` are not bit-reproducible across platforms. The fields split into two classes three orders apart: `profit` is the maximum itself and well conditioned; everything else is evaluated at the *argmax* of a flat maximum and so sqrt-amplified. **If you need a portable check of the solve, check `profit`.**