Using phylloptim from C++

The two umbrella headers, and why the extension matters

#include <phylloptim.hpp>    // the model. Plain C++, no R, no Rcpp.
#include <phylloptim.h>      // 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 <phylloptim/RcppR6_pre.hpp>, 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 <phylloptim.hpp> 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.

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

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

#include <phylloptim.hpp>
#include <cstdio>

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<double> psi_soil{2.0};     // positive magnitudes, MPa
  std::vector<double> 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++

<phylloptim/gradient.hpp>, 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.

#include <phylloptim.hpp>
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<g::Result> 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:

const double psi[2] = {2.7, 2.8};              // positive magnitudes, MPa
const std::vector<g::Result> 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:

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:

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, 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.