diff --git a/.gitignore b/.gitignore index e21d0379..3dca9d7f 100644 --- a/.gitignore +++ b/.gitignore @@ -23,6 +23,8 @@ Meta *.pdf *.png *.RDS +# ...except committed test fixtures (the AD-gradient regression baseline). +!tests/testthat/fixtures/gradient-baseline.rds # Overstorey guide preview figures — regenerated when the guide builds in the # Overstorey repo; not version-controlled here. diff --git a/Makefile b/Makefile index 7e896e47..3da6b833 100644 --- a/Makefile +++ b/Makefile @@ -32,6 +32,13 @@ roxygen: test: all Rscript -e 'library(methods); devtools::test()' +# AD tests self-skip under load_all (pkgload DLL). Run against the INSTALLED +# package: library(plant) loads the installed .so so the skip-gate passes; +# test_dir (not test_local) avoids re-loading via load_all; the namespace-parented +# env exposes unexported internals; parallel off so callr workers inherit the load. +test-ad: install + TESTTHAT_PARALLEL=false Rscript -e 'library(plant); e <- new.env(parent = asNamespace("plant")); testthat::test_dir("tests/testthat", filter = "ad", env = e)' + benchmark: Rscript scripts/benchmark.R @@ -49,4 +56,4 @@ check: build clean: rm -f src/*.o src/*.so src/*.o.tmp -.PHONY: all compile doc clean test attributes roxygen install build check +.PHONY: all compile doc clean test test-ad attributes roxygen install build check diff --git a/NAMESPACE b/NAMESPACE index 85789a24..b0486a28 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -43,6 +43,7 @@ export(TF24f_hyperpar) export(Weibull_Disturbance_Regime) export(add_mutant) export(add_strategies) +export(birth_rate_gradient) export(control) export(control_accurate) export(environment_type) @@ -52,6 +53,7 @@ export(export_patch_state) export(generate_strategy) export(grow_individual_to_height) export(grow_individual_to_size) +export(grow_individual_to_size_gradient) export(grow_individual_to_time) export(hyperpar) export(integrate_over_size_distribution) @@ -65,6 +67,7 @@ export(make_hyperpar) export(make_initial_state) export(mutant_parameters) export(node_schedule_times_default) +export(offspring_production_gradient) export(optimise_individual_rate_at_height_by_trait) export(optimise_individual_rate_at_size_by_trait) export(param_hyperpar) @@ -79,6 +82,8 @@ export(seq_log) export(seq_log_range) export(seq_range) export(set_initial_state) +export(stand_gradient) +export(stand_state_jacobian) export(strategy_list) export(tidy_individual) export(trait_matrix) diff --git a/R/RcppExports.R b/R/RcppExports.R index 9cfccc27..17cd8b88 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -89,6 +89,50 @@ Leaf__dprofit_droot_collar_psi <- function(obj_, opt_root_psi) { .Call('_plant_Leaf__dprofit_droot_collar_psi', PACKAGE = 'plant', obj_, opt_root_psi) } +Leaf__dprofit_dvcmax25 <- function(obj_, opt_root_psi) { + .Call('_plant_Leaf__dprofit_dvcmax25', PACKAGE = 'plant', obj_, opt_root_psi) +} + +Leaf__dprofit_dg1_TF24 <- function(obj_, opt_root_psi) { + .Call('_plant_Leaf__dprofit_dg1_TF24', PACKAGE = 'plant', obj_, opt_root_psi) +} + +Leaf__dprofit_dbeta2 <- function(obj_, opt_root_psi) { + .Call('_plant_Leaf__dprofit_dbeta2', PACKAGE = 'plant', obj_, opt_root_psi) +} + +Leaf__dprofit_dkmax <- function(obj_, opt_root_psi) { + .Call('_plant_Leaf__dprofit_dkmax', PACKAGE = 'plant', obj_, opt_root_psi) +} + +Leaf__dprofit_dEup <- function(obj_, opt_root_psi) { + .Call('_plant_Leaf__dprofit_dEup', PACKAGE = 'plant', obj_, opt_root_psi) +} + +Leaf__dprofit_db <- function(obj_, opt_root_psi) { + .Call('_plant_Leaf__dprofit_db', PACKAGE = 'plant', obj_, opt_root_psi) +} + +Leaf__dprofit_dc <- function(obj_, opt_root_psi) { + .Call('_plant_Leaf__dprofit_dc', PACKAGE = 'plant', obj_, opt_root_psi) +} + +Leaf__dprofit_djmax25 <- function(obj_, opt_root_psi) { + .Call('_plant_Leaf__dprofit_djmax25', PACKAGE = 'plant', obj_, opt_root_psi) +} + +Leaf__dprofit_da <- function(obj_, opt_root_psi) { + .Call('_plant_Leaf__dprofit_da', PACKAGE = 'plant', obj_, opt_root_psi) +} + +Leaf__dprofit_dcurv_elec <- function(obj_, opt_root_psi) { + .Call('_plant_Leaf__dprofit_dcurv_elec', PACKAGE = 'plant', obj_, opt_root_psi) +} + +Leaf__dprofit_dcurv_colim <- function(obj_, opt_root_psi) { + .Call('_plant_Leaf__dprofit_dcurv_colim', PACKAGE = 'plant', obj_, opt_root_psi) +} + Leaf__psi_stem_to_ci <- function(obj_, psi_stem, psi_upstream) { .Call('_plant_Leaf__psi_stem_to_ci', PACKAGE = 'plant', obj_, psi_stem, psi_upstream) } @@ -953,6 +997,10 @@ Individual___FF16__FF16_Env__net_mass_production_dt <- function(obj_, environmen .Call('_plant_Individual___FF16__FF16_Env__net_mass_production_dt', PACKAGE = 'plant', obj_, environment) } +Individual___FF16__FF16_Env__growth_rate_gradient_exact <- function(obj_, environment) { + .Call('_plant_Individual___FF16__FF16_Env__growth_rate_gradient_exact', PACKAGE = 'plant', obj_, environment) +} + Individual___FF16__FF16_Env__reset_mortality <- function(obj_) { invisible(.Call('_plant_Individual___FF16__FF16_Env__reset_mortality', PACKAGE = 'plant', obj_)) } @@ -1041,6 +1089,10 @@ Individual___TF24__TF24_Env__net_mass_production_dt <- function(obj_, environmen .Call('_plant_Individual___TF24__TF24_Env__net_mass_production_dt', PACKAGE = 'plant', obj_, environment) } +Individual___TF24__TF24_Env__growth_rate_gradient_exact <- function(obj_, environment) { + .Call('_plant_Individual___TF24__TF24_Env__growth_rate_gradient_exact', PACKAGE = 'plant', obj_, environment) +} + Individual___TF24__TF24_Env__reset_mortality <- function(obj_) { invisible(.Call('_plant_Individual___TF24__TF24_Env__reset_mortality', PACKAGE = 'plant', obj_)) } @@ -1129,6 +1181,10 @@ Individual___TF24f__TF24_Env__net_mass_production_dt <- function(obj_, environme .Call('_plant_Individual___TF24f__TF24_Env__net_mass_production_dt', PACKAGE = 'plant', obj_, environment) } +Individual___TF24f__TF24_Env__growth_rate_gradient_exact <- function(obj_, environment) { + .Call('_plant_Individual___TF24f__TF24_Env__growth_rate_gradient_exact', PACKAGE = 'plant', obj_, environment) +} + Individual___TF24f__TF24_Env__reset_mortality <- function(obj_) { invisible(.Call('_plant_Individual___TF24f__TF24_Env__reset_mortality', PACKAGE = 'plant', obj_)) } @@ -1217,6 +1273,10 @@ Individual___K93__K93_Env__net_mass_production_dt <- function(obj_, environment) .Call('_plant_Individual___K93__K93_Env__net_mass_production_dt', PACKAGE = 'plant', obj_, environment) } +Individual___K93__K93_Env__growth_rate_gradient_exact <- function(obj_, environment) { + .Call('_plant_Individual___K93__K93_Env__growth_rate_gradient_exact', PACKAGE = 'plant', obj_, environment) +} + Individual___K93__K93_Env__reset_mortality <- function(obj_) { invisible(.Call('_plant_Individual___K93__K93_Env__reset_mortality', PACKAGE = 'plant', obj_)) } @@ -2145,6 +2205,102 @@ Patch___FF16__FF16_Env__state__get <- function(obj_) { .Call('_plant_Patch___FF16__FF16_Env__state__get', PACKAGE = 'plant', obj_) } +Patch___FF16__FF16_Env__step_history__get <- function(obj_) { + .Call('_plant_Patch___FF16__FF16_Env__step_history__get', PACKAGE = 'plant', obj_) +} + +Patch___FF16__FF16_Env__step_history__set <- function(obj_, value) { + invisible(.Call('_plant_Patch___FF16__FF16_Env__step_history__set', PACKAGE = 'plant', obj_, value)) +} + +Patch___FF16__FF16_Env__environment_history__get <- function(obj_) { + .Call('_plant_Patch___FF16__FF16_Env__environment_history__get', PACKAGE = 'plant', obj_) +} + +Patch___FF16__FF16_Env__environment_history__set <- function(obj_, value) { + invisible(.Call('_plant_Patch___FF16__FF16_Env__environment_history__set', PACKAGE = 'plant', obj_, value)) +} + +Patch___FF16__FF16_Env__stand_height_history__get <- function(obj_) { + .Call('_plant_Patch___FF16__FF16_Env__stand_height_history__get', PACKAGE = 'plant', obj_) +} + +Patch___FF16__FF16_Env__stand_height_history__set <- function(obj_, value) { + invisible(.Call('_plant_Patch___FF16__FF16_Env__stand_height_history__set', PACKAGE = 'plant', obj_, value)) +} + +Patch___FF16__FF16_Env__stand_competition_history__get <- function(obj_) { + .Call('_plant_Patch___FF16__FF16_Env__stand_competition_history__get', PACKAGE = 'plant', obj_) +} + +Patch___FF16__FF16_Env__stand_competition_history__set <- function(obj_, value) { + invisible(.Call('_plant_Patch___FF16__FF16_Env__stand_competition_history__set', PACKAGE = 'plant', obj_, value)) +} + +Patch___FF16__FF16_Env__stand_height_stage_history__get <- function(obj_) { + .Call('_plant_Patch___FF16__FF16_Env__stand_height_stage_history__get', PACKAGE = 'plant', obj_) +} + +Patch___FF16__FF16_Env__stand_height_stage_history__set <- function(obj_, value) { + invisible(.Call('_plant_Patch___FF16__FF16_Env__stand_height_stage_history__set', PACKAGE = 'plant', obj_, value)) +} + +Patch___FF16__FF16_Env__stand_competition_stage_history__get <- function(obj_) { + .Call('_plant_Patch___FF16__FF16_Env__stand_competition_stage_history__get', PACKAGE = 'plant', obj_) +} + +Patch___FF16__FF16_Env__stand_competition_stage_history__set <- function(obj_, value) { + invisible(.Call('_plant_Patch___FF16__FF16_Env__stand_competition_stage_history__set', PACKAGE = 'plant', obj_, value)) +} + +Patch___FF16__FF16_Env__stand_newnode_height_stage_history__get <- function(obj_) { + .Call('_plant_Patch___FF16__FF16_Env__stand_newnode_height_stage_history__get', PACKAGE = 'plant', obj_) +} + +Patch___FF16__FF16_Env__stand_newnode_height_stage_history__set <- function(obj_, value) { + invisible(.Call('_plant_Patch___FF16__FF16_Env__stand_newnode_height_stage_history__set', PACKAGE = 'plant', obj_, value)) +} + +Patch___FF16__FF16_Env__stand_newnode_competition_stage_history__get <- function(obj_) { + .Call('_plant_Patch___FF16__FF16_Env__stand_newnode_competition_stage_history__get', PACKAGE = 'plant', obj_) +} + +Patch___FF16__FF16_Env__stand_newnode_competition_stage_history__set <- function(obj_, value) { + invisible(.Call('_plant_Patch___FF16__FF16_Env__stand_newnode_competition_stage_history__set', PACKAGE = 'plant', obj_, value)) +} + +Patch___FF16__FF16_Env__stand_height_stage_history_all__get <- function(obj_) { + .Call('_plant_Patch___FF16__FF16_Env__stand_height_stage_history_all__get', PACKAGE = 'plant', obj_) +} + +Patch___FF16__FF16_Env__stand_height_stage_history_all__set <- function(obj_, value) { + invisible(.Call('_plant_Patch___FF16__FF16_Env__stand_height_stage_history_all__set', PACKAGE = 'plant', obj_, value)) +} + +Patch___FF16__FF16_Env__stand_competition_stage_history_all__get <- function(obj_) { + .Call('_plant_Patch___FF16__FF16_Env__stand_competition_stage_history_all__get', PACKAGE = 'plant', obj_) +} + +Patch___FF16__FF16_Env__stand_competition_stage_history_all__set <- function(obj_, value) { + invisible(.Call('_plant_Patch___FF16__FF16_Env__stand_competition_stage_history_all__set', PACKAGE = 'plant', obj_, value)) +} + +Patch___FF16__FF16_Env__stand_newnode_height_stage_history_all__get <- function(obj_) { + .Call('_plant_Patch___FF16__FF16_Env__stand_newnode_height_stage_history_all__get', PACKAGE = 'plant', obj_) +} + +Patch___FF16__FF16_Env__stand_newnode_height_stage_history_all__set <- function(obj_, value) { + invisible(.Call('_plant_Patch___FF16__FF16_Env__stand_newnode_height_stage_history_all__set', PACKAGE = 'plant', obj_, value)) +} + +Patch___FF16__FF16_Env__stand_newnode_competition_stage_history_all__get <- function(obj_) { + .Call('_plant_Patch___FF16__FF16_Env__stand_newnode_competition_stage_history_all__get', PACKAGE = 'plant', obj_) +} + +Patch___FF16__FF16_Env__stand_newnode_competition_stage_history_all__set <- function(obj_, value) { + invisible(.Call('_plant_Patch___FF16__FF16_Env__stand_newnode_competition_stage_history_all__set', PACKAGE = 'plant', obj_, value)) +} + Patch___TF24__TF24_Env__ctor <- function(parameters, environment, control) { .Call('_plant_Patch___TF24__TF24_Env__ctor', PACKAGE = 'plant', parameters, environment, control) } @@ -2261,6 +2417,102 @@ Patch___TF24__TF24_Env__state__get <- function(obj_) { .Call('_plant_Patch___TF24__TF24_Env__state__get', PACKAGE = 'plant', obj_) } +Patch___TF24__TF24_Env__step_history__get <- function(obj_) { + .Call('_plant_Patch___TF24__TF24_Env__step_history__get', PACKAGE = 'plant', obj_) +} + +Patch___TF24__TF24_Env__step_history__set <- function(obj_, value) { + invisible(.Call('_plant_Patch___TF24__TF24_Env__step_history__set', PACKAGE = 'plant', obj_, value)) +} + +Patch___TF24__TF24_Env__environment_history__get <- function(obj_) { + .Call('_plant_Patch___TF24__TF24_Env__environment_history__get', PACKAGE = 'plant', obj_) +} + +Patch___TF24__TF24_Env__environment_history__set <- function(obj_, value) { + invisible(.Call('_plant_Patch___TF24__TF24_Env__environment_history__set', PACKAGE = 'plant', obj_, value)) +} + +Patch___TF24__TF24_Env__stand_height_history__get <- function(obj_) { + .Call('_plant_Patch___TF24__TF24_Env__stand_height_history__get', PACKAGE = 'plant', obj_) +} + +Patch___TF24__TF24_Env__stand_height_history__set <- function(obj_, value) { + invisible(.Call('_plant_Patch___TF24__TF24_Env__stand_height_history__set', PACKAGE = 'plant', obj_, value)) +} + +Patch___TF24__TF24_Env__stand_competition_history__get <- function(obj_) { + .Call('_plant_Patch___TF24__TF24_Env__stand_competition_history__get', PACKAGE = 'plant', obj_) +} + +Patch___TF24__TF24_Env__stand_competition_history__set <- function(obj_, value) { + invisible(.Call('_plant_Patch___TF24__TF24_Env__stand_competition_history__set', PACKAGE = 'plant', obj_, value)) +} + +Patch___TF24__TF24_Env__stand_height_stage_history__get <- function(obj_) { + .Call('_plant_Patch___TF24__TF24_Env__stand_height_stage_history__get', PACKAGE = 'plant', obj_) +} + +Patch___TF24__TF24_Env__stand_height_stage_history__set <- function(obj_, value) { + invisible(.Call('_plant_Patch___TF24__TF24_Env__stand_height_stage_history__set', PACKAGE = 'plant', obj_, value)) +} + +Patch___TF24__TF24_Env__stand_competition_stage_history__get <- function(obj_) { + .Call('_plant_Patch___TF24__TF24_Env__stand_competition_stage_history__get', PACKAGE = 'plant', obj_) +} + +Patch___TF24__TF24_Env__stand_competition_stage_history__set <- function(obj_, value) { + invisible(.Call('_plant_Patch___TF24__TF24_Env__stand_competition_stage_history__set', PACKAGE = 'plant', obj_, value)) +} + +Patch___TF24__TF24_Env__stand_newnode_height_stage_history__get <- function(obj_) { + .Call('_plant_Patch___TF24__TF24_Env__stand_newnode_height_stage_history__get', PACKAGE = 'plant', obj_) +} + +Patch___TF24__TF24_Env__stand_newnode_height_stage_history__set <- function(obj_, value) { + invisible(.Call('_plant_Patch___TF24__TF24_Env__stand_newnode_height_stage_history__set', PACKAGE = 'plant', obj_, value)) +} + +Patch___TF24__TF24_Env__stand_newnode_competition_stage_history__get <- function(obj_) { + .Call('_plant_Patch___TF24__TF24_Env__stand_newnode_competition_stage_history__get', PACKAGE = 'plant', obj_) +} + +Patch___TF24__TF24_Env__stand_newnode_competition_stage_history__set <- function(obj_, value) { + invisible(.Call('_plant_Patch___TF24__TF24_Env__stand_newnode_competition_stage_history__set', PACKAGE = 'plant', obj_, value)) +} + +Patch___TF24__TF24_Env__stand_height_stage_history_all__get <- function(obj_) { + .Call('_plant_Patch___TF24__TF24_Env__stand_height_stage_history_all__get', PACKAGE = 'plant', obj_) +} + +Patch___TF24__TF24_Env__stand_height_stage_history_all__set <- function(obj_, value) { + invisible(.Call('_plant_Patch___TF24__TF24_Env__stand_height_stage_history_all__set', PACKAGE = 'plant', obj_, value)) +} + +Patch___TF24__TF24_Env__stand_competition_stage_history_all__get <- function(obj_) { + .Call('_plant_Patch___TF24__TF24_Env__stand_competition_stage_history_all__get', PACKAGE = 'plant', obj_) +} + +Patch___TF24__TF24_Env__stand_competition_stage_history_all__set <- function(obj_, value) { + invisible(.Call('_plant_Patch___TF24__TF24_Env__stand_competition_stage_history_all__set', PACKAGE = 'plant', obj_, value)) +} + +Patch___TF24__TF24_Env__stand_newnode_height_stage_history_all__get <- function(obj_) { + .Call('_plant_Patch___TF24__TF24_Env__stand_newnode_height_stage_history_all__get', PACKAGE = 'plant', obj_) +} + +Patch___TF24__TF24_Env__stand_newnode_height_stage_history_all__set <- function(obj_, value) { + invisible(.Call('_plant_Patch___TF24__TF24_Env__stand_newnode_height_stage_history_all__set', PACKAGE = 'plant', obj_, value)) +} + +Patch___TF24__TF24_Env__stand_newnode_competition_stage_history_all__get <- function(obj_) { + .Call('_plant_Patch___TF24__TF24_Env__stand_newnode_competition_stage_history_all__get', PACKAGE = 'plant', obj_) +} + +Patch___TF24__TF24_Env__stand_newnode_competition_stage_history_all__set <- function(obj_, value) { + invisible(.Call('_plant_Patch___TF24__TF24_Env__stand_newnode_competition_stage_history_all__set', PACKAGE = 'plant', obj_, value)) +} + Patch___TF24f__TF24_Env__ctor <- function(parameters, environment, control) { .Call('_plant_Patch___TF24f__TF24_Env__ctor', PACKAGE = 'plant', parameters, environment, control) } @@ -2377,6 +2629,102 @@ Patch___TF24f__TF24_Env__state__get <- function(obj_) { .Call('_plant_Patch___TF24f__TF24_Env__state__get', PACKAGE = 'plant', obj_) } +Patch___TF24f__TF24_Env__step_history__get <- function(obj_) { + .Call('_plant_Patch___TF24f__TF24_Env__step_history__get', PACKAGE = 'plant', obj_) +} + +Patch___TF24f__TF24_Env__step_history__set <- function(obj_, value) { + invisible(.Call('_plant_Patch___TF24f__TF24_Env__step_history__set', PACKAGE = 'plant', obj_, value)) +} + +Patch___TF24f__TF24_Env__environment_history__get <- function(obj_) { + .Call('_plant_Patch___TF24f__TF24_Env__environment_history__get', PACKAGE = 'plant', obj_) +} + +Patch___TF24f__TF24_Env__environment_history__set <- function(obj_, value) { + invisible(.Call('_plant_Patch___TF24f__TF24_Env__environment_history__set', PACKAGE = 'plant', obj_, value)) +} + +Patch___TF24f__TF24_Env__stand_height_history__get <- function(obj_) { + .Call('_plant_Patch___TF24f__TF24_Env__stand_height_history__get', PACKAGE = 'plant', obj_) +} + +Patch___TF24f__TF24_Env__stand_height_history__set <- function(obj_, value) { + invisible(.Call('_plant_Patch___TF24f__TF24_Env__stand_height_history__set', PACKAGE = 'plant', obj_, value)) +} + +Patch___TF24f__TF24_Env__stand_competition_history__get <- function(obj_) { + .Call('_plant_Patch___TF24f__TF24_Env__stand_competition_history__get', PACKAGE = 'plant', obj_) +} + +Patch___TF24f__TF24_Env__stand_competition_history__set <- function(obj_, value) { + invisible(.Call('_plant_Patch___TF24f__TF24_Env__stand_competition_history__set', PACKAGE = 'plant', obj_, value)) +} + +Patch___TF24f__TF24_Env__stand_height_stage_history__get <- function(obj_) { + .Call('_plant_Patch___TF24f__TF24_Env__stand_height_stage_history__get', PACKAGE = 'plant', obj_) +} + +Patch___TF24f__TF24_Env__stand_height_stage_history__set <- function(obj_, value) { + invisible(.Call('_plant_Patch___TF24f__TF24_Env__stand_height_stage_history__set', PACKAGE = 'plant', obj_, value)) +} + +Patch___TF24f__TF24_Env__stand_competition_stage_history__get <- function(obj_) { + .Call('_plant_Patch___TF24f__TF24_Env__stand_competition_stage_history__get', PACKAGE = 'plant', obj_) +} + +Patch___TF24f__TF24_Env__stand_competition_stage_history__set <- function(obj_, value) { + invisible(.Call('_plant_Patch___TF24f__TF24_Env__stand_competition_stage_history__set', PACKAGE = 'plant', obj_, value)) +} + +Patch___TF24f__TF24_Env__stand_newnode_height_stage_history__get <- function(obj_) { + .Call('_plant_Patch___TF24f__TF24_Env__stand_newnode_height_stage_history__get', PACKAGE = 'plant', obj_) +} + +Patch___TF24f__TF24_Env__stand_newnode_height_stage_history__set <- function(obj_, value) { + invisible(.Call('_plant_Patch___TF24f__TF24_Env__stand_newnode_height_stage_history__set', PACKAGE = 'plant', obj_, value)) +} + +Patch___TF24f__TF24_Env__stand_newnode_competition_stage_history__get <- function(obj_) { + .Call('_plant_Patch___TF24f__TF24_Env__stand_newnode_competition_stage_history__get', PACKAGE = 'plant', obj_) +} + +Patch___TF24f__TF24_Env__stand_newnode_competition_stage_history__set <- function(obj_, value) { + invisible(.Call('_plant_Patch___TF24f__TF24_Env__stand_newnode_competition_stage_history__set', PACKAGE = 'plant', obj_, value)) +} + +Patch___TF24f__TF24_Env__stand_height_stage_history_all__get <- function(obj_) { + .Call('_plant_Patch___TF24f__TF24_Env__stand_height_stage_history_all__get', PACKAGE = 'plant', obj_) +} + +Patch___TF24f__TF24_Env__stand_height_stage_history_all__set <- function(obj_, value) { + invisible(.Call('_plant_Patch___TF24f__TF24_Env__stand_height_stage_history_all__set', PACKAGE = 'plant', obj_, value)) +} + +Patch___TF24f__TF24_Env__stand_competition_stage_history_all__get <- function(obj_) { + .Call('_plant_Patch___TF24f__TF24_Env__stand_competition_stage_history_all__get', PACKAGE = 'plant', obj_) +} + +Patch___TF24f__TF24_Env__stand_competition_stage_history_all__set <- function(obj_, value) { + invisible(.Call('_plant_Patch___TF24f__TF24_Env__stand_competition_stage_history_all__set', PACKAGE = 'plant', obj_, value)) +} + +Patch___TF24f__TF24_Env__stand_newnode_height_stage_history_all__get <- function(obj_) { + .Call('_plant_Patch___TF24f__TF24_Env__stand_newnode_height_stage_history_all__get', PACKAGE = 'plant', obj_) +} + +Patch___TF24f__TF24_Env__stand_newnode_height_stage_history_all__set <- function(obj_, value) { + invisible(.Call('_plant_Patch___TF24f__TF24_Env__stand_newnode_height_stage_history_all__set', PACKAGE = 'plant', obj_, value)) +} + +Patch___TF24f__TF24_Env__stand_newnode_competition_stage_history_all__get <- function(obj_) { + .Call('_plant_Patch___TF24f__TF24_Env__stand_newnode_competition_stage_history_all__get', PACKAGE = 'plant', obj_) +} + +Patch___TF24f__TF24_Env__stand_newnode_competition_stage_history_all__set <- function(obj_, value) { + invisible(.Call('_plant_Patch___TF24f__TF24_Env__stand_newnode_competition_stage_history_all__set', PACKAGE = 'plant', obj_, value)) +} + Patch___K93__K93_Env__ctor <- function(parameters, environment, control) { .Call('_plant_Patch___K93__K93_Env__ctor', PACKAGE = 'plant', parameters, environment, control) } @@ -2493,6 +2841,102 @@ Patch___K93__K93_Env__state__get <- function(obj_) { .Call('_plant_Patch___K93__K93_Env__state__get', PACKAGE = 'plant', obj_) } +Patch___K93__K93_Env__step_history__get <- function(obj_) { + .Call('_plant_Patch___K93__K93_Env__step_history__get', PACKAGE = 'plant', obj_) +} + +Patch___K93__K93_Env__step_history__set <- function(obj_, value) { + invisible(.Call('_plant_Patch___K93__K93_Env__step_history__set', PACKAGE = 'plant', obj_, value)) +} + +Patch___K93__K93_Env__environment_history__get <- function(obj_) { + .Call('_plant_Patch___K93__K93_Env__environment_history__get', PACKAGE = 'plant', obj_) +} + +Patch___K93__K93_Env__environment_history__set <- function(obj_, value) { + invisible(.Call('_plant_Patch___K93__K93_Env__environment_history__set', PACKAGE = 'plant', obj_, value)) +} + +Patch___K93__K93_Env__stand_height_history__get <- function(obj_) { + .Call('_plant_Patch___K93__K93_Env__stand_height_history__get', PACKAGE = 'plant', obj_) +} + +Patch___K93__K93_Env__stand_height_history__set <- function(obj_, value) { + invisible(.Call('_plant_Patch___K93__K93_Env__stand_height_history__set', PACKAGE = 'plant', obj_, value)) +} + +Patch___K93__K93_Env__stand_competition_history__get <- function(obj_) { + .Call('_plant_Patch___K93__K93_Env__stand_competition_history__get', PACKAGE = 'plant', obj_) +} + +Patch___K93__K93_Env__stand_competition_history__set <- function(obj_, value) { + invisible(.Call('_plant_Patch___K93__K93_Env__stand_competition_history__set', PACKAGE = 'plant', obj_, value)) +} + +Patch___K93__K93_Env__stand_height_stage_history__get <- function(obj_) { + .Call('_plant_Patch___K93__K93_Env__stand_height_stage_history__get', PACKAGE = 'plant', obj_) +} + +Patch___K93__K93_Env__stand_height_stage_history__set <- function(obj_, value) { + invisible(.Call('_plant_Patch___K93__K93_Env__stand_height_stage_history__set', PACKAGE = 'plant', obj_, value)) +} + +Patch___K93__K93_Env__stand_competition_stage_history__get <- function(obj_) { + .Call('_plant_Patch___K93__K93_Env__stand_competition_stage_history__get', PACKAGE = 'plant', obj_) +} + +Patch___K93__K93_Env__stand_competition_stage_history__set <- function(obj_, value) { + invisible(.Call('_plant_Patch___K93__K93_Env__stand_competition_stage_history__set', PACKAGE = 'plant', obj_, value)) +} + +Patch___K93__K93_Env__stand_newnode_height_stage_history__get <- function(obj_) { + .Call('_plant_Patch___K93__K93_Env__stand_newnode_height_stage_history__get', PACKAGE = 'plant', obj_) +} + +Patch___K93__K93_Env__stand_newnode_height_stage_history__set <- function(obj_, value) { + invisible(.Call('_plant_Patch___K93__K93_Env__stand_newnode_height_stage_history__set', PACKAGE = 'plant', obj_, value)) +} + +Patch___K93__K93_Env__stand_newnode_competition_stage_history__get <- function(obj_) { + .Call('_plant_Patch___K93__K93_Env__stand_newnode_competition_stage_history__get', PACKAGE = 'plant', obj_) +} + +Patch___K93__K93_Env__stand_newnode_competition_stage_history__set <- function(obj_, value) { + invisible(.Call('_plant_Patch___K93__K93_Env__stand_newnode_competition_stage_history__set', PACKAGE = 'plant', obj_, value)) +} + +Patch___K93__K93_Env__stand_height_stage_history_all__get <- function(obj_) { + .Call('_plant_Patch___K93__K93_Env__stand_height_stage_history_all__get', PACKAGE = 'plant', obj_) +} + +Patch___K93__K93_Env__stand_height_stage_history_all__set <- function(obj_, value) { + invisible(.Call('_plant_Patch___K93__K93_Env__stand_height_stage_history_all__set', PACKAGE = 'plant', obj_, value)) +} + +Patch___K93__K93_Env__stand_competition_stage_history_all__get <- function(obj_) { + .Call('_plant_Patch___K93__K93_Env__stand_competition_stage_history_all__get', PACKAGE = 'plant', obj_) +} + +Patch___K93__K93_Env__stand_competition_stage_history_all__set <- function(obj_, value) { + invisible(.Call('_plant_Patch___K93__K93_Env__stand_competition_stage_history_all__set', PACKAGE = 'plant', obj_, value)) +} + +Patch___K93__K93_Env__stand_newnode_height_stage_history_all__get <- function(obj_) { + .Call('_plant_Patch___K93__K93_Env__stand_newnode_height_stage_history_all__get', PACKAGE = 'plant', obj_) +} + +Patch___K93__K93_Env__stand_newnode_height_stage_history_all__set <- function(obj_, value) { + invisible(.Call('_plant_Patch___K93__K93_Env__stand_newnode_height_stage_history_all__set', PACKAGE = 'plant', obj_, value)) +} + +Patch___K93__K93_Env__stand_newnode_competition_stage_history_all__get <- function(obj_) { + .Call('_plant_Patch___K93__K93_Env__stand_newnode_competition_stage_history_all__get', PACKAGE = 'plant', obj_) +} + +Patch___K93__K93_Env__stand_newnode_competition_stage_history_all__set <- function(obj_, value) { + invisible(.Call('_plant_Patch___K93__K93_Env__stand_newnode_competition_stage_history_all__set', PACKAGE = 'plant', obj_, value)) +} + SCM___FF16__FF16_Env__ctor <- function(parameters, environment, control) { .Call('_plant_SCM___FF16__FF16_Env__ctor', PACKAGE = 'plant', parameters, environment, control) } @@ -4153,6 +4597,54 @@ TF24f_Strategy__ctor <- function() { .Call('_plant_TF24f_Strategy__ctor', PACKAGE = 'plant') } +ff16_offspring_production_gradient_impl <- function(pp, eh_list, sh, birth, ppsurv, ppsab, tw, traits) { + .Call('_plant_ff16_offspring_production_gradient_impl', PACKAGE = 'plant', pp, eh_list, sh, birth, ppsurv, ppsab, tw, traits) +} + +ff16_stand_gradient_impl <- function(pp, eh_list, sh, birth, ppsurv, ppsab, tw, traits, metrics, birth_rate, feedback, sh_h_list, sh_c_list, patch_area, al1_base, al2_base) { + .Call('_plant_ff16_stand_gradient_impl', PACKAGE = 'plant', pp, eh_list, sh, birth, ppsurv, ppsab, tw, traits, metrics, birth_rate, feedback, sh_h_list, sh_c_list, patch_area, al1_base, al2_base) +} + +ff16_stand_gradient_native <- function(scm_, pp, species, traits, metrics, birth_rate, feedback, sh_h_list, sh_c_list, patch_area, al1_base, al2_base) { + .Call('_plant_ff16_stand_gradient_native', PACKAGE = 'plant', scm_, pp, species, traits, metrics, birth_rate, feedback, sh_h_list, sh_c_list, patch_area, al1_base, al2_base) +} + +ff16_coupled_metrics_impl <- function(pp, eh_list, sh, birth, ppsurv, ppsab, tw, metrics, birth_rate, nn_h_list, nn_c_list, patch_area, active_birthenv = TRUE, birth_rate0 = -1.0) { + .Call('_plant_ff16_coupled_metrics_impl', PACKAGE = 'plant', pp, eh_list, sh, birth, ppsurv, ppsab, tw, metrics, birth_rate, nn_h_list, nn_c_list, patch_area, active_birthenv, birth_rate0) +} + +ff16_coupled_gradient_native <- function(scm_, pp, species, traits, metrics, birth_rate, patch_area, active_birthenv = TRUE) { + .Call('_plant_ff16_coupled_gradient_native', PACKAGE = 'plant', scm_, pp, species, traits, metrics, birth_rate, patch_area, active_birthenv) +} + +ff16_birth_rate_gradient_native <- function(scm_, pp, species, metrics, birth_rate, patch_area, active_birthenv = TRUE) { + .Call('_plant_ff16_birth_rate_gradient_native', PACKAGE = 'plant', scm_, pp, species, metrics, birth_rate, patch_area, active_birthenv) +} + +ff16_birth_rate_gradient_ms_native <- function(scm_, pp_list, metrics, birth_rate, patch_area, target, active_birthenv = TRUE) { + .Call('_plant_ff16_birth_rate_gradient_ms_native', PACKAGE = 'plant', scm_, pp_list, metrics, birth_rate, patch_area, target, active_birthenv) +} + +ff16_coupled_metrics_ms_impl <- function(pp_list, eh_list, sh, birth_list, metrics, birth_rate, nn_h_list, nn_c_list, patch_area, active_birthenv = TRUE) { + .Call('_plant_ff16_coupled_metrics_ms_impl', PACKAGE = 'plant', pp_list, eh_list, sh, birth_list, metrics, birth_rate, nn_h_list, nn_c_list, patch_area, active_birthenv) +} + +ff16_coupled_gradient_ms_impl <- function(pp_list, eh_list, sh, birth_list, traits, metrics, birth_rate, nn_h_list, nn_c_list, patch_area, target, active_birthenv = TRUE) { + .Call('_plant_ff16_coupled_gradient_ms_impl', PACKAGE = 'plant', pp_list, eh_list, sh, birth_list, traits, metrics, birth_rate, nn_h_list, nn_c_list, patch_area, target, active_birthenv) +} + +ff16_coupled_gradient_ms_native <- function(scm_, pp_list, traits, metrics, birth_rate, patch_area, target, active_birthenv = TRUE) { + .Call('_plant_ff16_coupled_gradient_ms_native', PACKAGE = 'plant', scm_, pp_list, traits, metrics, birth_rate, patch_area, target, active_birthenv) +} + +ff16_state_jacobian_impl <- function(pp, eh_list, sh, birth, ppsurv, ppsab, tw, traits) { + .Call('_plant_ff16_state_jacobian_impl', PACKAGE = 'plant', pp, eh_list, sh, birth, ppsurv, ppsab, tw, traits) +} + +ff16_grow_to_size_gradient_impl <- function(pp, env, y0v, sh, targets, sidx, traits, active_h0) { + .Call('_plant_ff16_grow_to_size_gradient_impl', PACKAGE = 'plant', pp, env, y0v, sh, targets, sidx, traits, active_h0) +} + node_schedule_default__Parameters___FF16__FF16_Env <- function(p) { .Call('_plant_node_schedule_default__Parameters___FF16__FF16_Env', PACKAGE = 'plant', p) } @@ -4161,6 +4653,14 @@ make_node_schedule__Parameters___FF16__FF16_Env <- function(p) { .Call('_plant_make_node_schedule__Parameters___FF16__FF16_Env', PACKAGE = 'plant', p) } +ff16_fecundity_dt_grad_ap1 <- function(height, light_E) { + .Call('_plant_ff16_fecundity_dt_grad_ap1', PACKAGE = 'plant', height, light_E) +} + +ff16_crown_top_fecundity_dt <- function(height, light_E, a_p1) { + .Call('_plant_ff16_crown_top_fecundity_dt', PACKAGE = 'plant', height, light_E, a_p1) +} + test_gradient_fd1 <- function(f, x, dx, direction, fx = NA_real_) { .Call('_plant_test_gradient_fd1', PACKAGE = 'plant', f, x, dx, direction, fx) } @@ -4232,6 +4732,18 @@ TF24f_strategy_expand_allometry <- function(s, height, area_heartwood, mass_hear .Call('_plant_TF24f_strategy_expand_allometry', PACKAGE = 'plant', s, height, area_heartwood, mass_heartwood) } +tf24_offspring_production_gradient_impl <- function(pp, eh_list, sh, birth, ppsurv, ppsab, tw, traits) { + .Call('_plant_tf24_offspring_production_gradient_impl', PACKAGE = 'plant', pp, eh_list, sh, birth, ppsurv, ppsab, tw, traits) +} + +tf24_offspring_production_gradient_native <- function(scm_, pp, species, birth_rate, traits) { + .Call('_plant_tf24_offspring_production_gradient_native', PACKAGE = 'plant', scm_, pp, species, birth_rate, traits) +} + +tf24_state_jacobian_impl <- function(pp, eh_list, sh, birth, ppsurv, ppsab, tw, traits, birth_rate) { + .Call('_plant_tf24_state_jacobian_impl', PACKAGE = 'plant', pp, eh_list, sh, birth, ppsurv, ppsab, tw, traits, birth_rate) +} + node_schedule_default__Parameters___TF24__TF24_Env <- function(p) { .Call('_plant_node_schedule_default__Parameters___TF24__TF24_Env', PACKAGE = 'plant', p) } @@ -4240,6 +4752,58 @@ make_node_schedule__Parameters___TF24__TF24_Env <- function(p) { .Call('_plant_make_node_schedule__Parameters___TF24__TF24_Env', PACKAGE = 'plant', p) } +tf24_crown_centre_fecundity_dt <- function(height, light_E, vcmax_25) { + .Call('_plant_tf24_crown_centre_fecundity_dt', PACKAGE = 'plant', height, light_E, vcmax_25) +} + +tf24_fecundity_dt_grad_vcmax <- function(height, light_E) { + .Call('_plant_tf24_fecundity_dt_grad_vcmax', PACKAGE = 'plant', height, light_E) +} + +tf24f_census_recon_native <- function(scm_, pp, species, birth_rate, k_acclim, use_ad_gradient, metrics, exact_ad_gprime) { + .Call('_plant_tf24f_census_recon_native', PACKAGE = 'plant', scm_, pp, species, birth_rate, k_acclim, use_ad_gradient, metrics, exact_ad_gprime) +} + +tf24f_census_gradient_ad_native <- function(scm_, pp, species, birth_rate, k_acclim, use_ad_gradient, traits, metrics, trait_rel_step) { + .Call('_plant_tf24f_census_gradient_ad_native', PACKAGE = 'plant', scm_, pp, species, birth_rate, k_acclim, use_ad_gradient, traits, metrics, trait_rel_step) +} + +tf24f_grow_to_size_gradient_impl <- function(pp, env, y0v, sh, targets, sidx, traits, k_acclim, use_ad_gradient, shading, gss_tol, trait_rel_step) { + .Call('_plant_tf24f_grow_to_size_gradient_impl', PACKAGE = 'plant', pp, env, y0v, sh, targets, sidx, traits, k_acclim, use_ad_gradient, shading, gss_tol, trait_rel_step) +} + +tf24f_coupled_metrics_impl <- function(pp, eh_list, sh, birth, birth_rate, k_acclim, use_ad_gradient, metrics, nn_h_list, nn_c_list, patch_area) { + .Call('_plant_tf24f_coupled_metrics_impl', PACKAGE = 'plant', pp, eh_list, sh, birth, birth_rate, k_acclim, use_ad_gradient, metrics, nn_h_list, nn_c_list, patch_area) +} + +tf24f_coupled_gradient_core <- function(pp, EH, sh, birth, birth_rate, k_acclim, use_ad_gradient, traits, metrics, NNH, NNC, patch_area, trait_rel_step) { + .Call('_plant_tf24f_coupled_gradient_core', PACKAGE = 'plant', pp, EH, sh, birth, birth_rate, k_acclim, use_ad_gradient, traits, metrics, NNH, NNC, patch_area, trait_rel_step) +} + +tf24f_coupled_gradient_native <- function(scm_, pp, species, birth_rate, k_acclim, use_ad_gradient, traits, metrics, patch_area, trait_rel_step) { + .Call('_plant_tf24f_coupled_gradient_native', PACKAGE = 'plant', scm_, pp, species, birth_rate, k_acclim, use_ad_gradient, traits, metrics, patch_area, trait_rel_step) +} + +tf24f_offspring_gradient_core <- function(pp, EH, sh, birth, ppsurv, ppsab, tw, k_acclim, use_ad_gradient, traits, trait_rel_step) { + .Call('_plant_tf24f_offspring_gradient_core', PACKAGE = 'plant', pp, EH, sh, birth, ppsurv, ppsab, tw, k_acclim, use_ad_gradient, traits, trait_rel_step) +} + +tf24f_offspring_gradient_impl <- function(pp, eh_list, sh, birth, ppsurv, ppsab, tw, k_acclim, use_ad_gradient, traits, trait_rel_step) { + .Call('_plant_tf24f_offspring_gradient_impl', PACKAGE = 'plant', pp, eh_list, sh, birth, ppsurv, ppsab, tw, k_acclim, use_ad_gradient, traits, trait_rel_step) +} + +tf24f_offspring_gradient_native <- function(scm_, pp, species, birth_rate, k_acclim, use_ad_gradient, traits, trait_rel_step) { + .Call('_plant_tf24f_offspring_gradient_native', PACKAGE = 'plant', scm_, pp, species, birth_rate, k_acclim, use_ad_gradient, traits, trait_rel_step) +} + +tf24f_coupled_gradient_ms_impl <- function(pp_list, eh_list, sh, birth_list, birth_rate, k_acclim, use_ad_gradient, traits, metrics, nn_h_list, nn_c_list, patch_area, target, trait_rel_step, gate_only) { + .Call('_plant_tf24f_coupled_gradient_ms_impl', PACKAGE = 'plant', pp_list, eh_list, sh, birth_list, birth_rate, k_acclim, use_ad_gradient, traits, metrics, nn_h_list, nn_c_list, patch_area, target, trait_rel_step, gate_only) +} + +tf24f_coupled_gradient_ms_native <- function(scm_, pp_list, birth_rate, k_acclim, use_ad_gradient, traits, metrics, patch_area, target, trait_rel_step, gate_only) { + .Call('_plant_tf24f_coupled_gradient_ms_native', PACKAGE = 'plant', scm_, pp_list, birth_rate, k_acclim, use_ad_gradient, traits, metrics, patch_area, target, trait_rel_step, gate_only) +} + node_schedule_default__Parameters___TF24f__TF24_Env <- function(p) { .Call('_plant_node_schedule_default__Parameters___TF24f__TF24_Env', PACKAGE = 'plant', p) } diff --git a/R/RcppR6.R b/R/RcppR6.R index be1bf993..0bade48f 100644 --- a/R/RcppR6.R +++ b/R/RcppR6.R @@ -1,6 +1,6 @@ ## Generated by RcppR6: do not edit by hand ## Version: 0.2.4 -## Hash: d66c262a7ae26c3ad786e6ef58df2d17 +## Hash: f82c59d833dcf16cfb8f7fb43e927f7d ##' @importFrom Rcpp evalCpp ##' @importFrom R6 R6Class @@ -120,6 +120,39 @@ check_type <- function(type, valid) { dprofit_droot_collar_psi = function(opt_root_psi) { Leaf__dprofit_droot_collar_psi(self, opt_root_psi) }, + dprofit_dvcmax25 = function(opt_root_psi) { + Leaf__dprofit_dvcmax25(self, opt_root_psi) + }, + dprofit_dg1_TF24 = function(opt_root_psi) { + Leaf__dprofit_dg1_TF24(self, opt_root_psi) + }, + dprofit_dbeta2 = function(opt_root_psi) { + Leaf__dprofit_dbeta2(self, opt_root_psi) + }, + dprofit_dkmax = function(opt_root_psi) { + Leaf__dprofit_dkmax(self, opt_root_psi) + }, + dprofit_dEup = function(opt_root_psi) { + Leaf__dprofit_dEup(self, opt_root_psi) + }, + dprofit_db = function(opt_root_psi) { + Leaf__dprofit_db(self, opt_root_psi) + }, + dprofit_dc = function(opt_root_psi) { + Leaf__dprofit_dc(self, opt_root_psi) + }, + dprofit_djmax25 = function(opt_root_psi) { + Leaf__dprofit_djmax25(self, opt_root_psi) + }, + dprofit_da = function(opt_root_psi) { + Leaf__dprofit_da(self, opt_root_psi) + }, + dprofit_dcurv_elec = function(opt_root_psi) { + Leaf__dprofit_dcurv_elec(self, opt_root_psi) + }, + dprofit_dcurv_colim = function(opt_root_psi) { + Leaf__dprofit_dcurv_colim(self, opt_root_psi) + }, psi_stem_to_ci = function(psi_stem, psi_upstream) { Leaf__psi_stem_to_ci(self, psi_stem, psi_upstream) }, @@ -1133,6 +1166,9 @@ Individual <- function(T, E) { net_mass_production_dt = function(environment) { Individual___FF16__FF16_Env__net_mass_production_dt(self, environment) }, + growth_rate_gradient_exact = function(environment) { + Individual___FF16__FF16_Env__growth_rate_gradient_exact(self, environment) + }, reset_mortality = function() { Individual___FF16__FF16_Env__reset_mortality(self) }, @@ -1249,6 +1285,9 @@ Individual <- function(T, E) { net_mass_production_dt = function(environment) { Individual___TF24__TF24_Env__net_mass_production_dt(self, environment) }, + growth_rate_gradient_exact = function(environment) { + Individual___TF24__TF24_Env__growth_rate_gradient_exact(self, environment) + }, reset_mortality = function() { Individual___TF24__TF24_Env__reset_mortality(self) }, @@ -1365,6 +1404,9 @@ Individual <- function(T, E) { net_mass_production_dt = function(environment) { Individual___TF24f__TF24_Env__net_mass_production_dt(self, environment) }, + growth_rate_gradient_exact = function(environment) { + Individual___TF24f__TF24_Env__growth_rate_gradient_exact(self, environment) + }, reset_mortality = function() { Individual___TF24f__TF24_Env__reset_mortality(self) }, @@ -1481,6 +1523,9 @@ Individual <- function(T, E) { net_mass_production_dt = function(environment) { Individual___K93__K93_Env__net_mass_production_dt(self, environment) }, + growth_rate_gradient_exact = function(environment) { + Individual___K93__K93_Env__growth_rate_gradient_exact(self, environment) + }, reset_mortality = function() { Individual___K93__K93_Env__reset_mortality(self) }, @@ -2877,6 +2922,90 @@ Patch <- function(T, E) { } else { stop("Patch$state is read-only") } + }, + step_history = function(value) { + if (missing(value)) { + Patch___FF16__FF16_Env__step_history__get(self) + } else { + Patch___FF16__FF16_Env__step_history__set(self, value) + } + }, + environment_history = function(value) { + if (missing(value)) { + Patch___FF16__FF16_Env__environment_history__get(self) + } else { + Patch___FF16__FF16_Env__environment_history__set(self, value) + } + }, + stand_height_history = function(value) { + if (missing(value)) { + Patch___FF16__FF16_Env__stand_height_history__get(self) + } else { + Patch___FF16__FF16_Env__stand_height_history__set(self, value) + } + }, + stand_competition_history = function(value) { + if (missing(value)) { + Patch___FF16__FF16_Env__stand_competition_history__get(self) + } else { + Patch___FF16__FF16_Env__stand_competition_history__set(self, value) + } + }, + stand_height_stage_history = function(value) { + if (missing(value)) { + Patch___FF16__FF16_Env__stand_height_stage_history__get(self) + } else { + Patch___FF16__FF16_Env__stand_height_stage_history__set(self, value) + } + }, + stand_competition_stage_history = function(value) { + if (missing(value)) { + Patch___FF16__FF16_Env__stand_competition_stage_history__get(self) + } else { + Patch___FF16__FF16_Env__stand_competition_stage_history__set(self, value) + } + }, + stand_newnode_height_stage_history = function(value) { + if (missing(value)) { + Patch___FF16__FF16_Env__stand_newnode_height_stage_history__get(self) + } else { + Patch___FF16__FF16_Env__stand_newnode_height_stage_history__set(self, value) + } + }, + stand_newnode_competition_stage_history = function(value) { + if (missing(value)) { + Patch___FF16__FF16_Env__stand_newnode_competition_stage_history__get(self) + } else { + Patch___FF16__FF16_Env__stand_newnode_competition_stage_history__set(self, value) + } + }, + stand_height_stage_history_all = function(value) { + if (missing(value)) { + Patch___FF16__FF16_Env__stand_height_stage_history_all__get(self) + } else { + Patch___FF16__FF16_Env__stand_height_stage_history_all__set(self, value) + } + }, + stand_competition_stage_history_all = function(value) { + if (missing(value)) { + Patch___FF16__FF16_Env__stand_competition_stage_history_all__get(self) + } else { + Patch___FF16__FF16_Env__stand_competition_stage_history_all__set(self, value) + } + }, + stand_newnode_height_stage_history_all = function(value) { + if (missing(value)) { + Patch___FF16__FF16_Env__stand_newnode_height_stage_history_all__get(self) + } else { + Patch___FF16__FF16_Env__stand_newnode_height_stage_history_all__set(self, value) + } + }, + stand_newnode_competition_stage_history_all = function(value) { + if (missing(value)) { + Patch___FF16__FF16_Env__stand_newnode_competition_stage_history_all__get(self) + } else { + Patch___FF16__FF16_Env__stand_newnode_competition_stage_history_all__set(self, value) + } })) @@ -3033,6 +3162,90 @@ Patch <- function(T, E) { } else { stop("Patch$state is read-only") } + }, + step_history = function(value) { + if (missing(value)) { + Patch___TF24__TF24_Env__step_history__get(self) + } else { + Patch___TF24__TF24_Env__step_history__set(self, value) + } + }, + environment_history = function(value) { + if (missing(value)) { + Patch___TF24__TF24_Env__environment_history__get(self) + } else { + Patch___TF24__TF24_Env__environment_history__set(self, value) + } + }, + stand_height_history = function(value) { + if (missing(value)) { + Patch___TF24__TF24_Env__stand_height_history__get(self) + } else { + Patch___TF24__TF24_Env__stand_height_history__set(self, value) + } + }, + stand_competition_history = function(value) { + if (missing(value)) { + Patch___TF24__TF24_Env__stand_competition_history__get(self) + } else { + Patch___TF24__TF24_Env__stand_competition_history__set(self, value) + } + }, + stand_height_stage_history = function(value) { + if (missing(value)) { + Patch___TF24__TF24_Env__stand_height_stage_history__get(self) + } else { + Patch___TF24__TF24_Env__stand_height_stage_history__set(self, value) + } + }, + stand_competition_stage_history = function(value) { + if (missing(value)) { + Patch___TF24__TF24_Env__stand_competition_stage_history__get(self) + } else { + Patch___TF24__TF24_Env__stand_competition_stage_history__set(self, value) + } + }, + stand_newnode_height_stage_history = function(value) { + if (missing(value)) { + Patch___TF24__TF24_Env__stand_newnode_height_stage_history__get(self) + } else { + Patch___TF24__TF24_Env__stand_newnode_height_stage_history__set(self, value) + } + }, + stand_newnode_competition_stage_history = function(value) { + if (missing(value)) { + Patch___TF24__TF24_Env__stand_newnode_competition_stage_history__get(self) + } else { + Patch___TF24__TF24_Env__stand_newnode_competition_stage_history__set(self, value) + } + }, + stand_height_stage_history_all = function(value) { + if (missing(value)) { + Patch___TF24__TF24_Env__stand_height_stage_history_all__get(self) + } else { + Patch___TF24__TF24_Env__stand_height_stage_history_all__set(self, value) + } + }, + stand_competition_stage_history_all = function(value) { + if (missing(value)) { + Patch___TF24__TF24_Env__stand_competition_stage_history_all__get(self) + } else { + Patch___TF24__TF24_Env__stand_competition_stage_history_all__set(self, value) + } + }, + stand_newnode_height_stage_history_all = function(value) { + if (missing(value)) { + Patch___TF24__TF24_Env__stand_newnode_height_stage_history_all__get(self) + } else { + Patch___TF24__TF24_Env__stand_newnode_height_stage_history_all__set(self, value) + } + }, + stand_newnode_competition_stage_history_all = function(value) { + if (missing(value)) { + Patch___TF24__TF24_Env__stand_newnode_competition_stage_history_all__get(self) + } else { + Patch___TF24__TF24_Env__stand_newnode_competition_stage_history_all__set(self, value) + } })) @@ -3189,6 +3402,90 @@ Patch <- function(T, E) { } else { stop("Patch$state is read-only") } + }, + step_history = function(value) { + if (missing(value)) { + Patch___TF24f__TF24_Env__step_history__get(self) + } else { + Patch___TF24f__TF24_Env__step_history__set(self, value) + } + }, + environment_history = function(value) { + if (missing(value)) { + Patch___TF24f__TF24_Env__environment_history__get(self) + } else { + Patch___TF24f__TF24_Env__environment_history__set(self, value) + } + }, + stand_height_history = function(value) { + if (missing(value)) { + Patch___TF24f__TF24_Env__stand_height_history__get(self) + } else { + Patch___TF24f__TF24_Env__stand_height_history__set(self, value) + } + }, + stand_competition_history = function(value) { + if (missing(value)) { + Patch___TF24f__TF24_Env__stand_competition_history__get(self) + } else { + Patch___TF24f__TF24_Env__stand_competition_history__set(self, value) + } + }, + stand_height_stage_history = function(value) { + if (missing(value)) { + Patch___TF24f__TF24_Env__stand_height_stage_history__get(self) + } else { + Patch___TF24f__TF24_Env__stand_height_stage_history__set(self, value) + } + }, + stand_competition_stage_history = function(value) { + if (missing(value)) { + Patch___TF24f__TF24_Env__stand_competition_stage_history__get(self) + } else { + Patch___TF24f__TF24_Env__stand_competition_stage_history__set(self, value) + } + }, + stand_newnode_height_stage_history = function(value) { + if (missing(value)) { + Patch___TF24f__TF24_Env__stand_newnode_height_stage_history__get(self) + } else { + Patch___TF24f__TF24_Env__stand_newnode_height_stage_history__set(self, value) + } + }, + stand_newnode_competition_stage_history = function(value) { + if (missing(value)) { + Patch___TF24f__TF24_Env__stand_newnode_competition_stage_history__get(self) + } else { + Patch___TF24f__TF24_Env__stand_newnode_competition_stage_history__set(self, value) + } + }, + stand_height_stage_history_all = function(value) { + if (missing(value)) { + Patch___TF24f__TF24_Env__stand_height_stage_history_all__get(self) + } else { + Patch___TF24f__TF24_Env__stand_height_stage_history_all__set(self, value) + } + }, + stand_competition_stage_history_all = function(value) { + if (missing(value)) { + Patch___TF24f__TF24_Env__stand_competition_stage_history_all__get(self) + } else { + Patch___TF24f__TF24_Env__stand_competition_stage_history_all__set(self, value) + } + }, + stand_newnode_height_stage_history_all = function(value) { + if (missing(value)) { + Patch___TF24f__TF24_Env__stand_newnode_height_stage_history_all__get(self) + } else { + Patch___TF24f__TF24_Env__stand_newnode_height_stage_history_all__set(self, value) + } + }, + stand_newnode_competition_stage_history_all = function(value) { + if (missing(value)) { + Patch___TF24f__TF24_Env__stand_newnode_competition_stage_history_all__get(self) + } else { + Patch___TF24f__TF24_Env__stand_newnode_competition_stage_history_all__set(self, value) + } })) @@ -3345,6 +3642,90 @@ Patch <- function(T, E) { } else { stop("Patch$state is read-only") } + }, + step_history = function(value) { + if (missing(value)) { + Patch___K93__K93_Env__step_history__get(self) + } else { + Patch___K93__K93_Env__step_history__set(self, value) + } + }, + environment_history = function(value) { + if (missing(value)) { + Patch___K93__K93_Env__environment_history__get(self) + } else { + Patch___K93__K93_Env__environment_history__set(self, value) + } + }, + stand_height_history = function(value) { + if (missing(value)) { + Patch___K93__K93_Env__stand_height_history__get(self) + } else { + Patch___K93__K93_Env__stand_height_history__set(self, value) + } + }, + stand_competition_history = function(value) { + if (missing(value)) { + Patch___K93__K93_Env__stand_competition_history__get(self) + } else { + Patch___K93__K93_Env__stand_competition_history__set(self, value) + } + }, + stand_height_stage_history = function(value) { + if (missing(value)) { + Patch___K93__K93_Env__stand_height_stage_history__get(self) + } else { + Patch___K93__K93_Env__stand_height_stage_history__set(self, value) + } + }, + stand_competition_stage_history = function(value) { + if (missing(value)) { + Patch___K93__K93_Env__stand_competition_stage_history__get(self) + } else { + Patch___K93__K93_Env__stand_competition_stage_history__set(self, value) + } + }, + stand_newnode_height_stage_history = function(value) { + if (missing(value)) { + Patch___K93__K93_Env__stand_newnode_height_stage_history__get(self) + } else { + Patch___K93__K93_Env__stand_newnode_height_stage_history__set(self, value) + } + }, + stand_newnode_competition_stage_history = function(value) { + if (missing(value)) { + Patch___K93__K93_Env__stand_newnode_competition_stage_history__get(self) + } else { + Patch___K93__K93_Env__stand_newnode_competition_stage_history__set(self, value) + } + }, + stand_height_stage_history_all = function(value) { + if (missing(value)) { + Patch___K93__K93_Env__stand_height_stage_history_all__get(self) + } else { + Patch___K93__K93_Env__stand_height_stage_history_all__set(self, value) + } + }, + stand_competition_stage_history_all = function(value) { + if (missing(value)) { + Patch___K93__K93_Env__stand_competition_stage_history_all__get(self) + } else { + Patch___K93__K93_Env__stand_competition_stage_history_all__set(self, value) + } + }, + stand_newnode_height_stage_history_all = function(value) { + if (missing(value)) { + Patch___K93__K93_Env__stand_newnode_height_stage_history_all__get(self) + } else { + Patch___K93__K93_Env__stand_newnode_height_stage_history_all__set(self, value) + } + }, + stand_newnode_competition_stage_history_all = function(value) { + if (missing(value)) { + Patch___K93__K93_Env__stand_newnode_competition_stage_history_all__get(self) + } else { + Patch___K93__K93_Env__stand_newnode_competition_stage_history_all__set(self, value) + } })) SCM <- function(T, E) { diff --git a/R/emergent_gradient.R b/R/emergent_gradient.R new file mode 100644 index 00000000..d039868e --- /dev/null +++ b/R/emergent_gradient.R @@ -0,0 +1,617 @@ +##' Reverse-mode trait gradient of an SCM's emergent \code{offspring_production} +##' (#472 scope B). +##' +##' A convenience wrapper over \code{\link{stand_gradient}} for the single +##' \code{offspring_production} metric, returning a named vector instead of a +##' one-row Jacobian. It works for \strong{any} strategy that \code{stand_gradient} +##' supports (FF16, TF24, TF24f) -- the dispatch and the reverse-mode replay live in +##' \code{stand_gradient}; this just selects the metric and reshapes. It is the +##' rare-mutant / invasion-fitness (selection) gradient: the resident light is held +##' frozen, so for a multi-species stand it is species \code{s}'s gradient against +##' the shared frozen canopy of all species. +##' +##' @title Reverse-mode gradient of emergent offspring_production +##' @param scm An \code{SCM} run with \code{control(save_RK45_cache = TRUE)} (FF16, +##' TF24, or TF24f). The cached schedule + resident light are read from its patch +##' natively; the SCM is not re-run. +##' @param traits Character vector of trait (parameter) names to differentiate. +##' \code{NULL} (default) uses the strategy's full production-trait set. +##' @param species Integer index of the species (cohort family) to differentiate. +##' Default \code{1}. +##' @param birth_rate The (constant) birth-rate driver used in the run. Recovered as +##' \code{offspring_production / net_reproduction_ratio} by default; pass it for a +##' time-varying driver. +##' @return A named numeric vector of trait derivatives, with attribute +##' \code{"offspring_production"} (the value reconstructed by the replay, which +##' should match \code{scm$offspring_production[[species]]}). +##' @seealso \code{\link{stand_gradient}}, the canonical metrics x traits engine. +##' @export +offspring_production_gradient <- function(scm, traits = NULL, species = 1L, + birth_rate = NULL) { + g <- stand_gradient(scm, metrics = "offspring_production", traits = traits, + species = species, birth_rate = birth_rate) + out <- stats::setNames(as.numeric(g$jacobian["offspring_production", ]), + colnames(g$jacobian)) + attr(out, "offspring_production") <- unname(g$values[["offspring_production"]]) + out +} + +# The 28 production-relevant FF16 trait (parameter) names the emergent gradients +# differentiate by default. +ff16_default_traits <- function() { + c("lma","rho","theta","a_b1","a_r1","eta_c","a_p1","a_p2","r_l","r_s", + "r_b","r_r","k_l","k_b","k_s","k_r","a_bio","a_y","a_l1","a_l2", + "a_f1","a_f2","hmat","omega","a_f3","d_I","a_dG1","a_dG2") +} + +# Harvest a run-with-cache FF16 SCM into the frozen pieces the two-pass replay +# consumes (the ResidentHarvest seam): the step schedule, the per-RK-stage resident +# environment, and species `species`'s cohort family (birth steps, node-spacing +# trapezoid weights, per-stage patch survival, survival-at-birth). Shared by every +# FF16 stand-gradient entry point so a new SCM variant only has to emit these. +ff16_harvest <- function(scm, species = 1L, birth_rate = NULL) { + types <- extract_RcppR6_template_types(scm$parameters, "Parameters") + if (!identical(types[[1]], "FF16")) { + stop("FF16 stand gradients are implemented for the FF16 strategy only") + } + if (species < 1L || species > length(scm$patch$species)) { + stop("species index out of range: stand has ", length(scm$patch$species), + " species") + } + # Cache the patch once: `scm$patch` rebuilds the whole RcppR6 patch object (every + # node + species) on each access, so repeatedly indexing it -- especially the + # per-stage pr_survival loop below -- is O(stand size) per call (~1600x slower). + patch <- scm$patch + sh <- patch$step_history + eh <- patch$environment_history + if (length(eh) < 1L) { + stop("No resident schedule cached: run the SCM with control(save_RK45_cache = TRUE)") + } + sp <- patch$species[[species]] + nt <- sp$node_times + pdens <- sp$patch_densities + ppsab <- sp$pr_patch_survival_at_birth + pp <- unlist(scm$parameters$strategies[[species]]$pars) + + if (is.null(birth_rate)) { + # Constant birth rate: offspring_production = birth_rate * net_reproduction_ratio. + birth_rate <- scm$offspring_production[[species]] / scm$net_reproduction_ratios[[species]] + } + + # Cohort birth steps (introductions land exactly on step times). + birth_step <- vapply(nt, function(t) which.min(abs(sh - t)) - 1L, integer(1)) + N <- length(eh) + # Node-spacing trapezoid weights so offspring_production == sum_i tw_i * offspring_i. + tcoef <- numeric(length(nt)); x <- nt; n <- length(x) + tcoef[1] <- 0.5 * (x[2] - x[1]); tcoef[n] <- 0.5 * (x[n] - x[n - 1]) + if (n > 2) tcoef[2:(n - 1)] <- 0.5 * (x[3:n] - x[1:(n - 2)]) + tw <- tcoef * pdens * pp[["S_D"]] * birth_rate + # pr_patch_survival at the exact Cash-Karp stage times sh[k] + ah[s]*h. + ah <- c(0, 0.2, 0.3, 0.6, 1.0, 0.875); hN <- diff(sh) + ppsurv <- matrix(0, N, 6) + for (k in seq_len(N)) for (s in 1:6) { + ppsurv[k, s] <- patch$pr_survival(sh[k] + ah[s] * hN[k]) + } + + # Per-RK-stage resident stand (species 0) for the RESIDENT feedback path (R0): + # [step][stage 0..5][cohort] heights + per-node competition effects, aligned 1:1 + # with environment_history. Empty unless the run cached them (older caches lack it). + sh_h <- patch$stand_height_stage_history + sh_c <- patch$stand_competition_stage_history + # Boundary new_node (height + competition effect) per RK stage: the trapezium tail + # term Species::compute_competition adds beyond `nodes`, needed by the COUPLED + # resident replay's per-stage canopy reconstruction at ground level. + nn_h <- patch$stand_newnode_height_stage_history + nn_c <- patch$stand_newnode_competition_stage_history + patch_area <- scm$parameters$patch_area + + list(pp = pp, eh = eh, sh = sh, birth_step = birth_step, ppsurv = ppsurv, + ppsab = ppsab, tw = tw, pdens = pdens, nt = nt, birth_rate = birth_rate, + sh_h = sh_h, sh_c = sh_c, nn_h = nn_h, nn_c = nn_c, patch_area = patch_area) +} + +# Harvest a run-with-cache FF16 SCM into the ALL-SPECIES arrays the multi-species +# coupled engine consumes (the cross-species resident Jacobian, #472 scope B R2): the +# shared schedule + joint env, per-species parameter vectors / cohort birth steps / +# constant birth rates, and the all-species per-RK-stage boundary harvest +# (stand_newnode_*_stage_history_all, [step][stage][species]). The joint stand light is +# reconstructed in C++ from each species' re-evolved cohorts; only these per-species +# pieces are needed from R. +ff16_harvest_ms <- function(scm) { + patch <- scm$patch + nsp <- length(scm$parameters$strategies) + sh <- patch$step_history + eh <- patch$environment_history + if (length(eh) < 1L) { + stop("No resident schedule cached: run the SCM with control(save_RK45_cache = TRUE)") + } + nn_h <- patch$stand_newnode_height_stage_history_all + nn_c <- patch$stand_newnode_competition_stage_history_all + if (length(nn_h) < 1L) { + stop("feedback = 'resident' on a multi-species stand needs the all-species ", + "per-RK-stage harvest; re-run the resident SCM on this plant version with ", + "control(save_RK45_cache = TRUE)") + } + list( + pp_list = lapply(seq_len(nsp), function(s) + unlist(scm$parameters$strategies[[s]]$pars)), + eh = eh, sh = sh, + birth_list = lapply(seq_len(nsp), function(s) + vapply(patch$species[[s]]$node_times, + function(t) which.min(abs(sh - t)) - 1L, integer(1))), + birth_rate = vapply(seq_len(nsp), function(s) + scm$offspring_production[[s]] / scm$net_reproduction_ratios[[s]], numeric(1)), + nn_h = nn_h, nn_c = nn_c, + patch_area = scm$parameters$patch_area, nsp = nsp) +} + +##' Reverse-mode derivative of an SCM's emergent census metrics w.r.t. the per-species +##' \emph{birth-rate driver} (#472 scope B, the birth-rate gradient AXIS). +##' +##' This differentiates a different axis than \code{\link{stand_gradient}}: not a trait, +##' but the constant birth-rate driver of a species. The use case is evolving a community +##' toward \strong{demographic equilibrium} -- birth_rate scales every cohort's +##' establishment density, so it re-shades the whole resident canopy, and its derivative +##' gives the Newton/gradient step that drives an emergent census target (e.g. a target +##' LAI / biomass, or self-replacement). +##' +##' Only the \strong{resident-coupled} reading needs a tape: with the canopy held frozen +##' every cohort's density is exactly linear in birth_rate, so the (rare-mutant) frozen +##' derivative is the trivial identity \eqn{\partial(\mathrm{metric})/\partial b = +##' \mathrm{metric}/b} (and \code{offspring_production} keeps that frozen reading even +##' under resident feedback -- the canopy a rare mutant invades is the resident's). So +##' this returns the resident TOTAL \eqn{\partial(\mathrm{census\ metric})/\partial b}: +##' birth_rate is registered as the sole input on the same coupled whole-stand replay the +##' resident trait gradient uses (one reverse sweep per metric). +##' +##' On a \strong{multi-species} stand this returns the CROSS-SPECIES total +##' \eqn{\partial(\mathrm{total\text{-}stand\ metric})/\partial b_s} for the chosen +##' \code{species} \eqn{s}: species \eqn{s}'s recruitment density re-shades the joint +##' canopy every species reads (the cross term the frozen reading sets to zero). It is +##' well-conditioned on a fixed node schedule and gated (with a clear error) for the stiff +##' refined-schedule case, like the cross-species trait gradient. +##' +##' FF16 only so far (its closed-form net keeps the coupled replay robust to long +##' horizons); TF24f's coupled replay is stiffness-gated (see +##' \code{stand_gradient(feedback = "resident")}). +##' +##' @title Resident-coupled birth-rate gradient of emergent census metrics (FF16) +##' @param scm An \code{SCM} run with \code{save_RK45_cache = TRUE} (FF16 strategy). +##' @param metrics Character vector, any of the census metrics \code{"LAI"}, +##' \code{"biomass"}, \code{"size_moment"} and -- single-species only -- +##' \code{"net_reproduction_ratio"} (the demographic-equilibrium / \eqn{R_0 = 1} axis: +##' \eqn{\partial R_0/\partial b}, with unit per-seed weights so \code{birth_rate} enters +##' only through the canopy -- the pure density feedback, i.e. the change in a +##' trait-identical mutant's fitness as the resident density moves; negative, which is +##' what makes an \eqn{R_0 = 1} equilibrium well-posed). (\code{offspring_production} is +##' excluded -- its birth-rate derivative is the trivial frozen identity above.) +##' @param species Integer species index; on a multi-species stand this selects which +##' species' birth-rate driver to differentiate (the cross-species total). Default \code{1}. +##' @param birth_rate The birth-rate driver to differentiate at; recovered from the run +##' by default (single-species only; the multi-species path recovers per species). +##' @return A list with \code{d_birth_rate} (named \code{d(metric)/d(birth_rate)}) and the +##' reconstructed \code{values} (plus \code{env_err}, the joint-canopy drift gauge, on a +##' multi-species stand). +##' @seealso \code{\link{stand_gradient}}, \code{\link{offspring_production_gradient}}. +##' @export +birth_rate_gradient <- function(scm, metrics = c("LAI", "biomass", "size_moment"), + species = 1L, birth_rate = NULL) { + strat <- extract_RcppR6_template_types(scm$parameters, "Parameters")[[1]] + if (!identical(strat, "FF16")) + stop("birth_rate_gradient is implemented for the FF16 strategy only so far (its ", + "closed-form net keeps the resident coupled replay robust); got ", strat) + bad <- setdiff(metrics, c("LAI", "biomass", "size_moment", "net_reproduction_ratio")) + if (length(bad)) + stop("birth_rate_gradient supports the census metrics (LAI, biomass, size_moment) ", + "and net_reproduction_ratio (the demographic-equilibrium / R0 = 1 axis). ", + "offspring_production's birth-rate derivative is the trivial frozen identity ", + "offspring_production / birth_rate. Got: ", paste(bad, collapse = ", ")) + nsp <- length(scm$parameters$strategies) + if (species < 1L || species > nsp) stop("species index out of range") + if ("net_reproduction_ratio" %in% metrics && nsp > 1L) + stop("net_reproduction_ratio is the single-species demographic-equilibrium axis; ", + "the cross-species birth_rate gradient supports the census metrics only") + patch_area <- scm$parameters$patch_area + if (nsp == 1L) { + pp <- unlist(scm$parameters$strategies[[species]]$pars) + return(ff16_birth_rate_gradient_native(scm, pp, as.integer(species - 1L), metrics, + if (is.null(birth_rate)) -1 else birth_rate, patch_area)) + } + # Multi-species: the CROSS-SPECIES total d(total-stand metric)/d(birth_rate of `species`) + # -- the joint canopy re-evolved from every species' cohorts, only the target species' + # birth_rate active (the cross term the frozen reading sets to zero). Well-conditioned on + # a fixed node schedule; gated for the stiff refined-schedule case like the MS trait path. + pp_list <- lapply(seq_len(nsp), function(s) unlist(scm$parameters$strategies[[s]]$pars)) + br_vec <- vapply(seq_len(nsp), function(s) + scm$offspring_production[[s]] / scm$net_reproduction_ratios[[s]], numeric(1)) + g <- ff16_birth_rate_gradient_ms_native(scm, pp_list, metrics, br_vec, patch_area, + as.integer(species)) + if (!all(is.finite(g$d_birth_rate)) || g$env_err > 1e-2) { + stop("the multi-species coupled resident re-evolution diverged on this node schedule ", + "(joint env drift = ", signif(g$env_err, 3), "); the cross-species birth_rate ", + "gradient needs a well-conditioned (fixed/uniform) node schedule -- re-run the ", + "resident SCM with refine_schedule = FALSE.") + } + g +} + +##' Reverse-mode trait gradient of an SCM's emergent stand metrics (#472 scope B, +##' the calibration-facing generic engine, FF16). +##' +##' The generic counterpart of \code{\link{offspring_production_gradient}}: given a +##' resident \code{SCM} run with \code{control(save_RK45_cache = TRUE)}, it returns a +##' \strong{metrics x traits Jacobian} \eqn{d(\mathrm{metric}_m)/d(\theta_k)} for a +##' set of emergent stand metrics, computed from ONE resident baseline. Every metric +##' is a weighted reduction over the replayed cohorts, \eqn{\mathrm{metric} = \sum_i +##' w_i\, f(\mathrm{state}_i)}; \code{offspring_production} is just one such entry -- +##' none is privileged. The engine records one forward replay onto a single adjoint +##' tape and takes one cheap reverse sweep \emph{per metric}, so M metrics cost a +##' replay plus M sweeps, not M replays. This is the calibration core: \code{plant} +##' returns the Jacobian; a downstream package composes likelihoods (data never +##' enters here), which is what lets many likelihood terms share one stand baseline. +##' +##' @title Reverse-mode Jacobian of emergent stand metrics (FF16) +##' @param scm An \code{SCM} run with \code{save_RK45_cache = TRUE} (FF16 strategy). +##' @param metrics Character vector of stand-metric names, any of +##' \code{"offspring_production"} (the seed-rain integral), \code{"LAI"} (leaf-area +##' index = the SCM's \code{compute_competition(0)}), \code{"biomass"} (the +##' size-distribution integral of live + heartwood mass) and \code{"size_moment"} +##' (the first moment of the size distribution, \eqn{\int n(h)\,h\,dh}). All are +##' symmetric reductions over the replayed cohorts; none is privileged. +##' @param traits Character vector of FF16 trait names. \code{NULL} (default) uses +##' all 28 production-relevant parameters. +##' @param species Integer index of the species (cohort family); see +##' \code{\link{offspring_production_gradient}}. Default \code{1}. +##' @param birth_rate The (constant) birth-rate driver; recovered from the run by +##' default. +##' @param feedback How the resident light responds to the trait. \code{"frozen"} +##' (default) holds the canopy fixed -- the rare-mutant / invasion-fitness gradient, +##' correct for \code{offspring_production}. \code{"resident"} is the resident TOTAL +##' gradient via the COUPLED whole-stand replay: all cohorts are re-evolved together +##' over the frozen schedule and the canopy light is reconstructed each RK stage from +##' the active stand (cohort heights AND densities respond to the trait, odelia #32 +##' active-knot spline). So EVERY trait feeds back -- a trait that changes growth or +##' mortality moves the canopy everyone reads -- and the feedback routinely dominates +##' and flips the sign of the census metrics (LAI / biomass / size-moment) relative +##' to the frozen reading. Applies to the census metrics; \code{offspring_production} +##' stays \code{"frozen"} (the invasion gradient) even under \code{"resident"}. FF16 +##' only so far (#472 scope B). \strong{Recommended use} (scope decision): for the +##' census metrics (LAI / biomass / size_moment) \code{"resident"} is the correct +##' stand-level total gradient and \code{"frozen"} the rare-mutant invasion gradient; +##' for \code{offspring_production} the two coincide (it is always the invasion +##' gradient). The default is left at \code{"frozen"} so the function is a safe, +##' backward-compatible no-feedback derivative unless feedback is asked for. +##' \strong{Multi-species:} on a stand with more than one species \code{"resident"} +##' returns the CROSS-SPECIES total gradient \eqn{d(\mathrm{total\ stand\ metric})/ +##' d(\theta_{\mathrm{species}})} -- the joint canopy is rebuilt from every species' +##' re-evolved cohorts and the census metrics are summed over species, so the +##' differentiated species' traits feed back through the canopy that every species +##' reads (the cross term \code{"frozen"} sets to zero). The joint re-evolution is +##' well-conditioned on a \emph{fixed} node schedule; on a strongly clustered +##' adaptively-refined schedule it can go stiff, and the function then raises a clear +##' error (gated by a cheap baseline pass) asking for a fixed/uniform schedule rather +##' than returning a diverged result. +##' @return A list with \code{$jacobian} (a metrics x traits matrix) and +##' \code{$values} (the reconstructed metric values, which should match the SCM's +##' emergent outputs). +##' @details Works for FF16, TF24 and TF24f residents (dispatched on the strategy). A +##' TF24/TF24f resident must have been run with \code{shading_model = "crown-centre"}. +##' FF16 supports all metrics (frozen + resident). TF24 currently supports +##' \code{"offspring_production"} (its census metrics need a leaf-optimisation +##' cross-sensitivity that is a follow-up). TF24f -- the fast-acclimation variant +##' whose tracked-collar leaf eval is analytic -- supports the CENSUS metrics +##' (\code{"LAI"} / \code{"biomass"} / \code{"size_moment"}) under \code{feedback = +##' "frozen"} via the reverse-mode AD census tape (the collar is carried as a taped +##' state with a curvature-linearised gradient-ascent rate); its resident (coupled) +##' census gradient and an offspring tape are follow-ups. \code{\link{stand_state_jacobian}} +##' works for FF16 and TF24. +##' @seealso \code{\link{offspring_production_gradient}}, \code{\link{birth_rate_gradient}}. +##' @export +stand_gradient <- function(scm, metrics = "offspring_production", traits = NULL, + species = 1L, birth_rate = NULL, + feedback = c("frozen", "resident")) { + # "resident_noanchor" is an undocumented validation mode (genuine recon value, for + # FD-checking R1); the public choices are "frozen"/"resident". + feedback <- if (length(feedback) > 1L) match.arg(feedback) else + match.arg(feedback, c("frozen", "resident", "resident_noanchor")) + is_resident <- feedback %in% c("resident", "resident_noanchor") + strat <- extract_RcppR6_template_types(scm$parameters, "Parameters")[[1]] + if (identical(strat, "FF16")) { + if (is.null(traits)) traits <- ff16_default_traits() + if (feedback == "frozen") { + # Frozen invasion gradient -> the FULLY native entry: the whole harvest (env + + # schedule + birth steps + weights + per-stage survival) is built in C++ from the + # live Patch, so the O(stand) R-side ff16_harvest (pr_survival loop) never runs. + pp <- unlist(scm$parameters$strategies[[species]]$pars) + return(ff16_stand_gradient_native(scm, pp, as.integer(species - 1L), traits, + metrics, if (is.null(birth_rate)) -1 else birth_rate, "frozen", + list(), list(), scm$parameters$patch_area, -1, -1)) + } + if (feedback != "resident") { + # The undocumented resident_noanchor validation mode -> the per-cohort frozen- + # canopy / leaf-area graft engine (needs the per-RK-stage stand harvest, so this + # rare path still uses the R harvest). + h <- ff16_harvest(scm, species, birth_rate) + if (is_resident && length(h$sh_h) < 1L) { + stop("feedback = 'resident' needs the per-RK-stage stand harvest; re-run the ", + "resident SCM on this plant version with control(save_RK45_cache = TRUE)") + } + return(ff16_stand_gradient_impl(h$pp, h$eh, h$sh, h$birth_step, h$ppsurv, + h$ppsab, h$tw, traits, metrics, h$birth_rate, feedback, + if (is.null(h$sh_h)) list() else h$sh_h, + if (is.null(h$sh_c)) list() else h$sh_c, h$patch_area, -1, -1)) + } + # feedback = "resident": the COUPLED whole-stand replay (every trait re-shades the + # canopy). Census metrics (LAI / biomass / size_moment) use the coupled total + # gradient; offspring_production stays the FROZEN invasion gradient (the canopy a + # rare mutant invades is the resident's, not co-moving with the mutant's trait). + # Fully native for single-species: cheap params from $parameters (no scm$patch + # rebuild, no ff16_harvest); the boundary-node guard lives in the C++ entry. + pp <- unlist(scm$parameters$strategies[[species]]$pars) + patch_area <- scm$parameters$patch_area + br <- if (is.null(birth_rate)) -1 else birth_rate + census_set <- c("LAI", "biomass", "size_moment") + bad <- setdiff(metrics, c(census_set, "offspring_production")) + if (length(bad)) stop("unknown stand metric: ", paste(bad, collapse = ", ")) + census <- metrics[metrics %in% census_set] + offsp <- metrics[metrics == "offspring_production"] + jac <- matrix(0, length(metrics), length(traits), + dimnames = list(metrics, traits)) + values <- stats::setNames(numeric(length(metrics)), metrics) + nsp <- length(scm$parameters$strategies) + if (length(census)) { + if (nsp == 1L) { + # Single-species coupled total gradient (canopy = this species' re-evolved stand). + # Fully native: env + harvest + boundary-node history from the live Patch. + gc <- ff16_coupled_gradient_native(scm, pp, as.integer(species - 1L), + traits, census, br, patch_area) + } else { + # Multi-species CROSS-SPECIES coupled Jacobian: the joint canopy is rebuilt from + # ALL species' re-evolved cohorts and the census metrics are TOTAL-stand sums, so + # this returns d(total-stand metric)/d(theta of species `species`) -- including + # the cross term whereby the differentiated species re-shades the canopy every + # species reads. The joint re-evolution is well-conditioned on a fixed node + # schedule but can go stiff (one cohort's log-density runs away) on a strongly + # clustered ADAPTIVELY-REFINED schedule; the native entry runs a cheap double R0 + # pass and returns env_err, which gates the result with a clear error (rather than + # a diverged Jacobian). Fully native: the joint env + schedule + birth steps + + # all-species boundary harvest come from the live Patch (cheap per-species pp + + # recovered birth rates from $parameters; no scm$patch rebuild, no Rcpp::as<> env). + pp_list <- lapply(seq_len(nsp), function(s) + unlist(scm$parameters$strategies[[s]]$pars)) + br_vec <- vapply(seq_len(nsp), function(s) + scm$offspring_production[[s]] / scm$net_reproduction_ratios[[s]], numeric(1)) + gc <- ff16_coupled_gradient_ms_native(scm, pp_list, traits, census, br_vec, + patch_area, as.integer(species)) + if (!all(is.finite(gc$values)) || gc$env_err > 1e-2) { + stop("the multi-species coupled resident re-evolution diverged on this node ", + "schedule (joint env drift = ", signif(gc$env_err, 3), "). The ", + "cross-species coupled gradient needs a well-conditioned node schedule; ", + "re-run the resident SCM with a fixed/uniform schedule (e.g. ", + "p$node_schedule_times <- list(seq(0, T, length.out = n), ...); ", + "run_scm(p, ..., refine_schedule = FALSE)) rather than an adaptively ", + "refined one.") + } + } + jac[census, ] <- gc$jacobian[census, , drop = FALSE] + values[census] <- gc$values[census] + } + if (length(offsp)) { + # offspring stays the FROZEN invasion gradient (fully native, like the frozen path). + go <- ff16_stand_gradient_native(scm, pp, as.integer(species - 1L), traits, + offsp, br, "frozen", list(), list(), patch_area, -1, -1) + jac[offsp, ] <- go$jacobian[offsp, , drop = FALSE] + values[offsp] <- go$values[offsp] + } + list(jacobian = jac, values = values) + } else if (identical(strat, "TF24")) { + if (is_resident) { + stop("feedback = 'resident' is implemented for FF16 only so far (R0-R1); ", + "TF24 resident light is R2") + } + # TF24 supports offspring_production only (census is a TF24 follow-up: its census + # number density needs a leaf-optimiser cross-sensitivity the linearised harvest + # does not differentiate faithfully -- TF24f, with its analytic tracked collar, + # gets census instead). Route to the native offspring entry (no R harvest). + bad <- setdiff(metrics, "offspring_production") + if (length(bad)) { + stop("stand_gradient for TF24 supports offspring_production only (census is a ", + "TF24 follow-up; use TF24f for census metrics). Got: ", + paste(bad, collapse = ", ")) + } + if (is.null(traits)) traits <- tf24_default_traits() + pp <- unlist(scm$parameters$strategies[[species]]$pars) + go <- tf24_offspring_production_gradient_native(scm, pp, as.integer(species - 1L), + if (is.null(birth_rate)) -1 else birth_rate, traits) + jac <- matrix(as.numeric(go), 1L, length(traits), + dimnames = list("offspring_production", traits)) + list(jacobian = jac, + values = stats::setNames(attr(go, "offspring_production"), "offspring_production")) + } else if (identical(strat, "TF24f")) { + # TF24f: the CENSUS metrics (LAI / biomass / size_moment) via the reverse-mode + # AD tape (#472 scope B). The tracked-collar leaf eval is analytic, so -- unlike + # TF24 -- the census number-density gradient is available (the collar is carried as + # a taped state with a curvature-linearised rate). feedback = "frozen" gives the + # rare-mutant / invasion census gradient (tf24f_census_gradient_ad, step 2); + # feedback = "resident" gives the coupled TOTAL stand gradient where every trait + # re-shades the canopy (tf24f_resident_census_gradient_ad, step 5). + # offspring_production is supported too (its own native tape, the frozen invasion + # gradient even under feedback = "resident"; wired below). + census_set <- c("LAI", "biomass", "size_moment") + bad <- setdiff(metrics, c(census_set, "offspring_production")) + if (length(bad)) { + stop("stand_gradient for TF24f supports offspring_production and the census ", + "metrics (", paste(census_set, collapse = ", "), "). Got: ", + paste(bad, collapse = ", ")) + } + if (is.null(traits)) traits <- tf24_default_traits() + census <- metrics[metrics %in% census_set] + offsp <- metrics[metrics == "offspring_production"] + jac <- matrix(0, length(metrics), length(traits), + dimnames = list(metrics, traits)) + values <- stats::setNames(numeric(length(metrics)), metrics) + if (length(census)) { + # feedback = "frozen" -> the rare-mutant / invasion census gradient (step 2); + # feedback = "resident" -> the coupled TOTAL stand gradient (step 5, single-species). + gc <- if (is_resident) { + if (length(scm$patch$species) > 1L) { + # Multi-species: the CROSS-SPECIES total gradient d(total-stand metric)/d(theta + # of species `species`) -- the joint canopy is rebuilt from every species' + # re-evolved cohorts, so the differentiated species feeds back through the canopy + # every species reads. Well-conditioned on a fixed node schedule (gated for stiff + # schedules inside the ms wrapper). + tf24f_resident_census_gradient_ms_ad(scm, metrics = census, traits = traits, + species = species) + } else { + tf24f_resident_census_gradient_ad(scm, metrics = census, traits = traits, + species = species, birth_rate = birth_rate) + } + } else { + tf24f_census_gradient_ad(scm, metrics = census, traits = traits, + species = species, birth_rate = birth_rate) + } + jac[census, ] <- gc$jacobian[census, , drop = FALSE] + values[census] <- gc$values[census] + } + if (length(offsp)) { + # offspring_production is always the FROZEN invasion gradient (the canopy a rare + # mutant invades is the resident's), even under feedback = "resident". + go <- tf24f_offspring_production_gradient(scm, traits = traits, species = species, + birth_rate = birth_rate) + jac["offspring_production", ] <- go$gradient + values["offspring_production"] <- go$value + } + list(jacobian = jac, values = values) + } else { + stop("stand_gradient is implemented for the FF16, TF24 and TF24f strategies only") + } +} + +##' Reverse-mode trait gradient of \code{\link{grow_individual_to_size}} (#472 scope +##' B, FF16): a single plant grown in a FIXED environment to target size(s), +##' differentiated w.r.t. traits. +##' +##' The per-plant, fixed-environment counterpart of \code{\link{stand_gradient}}: there +##' is no resident feedback (the environment is given), so the gradient is the exact +##' derivative of the grow-to-size solve. For each target size it returns the derivative +##' of the stopping TIME \eqn{t^*} (the time the plant reaches the target) and of every +##' ODE STATE component at \eqn{t^*}, w.r.t. all FF16 traits, as an optional output of +##' the same call you already make. This is the gradient a growth-rate or +##' time-to-size calibration / optimisation consumes (e.g. +##' \code{\link{optimise_individual_rate_at_size_by_trait}}). +##' +##' It is a two-pass method mirroring the SCM gradients. Pass 1 runs the ordinary +##' \code{\link{grow_individual_to_size}} (its adaptive Cash-Karp step schedule and +##' per-node trajectory are harvested). Pass 2 replays the demographic ODE over that +##' FROZEN schedule with the trait active, reading the fixed environment with the +##' default deep-crown assimilation, to a single partial final step landing on +##' \eqn{t^*}; one reverse sweep per state component gives \eqn{\partial +##' \mathrm{state}/\partial\theta} at fixed \eqn{t^*}. The stopping time itself responds +##' to the trait through the implicit function theorem on the stopping condition +##' \eqn{\mathrm{size}(t^*,\theta) = \mathrm{target}}: +##' \deqn{dt^*/d\theta = -\,(\partial\,\mathrm{size}/\partial\theta\,|_{t^*}) / +##' \dot{\mathrm{size}}(t^*),} +##' so the TOTAL derivative of each returned component \eqn{y_c} is \eqn{dy_c/d\theta = +##' \partial y_c/\partial\theta|_{t^*} + \dot y_c(t^*)\,dt^*/d\theta} (for the size +##' component itself the two terms cancel, as it is pinned to the target). The seedling +##' size \eqn{h_0} (which solves \eqn{\mathrm{mass}(h_0) = \mathrm{seed\ mass}}) carries +##' its own \eqn{dh_0/d\theta} by the same implicit-function step. +##' +##' @title Reverse-mode gradient of grow_individual_to_size (FF16 / TF24f) +##' @param individual An \code{Individual} object (FF16 or TF24f strategy), as passed to +##' \code{\link{grow_individual_to_size}}. For TF24f the trajectory carries the tracked +##' collar (\code{opt_root_psi_state}) as a 6th state and the gradient is dispatched to +##' the TF24f AD tape (same frozen-schedule replay + stopping-time IFT, with the collar +##' curvature-linearised); all other behaviour is identical. +##' @param sizes A vector of target sizes to grow the plant to (increasing). +##' @param size_name The size variable the targets refer to (one of the ODE state +##' names, e.g. \code{"height"}; FF16's monotonic size). +##' @param env An \code{Environment} object (the fixed environment). +##' @param traits Character vector of FF16 trait names; \code{NULL} (default) uses all +##' 28 production-relevant parameters. +##' @param time_max,warn Passed through to \code{\link{grow_individual_to_size}} for the +##' schedule-discovery pass. +##' @return A list with \code{$time} (the reconstructed \eqn{t^*} per size), +##' \code{$state} (a sizes x component matrix of the ODE state at \eqn{t^*}), +##' \code{$d_time} (a sizes x trait matrix \eqn{dt^*/d\theta}) and \code{$d_state} (a +##' sizes x component x trait array of TOTAL \eqn{d\,\mathrm{state}/d\theta}). The +##' reconstructed \code{$time}/\code{$state} match \code{grow_individual_to_size} to +##' the live \code{uniroot} tolerance. +##' @seealso \code{\link{grow_individual_to_size}}, \code{\link{stand_gradient}}. +##' @export +grow_individual_to_size_gradient <- function(individual, sizes, size_name, env, + traits = NULL, time_max = Inf, + warn = TRUE) { + # TF24f: the fast-acclimation variant carries the tracked collar as a 6th ODE state; + # its grow gradient reuses the same frozen-schedule replay + stopping-time IFT but with + # the curvature-linearised collar (build-order step 4). Dispatch to the TF24f AD tape. + if (grepl("^TF24f", individual$strategy_name)) { + return(tf24f_grow_individual_to_size_gradient_ad(individual, sizes, size_name, env, + traits, time_max, warn)) + } + if (!grepl("^FF16", individual$strategy_name)) + stop("grow_individual_to_size_gradient is implemented for the FF16 and TF24f ", + "strategies only") + if (is.unsorted(sizes) || length(sizes) == 0L) + stop("sizes must be non-empty and sorted") + sidx <- match(size_name, individual$ode_names) + if (is.na(sidx)) + stop("size_name must be one of the ODE state names: ", + paste(individual$ode_names, collapse = ", ")) + if (is.null(traits)) traits <- ff16_default_traits() + + # Pass 1: harvest the adaptive step schedule + initial state (the frozen schedule the + # replay reproduces). grow_individual_bracket returns the step times and per-node states. + brk <- grow_individual_bracket(individual, sizes, size_name, env, time_max, warn) + y0 <- stats::setNames(individual$ode_state, individual$ode_names) + pp <- unlist(individual$strategy$pars) + + res <- ff16_grow_to_size_gradient_impl(pp, env, y0, brk$time, as.numeric(sizes), + as.integer(sidx - 1L), traits, TRUE) + rownames(res$time) <- NULL + res$sizes <- sizes + res +} + +##' Per-cohort state x trait Jacobian of a resident SCM (#472 scope B, the +##' calibration-facing engine's escape hatch, FF16). +##' +##' The escape hatch for emergent metrics that are NOT a simple weighted reduction +##' (quantiles, ratios, bespoke statistics): rather than register a \code{(w, f)} +##' metric, this exposes \eqn{d(\mathrm{state}_{i,c})/d(\theta_k)} -- the derivative +##' of each replayed cohort's final demographic state component +##' \code{c} \eqn{\in} \{height, mortality, fecundity, area_heartwood, +##' mass_heartwood, offspring\} w.r.t. each trait. ANY smooth downstream metric over +##' the cohort states then has its gradient by the chain rule, with \code{plant} +##' never needing to know the metric -- the same boundary as "likelihoods live +##' downstream". Each cohort's final state is independent, so this tapes one cohort +##' at a time (one reverse sweep per state component). +##' +##' @title Per-cohort state x trait Jacobian (FF16) +##' @param scm An \code{SCM} run with \code{save_RK45_cache = TRUE} (FF16 strategy). +##' @param traits Character vector of FF16 trait names; \code{NULL} uses all 28. +##' @param species Integer species index (see \code{\link{stand_gradient}}). +##' @param birth_rate The (constant) birth-rate driver; recovered by default. +##' @return A list with \code{$states} (a cohort x component matrix of final-state +##' values) and \code{$jacobian} (a cohort x component x trait array). +##' @details Works for both FF16 and TF24 residents (dispatched on the strategy). A +##' TF24 resident must have been run with \code{shading_model = "crown-centre"}. +##' @seealso \code{\link{stand_gradient}}. +##' @export +stand_state_jacobian <- function(scm, traits = NULL, species = 1L, + birth_rate = NULL) { + strat <- extract_RcppR6_template_types(scm$parameters, "Parameters")[[1]] + if (identical(strat, "FF16")) { + if (is.null(traits)) traits <- ff16_default_traits() + h <- ff16_harvest(scm, species, birth_rate) + ff16_state_jacobian_impl(h$pp, h$eh, h$sh, h$birth_step, h$ppsurv, h$ppsab, h$tw, + traits) + } else if (identical(strat, "TF24")) { + if (is.null(traits)) traits <- tf24_default_traits() + h <- tf24_harvest(scm, species, birth_rate) + tf24_state_jacobian_impl(h$pp, h$eh, h$sh, h$birth_step, h$ppsurv, h$ppsab, h$tw, + traits, h$birth_rate) + } else { + stop("stand_state_jacobian is implemented for the FF16 and TF24 strategies only") + } +} diff --git a/R/tf24_emergent_gradient.R b/R/tf24_emergent_gradient.R new file mode 100644 index 00000000..c93a5d38 --- /dev/null +++ b/R/tf24_emergent_gradient.R @@ -0,0 +1,63 @@ +# TF24 emergent offspring_production gradient: there is no public tf24_*-named entry. +# offspring_production_gradient(scm) (R/emergent_gradient.R) dispatches on strategy and +# covers TF24 via stand_gradient()'s TF24 branch -> the native C++ entry +# tf24_offspring_production_gradient_native. (The single-strategy wrapper that used to +# live here was removed as redundant once stand_gradient unified the dispatch.) + +# The 27 net-production TF24 trait names (10 leaf + 17 mass-cascade) the emergent +# gradients differentiate by default. +tf24_default_traits <- function() { + c("vcmax_25","g1_TF24","beta2","K_s","b","c","jmax_25","a","curv_elec", + "curv_colim","lma","rho","a_b1","r_l","r_b","r_s","r_r","k_l","k_b", + "k_s","k_r","a_bio","a_y","a_l1","a_l2","theta","a_r1") +} + +# Harvest a run-with-cache TF24 SCM into the frozen pieces the two-pass replay +# consumes (the ResidentHarvest seam, TF24 mirror of ff16_harvest). Requires the +# resident to have been run with shading_model = "crown-centre" (the replay re-solves +# the crown-centre leaf optimisation). +tf24_harvest <- function(scm, species = 1L, birth_rate = NULL) { + types <- extract_RcppR6_template_types(scm$parameters, "Parameters") + if (!identical(types[[1]], "TF24")) { + stop("TF24 stand gradients are implemented for the TF24 strategy only") + } + if (species < 1L || species > length(scm$patch$species)) { + stop("species index out of range: stand has ", length(scm$patch$species), + " species") + } + # Cache the patch once: `scm$patch` rebuilds the whole RcppR6 patch object on each + # access, so the per-stage pr_survival loop below is O(stand size) per call otherwise. + patch <- scm$patch + sh <- patch$step_history + eh <- patch$environment_history + if (length(eh) < 1L) { + stop("No resident schedule cached: run the SCM with control(save_RK45_cache = TRUE)") + } + sp <- patch$species[[species]] + nt <- sp$node_times + pdens <- sp$patch_densities + ppsab <- sp$pr_patch_survival_at_birth + pp <- unlist(scm$parameters$strategies[[species]]$pars) + + if (is.null(birth_rate)) { + birth_rate <- scm$offspring_production[[species]] / scm$net_reproduction_ratios[[species]] + } + + # Cohort birth steps (introductions land on step times). + birth_step <- vapply(nt, function(t) which.min(abs(sh - t)) - 1L, integer(1)) + N <- length(eh) + # Node-spacing trapezoid weights so offspring_production == sum_i tw_i * offspring_i. + tcoef <- numeric(length(nt)); x <- nt; n <- length(x) + tcoef[1] <- 0.5 * (x[2] - x[1]); tcoef[n] <- 0.5 * (x[n] - x[n - 1]) + if (n > 2) tcoef[2:(n - 1)] <- 0.5 * (x[3:n] - x[1:(n - 2)]) + tw <- tcoef * pdens * pp[["S_D"]] * birth_rate + # pr_patch_survival at the exact Cash-Karp stage times sh[k] + ah[s]*h. + ah <- c(0, 0.2, 0.3, 0.6, 1.0, 0.875); hN <- diff(sh) + ppsurv <- matrix(0, N, 6) + for (k in seq_len(N)) for (s in 1:6) { + ppsurv[k, s] <- patch$pr_survival(sh[k] + ah[s] * hN[k]) + } + + list(pp = pp, eh = eh, sh = sh, birth_step = birth_step, ppsurv = ppsurv, + ppsab = ppsab, tw = tw, pdens = pdens, nt = nt, birth_rate = birth_rate) +} diff --git a/R/tf24f_emergent_gradient.R b/R/tf24f_emergent_gradient.R new file mode 100644 index 00000000..5adfdef4 --- /dev/null +++ b/R/tf24f_emergent_gradient.R @@ -0,0 +1,533 @@ +# TF24f stand-metric gradients (#472 scope B, build-order step 1). TF24f is the +# fast-acclimation TF24 variant whose optimal root-collar potential is a 6th ODE state +# tracked by gradient ascent (no per-step golden-section optimiser). TF24f -- not TF24 +# -- is the right target for the +# stand CENSUS gradients: its leaf evaluation at the tracked collar is analytic / IFT-able, +# so the census number density's growth-rate-gradient term needs no curvature harvest. The +# R0 GATE below is the first deliverable: a double-precision census reconstruction that +# proves the collar-state replay is faithful (the prerequisite §7 flags) before the +# reverse-mode tape (R1) is built. + +# Harvest a run-with-cache TF24f SCM into the pieces the census replay consumes (the +# TF24f mirror of tf24_harvest). Unlike TF24, no per-stage leaf-opt harvest is needed -- +# the replay evaluates the analytic leaf at the tracked collar live -- so this returns +# only the frozen schedule + env, the cohort birth steps, the parameter vector and the +# acclimation knobs (k_acclim / use_ad_gradient). Requires shading_model = "crown-centre". +tf24f_harvest <- function(scm, species = 1L, birth_rate = NULL) { + types <- extract_RcppR6_template_types(scm$parameters, "Parameters") + if (!identical(types[[1]], "TF24f")) { + stop("TF24f census gradients are implemented for the TF24f strategy only") + } + if (species < 1L || species > length(scm$patch$species)) { + stop("species index out of range: stand has ", length(scm$patch$species), + " species") + } + patch <- scm$patch + sh <- patch$step_history + eh <- patch$environment_history + if (length(eh) < 1L) { + stop("No resident schedule cached: run the SCM with control(save_RK45_cache = TRUE)") + } + sp <- patch$species[[species]] + nt <- sp$node_times + pdens <- sp$patch_densities + ppsab <- sp$pr_patch_survival_at_birth + strat <- scm$parameters$strategies[[species]] + pp <- unlist(strat$pars) + + if (is.null(birth_rate)) { + birth_rate <- scm$offspring_production[[species]] / scm$net_reproduction_ratios[[species]] + } + # Cohort birth steps (introductions land on step times). + birth_step <- vapply(nt, function(t) which.min(abs(sh - t)) - 1L, integer(1)) + N <- length(eh) + + # offspring_production weighting (the TF24f mirror of ff16_harvest / tf24_harvest): the + # node-spacing trapezoid weights so offspring_production == sum_i tw_i * offspring_i, + # and per-RK-stage patch survival at the exact Cash-Karp stage times. Only the offspring + # tape uses these (census / resident do not). + tcoef <- numeric(length(nt)); x <- nt; n <- length(x) + tcoef[1] <- 0.5 * (x[2] - x[1]); tcoef[n] <- 0.5 * (x[n] - x[n - 1]) + if (n > 2) tcoef[2:(n - 1)] <- 0.5 * (x[3:n] - x[1:(n - 2)]) + tw <- tcoef * pdens * pp[["S_D"]] * birth_rate + ah <- c(0, 0.2, 0.3, 0.6, 1.0, 0.875); hN <- diff(sh) + ppsurv <- matrix(0, N, 6) + for (k in seq_len(N)) for (s in 1:6) ppsurv[k, s] <- patch$pr_survival(sh[k] + ah[s] * hN[k]) + + # Boundary new_node (height + competition effect) per RK stage: the trapezium tail + # term Species::compute_competition adds beyond `nodes`, needed by the COUPLED + # resident replay's per-stage canopy reconstruction at ground level (the TF24f mirror + # of ff16_harvest's nn_h / nn_c). Empty on older caches. + nn_h <- patch$stand_newnode_height_stage_history + nn_c <- patch$stand_newnode_competition_stage_history + patch_area <- scm$parameters$patch_area + + list(pp = pp, eh = eh, sh = sh, birth_step = birth_step, birth_rate = birth_rate, + k_acclim = strat$k_acclim, use_ad_gradient = strat$use_ad_gradient, nt = nt, + tw = tw, ppsurv = ppsurv, ppsab = ppsab, + nn_h = nn_h, nn_c = nn_c, patch_area = patch_area) +} + +# R0 gate (internal): double-precision census reconstruction of a TF24f resident stand. +# Re-evolves every cohort's {5 demog, tracked collar, log_density} over the frozen +# schedule and returns the reconstructed per-cohort heights / collar / log-densities and +# the census metric values. Used to confirm the collar-state replay reproduces the SCM's +# stored stand (heights / log_densities / opt_root_psi_state) and that the LAI reduction +# matches compute_competition(0), before the reverse-mode tape (R1) is added. Not yet a +# public gradient entry point -- it returns the recon, not d(metric)/d(theta). +# Strategy guard for the native census entries (which take the live SCM via an +# RcppR6> cast). Mirrors the check tf24f_harvest used to provide, so a +# non-TF24f SCM still gets a clear message rather than a raw external-pointer cast error. +tf24f_require_strategy <- function(scm) { + types <- extract_RcppR6_template_types(scm$parameters, "Parameters") + if (!identical(types[[1]], "TF24f")) + stop("TF24f census gradients are implemented for the TF24f strategy only") +} + +tf24f_census_recon <- function(scm, metrics = c("LAI", "biomass", "size_moment"), + species = 1L, birth_rate = NULL) { + tf24f_require_strategy(scm) + # The reconstruction's growth-rate gradient g' must match whatever the resident + # SCM used (Node::growth_rate_gradient): the exact-AD path when the run set + # control(node_gradient_exact_ad = TRUE), else the backward finite difference. + exact_ad <- isTRUE(scm$parameters$strategies[[species]]$control$node_gradient_exact_ad) + # Fully native: env + birth steps from the live Patch (no Rcpp::as<> round-trip; the + # whole harvest is native, so the R-side tf24f_harvest -- whose ppsurv/tw loop is dead + # work for census -- is skipped). birth_rate < 0 recovers the rate natively. + strat <- scm$parameters$strategies[[species]] + tf24f_census_recon_native(scm, unlist(strat$pars), as.integer(species - 1L), + if (is.null(birth_rate)) -1 else birth_rate, + strat$k_acclim, strat$use_ad_gradient, metrics, exact_ad) +} + +# TF24f frozen census trait gradient (#472 scope B, build-order step 2 -- R1 GATE). +# d(census metric)/d(theta) for the FROZEN (rare-mutant / invasion) resident light, by a +# central finite difference over the R0 census reconstruction: for each trait, perturb the +# parameter vector by +/- a relative step, re-run the double-precision collar-state replay +# against the SAME frozen resident environment, and difference the metric. The resident +# schedule + per-RK-stage light are held fixed (the invasion gradient), so this is the +# faithful finite-difference reference the reverse-mode AD tape (the actual R1) must +# reproduce -- and a usable prototype gradient in its own right. Returns a metrics x traits +# Jacobian and the reconstructed metric values. (The AD tape will replace the per-trait +# replays with one reverse sweep; this gate fixes the target it must hit.) +tf24f_census_gradient_fd <- function(scm, metrics = c("LAI", "biomass", "size_moment"), + traits = NULL, species = 1L, birth_rate = NULL, + rel_step = 1e-5) { + if (is.null(traits)) traits <- tf24_default_traits() + h <- tf24f_harvest(scm, species, birth_rate) + exact_ad <- isTRUE(scm$parameters$strategies[[species]]$control$node_gradient_exact_ad) + # Native recon so the FD reference differentiates the SAME function the native AD + # tape does (faithful crown-sampled light, no Rcpp::as<> round-trip). birth steps are + # computed natively inside; h$birth_rate is the concrete recovered rate. + recon <- function(pp) tf24f_census_recon_native(scm, pp, as.integer(species - 1L), + h$birth_rate, h$k_acclim, h$use_ad_gradient, metrics, exact_ad)$values + + values <- recon(h$pp) + jac <- matrix(0, length(metrics), length(traits), + dimnames = list(metrics, traits)) + for (j in seq_along(traits)) { + tr <- traits[j] + if (!tr %in% names(h$pp)) stop("unknown TF24f trait: ", tr) + d <- rel_step * max(abs(h$pp[[tr]]), 1e-8) + pp_p <- h$pp; pp_p[[tr]] <- pp_p[[tr]] + d + pp_m <- h$pp; pp_m[[tr]] <- pp_m[[tr]] - d + jac[, j] <- (recon(pp_p) - recon(pp_m)) / (2 * d) + } + list(jacobian = jac, values = values) +} + +# TF24f frozen census trait gradient by reverse-mode AD (#472 scope B, build-order +# step 2 -- the R1 tape, the refine of tf24f_census_gradient_fd). One reverse sweep over +# the 7-state replay {5 demog, tracked collar, log_density} per metric, replacing the +# per-trait finite difference over the census reconstruction. The hard ingredient (the +# tracked collar is a strongly theta-dependent STATE that lags the optimum, so the +# envelope theorem does not zero it) is handled by a CURVATURE harvest: the collar is +# carried as a taped state with rate k_acclim * dprofit_dpsi, linearised with the second +# derivatives d2profit/dpsi2, d2profit/dpsi dh, d2profit/dpsi dtheta_k (all FD-harvested +# in the double discovery pass). The census g' = d(height_dt)/d(height) reproduces the +# SCM's backward-FD scheme by harvesting a second operating point at h - GEPS. Resident +# light is FROZEN (the rare-mutant / invasion gradient). Validated to ~1% (the recon +# noise floor) against tf24f_census_gradient_fd; see test-tf24f-census-gradient.R. +tf24f_census_gradient_ad <- function(scm, metrics = c("LAI", "biomass", "size_moment"), + traits = NULL, species = 1L, birth_rate = NULL, + trait_rel_step = 1e-5) { + tf24f_require_strategy(scm) + if (is.null(traits)) traits <- tf24_default_traits() + # Fully native: env + birth steps from the live Patch; tf24f_harvest skipped. + strat <- scm$parameters$strategies[[species]] + tf24f_census_gradient_ad_native(scm, unlist(strat$pars), as.integer(species - 1L), + if (is.null(birth_rate)) -1 else birth_rate, + strat$k_acclim, strat$use_ad_gradient, traits, metrics, + trait_rel_step) +} + +# TF24f individual grow-to-size trait gradient by reverse-mode AD (#472 scope B, +# build-order step 4 -- the AD refine of tf24f_grow_individual_to_size_gradient_fd). A +# single TF24f plant grown in a FIXED environment to target size(s); returns d(t*)/d(theta) +# and the TOTAL d(state at t*)/d(theta). No resident feedback / canopy / density, so the +# only machinery beyond the frozen-schedule replay is the stopping-time IFT. The tracked +# collar (opt_root_psi_state) is carried as the 6th replayed state with the SAME curvature- +# linearised gradient-ascent rate the census tape uses; it starts at the individual's birth +# value (theta-independent), so only the seedling height h0 carries an initial-condition +# derivative. Pass 1 (R) harvests the live grow's adaptive Cash-Karp schedule via +# grow_individual_bracket; pass 2 (C++) replays + sweeps. Validated against +# tf24f_grow_individual_to_size_gradient_fd; see test-tf24f-individual-gradient.R. +tf24f_grow_individual_to_size_gradient_ad <- function(individual, sizes, size_name, env, + traits = NULL, time_max = Inf, + warn = TRUE, trait_rel_step = 1e-5) { + if (!grepl("^TF24f", individual$strategy_name)) + stop("tf24f_grow_individual_to_size_gradient_ad is for the TF24f strategy only") + if (is.unsorted(sizes) || length(sizes) == 0L) + stop("sizes must be non-empty and sorted") + sidx <- match(size_name, individual$ode_names) + if (is.na(sidx)) + stop("size_name must be one of the ODE state names: ", + paste(individual$ode_names, collapse = ", ")) + if (is.null(traits)) traits <- tf24_default_traits() + + # Pass 1: harvest the adaptive step schedule (the frozen schedule the replay reproduces). + brk <- grow_individual_bracket(individual, sizes, size_name, env, time_max, warn) + y0 <- individual$ode_state + s <- individual$strategy + pp <- unlist(s$pars) + shading <- s$control$shading_model + if (is.null(shading) || !nzchar(shading)) shading <- "mean-light" # TF24's default + + res <- tf24f_grow_to_size_gradient_impl(pp, env, y0, brk$time, as.numeric(sizes), + as.integer(sidx - 1L), traits, s$k_acclim, + s$use_ad_gradient, shading, s$control$GSS_tol_abs, + trait_rel_step) + res$sizes <- sizes + res +} + +# TF24f individual grow-to-size trait gradient (#472 scope B, the "individuals" surface +# -- prototype). d(t*)/d(theta) and d(state at t*)/d(theta) for a single plant grown in a +# FIXED environment to target size(s), by central finite difference over +# grow_individual_to_size. No resident feedback (the env is given) and no canopy/density, +# so this is the lightest gradient surface; for TF24f the tracked collar is re-evolved +# inside each grow (it is one of the ODE states), so the FD captures the collar's response +# automatically -- which is exactly why an exact AD version is the heavier follow-up (the +# tracked collar is strongly theta-dependent). The +# FD here is the prototype + the reference that AD version must reproduce. Traits are +# perturbed on the (post-hyperpar) strategy parameters directly, matching the census FD. +tf24f_grow_individual_to_size_gradient_fd <- function(individual, sizes, size_name, env, + traits = NULL, time_max = Inf, + warn = FALSE, rel_step = 1e-5) { + if (!grepl("^TF24f", individual$strategy_name)) + stop("tf24f_grow_individual_to_size_gradient_fd is for the TF24f strategy only") + if (is.null(traits)) traits <- tf24_default_traits() + + grow <- function(ind) { + r <- grow_individual_to_size(ind, sizes, size_name, env, time_max, warn) + list(time = r$time, state = r$state) + } + perturb <- function(tr, delta) { + s <- individual$strategy + pars <- s$pars + if (!tr %in% names(pars)) stop("unknown TF24f trait: ", tr) + pars[[tr]] <- pars[[tr]] + delta + s$pars <- pars + grow(TF24f_Individual(s)) + } + + base <- grow(individual) + nS <- length(sizes); comp <- colnames(base$state); nC <- length(comp) + d_time <- matrix(0, nS, length(traits), dimnames = list(NULL, traits)) + d_state <- array(0, c(nS, nC, length(traits)), dimnames = list(NULL, comp, traits)) + s0 <- individual$strategy$pars + for (j in seq_along(traits)) { + tr <- traits[j] + d <- rel_step * max(abs(s0[[tr]]), 1e-8) + rp <- perturb(tr, d); rm <- perturb(tr, -d) + d_time[, j] <- (rp$time - rm$time) / (2 * d) + d_state[, , j] <- (rp$state - rm$state) / (2 * d) + } + list(sizes = sizes, time = base$time, state = base$state, + d_time = d_time, d_state = d_state) +} + +# TF24f offspring_production trait gradient by reverse-mode AD (#472 scope B, the +# offspring surface). d(offspring_production)/d(theta) for the seed-rain integral +# offspring_production = sum_i tw_i * offspring_i (survival-weighted lifetime fecundity), +# the FROZEN rare-mutant / invasion gradient. The TF24 offspring tape with the one TF24f +# difference the census tape already solves -- the tracked collar is carried as a taped +# state with the curvature-linearised gradient-ascent rate (the lag the envelope theorem +# does not zero), seeded at the birth optimum (IFT-injected). Returns a named gradient +# vector + the reconstructed value. +tf24f_offspring_production_gradient <- function(scm, traits = NULL, species = 1L, + birth_rate = NULL, trait_rel_step = 1e-5) { + tf24f_require_strategy(scm) + if (is.null(traits)) traits <- tf24_default_traits() + # Fully native (mirrors FF16's offspring path): env + offspring harvest from the + # live Patch; tf24f_harvest skipped. birth_rate < 0 recovers natively. + strat <- scm$parameters$strategies[[species]] + tf24f_offspring_gradient_native(scm, unlist(strat$pars), as.integer(species - 1L), + if (is.null(birth_rate)) -1 else birth_rate, + strat$k_acclim, strat$use_ad_gradient, traits, + trait_rel_step) +} + +# TF24f COUPLED census reconstruction -- the resident R0 gate (#472 scope B step 5, +# internal). A double-precision whole-stand re-evolution that reconstructs the ACTIVE +# canopy each RK stage and drives the real TF24f leaf at the reconstructed crown light. +# Returns the reconstructed TOTAL-stand census values and `env_err` (worst reconstructed +# vs SCM knot-light drift) -- the gauge that the coupled re-evolution reproduces the +# resident stand before the AD tape. Not a public gradient entry; used by the gate test. +tf24f_coupled_metrics <- function(scm, metrics = c("LAI", "size_moment"), + species = 1L, birth_rate = NULL) { + h <- tf24f_harvest(scm, species, birth_rate) + if (length(h$nn_h) < 1L) { + stop("the TF24f coupled (resident) census gradient needs the per-RK-stage ", + "boundary-node harvest; re-run the resident SCM with ", + "control(save_RK45_cache = TRUE)") + } + tf24f_coupled_metrics_impl(h$pp, h$eh, h$sh, h$birth_step, h$birth_rate, h$k_acclim, + h$use_ad_gradient, metrics, h$nn_h, h$nn_c, h$patch_area) +} + +# TF24f RESIDENT (coupled) census trait gradient by reverse-mode AD (#472 scope B, +# build-order step 5 -- the resident-feedback AD refine of tf24f_resident_census_gradient_fd). +# The TOTAL stand-level d(census metric)/d(theta): all cohorts are re-evolved TOGETHER over +# the frozen schedule and the canopy light each RK stage is reconstructed from the active +# stand (heights AND densities respond to theta), so EVERY trait re-shades the light every +# cohort reads -- the genuine resident feedback that routinely flips the sign relative to +# the frozen (rare-mutant) gradient. The new ingredient beyond the frozen census tape is the +# leaf's crown-centre LIGHT channel (dprofit_dL + the collar-rate light curvature, harvested +# by FD over a flat environment); on the tape the resident light correction is anchored so +# the baseline is exact at theta0 and the focal plant's own height->light slope stays in the +# frozen-census dprofit_dh (no double count). One reverse sweep per metric over the coupled +# whole-stand replay. Validated against tf24f_resident_census_gradient_fd (full SCM re-runs); +# see test-tf24f-census-gradient.R. Single-species; FIXED node schedule (TF24f is stiff). +tf24f_resident_census_gradient_ad <- function(scm, metrics = c("LAI", "size_moment"), + traits = NULL, species = 1L, + birth_rate = NULL, trait_rel_step = 1e-5) { + tf24f_require_strategy(scm) + if (is.null(traits)) traits <- tf24_default_traits() + # R0 gate (mirrors the multi-species path and FF16's coupled gate). The coupled double + # re-evolution must reproduce the resident stand before the linearised sensitivity tape + # is trusted. Past a short horizon the frozen-step replay drifts (joint env_err jumps + # ~1e-6 -> 1e-2 beyond patch lifetime ~4) and the coupled log_density<->canopy + # sensitivity runs away: TF24f's deeply-shaded leaf has a LARGE dprofit_dL, so the + # high-density shaded cohorts strongly amplify the canopy-reshaping feedback -- a + # stiffness the live SCM only tames by stepping ADAPTIVELY through it. The replay reuses + # the SCM's frozen step sizes, so the linearised sensitivity diverges (finite-but- + # astronomical, then NaN) with horizon. (FF16 stays robust to long horizons because its + # closed-form net has a bounded light response -- no leaf-solve light amplification.) + # Gate with a clear error rather than return NaN / garbage. feedback = "frozen" (the + # invasion gradient) is robust at all horizons. + diverged <- function(ee) + stop("the TF24f coupled (resident) census re-evolution is too stiff on this node ", + "schedule (joint env drift = ", signif(ee, 3), "). The single-species resident ", + "census gradient is reliable only on a short, finely-resolved FIXED schedule ", + "(patch lifetime ~4); use a shorter max_patch_lifetime or a finer fixed node ", + "schedule, or use feedback = 'frozen' (the invasion gradient, robust at all ", + "horizons).") + r0 <- tf24f_coupled_metrics(scm, metrics = metrics, species = species, + birth_rate = birth_rate) + if (!all(is.finite(r0$values)) || r0$env_err > 1e-2) diverged(r0$env_err) + # Fully native gradient: env + birth steps + boundary-node history from the live Patch + # (no tf24f_harvest, no Rcpp::as<> env). The boundary-node guard lives in the C++ entry. + strat <- scm$parameters$strategies[[species]] + g <- tf24f_coupled_gradient_native(scm, unlist(strat$pars), as.integer(species - 1L), + if (is.null(birth_rate)) -1 else birth_rate, + strat$k_acclim, strat$use_ad_gradient, traits, metrics, + scm$parameters$patch_area, trait_rel_step) + # Post-sweep guard: the linearised tape can still blow up (finite but astronomically + # large) where the gate passes but the sensitivity is stiff; reject rather than mislead. + if (!all(is.finite(g$jacobian)) || max(abs(g$jacobian)) > 1e12) diverged(r0$env_err) + g +} + +# Harvest a multi-species TF24f SCM into the all-species arrays the cross-species coupled +# engine consumes (#472 scope B, the cross-species resident Jacobian). The shared schedule +# + joint env, per-species parameter vectors / cohort birth steps / birth rates / TF24f +# acclimation knobs, and the all-species per-RK-stage boundary harvest +# (stand_newnode_*_stage_history_all, [step][stage][species]). The TF24f mirror of +# ff16_harvest_ms; the joint stand light is reconstructed in C++ from each species' +# re-evolved cohorts, so only these per-species pieces are needed from R. +tf24f_harvest_ms <- function(scm) { + patch <- scm$patch + nsp <- length(scm$parameters$strategies) + sh <- patch$step_history + eh <- patch$environment_history + if (length(eh) < 1L) { + stop("No resident schedule cached: run the SCM with control(save_RK45_cache = TRUE)") + } + nn_h <- patch$stand_newnode_height_stage_history_all + nn_c <- patch$stand_newnode_competition_stage_history_all + if (length(nn_h) < 1L) { + stop("feedback = 'resident' on a multi-species TF24f stand needs the all-species ", + "per-RK-stage harvest; re-run the resident SCM with control(save_RK45_cache = TRUE)") + } + list( + pp_list = lapply(seq_len(nsp), function(s) unlist(scm$parameters$strategies[[s]]$pars)), + eh = eh, sh = sh, + birth_list = lapply(seq_len(nsp), function(s) + vapply(patch$species[[s]]$node_times, + function(t) which.min(abs(sh - t)) - 1L, integer(1))), + birth_rate = vapply(seq_len(nsp), function(s) + scm$offspring_production[[s]] / scm$net_reproduction_ratios[[s]], numeric(1)), + k_acclim = vapply(seq_len(nsp), function(s) + scm$parameters$strategies[[s]]$k_acclim, numeric(1)), + use_ad_gradient = vapply(seq_len(nsp), function(s) + as.integer(isTRUE(scm$parameters$strategies[[s]]$use_ad_gradient)), integer(1)), + nn_h = nn_h, nn_c = nn_c, + patch_area = scm$parameters$patch_area, nsp = nsp) +} + +# TF24f CROSS-SPECIES resident (coupled) census trait gradient by reverse-mode AD (#472 +# scope B). On a multi-species stand, the TOTAL-stand d(census metric)/d(theta of species +# `species`): all species' cohorts are re-evolved together over the frozen schedule and +# the canopy light each RK stage is the JOINT reconstruction (sum over species of each +# species' trapezium); only the target species' traits are registered as tape inputs, so +# one reverse sweep per metric gives the cross-species total -- the target's traits re-shade +# the joint canopy that EVERY species reads (the cross term the frozen gradient zeroes). +# A cheap double R0 pass gates a diverged (stiff) node schedule before the sweep. Unlike +# FF16 (closed-form leaf), the TF24f cross-species tape's linearised log_density / g' +# derivative is sensitive to the JOINT-canopy reshaping, and amplifies when the species +# are widely separated or the schedule is coarse -- so the gate is tighter (joint-env drift +# > 1e-3, vs FF16's 1e-2) and a post-sweep finiteness / magnitude guard catches any +# residual blow-up. The cross-species gradient therefore needs a WELL-CONDITIONED stand: +# closely-spaced species on a fine FIXED node schedule (where it matches the FD reference; +# see tf24f_resident_census_gradient_ms_fd). Validated on such a stand. +tf24f_resident_census_gradient_ms_ad <- function(scm, metrics = c("LAI", "size_moment"), + traits = NULL, species = 1L, + trait_rel_step = 1e-5) { + tf24f_require_strategy(scm) + if (is.null(traits)) traits <- tf24_default_traits() + # Fully native (mirrors ff16_coupled_gradient_ms_native): the joint env + step schedule + # + per-species birth steps + all-species boundary harvest come from the live Patch; only + # the cheap per-species scalars (pp / recovered birth rates / acclimation knobs) are read + # from $parameters in R. No tf24f_harvest_ms (no Rcpp::as<> env, no O(stand) patch rebuild). + nsp <- length(scm$parameters$strategies) + if (species < 1L || species > nsp) stop("target species out of range") + pp_list <- lapply(seq_len(nsp), function(s) unlist(scm$parameters$strategies[[s]]$pars)) + br_vec <- vapply(seq_len(nsp), function(s) + scm$offspring_production[[s]] / scm$net_reproduction_ratios[[s]], numeric(1)) + k_acclim <- vapply(seq_len(nsp), function(s) + scm$parameters$strategies[[s]]$k_acclim, numeric(1)) + use_ad <- vapply(seq_len(nsp), function(s) + as.integer(isTRUE(scm$parameters$strategies[[s]]$use_ad_gradient)), integer(1)) + patch_area <- scm$parameters$patch_area + diverged <- function(ee) + stop("the multi-species TF24f coupled re-evolution is too stiff on this node ", + "schedule (joint env drift = ", signif(ee, 3), "). The cross-species gradient ", + "needs a well-conditioned stand: closely-spaced species on a fine FIXED node ", + "schedule (refine_schedule = FALSE). Re-run the resident SCM accordingly.") + r0 <- tf24f_coupled_gradient_ms_native(scm, pp_list, br_vec, k_acclim, use_ad, traits, + metrics, patch_area, as.integer(species), trait_rel_step, TRUE) + if (!all(is.finite(r0$values)) || r0$env_err > 1e-3) diverged(r0$env_err) + g <- tf24f_coupled_gradient_ms_native(scm, pp_list, br_vec, k_acclim, use_ad, traits, + metrics, patch_area, as.integer(species), trait_rel_step, FALSE) + # Post-sweep guard: the linearised tape can still blow up (finite but astronomically + # large) where the gate passes but the derivative is stiff; reject rather than mislead. + if (!all(is.finite(g$jacobian)) || max(abs(g$jacobian)) > 1e12) diverged(r0$env_err) + g +} + +# TF24f multi-species CROSS-SPECIES resident census FD reference (ground truth). Perturbs +# one species' (post-hyperpar) trait, re-runs the full SCM on the SAME fixed node schedule, +# and differences the TOTAL-stand metric (LAI = compute_competition(0); size_moment summed +# over species). Slow (one SCM solve per trait per side); keep `traits` small. +tf24f_resident_census_gradient_ms_fd <- function(p, env, ctrl, + metrics = c("LAI", "size_moment"), + traits = NULL, species = 1L, + rel_step = 1e-4) { + if (is.null(traits)) traits <- tf24_default_traits() + stand_metrics <- function(scm) { + patch <- scm$patch + out <- c(LAI = patch$compute_competition(0)) + sm <- 0 + for (s in seq_along(patch$species)) { + sp <- patch$species[[s]] + h <- sp$heights; dens <- exp(sp$log_densities) + ord <- order(h, decreasing = TRUE); h <- h[ord]; dens <- dens[ord] + phi <- dens * h + if (length(h) > 1) sm <- sm + sum(0.5 * (h[-length(h)] - h[-1]) * (phi[-length(h)] + phi[-1])) + } + out["size_moment"] <- sm + out[metrics] + } + run_with <- function(pp_override) { + pmod <- p + if (!is.null(pp_override)) { + s <- pmod$strategies[[species]]; s$pars <- pp_override; pmod$strategies[[species]] <- s + } + run_scm(pmod, env, ctrl, refine_schedule = FALSE) + } + values <- stand_metrics(run_with(NULL)) + pars0 <- p$strategies[[species]]$pars + jac <- matrix(0, length(metrics), length(traits), dimnames = list(metrics, traits)) + for (j in seq_along(traits)) { + tr <- traits[j] + if (!tr %in% names(pars0)) stop("unknown TF24f trait: ", tr) + d <- rel_step * max(abs(pars0[[tr]]), 1e-8) + pp_p <- pars0; pp_p[[tr]] <- pp_p[[tr]] + d + pp_m <- pars0; pp_m[[tr]] <- pp_m[[tr]] - d + jac[, j] <- (stand_metrics(run_with(pp_p)) - stand_metrics(run_with(pp_m))) / (2 * d) + } + list(jacobian = jac, values = values) +} + +# TF24f RESIDENT (coupled) census trait gradient (#472 scope B, the resident-feedback +# surface -- prototype). The TOTAL stand-level d(census metric)/d(theta): unlike the +# frozen (rare-mutant) gradient, every cohort's height + density feeds back through the +# canopy light that the whole stand reads, so the feedback routinely dominates and can +# flip the sign relative to the frozen reading. Computed as a central finite difference +# over the FULL SCM: perturb a (post-hyperpar) strategy parameter, re-run run_scm on the +# SAME fixed node schedule, and difference the realised stand metric. This is the ground- +# truth resident gradient (the SCM responds in full) and the reference the coupled AD +# engine must reproduce; it is slow (one SCM solve per trait per side), so keep `traits` +# small. LAI is read from the realised patch (compute_competition(0)); size_moment is the +# size-distribution first moment Sum density_i * height_i (trapezium over the stand). +tf24f_resident_census_gradient_fd <- function(p, env, ctrl, + metrics = c("LAI", "size_moment"), + traits = NULL, species = 1L, + rel_step = 1e-4) { + if (is.null(traits)) traits <- tf24_default_traits() + metric_set <- c("LAI", "size_moment") + bad <- setdiff(metrics, metric_set) + if (length(bad)) stop("unsupported resident metric(s): ", paste(bad, collapse = ", ")) + + # Realised stand metrics from a completed SCM (the coupled, self-shaded stand). + stand_metrics <- function(scm) { + patch <- scm$patch + out <- c(LAI = patch$compute_competition(0)) + sp <- patch$species[[species]] + h <- sp$heights; dens <- exp(sp$log_densities) + ord <- order(h, decreasing = TRUE) + h <- h[ord]; dens <- dens[ord] + # size_moment = trapezium of density*height over descending heights down to h0. + phi <- dens * h + sm <- 0 + if (length(h) > 1) sm <- sum(0.5 * (h[-length(h)] - h[-1]) * (phi[-length(h)] + phi[-1])) + out["size_moment"] <- sm + out[metrics] + } + run_with <- function(pp_override) { + pmod <- p + if (!is.null(pp_override)) { + s <- pmod$strategies[[species]] + s$pars <- pp_override + pmod$strategies[[species]] <- s + } + run_scm(pmod, env, ctrl, refine_schedule = FALSE) + } + + base_scm <- run_with(NULL) + values <- stand_metrics(base_scm) + pars0 <- p$strategies[[species]]$pars + jac <- matrix(0, length(metrics), length(traits), + dimnames = list(metrics, traits)) + for (j in seq_along(traits)) { + tr <- traits[j] + if (!tr %in% names(pars0)) stop("unknown TF24f trait: ", tr) + d <- rel_step * max(abs(pars0[[tr]]), 1e-8) + pp_p <- pars0; pp_p[[tr]] <- pp_p[[tr]] + d + pp_m <- pars0; pp_m[[tr]] <- pp_m[[tr]] - d + mp <- stand_metrics(run_with(pp_p)); mm <- stand_metrics(run_with(pp_m)) + jac[, j] <- (mp - mm) / (2 * d) + } + list(jacobian = jac, values = values) +} diff --git a/inst/RcppR6_classes.yml b/inst/RcppR6_classes.yml index 235c07a1..d2fbfe8f 100644 --- a/inst/RcppR6_classes.yml +++ b/inst/RcppR6_classes.yml @@ -174,6 +174,50 @@ Leaf: return_type: double args: [opt_root_psi: double] + dprofit_dvcmax25: + return_type: double + args: [opt_root_psi: double] + + dprofit_dg1_TF24: + return_type: double + args: [opt_root_psi: double] + + dprofit_dbeta2: + return_type: double + args: [opt_root_psi: double] + + dprofit_dkmax: + return_type: double + args: [opt_root_psi: double] + + dprofit_dEup: + return_type: double + args: [opt_root_psi: double] + + dprofit_db: + return_type: double + args: [opt_root_psi: double] + + dprofit_dc: + return_type: double + args: [opt_root_psi: double] + + dprofit_djmax25: + return_type: double + args: [opt_root_psi: double] + + dprofit_da: + return_type: double + args: [opt_root_psi: double] + + dprofit_dcurv_elec: + return_type: double + args: [opt_root_psi: double] + + dprofit_dcurv_colim: + return_type: double + args: [opt_root_psi: double] + psi_stem_to_ci: return_type: double args: [psi_stem: double, psi_upstream: double] @@ -305,6 +349,7 @@ Control: - node_gradient_direction: int - node_gradient_richardson: bool - node_gradient_richardson_depth: size_t + - node_gradient_exact_ad: bool - ode_step_size_initial: double - ode_step_size_min: double - ode_step_size_max: double @@ -439,6 +484,9 @@ Individual: net_mass_production_dt: args: [environment: E] return_type: double + growth_rate_gradient_exact: + args: [environment: E] + return_type: double reset_mortality: return_type: void resource_compensation_point: @@ -641,6 +689,33 @@ Patch: ode_aux: {type: "std::vector", access: function, name_cpp: "odelia::ode::r_ode_aux"} node_ode_size: {type: size_t, access: member} state: {type: "Rcpp::List", access: member, name_cpp: r_get_state} + # Resident environment trajectory cached during a save_RK45_cache run (the + # mutant-replay landscape). step_history holds the ODE step times + # {0, t_1, ...}; environment_history[n] holds the 6 frozen per-RK-stage + # environments for the step advancing step_history[n] -> step_history[n+1]. + # Exposed so the two-pass AD emergent-gradient driver (#472 scope B) can + # harvest the frozen resident light schedule without re-running in C++. + step_history: {type: "std::vector", access: field} + environment_history: {type: "std::vector >", access: field} + # Per-ODE-step resident stand (species 0) for the active-knot self-shading + # reconstruction: node heights and per-node competition effects per step. + stand_height_history: {type: "std::vector >", access: field} + stand_competition_history: {type: "std::vector >", access: field} + # Per-RK-STAGE resident stand (species 0): [step][stage 0..5][cohort], aligned + # 1:1 with environment_history's 6 per-stage environments. Faithful harvest for + # the RESIDENT total-gradient reconstruction (#472 scope B, R0). + stand_height_stage_history: {type: "std::vector > >", access: field} + stand_competition_stage_history: {type: "std::vector > >", access: field} + # Boundary (new_node) per-RK-stage state (species 0): height + competition effect, + # the tail term Species::compute_competition adds beyond the `nodes` loop. [step][stage]. + stand_newnode_height_stage_history: {type: "std::vector >", access: field} + stand_newnode_competition_stage_history: {type: "std::vector >", access: field} + # ALL-SPECIES per-RK-stage harvest (#472 scope B, R2 cross-species coupled): the + # stand is [step][stage][species][cohort], the boundary node [step][stage][species]. + stand_height_stage_history_all: {type: "std::vector > > > ", access: field} + stand_competition_stage_history_all: {type: "std::vector > > > ", access: field} + stand_newnode_height_stage_history_all: {type: "std::vector > >", access: field} + stand_newnode_competition_stage_history_all: {type: "std::vector > >", access: field} methods: introduce_new_node: return_type: void diff --git a/inst/include/plant/RcppR6_post.hpp b/inst/include/plant/RcppR6_post.hpp index 1356117a..09ad2b60 100644 --- a/inst/include/plant/RcppR6_post.hpp +++ b/inst/include/plant/RcppR6_post.hpp @@ -323,6 +323,7 @@ template <> inline SEXP wrap(const plant::Control& x) { ret["node_gradient_direction"] = Rcpp::wrap(x.node_gradient_direction); ret["node_gradient_richardson"] = Rcpp::wrap(x.node_gradient_richardson); ret["node_gradient_richardson_depth"] = Rcpp::wrap(x.node_gradient_richardson_depth); + ret["node_gradient_exact_ad"] = Rcpp::wrap(x.node_gradient_exact_ad); ret["ode_step_size_initial"] = Rcpp::wrap(x.ode_step_size_initial); ret["ode_step_size_min"] = Rcpp::wrap(x.ode_step_size_min); ret["ode_step_size_max"] = Rcpp::wrap(x.ode_step_size_max); @@ -371,6 +372,8 @@ template <> inline plant::Control as(SEXP x) { ret.node_gradient_richardson = Rcpp::as(xl["node_gradient_richardson"]); // ret.node_gradient_richardson_depth = Rcpp::as(xl["node_gradient_richardson_depth"]); ret.node_gradient_richardson_depth = Rcpp::as(xl["node_gradient_richardson_depth"]); + // ret.node_gradient_exact_ad = Rcpp::as(xl["node_gradient_exact_ad"]); + ret.node_gradient_exact_ad = Rcpp::as(xl["node_gradient_exact_ad"]); // ret.ode_step_size_initial = Rcpp::as(xl["ode_step_size_initial"]); ret.ode_step_size_initial = Rcpp::as(xl["ode_step_size_initial"]); // ret.ode_step_size_min = Rcpp::as(xl["ode_step_size_min"]); diff --git a/inst/include/plant/control.h b/inst/include/plant/control.h index 95d5fdfe..9fdd641e 100644 --- a/inst/include/plant/control.h +++ b/inst/include/plant/control.h @@ -59,6 +59,11 @@ struct Control { int node_gradient_direction; bool node_gradient_richardson; size_t node_gradient_richardson_depth; + // Use the strategy's exact AD growth-rate gradient in Node::growth_rate_gradient + // when it provides one (#537 A1), instead of the finite difference. Default + // false (the FD path is unchanged); strategies without an AD gradient fall + // back to FD regardless. + bool node_gradient_exact_ad; double ode_step_size_initial; double ode_step_size_min; diff --git a/inst/include/plant/gradient/coupled_canopy.h b/inst/include/plant/gradient/coupled_canopy.h new file mode 100644 index 00000000..ae392f13 --- /dev/null +++ b/inst/include/plant/gradient/coupled_canopy.h @@ -0,0 +1,56 @@ +// -*-c++-*- +#ifndef PLANT_GRADIENT_COUPLED_CANOPY_H_ +#define PLANT_GRADIENT_COUPLED_CANOPY_H_ + +// Shared coupled-resident canopy kernel (#472 scope B, refactor+optimize phase). The +// Yokozawa light-competition trapezium is identical across FF16 and TF24/TF24f -- in +// all of them area_leaf is pure allometry and competition is density*k_I*area_leaf*Q +// with Q = (1-(z/h)^eta)^2; the leaf optimiser (TF24) affects only the demographic +// RATES, not the light-field geometry. So the +// per-RK-stage canopy reconstruction shared this code by hand-copy (ff16's +// coupled_comp_at == tf24f's tf24f_comp_at); this is the single templated source. + +#include +#include +#include +#include // active-type value extraction for the crown-height check + +namespace plant { +namespace gradient { + +// Double value of a scalar that may be a plain double or an XAD active type. +inline double scalar_value(double v) { return v; } +template double scalar_value(const T& v) { return xad::value(v); } + +// Canopy competition at height z: the descending-height trapezium of geff_i * Q(z/h_i), +// Q = (1-(z/h)^eta)^2, over the active stand (h sorted descending, geff the per-cohort +// weight density*k_I*area_leaf). Cohorts whose crown is below z contribute nothing. +// Returned UN-divided by patch area (the caller applies Beer's law / area). Active in +// the heights/weights, so every trait that moves a height or density re-shades the +// canopy -- the resident coupling. Bit-for-bit the hand-rolled comp_at it replaces. +template +S canopy_comp_at(double z, const std::vector& h, const std::vector& geff, + double eta) { + using std::pow; + const std::size_t n = h.size(); + if (n < 2) return S(0.0); + auto g = [&](std::size_t i) -> S { + if (z >= scalar_value(h[i])) return S(0.0); // no leaf area above the crown + const S u = S(z) / h[i]; + const S om = S(1.0) - pow(u, S(eta)); + return geff[i] * (om * om); + }; + S comp = S(0.0); + S gp = g(0); S hp = h[0]; + for (std::size_t i = 1; i < n; ++i) { + S gi = g(i); + comp = comp + (hp - h[i]) * (gp + gi); + hp = h[i]; gp = gi; + } + return S(0.5) * comp; +} + +} // namespace gradient +} // namespace plant + +#endif diff --git a/inst/include/plant/gradient/scm_harvest.h b/inst/include/plant/gradient/scm_harvest.h new file mode 100644 index 00000000..89a516bc --- /dev/null +++ b/inst/include/plant/gradient/scm_harvest.h @@ -0,0 +1,122 @@ +// -*-c++-*- +#ifndef PLANT_GRADIENT_SCM_HARVEST_H_ +#define PLANT_GRADIENT_SCM_HARVEST_H_ + +// Strategy-agnostic helpers shared by the emergent-gradient engines (#472 scope B, +// refactor+optimize phase). These collapse boilerplate that was hand-copied across +// ff16_emergent.cpp / tf24_emergent.cpp / tf24f_emergent.cpp once the harvest moved +// into C++ (the native-SCM read + birth-step computation), plus the census +// height-trapezium reduction. Templated so the same code serves every strategy / +// scalar type; header-only so each engine TU instantiates what it needs. + +#include +#include +#include + +namespace plant { +namespace gradient { + +// Recover the constant birth rate from a run if not supplied (given < 0): the SCM's +// offspring_production / net_reproduction_ratio for the species (both per-species). +// `scm` is anything with operator-> to the SCM (a raw SCM* or an RcppR6 wrapper). +template +double recover_birth_rate(Scm& scm, std::size_t species, double given) { + if (given >= 0.0) return given; + return scm->offspring_production()[species] / + scm->net_reproduction_ratios()[species]; +} + +// Per-cohort birth steps: the (0-based) step time nearest each node introduction time. +// Mirrors the which.min(abs(sh - node_time)) the R harvests used. +template +std::vector birth_steps(const Patch& patch, std::size_t species) { + const std::vector& sh = patch.step_history; + const std::vector nt = patch.at_species(species).node_times(); + std::vector birth(nt.size()); + for (std::size_t i = 0; i < nt.size(); ++i) { + std::size_t best = 0; + double bd = std::abs(sh[0] - nt[i]); + for (std::size_t k = 1; k < sh.size(); ++k) { + const double d = std::abs(sh[k] - nt[i]); + if (d < bd) { bd = d; best = k; } + } + birth[i] = static_cast(best); + } + return birth; +} + +// The offspring-tape harvest computed natively from a live Patch: the per-cohort +// birth steps, node-spacing trapezoid weights (so offspring_production = sum_i tw_i * +// offspring_i), per-RK-stage patch survival, and survival-at-birth. Strategy-agnostic +// (every Patch exposes step_history / at_species / r_pr_survival / environment_history); +// shared by FF16's build_frozen_scm and the TF24/TF24f native offspring entries so the +// O(stand) survival loop runs once in C++ instead of being rebuilt in R. ppsurv is +// returned row-major [N*6] (caller wraps into the engine's matrix form). +struct OffspringWeights { + std::size_t N; // = environment_history.size() (Cash-Karp steps) + std::vector birth; // per-cohort birth step (0-based) + std::vector tw; // trapezoid node weight * patch_density * S_D * birth_rate + std::vector ppsab; // pr_patch_survival_at_birth (per cohort) + std::vector ppsurv; // [N*6] row-major: pr_survival at sh[k] + ah[s]*h_k +}; +template +OffspringWeights offspring_weights(const Patch& patch, std::size_t species, + double S_D, double birth_rate) { + OffspringWeights w; + const std::vector& sh = patch.step_history; + const auto& sp = patch.at_species(species); + const std::vector nt = sp.node_times(); + const std::vector pdens = sp.r_patch_densities(); + w.ppsab = sp.r_pr_patch_survival_at_birth(); + const std::size_t nC = nt.size(); + w.birth = birth_steps(patch, species); + + std::vector tcoef(nC, 0.0); + if (nC >= 2) { + tcoef[0] = 0.5 * (nt[1] - nt[0]); + tcoef[nC - 1] = 0.5 * (nt[nC - 1] - nt[nC - 2]); + for (std::size_t i = 1; i + 1 < nC; ++i) tcoef[i] = 0.5 * (nt[i + 1] - nt[i - 1]); + } + w.tw.resize(nC); + for (std::size_t i = 0; i < nC; ++i) w.tw[i] = tcoef[i] * pdens[i] * S_D * birth_rate; + + w.N = patch.environment_history.size(); + const double ah[6] = {0.0, 0.2, 0.3, 0.6, 1.0, 0.875}; + w.ppsurv.assign(w.N * 6, 0.0); + for (std::size_t k = 0; k < w.N; ++k) { + const double hN = sh[k + 1] - sh[k]; + for (int st = 0; st < 6; ++st) + w.ppsurv[k * 6 + st] = patch.r_pr_survival(sh[k] + ah[st] * hN); + } + return w; +} + +// Census height-trapezium over the replayed cohorts, descending height + the +// pending-seed tail term (mirrors Species::compute_competition / Patch census +// integral). Generic over the cohort container: `geth(i)` / `getdens(i)` / `getmhw(i)` +// extract the per-cohort height / number density / heartwood mass at index i, and +// `psi(h, dens, mhw)` is the metric kernel. `ord` indexes cohorts in descending +// height; (h0, dens_new) are the pending-seed height / density for the ground tail. +// Identical arithmetic to the hand-rolled reductions it replaces (bit-for-bit). +template +S census_trapezium(std::size_t nC, const std::vector& ord, S h0, S dens_new, + GetH geth, GetDens getdens, GetMhw getmhw, Psi psi) { + std::vector phi(nC); + for (std::size_t i = 0; i < nC; ++i) phi[i] = psi(geth(i), getdens(i), getmhw(i)); + S J = S(0.0); + for (std::size_t j = 0; j + 1 < nC; ++j) { + const std::size_t a = ord[j], b = ord[j + 1]; + J = J + S(0.5) * (geth(a) - geth(b)) * (phi[a] + phi[b]); + } + if (nC > 0) { + const std::size_t last = ord[nC - 1]; + S phi_new = psi(h0, dens_new, S(0.0)); + J = J + S(0.5) * (geth(last) - h0) * (phi[last] + phi_new); + } + return J; +} + +} // namespace gradient +} // namespace plant + +#endif diff --git a/inst/include/plant/individual.h b/inst/include/plant/individual.h index 94939ee0..2af181db 100644 --- a/inst/include/plant/individual.h +++ b/inst/include/plant/individual.h @@ -11,7 +11,13 @@ namespace plant { -template class Individual { +// Templated on the scalar type S (#472 scope B / #537, Milestone C) so a plant's +// ODE state can be an AD active type for reverse-mode gradients. S defaults to +// double, so every existing `Individual` is unchanged and bit-identical; +// only AD paths instantiate Individual (and then only the members +// they use are compiled -- the double-only ODE-iterator methods stay uncompiled +// for the AD instantiation until the ODE-state boundary is wired). +template class Individual { public: typedef T strategy_type; typedef E environment_type; @@ -32,39 +38,39 @@ template class Individual { } // useage: state(HEIGHT_INDEX) - double state(std::string name) const { + S state(std::string name) const { return vars.state(strategy->state_index.at(name)); } - double state(int i) const { return vars.state(i); } - + S state(int i) const { return vars.state(i); } + // useage:_rate("area_heartwood") - double rate(std::string name) const { + S rate(std::string name) const { return vars.rate(strategy->state_index.at(name)); } - double rate(int i) const { return vars.rate(i); } + S rate(int i) const { return vars.rate(i); } // useage: set_state("height", 2.0) - void set_state(std::string name, double v) { + void set_state(std::string name, S v) { int i = strategy->state_index.at(name); vars.set_state(i, v); strategy->update_dependent_aux(i, vars); } - void set_state(int i, double v) { + void set_state(int i, S v) { vars.set_state(i, v); strategy->update_dependent_aux(i, vars); } // aux vars by name and index - double aux(std::string name) const { + S aux(std::string name) const { return vars.aux(strategy->aux_index.at(name)); } - double aux(int i) const { return vars.aux(i); } + S aux(int i) const { return vars.aux(i); } // set # consumable resources based on env. variables void resize_consumption_rates(int i) { vars.resize_consumption_rates(i); } - double consumption_rate(int i) const { return vars.consumption_rate(i); } + S consumption_rate(int i) const { return vars.consumption_rate(i); } double compute_competition(double z) const { return strategy->compute_competition(z, vars); @@ -135,6 +141,14 @@ template class Individual { void reset_mortality() { set_state("mortality", 0.0); } + // Exact d(growth rate)/d(height) at the current height, delegated to the + // strategy's AD gradient (#537 A1); returns NA if the strategy provides none, + // so Node::growth_rate_gradient can fall back to finite differences. + double growth_rate_gradient_exact(const environment_type& environment) const { + return strategy->growth_rate_gradient_height_ad(vars.state(HEIGHT_INDEX), + environment); + } + double growth_rate_given_height(double height, const environment_type& environment) { // Called repeatedly from the finite-difference gradient (Node:: // growth_rate_gradient), so address height by integer slot rather than the @@ -172,12 +186,12 @@ template class Individual { // ! External R code depends on knowing r internals for like growing plant to // ! height or something - Internals r_internals() const { return vars; } + basic_internals r_internals() const { return vars; } const Control &control() const { return strategy->control; } private: strategy_type_ptr strategy; - Internals vars; + basic_internals vars; }; template Individual make_individual(T s) { diff --git a/inst/include/plant/internals.h b/inst/include/plant/internals.h index 1a7389e1..e24359a0 100644 --- a/inst/include/plant/internals.h +++ b/inst/include/plant/internals.h @@ -14,17 +14,24 @@ // TODO(#483): extra_state bounds, upper and lower limits namespace plant { -class Internals { +// Templated on the scalar type S so the per-plant state container can hold an +// AD active type for reverse-mode gradients (#472 scope B / #537, Milestone C). +// The `Internals` alias below pins S = double, so every existing use across the +// package keeps compiling and stays bit-identical; only AD paths instantiate +// basic_internals. NA_REAL initialisers are wrapped in S(...) so they +// also work when S is an AD type. +template +class basic_internals { public: - Internals(size_t s_size=0, size_t a_size=0, size_t r_size=0) + basic_internals(size_t s_size=0, size_t a_size=0, size_t r_size=0) : state_size(s_size), aux_size(a_size), resource_size(r_size), - states(s_size, 0.0), - rates(s_size, NA_REAL) , - auxs(a_size, 0.0), - consumption_rates(r_size, NA_REAL) + states(s_size, S(0.0)), + rates(s_size, S(NA_REAL)) , + auxs(a_size, S(0.0)), + consumption_rates(r_size, S(NA_REAL)) {} size_t state_size; size_t aux_size; @@ -32,35 +39,39 @@ class Internals { // Perhaps make these private so the () overloads below have some use - std::vector states; - std::vector rates; - std::vector auxs; - std::vector consumption_rates; // not quite as pithy + std::vector states; + std::vector rates; + std::vector auxs; + std::vector consumption_rates; // not quite as pithy - double state(int i) const { return states[i]; } - double rate(int i) const { return rates[i]; } - double aux(int i) const { return auxs[i]; } - double consumption_rate(int i) const { return consumption_rates[i]; } + S state(int i) const { return states[i]; } + S rate(int i) const { return rates[i]; } + S aux(int i) const { return auxs[i]; } + S consumption_rate(int i) const { return consumption_rates[i]; } - void set_state(int i, double v) { states[i] = v; } - void set_rate(int i, double v) { rates[i] = v; } - void set_aux(int i, double v) { auxs[i] = v; } - void set_consumption_rate(int i, double v) { consumption_rates[i] = v; } + void set_state(int i, S v) { states[i] = v; } + void set_rate(int i, S v) { rates[i] = v; } + void set_aux(int i, S v) { auxs[i] = v; } + void set_consumption_rate(int i, S v) { consumption_rates[i] = v; } void resize(size_t new_size, size_t new_aux_size) { state_size = new_size; aux_size = new_aux_size; - states.resize(new_size, 0.0); - rates.resize(new_size, NA_REAL); - auxs.resize(new_aux_size, 0.0); + states.resize(new_size, S(0.0)); + rates.resize(new_size, S(NA_REAL)); + auxs.resize(new_aux_size, S(0.0)); } void resize_consumption_rates(size_t new_resource_size) { resource_size = new_resource_size; - consumption_rates.resize(new_resource_size, NA_REAL); + consumption_rates.resize(new_resource_size, S(NA_REAL)); } }; +// Default state container used everywhere in the package (bit-identical to the +// pre-templating concrete class). +using Internals = basic_internals; + } // namespace plant #endif diff --git a/inst/include/plant/leaf_model.h b/inst/include/plant/leaf_model.h index eee748fd..31986a52 100644 --- a/inst/include/plant/leaf_model.h +++ b/inst/include/plant/leaf_model.h @@ -332,6 +332,50 @@ class Leaf { // the noisy finite-difference gradient. Assumes prepare_collar_solve setup has // run (psi_soil_inverted_ etc.), as evaluate_root_collar_psi does. double dprofit_droot_collar_psi(double opt_root_psi); + // Exact d(profit*)/d(vcmax_25) at the optimised operating point (#472 scope B / + // Phase F): the envelope theorem fixes the optimal collar potential, the IFT + // handles the psi_stem->ci root-find. The first TF24 trait gradient; the pattern + // the TF24 net-production kernel reuses for every leaf trait. + double dprofit_dvcmax25(double opt_root_psi); + // Exact d(profit*)/d(hydraulic trait) at the optimised operating point (#472 + // scope B / Phase F1-full). g1_TF24 and beta2 enter only the hydraulic cost + // (no transport / assimilation change), so the envelope theorem reduces each + // to minus the explicit cost derivative (forward-mode AD of the templated + // cost). The harder hydraulic traits (b, c, K_s) -- which also move psi_stem + // and ci -- follow the dprofit_droot_collar_psi transport+IFT pattern. + double dprofit_dg1_TF24(double opt_root_psi); + double dprofit_dbeta2(double opt_root_psi); + // d(profit*)/d(leaf_specific_conductance_max_): a TRANSPORT trait (moves + // psi_stem and ci, not the cost explicitly). The TF24 trait K_s scales k_max + // linearly (k_max = K_s*theta/(h*eta_c)), so the strategy chains by k_max/K_s. + double dprofit_dkmax(double opt_root_psi); + // d(profit*)/d(E_up_): sensitivity to the soil->collar water uptake. Used by + // the mass-cascade trait a_r1 (root mass per leaf area), which scales every + // root resistance by 1/a_r1 hence E_up_ linearly (d E_up_/d a_r1 = E_up_/a_r1). + double dprofit_dEup(double opt_root_psi); + // d(profit*)/d(b) and d(profit*)/d(c): the xylem vulnerability shape traits + // (prop_cond = exp(-(psi/b)^c)). They reshape the transpiration spline (so + // psi_stem moves) and enter the cost explicitly; ci/benefit are frozen because + // the operating-point transpiration equals the (root-vulnerability) uptake + // E_up_. dprofit_dbc(.., wrt_b) is the shared core. dtranspiration_integral_ + // dtrait is dS/d(trait) for the cumulative transpiration curve S. + double dprofit_db(double opt_root_psi); + double dprofit_dc(double opt_root_psi); + double dprofit_dbc(double opt_root_psi, bool wrt_b); + double dtranspiration_integral_dtrait(double x, bool wrt_b); + // d(profit*)/d(photosynthesis trait): jmax_25, a (quantum yield) and the two + // curvature factors affect only assimilation (vcmax-like), so the envelope + + // IFT pattern of dprofit_dvcmax25 applies. dprofit_dphoto(.., which) is the + // shared core (which: 0=jmax_25, 1=a, 2=curv_elec, 3=curv_colim, 4=PPFD). + double dprofit_djmax25(double opt_root_psi); + double dprofit_da(double opt_root_psi); + double dprofit_dcurv_elec(double opt_root_psi); + double dprofit_dcurv_colim(double opt_root_psi); + // d(profit*)/d(absorbed radiation PPFD) -- the light channel of d(profit)/ + // d(height) for the exact-AD growth-rate gradient (#472 scope B). PPFD enters + // only through electron_transport, so it reuses the dprofit_dphoto core. + double dprofit_dPPFD(double opt_root_psi); + double dprofit_dphoto(double opt_root_psi, int which); // Analytic d(E_up_)/d(collar potential) for the soil->root-collar uptake // (kg H2O m^-2 s^-1 per MPa of signed collar potential P_x_r), mirroring the // general branch of E_from_Soil_to_Root_Collar layer by layer. The integral's diff --git a/inst/include/plant/models/ff16_environment.h b/inst/include/plant/models/ff16_environment.h index ce795d1b..d2be963a 100644 --- a/inst/include/plant/models/ff16_environment.h +++ b/inst/include/plant/models/ff16_environment.h @@ -7,6 +7,7 @@ #include #include // ShadingModel, shading_model_from_string #include // std::log, std::exp, std::floor (PPA stepping) +#include // std::numeric_limits (AD-deriv NaN signal) using namespace Rcpp; @@ -81,6 +82,16 @@ class FF16_Environment : public Environment { return step_light(light_availability.get_value_at_height(height, cap)); } + // Analytic d(light)/d(height) for the SMOOTH models (deep-crown/crown-centre), + // where step_light is the identity so the derivative is just the resource + // spline's. Returns NaN for the PPA stepped profile (non-smooth) so AD callers + // fall back. Enables exact AD gradients through the light environment + // (#472 scope B / #537). + double get_environment_deriv_at_height(double height) const { + if (light_profile_stepped) return std::numeric_limits::quiet_NaN(); + return light_availability.get_value_deriv_at_height(height); + } + // Discretise a smooth light value into PPA canopy layers. For the smooth // models this is a single predicted branch returning the input unchanged, so // it adds no measurable cost to deep-crown/crown-centre. For PPA it maps the diff --git a/inst/include/plant/models/ff16_production_kernel.h b/inst/include/plant/models/ff16_production_kernel.h index 702956c3..5c697cc5 100644 --- a/inst/include/plant/models/ff16_production_kernel.h +++ b/inst/include/plant/models/ff16_production_kernel.h @@ -62,6 +62,13 @@ struct FF16ProdPars { S r_l, r_s, r_b, r_r; S k_l, k_b, k_s, k_r; S a_bio, a_y; + // Allometry + allocation parameters for the height-growth rate (Milestone C). + S a_l1, a_l2; // height <-> leaf-area allometry [eqn 2/3] + S a_f1, a_f2, hmat; // reproduction-allocation logistic [eqn 16] + // Demographic rate parameters for the full compute_rates fill (Milestone C): + // fecundity [eqn 17] and mortality [eqn 21]. + S omega, a_f3; // seed mass + accessory reproduction cost (fecundity_dt) + S d_I, a_dG1, a_dG2; // mortality: growth-independent + growth-dependent }; // Whole single-plant net production under the CROWN-TOP assimilation variant (a @@ -120,6 +127,491 @@ S ff16_assimilation_deep_crown_replay(S a_p1, S a_p2, S area_leaf, return area_leaf * A; } +// Total live mass given height (#472 scope B): the same leaf+sapwood+bark+root mass +// cascade as ff16_net_from_components, as a function of height. Mirrors +// FF16_Strategy::mass_live_given_height. The seedling height_0 solves +// mass_live_given_height(h0) = omega, so this is the residual whose root is height_0; +// differentiating it gives d(height_0)/d(trait) by the implicit function theorem +// d(h0)/d(theta) = -(d mass_live/d theta - [theta==omega]) / (d mass_live/d height) +// at h0 -- the seedling-size response of the emergent gradient (the #539 IFT pattern). +template +S ff16_mass_live_given_height(const FF16ProdPars& p, S height) { + const S area_leaf = ff16_area_leaf(p.a_l1, p.a_l2, height); + const S area_sapwood = area_leaf * p.theta; + const S mass_leaf = area_leaf * p.lma; + const S mass_sapwood = area_sapwood * height * p.eta_c * p.rho; + const S mass_bark = p.a_b1 * area_sapwood * height * p.eta_c * p.rho; + const S mass_root = p.a_r1 * area_leaf; + return mass_leaf + mass_sapwood + mass_bark + mass_root; +} + +// Establishment probability (the recruitment filter), scalar-templated as a +// function of the SEEDLING net production (#472 scope B). Mirrors +// FF16_Strategy::establishment_probability: +// pr_estab = decay_over_time / ((a_d0 * area_leaf_0 / net0)^2 + 1) (net0 > 0), +// where net0 is the seedling's net production in the birth environment and +// decay_over_time = exp(-recruitment_decay * birth_time). recruitment_decay, the +// birth time and a_d0 are not physiology traits, so they fold to doubles; the trait +// dependence enters through net0 (and area_leaf_0). A node's initial mortality state +// is -log(pr_estab), so seeding the taped replay with -log(ff16_establishment_ +// probability(...)) -- rather than a frozen constant -- makes the recruitment filter +// part of the emergent gradient. net0 > 0 is a frozen pass-1 sign. +template +S ff16_establishment_probability(S area_leaf_0, S net0, double a_d0, + double decay_over_time) { + const S tmp = a_d0 * area_leaf_0 / net0; + return decay_over_time / (tmp * tmp + 1.0); +} + +// --------------------------------------------------------------------------- +// Height-growth rate pieces (#472 scope B, Milestone C). Mirror +// FF16_Strategy::{fraction_allocation_growth, dheight_darea_leaf, +// dmass_*_darea_leaf, darea_leaf_dmass_live} and the dheight/dt assembly in +// compute_rates. Elementary, so dheight/dt is differentiable w.r.t. height +// (A1 -- the exact growth-rate gradient now done by finite difference in +// Node::growth_rate_gradient) and w.r.t. traits. +// --------------------------------------------------------------------------- + +// [eqn 16] Fraction of production allocated to reproduction (logistic in height). +template +S ff16_fraction_allocation_reproduction(S a_f1, S a_f2, S hmat, S height) { + using std::exp; + return a_f1 / (1.0 + exp(a_f2 * (1.0 - height / hmat))); +} + +// [eqn 16] Fraction of production allocated to growth = 1 - reproduction. +template +S ff16_fraction_allocation_growth(S a_f1, S a_f2, S hmat, S height) { + return 1.0 - ff16_fraction_allocation_reproduction(a_f1, a_f2, hmat, height); +} + +// d(height)/d(area_leaf): derivative of the [eqn 2] allometry. +template +S ff16_dheight_darea_leaf(S a_l1, S a_l2, S area_leaf) { + using std::pow; + return a_l1 * a_l2 * pow(area_leaf, a_l2 - 1.0); +} + +// d(area_leaf)/d(mass_live): reciprocal of the summed per-component mass +// derivatives (leaf + sapwood + bark + root) w.r.t. area_leaf. +template +S ff16_darea_leaf_dmass_live(const FF16ProdPars& p, S area_leaf) { + using std::pow; + const S dmass_leaf = p.lma; // d(area_leaf*lma) + const S dmass_sapwood = p.rho * p.eta_c * p.a_l1 * p.theta * + (p.a_l2 + 1.0) * pow(area_leaf, p.a_l2); + const S dmass_bark = p.a_b1 * dmass_sapwood; + const S dmass_root = p.a_r1; + return 1.0 / (dmass_leaf + dmass_sapwood + dmass_bark + dmass_root); +} + +// dheight/dt given net production (the compute_rates growth assembly): returns +// 0 when net is non-positive (the growth clamp), else dheight_darea_leaf * +// area_leaf_dt. Shared by every assimilation variant. +template +S ff16_height_dt_from_net(const FF16ProdPars& p, S height, S area_leaf, S net) { + if (net <= 0.0) return S(0.0); + const S frac_growth = ff16_fraction_allocation_growth(p.a_f1, p.a_f2, p.hmat, height); + const S darea_dmass = ff16_darea_leaf_dmass_live(p, area_leaf); + const S area_leaf_dt = net * frac_growth * darea_dmass; + return ff16_dheight_darea_leaf(p.a_l1, p.a_l2, area_leaf) * area_leaf_dt; +} + +// dheight/dt for a plant of the given height under the CROWN-TOP assimilation +// variant, in light light_E. area_leaf is derived from height so the gradient +// w.r.t. height flows through the whole chain. +template +S ff16_height_dt_crown_top(const FF16ProdPars& p, S height, S light_E) { + const S area_leaf = ff16_area_leaf(p.a_l1, p.a_l2, height); + const S net = ff16_net_mass_production_crown_top(p, height, area_leaf, light_E); + return ff16_height_dt_from_net(p, height, area_leaf, net); +} + +// The five ODE state rates FF16_Strategy::compute_rates writes, plus the net +// production aux. Scalar-templated (#472 scope B, Milestone C) so the whole +// demographic rate fill differentiates w.r.t. a trait by reverse-mode AD. +template +struct FF16Rates { + S net_mass_production_dt; + S height_dt; + S fecundity_dt; + S area_heartwood_dt; + S mass_heartwood_dt; + S mortality_dt; +}; + +// Full FF16 compute_rates fill for the CROWN-TOP assimilation variant (single +// light evaluation light_E). Mirrors FF16_Strategy::compute_rates EXACTLY: the +// net>0 growth clamp gates the growth/fecundity/heartwood rates, and mortality +// is the [eqn 21] growth-independent + growth-dependent sum (productivity = +// net/area_leaf). `mortality_finite` is the frozen util::is_finite(cumulative +// mortality) branch -- a pass-1 (double) control-flow decision, passed in so the +// taped replay is branch-free (it never differentiates the is_finite test). +// Deep-crown differs only in how `net` is formed (the frozen-replay crown +// integral, ff16_assimilation_deep_crown_replay -> ff16_net_from_components), +// so a deep-crown fill reuses everything below by substituting that `net`. +// The rate fill SHARED by every assimilation variant: given the net production +// (however it was formed -- single-light crown-top, or the deep-crown crown +// integral) and the area_leaf, write the five ODE rates. This is the part of +// FF16_Strategy::compute_rates downstream of `net`; the net>0 growth clamp gates +// growth/fecundity/heartwood, mortality is the [eqn 21] sum (productivity = +// net/area_leaf). `mortality_finite` is the frozen is_finite branch (a pass-1 +// control-flow decision, so the taped replay never differentiates the test). +template +FF16Rates ff16_compute_rates_from_net(const FF16ProdPars& p, S height, + S area_leaf, S net, + bool mortality_finite) { + FF16Rates r; + r.net_mass_production_dt = net; + if (net > 0.0) { + const S frac_repro = ff16_fraction_allocation_reproduction(p.a_f1, p.a_f2, + p.hmat, height); + r.height_dt = ff16_height_dt_from_net(p, height, area_leaf, net); + r.fecundity_dt = net * frac_repro / (p.omega + p.a_f3); + const S area_sapwood = area_leaf * p.theta; // [eqn 4] + r.area_heartwood_dt = p.k_s * area_sapwood; // turnover of sapwood area + const S mass_sapwood = area_sapwood * height * p.eta_c * p.rho; + r.mass_heartwood_dt = p.k_s * mass_sapwood; // turnover_sapwood(mass) + } else { + r.height_dt = S(0.0); r.fecundity_dt = S(0.0); + r.area_heartwood_dt = S(0.0); r.mass_heartwood_dt = S(0.0); + } + // [eqn 21] instantaneous mortality rate; productivity_area = net / area_leaf. + using std::exp; + if (mortality_finite) { + const S productivity_area = net / area_leaf; + r.mortality_dt = p.d_I + p.a_dG1 * exp(-p.a_dG2 * productivity_area); + } else { + r.mortality_dt = S(0.0); + } + return r; +} + +template +FF16Rates ff16_compute_rates_crown_top(const FF16ProdPars& p, S height, + S light_E, bool mortality_finite) { + const S area_leaf = ff16_area_leaf(p.a_l1, p.a_l2, height); + const S net = ff16_net_mass_production_crown_top(p, height, area_leaf, light_E); + return ff16_compute_rates_from_net(p, height, area_leaf, net, mortality_finite); +} + +// [eqn] Yokozawa leaf-area density q(z,H) = 2 eta (1 - u^eta) u^eta / z, +// u = z/H. Mirrors CanopyShape::q exactly (eta is a fixed double; the gradient +// flows through the active u and z). The deep-crown crown integral weights the +// per-depth assimilation by this. +template +S ff16_canopy_q(double eta, S u, S z) { + using std::pow; + const S u_eta = pow(u, eta); + return 2.0 * eta * (1.0 - u_eta) * u_eta / z; +} + +// Single-plant height TRAJECTORY: integrate dheight/dt from h0 over n fixed RK4 +// steps to age t_end, in a fixed crown-top light light_E (#472 scope B). This is +// the bridge from instantaneous-rate gradients to time-integrated emergent +// outputs: the whole trajectory is scalar-templated, so reverse/forward AD gives +// d(height at age t_end)/d(trait) -- a calibration gradient through the growth +// ODE. A fixed step schedule is exactly the frozen-schedule formulation +// end-to-end AD needs (the adaptive stepper is a pass-1 discovery, replayed). +template +S ff16_grow_height(const FF16ProdPars& p, S h0, S light_E, + double t_end, int n_steps) { + S h = h0; + const double dt = t_end / n_steps; + for (int i = 0; i < n_steps; ++i) { + // Materialise each stage as S: with XAD expression templates `h + c*k` + // is an expression type, not S, which would break template deduction of + // ff16_height_dt_crown_top's scalar. + const S k1 = ff16_height_dt_crown_top(p, h, light_E); + const S h2 = h + S(0.5 * dt) * k1; + const S k2 = ff16_height_dt_crown_top(p, h2, light_E); + const S h3 = h + S(0.5 * dt) * k2; + const S k3 = ff16_height_dt_crown_top(p, h3, light_E); + const S h4 = h + S(dt) * k3; + const S k4 = ff16_height_dt_crown_top(p, h4, light_E); + h = h + S(dt / 6.0) * (k1 + S(2.0) * k2 + S(2.0) * k3 + k4); + } + return h; +} + +// The five FF16 ODE states, as a value the trajectory integrator carries. +template +struct FF16State { + S height, mortality, fecundity, area_heartwood, mass_heartwood; +}; + +// Single-plant FULL-STATE demographic TRAJECTORY (#472 scope B, Milestone C): +// integrate all five FF16 states from y0 over n_steps fixed RK4 steps to age +// t_end, in a fixed crown-top light light_E, using ff16_compute_rates_crown_top. +// Generalises ff16_grow_height (height only) to the demographic vector, so +// reverse AD gives d(any emergent state at age t_end)/d(trait) -- e.g. lifetime +// fecundity (cumulative offspring) sensitivity, a calibration target through the +// whole demographic ODE. mortality_finite is the frozen pass-1 branch (see the +// rate kernel). Each RK4 stage is materialised as S so XAD expression templates +// don't break scalar deduction (same caveat as ff16_grow_height). +template +FF16State ff16_grow_demography(const FF16ProdPars& p, FF16State y, + S light_E, double t_end, int n_steps, + bool mortality_finite) { + const double dt = t_end / n_steps; + auto deriv = [&](const FF16State& s) -> FF16State { + const FF16Rates r = + ff16_compute_rates_crown_top(p, s.height, light_E, mortality_finite); + return FF16State{r.height_dt, r.mortality_dt, r.fecundity_dt, + r.area_heartwood_dt, r.mass_heartwood_dt}; + }; + auto axpy = [](const FF16State& a, S c, const FF16State& k) -> FF16State { + return FF16State{a.height + c * k.height, a.mortality + c * k.mortality, + a.fecundity + c * k.fecundity, + a.area_heartwood + c * k.area_heartwood, + a.mass_heartwood + c * k.mass_heartwood}; + }; + for (int i = 0; i < n_steps; ++i) { + const FF16State k1 = deriv(y); + const FF16State k2 = deriv(axpy(y, S(0.5 * dt), k1)); + const FF16State k3 = deriv(axpy(y, S(0.5 * dt), k2)); + const FF16State k4 = deriv(axpy(y, S(dt), k3)); + const S c = S(dt / 6.0); + y.height = y.height + c * (k1.height + S(2.0) * k2.height + S(2.0) * k3.height + k4.height); + y.mortality = y.mortality + c * (k1.mortality + S(2.0) * k2.mortality + S(2.0) * k3.mortality + k4.mortality); + y.fecundity = y.fecundity + c * (k1.fecundity + S(2.0) * k2.fecundity + S(2.0) * k3.fecundity + k4.fecundity); + y.area_heartwood = y.area_heartwood + c * (k1.area_heartwood + S(2.0) * k2.area_heartwood + S(2.0) * k3.area_heartwood + k4.area_heartwood); + y.mass_heartwood = y.mass_heartwood + c * (k1.mass_heartwood + S(2.0) * k2.mass_heartwood + S(2.0) * k3.mass_heartwood + k4.mass_heartwood); + } + return y; +} + +// Frozen-schedule forward-Euler replay of ONE cohort's full demographic state +// over a per-step crown-light schedule (#472 scope B, Milestone C -- the two-pass +// SCM replay primitive). `light[k]` is the (frozen, pass-1 double) light the +// cohort's crown reads at global replay step k; the cohort is born at step0 and +// integrated to the end of the schedule with forward Euler at fixed dt. Euler +// (not RK4) is deliberate: it reproduces the SCM's control.fixed_time_step +// integration EXACTLY, so replaying a fixed_time_step resident run is faithful. +// Templated on S, so reverse AD over a weighted sum of cohort outcomes +// J(theta) = sum_i w_i * f(replay_i) (w_i, light frozen from pass 1) +// gives d(emergent stand output)/d(trait), holding the resident light schedule +// fixed -- the legitimate "resident-light-frozen" gradient. The full resident +// self-shading gradient additionally makes light active (odelia #32 active-query +// spline / ff16_assimilation_deep_crown_replay), deferred. +template +FF16State ff16_replay_cohort(const FF16ProdPars& p, FF16State y, + double dt, const std::vector& light, + std::size_t step0, bool mortality_finite) { + for (std::size_t k = step0; k < light.size(); ++k) { + const FF16Rates r = + ff16_compute_rates_crown_top(p, y.height, S(light[k]), mortality_finite); + y.height = y.height + S(dt) * r.height_dt; + y.mortality = y.mortality + S(dt) * r.mortality_dt; + y.fecundity = y.fecundity + S(dt) * r.fecundity_dt; + y.area_heartwood = y.area_heartwood + S(dt) * r.area_heartwood_dt; + y.mass_heartwood = y.mass_heartwood + S(dt) * r.mass_heartwood_dt; + } + return y; +} + +// As ff16_replay_cohort, but the cohort reads its crown light ACTIVELY from a +// frozen resident profile at each step (#472 scope B, Milestone C). `crown_light` +// is a caller-supplied callable S -> S returning the light at the cohort's crown +// for its current (active) height; in the AD context the caller seeds it from the +// resident profile's value + slope (FF16_Environment::get_environment_at_height / +// get_environment_deriv_at_height) so d(light)/d(height) flows -- the within-cohort +// self-shading feedback (a taller cohort reads higher in the canopy -> more light +// -> faster growth) that the frozen per-step light of the plain overload omits. +// The profile KNOTS stay frozen double (resident held fixed); making them active +// is the full self-shading gradient (odelia #32 active-knot spline). LightFn is a +// template so XAD never enters this header (mirrors ff16_assimilation_deep_crown_replay). +template +FF16State ff16_replay_cohort_active_light(const FF16ProdPars& p, FF16State y, + double dt, LightFn&& crown_light, + int n_steps, bool mortality_finite) { + for (int k = 0; k < n_steps; ++k) { + const S light_E = crown_light(y.height); + const FF16Rates r = + ff16_compute_rates_crown_top(p, y.height, light_E, mortality_finite); + y.height = y.height + S(dt) * r.height_dt; + y.mortality = y.mortality + S(dt) * r.mortality_dt; + y.fecundity = y.fecundity + S(dt) * r.fecundity_dt; + y.area_heartwood = y.area_heartwood + S(dt) * r.area_heartwood_dt; + y.mass_heartwood = y.mass_heartwood + S(dt) * r.mass_heartwood_dt; + } + return y; +} + +// PRODUCTION two-pass replay primitive: a single cohort's full demographic state +// integrated with the SAME adaptive Cash-Karp RKCK scheme the live SCM used +// (#472 scope B, Milestone C -- the FAITHFUL replacement for the Euler +// ff16_replay_cohort). Forward Euler mirrors only the non-default +// control.fixed_time_step path; the real SCM integrates with odelia's embedded +// 4/5 RKCK (ode_step.hpp), and the mutant-fitness replay (run_mutant -> +// advance_fixed -> step_to) re-uses that SAME stepper over the resident's pinned +// step times, swapping in a FROZEN per-RK-stage environment (environment_history +// [step][stage], 6 stages/step). This kernel is that path lifted to the scalar S. +// +// step_h[n] = the resident's actual adaptive step size for global step n (= +// step_history[n+1]-step_history[n]); the cohort is integrated from global step +// `step0` (its birth step) to the end of the schedule. The per-stage crown light +// is supplied by `crown_light(n, stage, height) -> S` so XAD/odelia stay out of +// this header (mirrors ff16_replay_cohort_active_light). The caller wires it to +// the FROZEN resident env: in the AD context it seeds value + slope from +// FF16_Environment::get_environment_at_height / get_environment_deriv_at_height +// so d(light)/d(height) flows -- the mutant-through-frozen-canopy feedback (a +// taller focal cohort reads higher in the resident profile). Stage codes: +// stage 0 -> k1 env = the env at the step START (environment_history[n-1][5], +// or the birth env for n==step0); recomputing k1 against +// it is numerically identical to the solver's FSAL reuse; +// stage 1..5 -> the envs for the k2..k6 derivs (environment_history[n][0..4]). +// The 6th cached env (environment_history[n][5], the solver's dydt_out stage) is +// re-used as the next step's stage-0 env, exactly as first_same_as_last does. The +// y-update uses k1,k3,k4,k6 (c2==c5==0), matching ode_step.hpp::step line-for-line. +// Generic Cash-Karp RKCK driver over the FROZEN resident schedule, shared by the +// demographic replay (ff16_replay_cohort_rkck) and the lifetime-offspring replay +// (ff16_replay_cohort_offspring_rkck). The integrator logic -- the GSL Cash-Karp +// tableau, the FSAL stage-0 reuse, the c2==c5==0 final sum -- lives here ONCE. +// Callers supply the state type and its two operations: +// deriv(state, n, stage) -> State : the state-derivative at RK `stage` (0..5) of +// global step n. stage 0 is the FSAL k1 (evaluated at the step START); +// stages 1..5 are the k2..k6 derivs. Callers map (n, stage) to the frozen +// per-RK-stage resident environment (see ff16_replay_cohort_rkck). +// axpy(a, c, k) -> State : a + c*k with c a double RK coefficient; each component +// MUST be materialised into the scalar type in the returned State (brace- +// init), so XAD expression templates never escape with dangling references. +// The y-update is the 5th-order sum y += h*(c1 k1 + c3 k3 + c4 k4 + c6 k6) written +// as an axpy chain (c2==c5==0), matching odelia ode_step.hpp::step. +template +State ff16_cashkarp_replay(State y, const std::vector& step_h, + std::size_t step0, DerivFn&& deriv, AxpyFn&& axpy) { + // Cash-Karp coefficients, identical to odelia::ode::Step (from GSL). + const double b21 = 1.0 / 5.0; + const double b3[2] = {3.0 / 40.0, 9.0 / 40.0}; + const double b4[3] = {0.3, -0.9, 1.2}; + const double b5[4] = {-11.0 / 54.0, 2.5, -70.0 / 27.0, 35.0 / 27.0}; + const double b6[5] = {1631.0 / 55296.0, 175.0 / 512.0, 575.0 / 13824.0, + 44275.0 / 110592.0, 253.0 / 4096.0}; + const double c1 = 37.0 / 378.0, c3 = 250.0 / 621.0, + c4 = 125.0 / 594.0, c6 = 512.0 / 1771.0; + + for (std::size_t n = step0; n < step_h.size(); ++n) { + const double h = step_h[n]; + const State k1 = deriv(y, n, 0); + const State k2 = deriv(axpy(y, b21 * h, k1), n, 1); + State y3 = axpy(y, h * b3[0], k1); y3 = axpy(y3, h * b3[1], k2); + const State k3 = deriv(y3, n, 2); + State y4 = axpy(y, h * b4[0], k1); + y4 = axpy(y4, h * b4[1], k2); y4 = axpy(y4, h * b4[2], k3); + const State k4 = deriv(y4, n, 3); + State y5 = axpy(y, h * b5[0], k1); + y5 = axpy(y5, h * b5[1], k2); y5 = axpy(y5, h * b5[2], k3); y5 = axpy(y5, h * b5[3], k4); + const State k5 = deriv(y5, n, 4); + State y6 = axpy(y, h * b6[0], k1); + y6 = axpy(y6, h * b6[1], k2); y6 = axpy(y6, h * b6[2], k3); + y6 = axpy(y6, h * b6[3], k4); y6 = axpy(y6, h * b6[4], k5); + const State k6 = deriv(y6, n, 5); + y = axpy(axpy(axpy(axpy(y, h * c1, k1), h * c3, k3), h * c4, k4), h * c6, k6); + } + return y; +} + +template +FF16State ff16_replay_cohort_rkck(const FF16ProdPars& p, FF16State y, + const std::vector& step_h, + std::size_t step0, + StageLightFn&& crown_light, + bool mortality_finite) { + // 5-state FF16 derivative at a trial state, reading the frozen stage env. + auto deriv = [&](const FF16State& s, std::size_t n, int stage) -> FF16State { + const S light_E = crown_light(n, stage, s.height); + const FF16Rates r = + ff16_compute_rates_crown_top(p, s.height, light_E, mortality_finite); + return FF16State{r.height_dt, r.mortality_dt, r.fecundity_dt, + r.area_heartwood_dt, r.mass_heartwood_dt}; + }; + // a + c*k, each component materialised as S in the brace-init. + auto axpy = [](const FF16State& a, double c, const FF16State& k) -> FF16State { + return FF16State{a.height + c * k.height, a.mortality + c * k.mortality, + a.fecundity + c * k.fecundity, + a.area_heartwood + c * k.area_heartwood, + a.mass_heartwood + c * k.mass_heartwood}; + }; + return ff16_cashkarp_replay(y, step_h, step0, deriv, axpy); +} + +// State for the lifetime-offspring replay: the 5 FF16 states + the cumulative +// survival-weighted offspring (the SCM's offspring_produced_survival_weighted). +template +struct FF16LifeState { + FF16State demog; + S offspring; +}; + +// Lifetime survival-weighted offspring replay (#472 scope B): augments the +// demographic replay with a 6th accumulator mirroring Node::compute_rates, +// d(offspring)/dt = fecundity_dt * exp(-mortality) * surv_weight(n, stage), +// integrated with the SAME Cash-Karp driver. `surv_weight(n, stage) -> double` +// supplies the FROZEN pr_patch_survival(t_stage)/pr_patch_survival_at_birth at the +// step's RK stage time; set y.demog.mortality = -log(establishment_probability) at +// birth (the node's initial condition) before calling. The stand's emergent +// offspring_production is then the node-spacing trapezium of +// offspring * patch_density_at_birth * S_D * birth_rate +// over the cohorts (a frozen, linear post-weighting), so one reverse sweep of that +// sum gives d(offspring_production)/d(trait). crown_light/surv_weight are callables +// so XAD/odelia stay out of this header. +template +FF16LifeState ff16_replay_cohort_offspring_rkck( + const FF16ProdPars& p, FF16LifeState y, const std::vector& step_h, + std::size_t step0, StageLightFn&& crown_light, SurvFn&& surv_weight, + bool mortality_finite) { + using std::exp; // XAD provides exp for active types via ADL + auto deriv = [&](const FF16LifeState& s, std::size_t n, int stage) -> FF16LifeState { + const S light_E = crown_light(n, stage, s.demog.height); + const FF16Rates r = + ff16_compute_rates_crown_top(p, s.demog.height, light_E, mortality_finite); + const S off_dt = r.fecundity_dt * exp(-s.demog.mortality) * S(surv_weight(n, stage)); + return FF16LifeState{FF16State{r.height_dt, r.mortality_dt, r.fecundity_dt, + r.area_heartwood_dt, r.mass_heartwood_dt}, + off_dt}; + }; + auto axpy = [](const FF16LifeState& a, double c, const FF16LifeState& k) -> FF16LifeState { + return FF16LifeState{ + FF16State{a.demog.height + c * k.demog.height, + a.demog.mortality + c * k.demog.mortality, + a.demog.fecundity + c * k.demog.fecundity, + a.demog.area_heartwood + c * k.demog.area_heartwood, + a.demog.mass_heartwood + c * k.demog.mass_heartwood}, + a.offspring + c * k.offspring}; + }; + return ff16_cashkarp_replay(y, step_h, step0, deriv, axpy); +} + +// Resident light availability E(z) = exp( - sum_i density_i * k_I * area_leaf_i * +// Q(z/h_i) ) at height z from a FROZEN stand (heights/densities are pass-1 +// doubles), ACTIVE in the traits through each cohort's area_leaf [eqn 2] +// (#472 scope B, Milestone C -- the resident self-shading coupling). Q is the +// deep/Yokozawa leaf-area-above (1 - u^eta)^2 with eta a fixed double (other +// shading variants swap Q); contributions vanish above each plant's top. Beer's +// law E = exp(-projected leaf area), matching FF16_Environment::compute_environment. +// Evaluated at FROZEN knot positions z_k to fill an active-VALUE light spline +// (odelia basic_interpolator), this is what makes a resident emergent output +// differentiable w.r.t. a trait THROUGH the self-shaded light profile -- the +// full self-shading gradient (vs the frozen-knot active-query of +// ff16_replay_cohort_active_light). Heights frozen here (a fixed stand census); +// coupling growth back in is the live two-pass replay. +template +S ff16_resident_light_at(double z, S a_l1, S a_l2, double k_I, double eta, + const std::vector& height, + const std::vector& density) { + using std::pow; using std::exp; + S L = S(0.0); + for (std::size_t i = 0; i < height.size(); ++i) { + if (z >= height[i]) continue; // no leaf area above the plant's crown + const double u = z / height[i]; + const double one_minus = 1.0 - pow(u, eta); + const double Q = one_minus * one_minus; // Yokozawa leaf-area-above + L += S(density[i] * k_I) * ff16_area_leaf(a_l1, a_l2, S(height[i])) * S(Q); + } + return exp(-L); +} + } // namespace plant #endif diff --git a/inst/include/plant/models/ff16_strategy.h b/inst/include/plant/models/ff16_strategy.h index 29e00a62..5b500f72 100644 --- a/inst/include/plant/models/ff16_strategy.h +++ b/inst/include/plant/models/ff16_strategy.h @@ -15,25 +15,29 @@ namespace plant { // (so R access is `s$pars$lma`). Derived/precomputed quantities (eta_c, // height_0, canopy_shape, ...) are NOT here -- they are outputs of // prepare_strategy() and stay as plain members on the strategy. -struct FF16_Pars { +// Templated on the scalar S (#472 scope B / #537, Milestone C) so traits can be +// AD active types for reverse-mode calibration. `FF16_Pars` (alias below) pins +// S = double, leaving the R/RcppR6 interface and every existing use unchanged. +template +struct basic_FF16_Pars { // * Core traits - double lma = 0.1978791; // Leaf mass per area [kg / m2] - double rho = 608.0; // Wood density [kg/m3] - double hmat = 16.5958691; // Height at maturation [m] - double omega = 3.8e-5; // Seed mass [kg] + S lma = 0.1978791; // Leaf mass per area [kg / m2] + S rho = 608.0; // Wood density [kg/m3] + S hmat = 16.5958691; // Height at maturation [m] + S omega = 3.8e-5; // Seed mass [kg] // * Individual allometry // Canopy shape parameter - double eta = 12.0; // [dimensionless] + S eta = 12.0; // [dimensionless] // Sapwood area per leaf area // Ratio sapwood area area to leaf area - double theta = 1.0/4669; // [dimensionless] + S theta = 1.0/4669; // [dimensionless] // Height - leaf mass scaling - double a_l1 = 5.44; // height with 1m2 leaf [m] - double a_l2 = 0.306; // dimensionless scaling of height with leaf area + S a_l1 = 5.44; // height with 1m2 leaf [m] + S a_l2 = 0.306; // dimensionless scaling of height with leaf area // Root mass per leaf area - double a_r1 = 0.07; //[kg / m] + S a_r1 = 0.07; //[kg / m] // Ratio of bark area : sapwood area - double a_b1 = 0.17; // [dimensionless] + S a_b1 = 0.17; // [dimensionless] // * Production // Ratio of leaf dark respiration to leaf mass [mol CO2 / yr / kg] @@ -41,63 +45,66 @@ struct FF16_Pars { // / [kg(leaf) / m2 ] | / (0.1978791) | lma // Hard coded in value of lma here so that this value doesn't change // if that trait changes above. - double r_l = 39.27 / 0.1978791; + S r_l = 39.27 / 0.1978791; // Root respiration per mass [mol CO2 / yr / kg] - double r_r = 217.0; + S r_r = 217.0; // Sapwood respiration per stem mass [mol CO2 / yr / kg] // = respiration per volume [mol CO2 / m3 / yr] // / wood density [kg/m3] - double r_s = 4012.0 / 608.0; + S r_s = 4012.0 / 608.0; // Bark respiration per stem mass // assumed to be twice rate of sapwood // (NOTE that there is a re-parametrisation here relative to the paper // -- r_b is defined (new) as 2*r_s, whereas the paper assumes a // fixed multiplication by 2) - double r_b = 2.0 * r_s; + S r_b = 2.0 * r_s; // Carbon conversion parameter - double a_y = 0.7; + S a_y = 0.7; // Constant converting assimilated CO2 to dry mass [kg / mol] // (12E-3 / 0.49) - double a_bio = 2.45e-2; + S a_bio = 2.45e-2; // Leaf turnover [/yr] - double k_l = 0.4565855; + S k_l = 0.4565855; // Bark turnover [/yr] - double k_b = 0.2; + S k_b = 0.2; // Sapwood turnover [/yr] - double k_s = 0.2; + S k_s = 0.2; // Root turnover [/yr] - double k_r = 1.0; + S k_r = 1.0; // Parameters of the hyperbola for annual LRC - double a_p1 = 151.177775377968; // [mol CO2 / yr / m2] - double a_p2 = 0.204716166503633; // [dimensionless] + S a_p1 = 151.177775377968; // [mol CO2 / yr / m2] + S a_p2 = 0.204716166503633; // [dimensionless] // * Seed production // Accessory cost of reproduction - double a_f3 = 3.0 * 3.8e-5; // [kg per seed] + S a_f3 = 3.0 * 3.8e-5; // [kg per seed] // Maximum allocation to reproduction - double a_f1 = 1.0; //[dimensionless] + S a_f1 = 1.0; //[dimensionless] // Size range across which individuals mature - double a_f2 = 50; // [dimensionless] + S a_f2 = 50; // [dimensionless] // * Mortality parameters // Probability of survival during dispersal - double S_D = 0.25; // [dimensionless] + S S_D = 0.25; // [dimensionless] // Parameter for seedling survival - double a_d0 = 0.1; //[kg / yr / m2] + S a_d0 = 0.1; //[kg / yr / m2] // Baseline for intrinsic mortality - double d_I = 0.01; // [ / yr] + S d_I = 0.01; // [ / yr] // Baseline rate for growth-related mortality - double a_dG1 = 5.5; // [ / yr] + S a_dG1 = 5.5; // [ / yr] // Risk coefficient for dry mass production (per area) - double a_dG2 = 20.0;// [yr m2 / kg ] + S a_dG2 = 20.0;// [yr m2 / kg ] // Germination - double recruitment_decay = 0.0; + S recruitment_decay = 0.0; // * Light capture parameters - double k_I = 0.5; + S k_I = 0.5; }; +// Default parameter set used by FF16_Strategy and the R/RcppR6 interface. +using FF16_Pars = basic_FF16_Pars; + class FF16_Strategy: public Strategy { public: typedef std::shared_ptr ptr; @@ -167,6 +174,15 @@ class FF16_Strategy: public Strategy { // reference-comparison test. return ff16_area_leaf(pars.a_l1, pars.a_l2, height); } + // Scalar-templated overload (#472 scope B / #537, Milestone C): area_leaf with + // an AD-active height (the live ODE state), keeping the allometry pars double + // (S(pars.*) lifts them). The non-template double overload above still wins for + // a double argument, so the existing hot path is unchanged; only ad heights + // (e.g. update_dependent_aux on Individual<...,ad>) select this. + template + S area_leaf(S height) const { + return ff16_area_leaf(S(pars.a_l1), S(pars.a_l2), height); + } // [eqn 1] mass_leaf (inverse of [eqn 2]) double mass_leaf(double area_leaf) const; @@ -204,11 +220,18 @@ class FF16_Strategy: public Strategy { // Inline (header): called per state-set / ODE-state update from templated // Individual code, so inlining avoids a cross-TU call (no LTO build) // and lets the now-inline area_leaf fold in. - void update_dependent_aux(const int index, Internals& vars) { + // Scalar-templated on the Internals' value type S (#472 scope B / #537, + // Milestone C) so a Individual<...,ad> can be constructed and have its height + // state set: the dependent aux (competition_effect = area_leaf, height_inverse) + // then carry the active scalar. The non-template path is gone, but S is + // deduced as double for every existing caller (Individual<...,double>), so the + // double codegen and the FF16 reference test are unchanged. + template + void update_dependent_aux(const int index, basic_internals& vars) { if (index == HEIGHT_INDEX) { - double height = vars.state(HEIGHT_INDEX); + S height = vars.state(HEIGHT_INDEX); vars.set_aux(COMPETITION_EFFECT_AUX_INDEX, area_leaf(height)); - vars.set_aux(HEIGHT_INVERSE_AUX_INDEX, 1.0 / height); + vars.set_aux(HEIGHT_INVERSE_AUX_INDEX, S(1.0) / height); } } @@ -388,6 +411,37 @@ class FF16_Strategy: public Strategy { // Set constants within FF16_Strategy void prepare_strategy(); + // The net-production kernel's parameter set, gathered from this (prepared) + // strategy's pars + derived eta_c (#472 scope B, Milestone C). Bridges a live, + // prepared FF16_Strategy to the scalar-templated AD kernel: lift the result to + // FF16ProdPars (registering the trait of interest as a tape input) to + // get reverse-mode trait gradients of net production from the real model + // configuration rather than hand-supplied numbers. + FF16ProdPars prod_pars() const { + FF16ProdPars p; + p.lma = pars.lma; p.rho = pars.rho; p.theta = pars.theta; + p.a_b1 = pars.a_b1; p.a_r1 = pars.a_r1; p.eta_c = eta_c; + p.a_p1 = pars.a_p1; p.a_p2 = pars.a_p2; + p.r_l = pars.r_l; p.r_s = pars.r_s; p.r_b = pars.r_b; p.r_r = pars.r_r; + p.k_l = pars.k_l; p.k_b = pars.k_b; p.k_s = pars.k_s; p.k_r = pars.k_r; + p.a_bio = pars.a_bio; p.a_y = pars.a_y; + p.a_l1 = pars.a_l1; p.a_l2 = pars.a_l2; + p.a_f1 = pars.a_f1; p.a_f2 = pars.a_f2; p.hmat = pars.hmat; + // Demographic rate params for the full ff16_compute_rates_* fill (Milestone C). + p.omega = pars.omega; p.a_f3 = pars.a_f3; + p.d_I = pars.d_I; p.a_dG1 = pars.a_dG1; p.a_dG2 = pars.a_dG2; + return p; + } + + // Exact d(dheight/dt)/d(height) at the given height in environment `env`, via + // forward-mode AD over the scalar-templated growth kernel (#537 A1; the + // gradient Node::growth_rate_gradient currently obtains by finite difference). + // Crown-top assimilation: the light the crown reads is taken at the operating + // height, so in a fixed environment it is exact. Defined in the .cpp (XAD + // include), mirroring Leaf::dprofit_droot_collar_psi. + double growth_rate_gradient_height_ad(double height, + const FF16_Environment& env); + // Birth height of a (germinated) seed. Strategy-agnostic accessor used by // the templated Individual; here height_0 is derived in prepare_strategy(). double initial_height() const { return height_0; } diff --git a/inst/include/plant/models/tf24_production_kernel.h b/inst/include/plant/models/tf24_production_kernel.h new file mode 100644 index 00000000..ac95cc2b --- /dev/null +++ b/inst/include/plant/models/tf24_production_kernel.h @@ -0,0 +1,240 @@ +// -*-c++-*- +#ifndef PLANT_PLANT_TF24_PRODUCTION_KERNEL_H_ +#define PLANT_PLANT_TF24_PRODUCTION_KERNEL_H_ + +#include // std::pow / std::exp; XAD provides these for active types via ADL + +// Scalar-templated core of the TF24 net-mass-production + demographic rate chain +// (#472 scope B / traitecoevo/plant#537, Phase F1-full). These pieces are the +// TF24 analogues of ff16_production_kernel.h. The MASS CASCADE, respiration, +// turnover and the whole demographic rate fill (growth / fecundity / heartwood / +// mortality downstream of net production) are ALGEBRAICALLY IDENTICAL to FF16; +// only the assimilation term differs -- TF24 forms it from the OPTIMISED leaf +// profit (assim = profit * area_leaf * conv) rather than the FF16 light-response +// hyperbola. They are the SINGLE SOURCE OF TRUTH: TF24_Strategy's double rate +// methods delegate to them (so the existing TF24 reference test validates +// faithfulness, bit-identical), and the AD calibration path instantiates them +// with an active scalar so net / the demographic rates differentiate w.r.t. +// traits by forward- or reverse-mode AD with no special handling. +// +// The leaf optimisation itself (profit*, a max over collar potential nesting a +// psi_stem->ci root-find) is NOT in this header: it stays in the leaf submodel, +// and its trait sensitivities enter the AD path as the Leaf::dprofit_d* numbers +// injected first-order into an active `profit` (the #539 IFT / FF16 height_0 +// injection pattern). This kernel takes `profit` (or `assimilation`) as a given +// and carries it -- and every mass-cascade trait -- through to the five ODE rates. + +namespace plant { + +// Unit conversion folding leaf-area assimilation (umol CO2 m^-2 s^-1) to +// canopy-level yearly assimilation (mol yr^-1): 60*60 s/h * 12 h/day daylight * +// 365 d/yr / 1e6 umol/mol. Recurs in TF24_Strategy::net_mass_production_dt. +constexpr double tf24_assimilation_conv = 60.0 * 60.0 * 12.0 * 365.0 / 1e6; + +// [eqn 2] Leaf area as a function of height (inverse of [eqn 3]). Carries the +// allometric traits a_l1, a_l2 -- the entry point through which an allometric +// trait reaches both the mass cascade and area_leaf itself. +template +S tf24_area_leaf(S a_l1, S a_l2, S height) { + using std::pow; + return pow(height / a_l1, 1.0 / a_l2); +} + +// [eqn 13] Total maintenance respiration (linear in the mass cascade). +template +S tf24_respiration(S mass_leaf, S mass_sapwood, S mass_bark, S mass_root, + S r_l, S r_s, S r_b, S r_r) { + return r_l * mass_leaf + r_b * mass_bark + r_s * mass_sapwood + r_r * mass_root; +} + +// [eqn 14] Total turnover. +template +S tf24_turnover(S mass_leaf, S mass_bark, S mass_sapwood, S mass_root, + S k_l, S k_b, S k_s, S k_r) { + return k_l * mass_leaf + k_b * mass_bark + k_s * mass_sapwood + k_r * mass_root; +} + +// [eqn 15] Net production from assimilation/respiration/turnover. +template +S tf24_net_production_A(S a_bio, S a_y, S assimilation, S respiration, S turnover) { + return a_bio * a_y * (assimilation - respiration) - turnover; +} + +// The TF24 parameters the net-production chain + demographic rate fill read, plus +// the prepare_strategy()-derived eta_c. The hot double path uses the per-piece +// functions directly with pars.* members; the all-in-one functions below (the AD +// calibration entry points) read these. NOTE: unlike FF16ProdPars there are no +// a_p1/a_p2 fields -- TF24 assimilation comes from the optimised leaf profit, not +// the light-response hyperbola. +template +struct TF24ProdPars { + S lma, rho, theta, a_b1, a_r1, eta_c; + S r_l, r_s, r_b, r_r; + S k_l, k_b, k_s, k_r; + S a_bio, a_y; + // Allometry + reproduction-allocation parameters for the growth/fecundity rates. + S a_l1, a_l2; // height <-> leaf-area allometry [eqn 2/3] + S a_f1, a_f2, hmat; // reproduction-allocation logistic [eqn 16] + S omega, a_f3; // seed mass + accessory reproduction cost (fecundity_dt) [eqn 17] + S d_I, a_dG1, a_dG2; // mortality: growth-independent + growth-dependent [eqn 21] +}; + +// Mass cascade -> respiration/turnover -> net production, GIVEN the assimilation +// rate. Shared by every assimilation variant (crown-centre, mean-light, +// deep-crown), which differ only in how `assimilation` is formed. Mirrors +// TF24_Strategy::net_mass_production_dt's tail exactly (FF16-identical algebra). +template +S tf24_net_from_components(const TF24ProdPars& p, S height, S area_leaf, + S assimilation) { + const S mass_leaf = area_leaf * p.lma; + const S area_sapwood = area_leaf * p.theta; + const S mass_sapwood = area_sapwood * height * p.eta_c * p.rho; + const S area_bark = p.a_b1 * area_leaf * p.theta; + const S mass_bark = area_bark * height * p.eta_c * p.rho; + const S mass_root = p.a_r1 * area_leaf; + + const S respiration = tf24_respiration(mass_leaf, mass_sapwood, mass_bark, mass_root, + p.r_l, p.r_s, p.r_b, p.r_r); + const S turnover = tf24_turnover(mass_leaf, mass_bark, mass_sapwood, mass_root, + p.k_l, p.k_b, p.k_s, p.k_r); + return tf24_net_production_A(p.a_bio, p.a_y, assimilation, respiration, turnover); +} + +// Net production from the OPTIMISED leaf profit: assim = profit * area_leaf * conv, +// then the shared mass cascade. `profit` is the leaf submodel's profit_ (a max over +// collar potential); in the AD path it is an active scalar holding the injected +// Leaf::dprofit_d* sensitivities. Mirrors net_mass_production_dt's last lines. +template +S tf24_net_mass_production(const TF24ProdPars& p, S height, S area_leaf, S profit) { + const S assimilation = profit * area_leaf * S(tf24_assimilation_conv); + return tf24_net_from_components(p, height, area_leaf, assimilation); +} + +// --------------------------------------------------------------------------- +// Demographic rate pieces (the rates downstream of net production). Mirror +// TF24_Strategy::{fraction_allocation_reproduction, fraction_allocation_growth, +// fecundity_dt, dheight_darea_leaf, dmass_*_darea_leaf, darea_leaf_dmass_live, +// mortality_growth_independent_dt, mortality_growth_dependent_dt} and the +// dheight/dt + heartwood assembly in compute_rates. Elementary arithmetic, so the +// whole demographic rate fill differentiates w.r.t. a trait. +// --------------------------------------------------------------------------- + +// [eqn 16] Fraction of production allocated to reproduction (logistic in height). +template +S tf24_fraction_allocation_reproduction(S a_f1, S a_f2, S hmat, S height) { + using std::exp; + return a_f1 / (1.0 + exp(a_f2 * (1.0 - height / hmat))); +} + +// [eqn 16] Fraction of production allocated to growth = 1 - reproduction. +template +S tf24_fraction_allocation_growth(S a_f1, S a_f2, S hmat, S height) { + return 1.0 - tf24_fraction_allocation_reproduction(a_f1, a_f2, hmat, height); +} + +// [eqn 17] Rate of offspring production. +template +S tf24_fecundity_dt(S net, S fraction_allocation_reproduction, S omega, S a_f3) { + return net * fraction_allocation_reproduction / (omega + a_f3); +} + +// d(height)/d(area_leaf): derivative of the [eqn 2] allometry. +template +S tf24_dheight_darea_leaf(S a_l1, S a_l2, S area_leaf) { + using std::pow; + return a_l1 * a_l2 * pow(area_leaf, a_l2 - 1.0); +} + +// d(area_leaf)/d(mass_live): reciprocal of the summed per-component mass +// derivatives (leaf + sapwood + bark + root) w.r.t. area_leaf. Matches the +// TF24_Strategy::dmass_*_darea_leaf sum exactly. +template +S tf24_darea_leaf_dmass_live(const TF24ProdPars& p, S area_leaf) { + using std::pow; + const S dmass_leaf = p.lma; // d(area_leaf*lma) + const S dmass_sapwood = p.rho * p.eta_c * p.a_l1 * p.theta * + (p.a_l2 + 1.0) * pow(area_leaf, p.a_l2); + const S dmass_bark = p.a_b1 * dmass_sapwood; + const S dmass_root = p.a_r1; + return 1.0 / (dmass_leaf + dmass_sapwood + dmass_bark + dmass_root); +} + +// dheight/dt given net production (the compute_rates growth assembly): returns 0 +// when net is non-positive (the growth clamp), else dheight_darea_leaf * +// area_leaf_dt with area_leaf_dt = net * frac_growth * darea_leaf_dmass_live. +template +S tf24_height_dt_from_net(const TF24ProdPars& p, S height, S area_leaf, S net) { + if (net <= 0.0) return S(0.0); + const S frac_growth = tf24_fraction_allocation_growth(p.a_f1, p.a_f2, p.hmat, height); + const S darea_dmass = tf24_darea_leaf_dmass_live(p, area_leaf); + const S area_leaf_dt = net * frac_growth * darea_dmass; + return tf24_dheight_darea_leaf(p.a_l1, p.a_l2, area_leaf) * area_leaf_dt; +} + +// [eqn 21] growth-independent mortality (intrinsic baseline). +template +S tf24_mortality_growth_independent_dt(S d_I) { return d_I; } + +// [eqn 21] growth-dependent mortality; productivity_area = net / area_leaf. +template +S tf24_mortality_growth_dependent_dt(S a_dG1, S a_dG2, S productivity_area) { + using std::exp; + return a_dG1 * exp(-a_dG2 * productivity_area); +} + +// The five ODE state rates TF24_Strategy::compute_rates writes, plus the net +// production aux. Scalar-templated (#472 scope B) so the whole demographic rate +// fill differentiates w.r.t. a trait. +template +struct TF24Rates { + S net_mass_production_dt; + S height_dt; + S fecundity_dt; + S area_heartwood_dt; + S mass_heartwood_dt; + S mortality_dt; +}; + +// Full TF24 compute_rates fill GIVEN the net production (however `net` was +// formed -- crown-centre / mean-light / deep-crown). Mirrors +// TF24_Strategy::compute_rates downstream of net EXACTLY: the net>0 growth clamp +// gates the growth/fecundity/heartwood rates, and mortality is the [eqn 21] +// growth-independent + growth-dependent sum (productivity = net/area_leaf). +// `mortality_finite` is the frozen util::is_finite(cumulative mortality) branch -- +// a pass-1 (double) control-flow decision, passed in so the taped replay is +// branch-free (it never differentiates the is_finite test). This is the part of +// compute_rates the demographic kernel owns; the leaf optimisation upstream +// supplies `net`. +template +TF24Rates tf24_compute_rates_from_net(const TF24ProdPars& p, S height, + S area_leaf, S net, + bool mortality_finite) { + TF24Rates r; + r.net_mass_production_dt = net; + if (net > 0.0) { + const S frac_repro = tf24_fraction_allocation_reproduction(p.a_f1, p.a_f2, + p.hmat, height); + r.height_dt = tf24_height_dt_from_net(p, height, area_leaf, net); + r.fecundity_dt = tf24_fecundity_dt(net, frac_repro, p.omega, p.a_f3); + const S area_sapwood = area_leaf * p.theta; // [eqn 4] + r.area_heartwood_dt = p.k_s * area_sapwood; // turnover of sapwood area + const S mass_sapwood = area_sapwood * height * p.eta_c * p.rho; + r.mass_heartwood_dt = p.k_s * mass_sapwood; // turnover_sapwood(mass) + } else { + r.height_dt = S(0.0); r.fecundity_dt = S(0.0); + r.area_heartwood_dt = S(0.0); r.mass_heartwood_dt = S(0.0); + } + // [eqn 21] instantaneous mortality rate; productivity_area = net / area_leaf. + if (mortality_finite) { + const S productivity_area = net / area_leaf; + r.mortality_dt = tf24_mortality_growth_independent_dt(p.d_I) + + tf24_mortality_growth_dependent_dt(p.a_dG1, p.a_dG2, productivity_area); + } else { + r.mortality_dt = S(0.0); + } + return r; +} + +} // namespace plant + +#endif diff --git a/inst/include/plant/models/tf24_strategy.h b/inst/include/plant/models/tf24_strategy.h index 945ba543..28f82691 100644 --- a/inst/include/plant/models/tf24_strategy.h +++ b/inst/include/plant/models/tf24_strategy.h @@ -8,6 +8,7 @@ #include #include #include // ShadingModel +#include // scalar-templated net + rate kernel namespace plant { @@ -320,6 +321,25 @@ class TF24_Strategy: public Strategy { // the templated Individual; here height_0 is derived in prepare_strategy(). double initial_height() const { return height_0; } + // Gather the net-production + demographic-rate kernel parameters from the live + // strategy configuration (pars + the prepare_strategy()-derived eta_c). The + // double path uses the per-piece kernel functions directly; this is the AD + // calibration entry point (lift to -> the demographic rates differentiate + // w.r.t. traits). Mirrors FF16_Strategy::prod_pars (#472 scope B, Phase F1). + TF24ProdPars prod_pars() const { + TF24ProdPars p; + p.lma = pars.lma; p.rho = pars.rho; p.theta = pars.theta; + p.a_b1 = pars.a_b1; p.a_r1 = pars.a_r1; p.eta_c = eta_c; + p.r_l = pars.r_l; p.r_s = pars.r_s; p.r_b = pars.r_b; p.r_r = pars.r_r; + p.k_l = pars.k_l; p.k_b = pars.k_b; p.k_s = pars.k_s; p.k_r = pars.k_r; + p.a_bio = pars.a_bio; p.a_y = pars.a_y; + p.a_l1 = pars.a_l1; p.a_l2 = pars.a_l2; + p.a_f1 = pars.a_f1; p.a_f2 = pars.a_f2; p.hmat = pars.hmat; + p.omega = pars.omega; p.a_f3 = pars.a_f3; + p.d_I = pars.d_I; p.a_dG1 = pars.a_dG1; p.a_dG2 = pars.a_dG2; + return p; + } + // Crown shading model, resolved once from control.shading_model in // prepare_strategy(). TF24 supports deep-crown, mean-light (its default) // and crown-centre; PPA is not available for TF24. diff --git a/inst/include/plant/models/tf24f_strategy.h b/inst/include/plant/models/tf24f_strategy.h index da53eadc..903e7f8a 100644 --- a/inst/include/plant/models/tf24f_strategy.h +++ b/inst/include/plant/models/tf24f_strategy.h @@ -43,6 +43,18 @@ class TF24f_Strategy : public TF24_Strategy { void compute_rates(const TF24_Environment& environment, Internals& vars); + // Exact d(dheight/dt)/d(height) at the tracked collar (#472 scope B / #537 A1), + // so run_scm's log_density (and the census stand gradients) use an analytic + // growth-rate gradient instead of Node::growth_rate_gradient's 1e-6 backward FD + // -- which, on TF24's hydraulic-leaf height_dt, amplifies leaf-solve noise. At + // the fixed tracked collar, height enters profit through kmax (= K_s*theta/ + // (h*eta_c)), the crown-centre light PPFD = k_I*L(h*eta_c)*PPFD_top, and (weakly, + // below the rooting-depth clamp) the soil uptake E_up; the demographic kernel + // adds the explicit-height allometry/cascade terms. Returns NA_REAL on an + // infeasible leaf / non-smooth light so Node falls back to the finite difference. + double growth_rate_gradient_height_ad(double height, + const TF24_Environment& environment); + // Override the leaf solve: instead of optimising the root-collar psi, evaluate // the leaf at the tracked state and finite-difference the profit gradient // (left in dprofit_dpsi_ for compute_rates to turn into the state's rate). diff --git a/inst/include/plant/node.h b/inst/include/plant/node.h index 12f3c612..11b54dfa 100644 --- a/inst/include/plant/node.h +++ b/inst/include/plant/node.h @@ -10,12 +10,20 @@ namespace plant { -template +// Templated on the scalar type S (#472 scope B / #537, Milestone C) to mirror +// Individual: S defaults to double so every existing `Node` is +// `Node` and bit-identical. The *demographic* bookkeeping +// (log_density, density, fecundity, offspring) stays double -- a Node<...,ad> +// is intentionally MIXED: only the individual's physiological state carries the +// active scalar, and (like Individual) members are compiled per-member-on-use, +// so the double-bound ODE-iterator / demographic methods stay uncompiled for ad +// until the ODE-state boundary is wired. +template class Node { public: typedef T strategy_type; typedef E environment_type; - typedef Individual individual_type; + typedef Individual individual_type; typedef typename strategy_type::ptr strategy_type_ptr; Node(strategy_type_ptr s); @@ -25,8 +33,8 @@ class Node { // Wrapper to growth_rate_gradient for testing double r_growth_rate_gradient(const environment_type& environment); - double height() const {return individual.state(HEIGHT_INDEX);} - double compute_competition(double z) const; + S height() const {return individual.state(HEIGHT_INDEX);} + S compute_competition(double z) const; double fecundity() const {return offspring_produced_survival_weighted;} // Bookkeeping recorded at the moment the node is introduced, so that @@ -91,7 +99,7 @@ class Node { individual.resize_consumption_rates(i); } - double consumption_rate(int i) const { + S consumption_rate(int i) const { return individual.consumption_rate(i) * density; } @@ -113,8 +121,8 @@ class Node { double patch_density_at_birth; }; -template -Node::Node(strategy_type_ptr s) +template +Node::Node(strategy_type_ptr s) : individual(s), log_density(-std::numeric_limits::infinity()), log_density_dt(0), @@ -125,8 +133,8 @@ Node::Node(strategy_type_ptr s) patch_density_at_birth(0) { } -template -void Node::compute_rates(const environment_type& environment, +template +void Node::compute_rates(const environment_type& environment, double pr_patch_survival) { individual.compute_rates(environment); @@ -157,8 +165,8 @@ void Node::compute_rates(const environment_type& environment, // // NOTE: The initial condition for log_density is also a bit tricky, and // defined on p 7 at the moment. -template -void Node::compute_initial_conditions(const environment_type& environment, +template +void Node::compute_initial_conditions(const environment_type& environment, double pr_patch_survival, double birth_rate) { pr_patch_survival_at_birth = pr_patch_survival; // Seed strategy-specific initial states (e.g. TF24f's tracked psi at its @@ -185,8 +193,16 @@ void Node::compute_initial_conditions(const environment_type& environment, // likely. } -template -double Node::growth_rate_gradient(const environment_type& environment) const { +template +double Node::growth_rate_gradient(const environment_type& environment) const { + // Exact AD growth-rate gradient when enabled and the strategy provides one + // (#537 A1). Returns NA otherwise, falling through to the finite difference. + if (individual.control().node_gradient_exact_ad) { + const double g = individual.growth_rate_gradient_exact(environment); + if (util::is_finite(g)) { + return g; + } + } // Finite-differencing the growth rate needs a mutable Individual to perturb // height on, but it must not disturb this node's already-computed state and // rates. Rather than copy-construct a fresh Individual (and its four @@ -217,8 +233,8 @@ double Node::growth_rate_gradient(const environment_type& environment) cons } // Wrapper to growth_rate_gradient for testing -template -double Node::r_growth_rate_gradient(const environment_type& environment) { +template +double Node::r_growth_rate_gradient(const environment_type& environment) { // We need to compute the physiological variables here, first, so // that reusing intervals works as expected. This would ordinarily // be taken care of because of the calling order of @@ -227,15 +243,15 @@ double Node::r_growth_rate_gradient(const environment_type& environment) { return growth_rate_gradient(environment); } -template -double Node::compute_competition(double height_) const { +template +S Node::compute_competition(double height_) const { return density * individual.compute_competition(height_); } // ODE interface -- note that the don't care about time in the node; // only Patch and above does. -template -odelia::ode::const_iterator Node::set_ode_state(odelia::ode::const_iterator it) { +template +odelia::ode::const_iterator Node::set_ode_state(odelia::ode::const_iterator it) { for (size_t i = 0; i < individual.ode_size(); i++) { individual.set_state(i, *it++); } @@ -243,8 +259,8 @@ odelia::ode::const_iterator Node::set_ode_state(odelia::ode::const_iterator set_log_density(*it++); return it; } -template -odelia::ode::iterator Node::ode_state(odelia::ode::iterator it) const { +template +odelia::ode::iterator Node::ode_state(odelia::ode::iterator it) const { for (size_t i = 0; i < individual.ode_size(); i++) { *it++ = individual.state(i); } @@ -252,8 +268,8 @@ odelia::ode::iterator Node::ode_state(odelia::ode::iterator it) const { *it++ = log_density; return it; } -template -odelia::ode::iterator Node::ode_rates(odelia::ode::iterator it) const { +template +odelia::ode::iterator Node::ode_rates(odelia::ode::iterator it) const { for (size_t i = 0; i < individual.ode_size(); i++) { *it++ = individual.rate(i); } @@ -262,8 +278,8 @@ odelia::ode::iterator Node::ode_rates(odelia::ode::iterator it) const { return it; } -template -odelia::ode::iterator Node::ode_aux(odelia::ode::iterator it) const { +template +odelia::ode::iterator Node::ode_aux(odelia::ode::iterator it) const { for (size_t i = 0; i < individual.aux_size(); i++) { *it++ = individual.aux(i); } diff --git a/inst/include/plant/patch.h b/inst/include/plant/patch.h index b15485eb..fbfb4ada 100644 --- a/inst/include/plant/patch.h +++ b/inst/include/plant/patch.h @@ -16,16 +16,25 @@ using namespace Rcpp; namespace plant { -template +// Templated on the scalar type S (#472 scope B / #537, Milestone C) to mirror +// Species / Node: S defaults to double so every existing +// Patch is Patch and bit-identical. The species storage is +// Species, so a Patch<...,ad> holds ad-typed individual state. NOTE the +// environment member stays type E (double): a Patch<...,ad> is MIXED -- ad +// species over a frozen-double resident environment, matching the two-pass +// replay design (freeze the resident light schedule, replay the node ODE with +// the active scalar). An ad-valued resident light spline (self-shading through +// traits) is a later piece -- the odelia AD interpolator (odelia PR #32). +template class Patch { public: using value_type = double; typedef T strategy_type; typedef E environment_type; - typedef Individual individual_type; - typedef Node node_type; - typedef Species species_type; + typedef Individual individual_type; + typedef Node node_type; + typedef Species species_type; typedef Parameters parameters_type; Patch(parameters_type p, environment_type e, plant::Control c); @@ -144,6 +153,70 @@ class Patch { std::vector> environment_history; std::vector environment_cache; + // Per-ODE-step resident STAND state for species 0, captured alongside + // environment_history during a save_RK45_cache run: each entry is the node + // heights / per-node competition effects (node.compute_competition(0) = + // k_I * area_leaf, so resident competition(z) = trapezium_i(ce_i * Q(z/h_i))). + // Exposed so the active-knot self-shading AD driver (#472 scope B) can + // reconstruct the resident light DIFFERENTIABLY per step (light responds to an + // allometric trait via area_leaf) instead of holding the cached env frozen. + // Single-species for now (the FF16 self-shading demo); multi-species is additive. + std::vector> stand_height_history; + std::vector> stand_competition_history; + + // Per-RK-STAGE resident stand state for species 0, captured alongside the 6 + // per-stage environments in environment_history (one inner vector per Cash-Karp + // stage, mirroring environment_cache exactly). This is the faithful-to-the-SCM + // harvest the RESIDENT total-gradient build (#472 scope B, R0) needs: at each + // RK stage the active light is reconstructed from the stand that PRODUCED that + // stage's frozen env (heights h_i, per-node competition effect ce_i), so the + // resident feedback derivative tracks the SCM stage-by-stage rather than holding + // a step-start census across stages (the cheaper per-step C-28 reconstruction). + // [step][stage 0..5][cohort]. Single-species for now (the FF16 light demo). + std::vector>> stand_height_stage_history; + std::vector>> stand_competition_stage_history; + // Per-step working caches (6 stages), filled by cache_RK45_step and flushed to + // the *_stage_history by cache_ode_step, exactly as environment_cache is. + std::vector> stand_height_stage_cache; + std::vector> stand_competition_stage_cache; + + // Boundary (new_node) per-RK-stage state for species 0: its height + competition + // effect. r_compute_competition_effect_by_nodes / r_heights iterate `nodes` only + // and OMIT the boundary new_node, whose tail term Species::compute_competition + // DOES add (it is the bottom of the last trapezium segment). Without it the + // coupled-replay env reconstruction (#472 scope B, R0) is wrong at ground level + // (z < the seedling top), where the boundary node's projected leaf area lands. + // [step][stage 0..5], a single scalar per stage (species 0). Empty for an empty + // species (a NaN-free placeholder is pushed so indices stay aligned). + std::vector> stand_newnode_height_stage_history; + std::vector> stand_newnode_competition_stage_history; + std::vector stand_newnode_height_stage_cache; + std::vector stand_newnode_competition_stage_cache; + // Boundary-node state SNAPSHOT taken at the moment compute_environment builds the + // light spline (species 0). compute_rates() then calls new_node.compute_initial_ + // conditions and MUTATES the boundary node, so querying it at cache time (after + // derivs = compute_environment + compute_rates) reads the wrong state. These hold + // the new_node baked into the cached env. (#472 scope B, R0 reconstruction.) + double newnode_height_env_snapshot = 0.0; + double newnode_competition_env_snapshot = 0.0; + + // ALL-SPECIES per-RK-stage harvest (#472 scope B, R2 -- the cross-species coupled + // Jacobian). The species-0 fields above drive the single-species coupled replay; + // these add the species dimension so the JOINT canopy can be reconstructed from + // every species' re-evolved cohorts. [step][stage 0..5][species][cohort] for the + // stand, [step][stage][species] for the boundary node. Additive: the single-species + // path is untouched. Filled only on save_RK45_cache runs. + std::vector>>> stand_height_stage_history_all; + std::vector>>> stand_competition_stage_history_all; + std::vector>> stand_newnode_height_stage_history_all; + std::vector>> stand_newnode_competition_stage_history_all; + std::vector>> stand_height_stage_cache_all; + std::vector>> stand_competition_stage_cache_all; + std::vector> stand_newnode_height_stage_cache_all; + std::vector> stand_newnode_competition_stage_cache_all; + std::vector newnode_height_env_snapshot_all; + std::vector newnode_competition_env_snapshot_all; + void cache_ode_step(); void cache_RK45_step(int step); void load_ode_step(); @@ -191,8 +264,8 @@ class Patch { std::vector> competition_error_by_node; }; -template -Patch::Patch(parameters_type p, environment_type e, Control c) +template +Patch::Patch(parameters_type p, environment_type e, Control c) : parameters(p), area(p.patch_area), environment(e), @@ -216,24 +289,24 @@ Patch::Patch(parameters_type p, environment_type e, Control c) reset(); } -template -void Patch::overwrite_strategies(std::vector strategies) { +template +void Patch::overwrite_strategies(std::vector strategies) { species.clear(); add_strategies(strategies); } -template -void Patch::add_strategies(std::vector strategies) { +template +void Patch::add_strategies(std::vector strategies) { for (auto i = 0; i < strategies.size(); ++i) { auto s = strategies[i]; s.control = control; // Overwrite to take the patch control object - auto spec = Species(s); + auto spec = Species(s); species.push_back(spec); } } -template -void Patch::set_mutant() { +template +void Patch::set_mutant() { if (environment_history.empty()) { util::stop("Run a resident first to generate a competitve landscape"); } @@ -244,8 +317,8 @@ void Patch::set_mutant() { idx = 0; } -template -void Patch::reset() { +template +void Patch::reset() { for (auto& s : species) { s.clear(); // allocate variables for tracking resource consumption @@ -284,8 +357,8 @@ void Patch::reset() { // double-arg set_ode_state) so the first environment build is a full // compute_environment(false): a rescale of the not-yet-built light spline would // read uninitialised grid state. -template -void Patch::set_initial_state() { +template +void Patch::set_initial_state() { const size_t n_species = species.size(); util::check_length(parameters.n_initial_cohorts.size(), n_species); @@ -332,8 +405,8 @@ void Patch::set_initial_state() { compute_rates(); } -template -void Patch::check_initial_density_rates() const { +template +void Patch::check_initial_density_rates() const { for (const auto& s : species) { std::vector rates = s.r_log_density_rates(); if (std::any_of(rates.begin(), rates.end(), @@ -345,8 +418,8 @@ void Patch::check_initial_density_rates() const { } } -template -double Patch::height_max() const { +template +double Patch::height_max() const { double ret = 0.0; for (size_t i = 0; i < species.size(); ++i) { if (!is_mutant_run) { @@ -356,8 +429,8 @@ double Patch::height_max() const { return ret; } -template -double Patch::compute_competition(double height) const { +template +double Patch::compute_competition(double height) const { double tot = 0.0; for (size_t i = 0; i < species.size(); ++i) { if (!is_mutant_run) { @@ -367,15 +440,15 @@ double Patch::compute_competition(double height) const { return tot; } -template -std::vector Patch::r_compute_competition_effect_error_by_node_for_species_i(size_t species_index) const { +template +std::vector Patch::r_compute_competition_effect_error_by_node_for_species_i(size_t species_index) const { const double tot_competition_effect = compute_competition(0.0); return species[species_index].r_compute_competition_effect_by_nodes_error(tot_competition_effect); } // Integrate over lifetime fitness of individual nodes, scaled per node. -template -double Patch::net_reproduction_ratio_for_species( +template +double Patch::net_reproduction_ratio_for_species( size_t species_index, std::vector const& scalars) const { auto net_prod = species[species_index].net_reproduction_ratio_by_node_weighted(); auto const times = species[species_index].node_times(); @@ -387,8 +460,8 @@ double Patch::net_reproduction_ratio_for_species( } // Offspring production, equal to overall fitness scaled by the birth rate. -template -std::vector Patch::offspring_production() const { +template +std::vector Patch::offspring_production() const { auto ret = std::vector(species.size()); for (size_t i = 0; i < species.size(); ++i) { // scale by birth rate function over time @@ -403,8 +476,8 @@ std::vector Patch::offspring_production() const { } // Overall fitness (no scaling, ie scalars set to 1.0). -template -std::vector Patch::net_reproduction_ratios() const { +template +std::vector Patch::net_reproduction_ratios() const { auto ret = std::vector(species.size()); for (size_t i = 0; i < species.size(); ++i) { auto scalars = std::vector(species[i].size(), 1.0); @@ -414,8 +487,8 @@ std::vector Patch::net_reproduction_ratios() const { } // Sum up all offspring produced. -template -double Patch::total_offspring_production() const { +template +double Patch::total_offspring_production() const { double total = 0.0; std::vector offspring = offspring_production(); for (size_t i = 0; i < species.size(); ++i) { @@ -425,8 +498,8 @@ double Patch::total_offspring_production() const { } // Check integration errors for each species' reproduction integral. -template -std::vector> Patch::net_reproduction_ratio_errors() const { +template +std::vector> Patch::net_reproduction_ratio_errors() const { std::vector> ret; double total_offspring = total_offspring_production(); for (size_t i = 0; i < species.size(); ++i) { @@ -440,8 +513,8 @@ std::vector> Patch::net_reproduction_ratio_errors() con // Sample the competition error for each species introduced this step and fold // it into the running per-node max (ignoring NA, matching na.rm=TRUE in R). -template -void Patch::collect_competition_errors(const std::vector& added) { +template +void Patch::collect_competition_errors(const std::vector& added) { for (size_t idx : added) { std::vector v = r_compute_competition_effect_error_by_node_for_species_i(idx); @@ -460,8 +533,8 @@ void Patch::collect_competition_errors(const std::vector& added) { // Combine the competition error (sampled during the run) with the reproduction // error (computed now) into a single per-node error vector per species. An // all-NA node yields -Inf, matching apply(rbind(...), 2, max, na.rm=TRUE) in R. -template -std::vector> Patch::refinement_error_by_node() const { +template +std::vector> Patch::refinement_error_by_node() const { std::vector> repro = net_reproduction_ratio_errors(); std::vector> ret(species.size()); for (size_t i = 0; i < species.size(); ++i) { @@ -484,20 +557,35 @@ std::vector> Patch::refinement_error_by_node() const { // Pre-compute environment, as shaped by residents // Creates splines of resource availability -template -void Patch::compute_environment(bool rescale) { +template +void Patch::compute_environment(bool rescale) { // Define an anonymous function to use in creation of environment auto f = [&](double x) -> double { return compute_competition(x); }; if (size() > 0 & !is_mutant_run) { environment.compute_environment(f, height_max(), rescale); + // Snapshot the boundary node baked into THIS spline, before compute_rates() + // mutates it (species 0; the per-RK-stage harvest reads these in cache). + if (!species.empty()) { + const auto& nn = species[0].r_new_node(); + newnode_height_env_snapshot = nn.height(); + newnode_competition_env_snapshot = nn.compute_competition(0.0); + } + // Per-species boundary-node snapshot for the all-species coupled harvest (R2). + newnode_height_env_snapshot_all.assign(species.size(), 0.0); + newnode_competition_env_snapshot_all.assign(species.size(), 0.0); + for (size_t k = 0; k < species.size(); ++k) { + const auto& nnk = species[k].r_new_node(); + newnode_height_env_snapshot_all[k] = nnk.height(); + newnode_competition_env_snapshot_all[k] = nnk.compute_competition(0.0); + } } } -template -void Patch::compute_rates() { +template +void Patch::compute_rates() { // Computes rates of change for the patch, including all the component species // While the patch has an `environment`, the rates here are calculated from @@ -534,16 +622,16 @@ void Patch::compute_rates() { // TODO(#478): We should only be recomputing the light environment for the // points that are below the height of the seedling -- not the entire // light environment; probably worth just doing a rescale there? -template -void Patch::introduce_new_node(size_t species_index) { +template +void Patch::introduce_new_node(size_t species_index) { species[species_index].introduce_new_node(); compute_environment(false); } -template -void Patch::introduce_new_nodes(const std::vector& species_index) { +template +void Patch::introduce_new_nodes(const std::vector& species_index) { // Record introduction time and patch-age density on each node as it is // introduced, so lifetime-fitness calcs need not look these up later. const double t = time(); @@ -555,8 +643,8 @@ void Patch::introduce_new_nodes(const std::vector& species_index) { compute_environment(false); } -template -void Patch::r_set_time(double time) { +template +void Patch::r_set_time(double time) { environment.time = time; } @@ -564,8 +652,8 @@ void Patch::r_set_time(double time) { // time: time // state: vector of ode state; we'll pass an iterator with that in // n: number of *individuals* of each species -template -void Patch::r_set_state(double time, +template +void Patch::r_set_state(double time, const std::vector& state, const std::vector& n, const std::vector& light_availability) { @@ -583,26 +671,26 @@ void Patch::r_set_state(double time, } // ODE interface -template -size_t Patch::ode_size() const { +template +size_t Patch::ode_size() const { return odelia::ode::ode_size(species.begin(), species.end()) + environment.ode_size(); } -template -size_t Patch::aux_size() const { +template +size_t Patch::aux_size() const { // TODO(#478): Is this useful for environment vectors? // no use for auxiliary environment variables (yet) return odelia::ode::aux_size(species.begin(), species.end());// + environment.ode_size(); } -template -double Patch::ode_time() const { +template +double Patch::ode_time() const { return time(); } // First set_ode_state function is for resident runs. Second is for mutant runs -template -odelia::ode::const_iterator Patch::set_ode_state(odelia::ode::const_iterator it, +template +odelia::ode::const_iterator Patch::set_ode_state(odelia::ode::const_iterator it, double time) { // Set ode states @@ -624,8 +712,8 @@ odelia::ode::const_iterator Patch::set_ode_state(odelia::ode::const_iterato // used for mutant runs // -- differs from above in that an index is passed in as argument // -- environments are loaded from ODE history, instead of being calculated -template -odelia::ode::const_iterator Patch::set_ode_state(odelia::ode::const_iterator it, +template +odelia::ode::const_iterator Patch::set_ode_state(odelia::ode::const_iterator it, int index) { it = odelia::ode::set_ode_state(species.begin(), species.end(), it); @@ -644,29 +732,87 @@ odelia::ode::const_iterator Patch::set_ode_state(odelia::ode::const_iterato // called from ode_solver->cache // saves cached set of environments(6) from each ODE step to the step history -template -void Patch::cache_ode_step() { - if(save_RK45_cache) { +template +void Patch::cache_ode_step() { + if(save_RK45_cache) { step_history.push_back(time()); environment_history.push_back(environment_cache); + // Capture the species-0 stand (heights + per-node competition effect) at this + // step boundary, for the active-knot resident-light reconstruction. + if (!species.empty()) { + stand_height_history.push_back(species[0].r_heights()); + stand_competition_history.push_back(species[0].r_compute_competition_effect_by_nodes()); + } + // Flush the per-RK-stage stand caches accumulated over this step's 6 stages. + stand_height_stage_history.push_back(stand_height_stage_cache); + stand_competition_stage_history.push_back(stand_competition_stage_cache); + stand_newnode_height_stage_history.push_back(stand_newnode_height_stage_cache); + stand_newnode_competition_stage_history.push_back(stand_newnode_competition_stage_cache); + // All-species harvest (R2). + stand_height_stage_history_all.push_back(stand_height_stage_cache_all); + stand_competition_stage_history_all.push_back(stand_competition_stage_cache_all); + stand_newnode_height_stage_history_all.push_back(stand_newnode_height_stage_cache_all); + stand_newnode_competition_stage_history_all.push_back(stand_newnode_competition_stage_cache_all); } } // called from ode_step->cache // saves environment at each RK45 step to the environment cache -template -void Patch::cache_RK45_step(int step) { - if(save_RK45_cache) { +template +void Patch::cache_RK45_step(int step) { + if(save_RK45_cache) { if(step == 0) { environment_cache.clear(); + stand_height_stage_cache.clear(); + stand_competition_stage_cache.clear(); + stand_newnode_height_stage_cache.clear(); + stand_newnode_competition_stage_cache.clear(); + stand_height_stage_cache_all.clear(); + stand_competition_stage_cache_all.clear(); + stand_newnode_height_stage_cache_all.clear(); + stand_newnode_competition_stage_cache_all.clear(); } environment_cache.push_back(environment); + // All-species stand for THIS stage: [species][cohort] heights + per-node effect, + // and the per-species boundary node. Mirrors the species-0 capture below. + { + std::vector> h_all(species.size()), c_all(species.size()); + std::vector nnh_all(species.size(), 0.0), nnc_all(species.size(), 0.0); + for (size_t k = 0; k < species.size(); ++k) { + h_all[k] = species[k].r_heights(); + c_all[k] = species[k].r_compute_competition_effect_by_nodes(); + if (k < newnode_height_env_snapshot_all.size()) { + nnh_all[k] = newnode_height_env_snapshot_all[k]; + nnc_all[k] = newnode_competition_env_snapshot_all[k]; + } + } + stand_height_stage_cache_all.push_back(h_all); + stand_competition_stage_cache_all.push_back(c_all); + stand_newnode_height_stage_cache_all.push_back(nnh_all); + stand_newnode_competition_stage_cache_all.push_back(nnc_all); + } + // Capture the species-0 stand that produced THIS stage's environment (the ODE + // state was just set to the stage trial point and compute_environment ran in + // derivs, immediately before this cache call), aligned 1:1 with environment_cache. + if (!species.empty()) { + stand_height_stage_cache.push_back(species[0].r_heights()); + stand_competition_stage_cache.push_back(species[0].r_compute_competition_effect_by_nodes()); + // Boundary node tail term (compute_competition adds it beyond the `nodes` loop), + // snapshotted in compute_environment before compute_rates mutated it. + stand_newnode_height_stage_cache.push_back(newnode_height_env_snapshot); + stand_newnode_competition_stage_cache.push_back(newnode_competition_env_snapshot); + } else { + stand_height_stage_cache.emplace_back(); + stand_competition_stage_cache.emplace_back(); + stand_newnode_height_stage_cache.push_back(0.0); + stand_newnode_competition_stage_cache.push_back(0.0); + } } } // called from ode_solver->load, only gets called for mutant runs -template -void Patch::load_ode_step() { +template +void Patch::load_ode_step() { if (use_cached_environment) { // Minor optimization to check the current and next index before doing a search, as the most common case is that the ODE solver is stepping through the cached environments in order. If the call sequence was not strictly sequential, we fallback to a search through the step history to find the correct environment. @@ -694,15 +840,15 @@ void Patch::load_ode_step() { } } -template -odelia::ode::iterator Patch::ode_state(odelia::ode::iterator it) const { +template +odelia::ode::iterator Patch::ode_state(odelia::ode::iterator it) const { it = odelia::ode::ode_state(species.begin(), species.end(), it); it = environment.ode_state(it); return it; } -template -Rcpp::List Patch::r_get_state() const +template +Rcpp::List Patch::r_get_state() const { // Aseemble commkunity state, icnluding auxiallry variables @@ -718,15 +864,15 @@ Rcpp::List Patch::r_get_state() const _["env"] = environment.r_get_state()); } -template -odelia::ode::iterator Patch::ode_rates(odelia::ode::iterator it) const { +template +odelia::ode::iterator Patch::ode_rates(odelia::ode::iterator it) const { it = odelia::ode::ode_rates(species.begin(), species.end(), it); it = environment.ode_rates(it); return it; } -template -odelia::ode::iterator Patch::ode_aux(odelia::ode::iterator it) const { +template +odelia::ode::iterator Patch::ode_aux(odelia::ode::iterator it) const { it = odelia::ode::ode_aux(species.begin(), species.end(), it); return it; } diff --git a/inst/include/plant/qk.h b/inst/include/plant/qk.h index 950f8fc0..1c15f6de 100644 --- a/inst/include/plant/qk.h +++ b/inst/include/plant/qk.h @@ -26,6 +26,30 @@ class QK { template double integrate(Function f, double a, double b); + // Scalar-templated Kronrod estimate for AD (#472 scope B / #537): the same + // fixed Gauss-Kronrod rule (xgk/wgk constants), but the bounds, abscissae and + // accumulator are the integrand's scalar type S, so an AD active bound (e.g. a + // plant height) propagates through the MOVING nodes -- which a frozen-node + // replay would miss. Returns only the Kronrod result (no error/abs/asc + // estimate, which the gradient does not need). Stateless (no last_* writes). + template + S integrate_ad(Function f, S a, S b) const { + const S center = 0.5 * (a + b); + const S half_length = 0.5 * (b - a); + S result_kronrod = f(center) * wgk[n - 1]; + for (size_t j = 0; j < (n - 1) / 2; j++) { + const size_t jtw = j * 2 + 1; + const S abscissa = half_length * xgk[jtw]; + result_kronrod += wgk[jtw] * (f(center - abscissa) + f(center + abscissa)); + } + for (size_t j = 0; j < n / 2; j++) { + const size_t jtwm1 = j * 2; + const S abscissa = half_length * xgk[jtwm1]; + result_kronrod += wgk[jtwm1] * (f(center - abscissa) + f(center + abscissa)); + } + return result_kronrod * half_length; + } + // These two provide very low level access to the integration // routines. std::vector integrate_vector_x(double a, double b) const; diff --git a/inst/include/plant/resource_spline.h b/inst/include/plant/resource_spline.h index 1e404512..7d0aaad5 100644 --- a/inst/include/plant/resource_spline.h +++ b/inst/include/plant/resource_spline.h @@ -88,6 +88,18 @@ class ResourceSpline { return height <= cap ? std::max(0.0, spline(height)) : 1.0; } + // Analytic d(value)/d(height), consistent with get_value_at_height: 0 above + // the cap (constant open value) and where the #253 floor clamps a negative + // undershoot, else the spline's analytic derivative. Enables exact AD + // gradients through the resource environment (#472 scope B / #537). + double get_value_deriv_at_height(double height) const { + return get_value_deriv_at_height(height, spline.max()); + } + double get_value_deriv_at_height(double height, double cap) const { + if (height > cap) return 0.0; + return spline(height) > 0.0 ? spline.deriv(height) : 0.0; + } + virtual void r_init_interpolators(const std::vector& state) { // See issue #144; this is important as we have to at least refine // the light environment, but doing this is better because it means diff --git a/inst/include/plant/species.h b/inst/include/plant/species.h index ff83d6fb..e38fe023 100644 --- a/inst/include/plant/species.h +++ b/inst/include/plant/species.h @@ -18,14 +18,22 @@ namespace plant { // survival-weighted rates, lifetime fitness, schedule-refinement error) stays // here. -template -class Species : public SpeciesBase, T, E, Node> { - typedef SpeciesBase, T, E, Node> base_type; +// Templated on the scalar type S (#472 scope B / #537, Milestone C) to mirror +// Node / Individual: S defaults to double so every existing +// Species is Species and bit-identical. The node storage is +// Node (the element passed to SpeciesBase), so a Species<...,ad> holds +// ad-typed individual state. The R-facing accessors (r_heights, r_get_state, +// the std::vector reductions) stay double-returning and, like the rest +// of the hierarchy, compile per-member-on-use -- for the ad instantiation they +// stay uncompiled, only the trait-carrying compute_rates/competition path is. +template +class Species : public SpeciesBase, T, E, Node> { + typedef SpeciesBase, T, E, Node> base_type; public: typedef T strategy_type; typedef E environment_type; - typedef Individual individual_type; - typedef Node node_type; + typedef Individual individual_type; + typedef Node node_type; typedef typename strategy_type::ptr strategy_type_ptr; Species(strategy_type s); @@ -129,26 +137,26 @@ class Species : public SpeciesBase, T, E, Node> { typedef typename std::vector::const_iterator nodes_const_iterator; }; -template -Species::Species(strategy_type s) +template +Species::Species(strategy_type s) : base_type(s), new_node(this->strategy) { } -template -size_t Species::size() const { +template +size_t Species::size() const { return nodes.size(); } -template -void Species::clear() { +template +void Species::clear() { nodes.clear(); // Reset the new_node to a blank new_node, too. new_node = node_type(strategy); } -template -void Species::introduce_new_node() { +template +void Species::introduce_new_node() { // new_node already holds the initial conditions computed against the current // environment by the most recent compute_rates() call (see compute_rates -> // new_node.compute_initial_conditions above), and the member is refreshed @@ -163,8 +171,8 @@ void Species::introduce_new_node() { // seed of the species. Otherwise we return the height of the largest // individual (always the first in the list) which will be at least // tall as a seed. -template -double Species::height_max() const { +template +double Species::height_max() const { return nodes.empty() ? new_node.height() : nodes.front().height(); } @@ -193,8 +201,8 @@ double Species::height_max() const { // single node (needed to be the second half of the trapezium) and // also needed if the last looked at plant was still contributing to // the integral). -template -double Species::compute_competition(double height) const { +template +double Species::compute_competition(double height) const { if (size() == 0 || height_max() < height) { return 0.0; } @@ -228,24 +236,24 @@ double Species::compute_competition(double height) const { // NOTE: We should probably prefer to rescale when this is called // through the ode stepper. -template -void Species::compute_rates(const E& environment, double pr_patch_survival, double birth_rate) { +template +void Species::compute_rates(const E& environment, double pr_patch_survival, double birth_rate) { for (auto& c : nodes) { c.compute_rates(environment, pr_patch_survival); } new_node.compute_initial_conditions(environment, pr_patch_survival, birth_rate); } -template -void Species::introduce_new_node(double time, double patch_density) { +template +void Species::introduce_new_node(double time, double patch_density) { // Stamp the pushed copy (not new_node) so the member stays pristine for // the no-arg introduction paths. nodes.push_back(new_node); nodes.back().set_introduction(time, patch_density); } -template -std::vector Species::net_reproduction_ratio_by_node() const { +template +std::vector Species::net_reproduction_ratio_by_node() const { std::vector ret; ret.reserve(size()); for (auto& c : nodes) { @@ -254,8 +262,8 @@ std::vector Species::net_reproduction_ratio_by_node() const { return ret; } -template -std::vector Species::net_reproduction_ratio_by_node_weighted() const { +template +std::vector Species::net_reproduction_ratio_by_node_weighted() const { std::vector ret; ret.reserve(size()); for (auto& c : nodes) { @@ -264,8 +272,8 @@ std::vector Species::net_reproduction_ratio_by_node_weighted() cons return ret; } -template -std::vector Species::node_times() const { +template +std::vector Species::node_times() const { std::vector ret; ret.reserve(size()); for (auto& c : nodes) { @@ -274,13 +282,13 @@ std::vector Species::node_times() const { return ret; } -template -void Species::resize_consumption_rates(int r) { +template +void Species::resize_consumption_rates(int r) { new_node.resize_consumption_rates(r); } -template -double Species::consumption_rate(int i) const { +template +double Species::consumption_rate(int i) const { // can't determine density for one node if(size() < 2) { return 0.0; @@ -290,8 +298,8 @@ double Species::consumption_rate(int i) const { } } -template -std::vector Species::consumption_rate_by_node_rev(int i) const { +template +std::vector Species::consumption_rate_by_node_rev(int i) const { std::vector ret; ret.reserve(size()); for(auto it = nodes.rbegin(); it != nodes.rend(); ++it) { @@ -301,18 +309,18 @@ std::vector Species::consumption_rate_by_node_rev(int i) const { } // bit clunky... -template -size_t Species::aux_size() const { +template +size_t Species::aux_size() const { return size() * strategy->aux_size(); } -template -odelia::ode::iterator Species::ode_aux(odelia::ode::iterator it) const { +template +odelia::ode::iterator Species::ode_aux(odelia::ode::iterator it) const { return odelia::ode::ode_aux(nodes.begin(), nodes.end(), it); } -template -Rcpp::NumericMatrix Species::r_get_state() const { +template +Rcpp::NumericMatrix Species::r_get_state() const { size_t ode_size = node_type::ode_size(), n_nodes = size(); size_t aux_size = strategy->aux_size(); @@ -340,8 +348,8 @@ Rcpp::NumericMatrix Species::r_get_state() const { return ret; } -template -std::vector Species::r_heights() const { +template +std::vector Species::r_heights() const { std::vector ret; ret.reserve(size()); for (nodes_const_iterator it = nodes.begin(); @@ -351,8 +359,8 @@ std::vector Species::r_heights() const { return ret; } -template -std::vector Species::r_heights_rev() const { +template +std::vector Species::r_heights_rev() const { std::vector ret; ret.reserve(size()); for (nodes_const_iterator it = nodes.begin(); @@ -363,8 +371,8 @@ std::vector Species::r_heights_rev() const { return ret; } -template -void Species::r_set_heights(std::vector heights) { +template +void Species::r_set_heights(std::vector heights) { util::check_length(heights.size(), size()); if (!util::is_decreasing(heights.begin(), heights.end())) { util::stop("height must be decreasing (ties allowed)"); @@ -375,8 +383,8 @@ void Species::r_set_heights(std::vector heights) { } } -template -std::vector Species::r_compute_competition_effect_by_nodes() const { +template +std::vector Species::r_compute_competition_effect_by_nodes() const { std::vector ret; ret.reserve(size()); for (auto& c : nodes) { @@ -385,13 +393,13 @@ std::vector Species::r_compute_competition_effect_by_nodes() const return ret; } -template -std::vector Species::r_compute_competition_effect_by_nodes_error(double scal) const { +template +std::vector Species::r_compute_competition_effect_by_nodes_error(double scal) const { return util::local_error_integration(r_heights(), r_compute_competition_effect_by_nodes(), scal); } -template -std::vector Species::r_log_densities() const { +template +std::vector Species::r_log_densities() const { std::vector ret; ret.reserve(size()); for (nodes_const_iterator it = nodes.begin(); @@ -401,8 +409,8 @@ std::vector Species::r_log_densities() const { return ret; } -template -std::vector Species::r_log_density_rates() const { +template +std::vector Species::r_log_density_rates() const { std::vector ret; ret.reserve(size()); for (nodes_const_iterator it = nodes.begin(); it != nodes.end(); ++it) { @@ -411,8 +419,8 @@ std::vector Species::r_log_density_rates() const { return ret; } -template -std::vector Species::r_patch_densities() const { +template +std::vector Species::r_patch_densities() const { std::vector ret; ret.reserve(size()); for (nodes_const_iterator it = nodes.begin(); it != nodes.end(); ++it) { @@ -421,8 +429,8 @@ std::vector Species::r_patch_densities() const { return ret; } -template -std::vector Species::r_pr_patch_survival_at_birth() const { +template +std::vector Species::r_pr_patch_survival_at_birth() const { std::vector ret; ret.reserve(size()); for (nodes_const_iterator it = nodes.begin(); it != nodes.end(); ++it) { @@ -431,8 +439,8 @@ std::vector Species::r_pr_patch_survival_at_birth() const { return ret; } -template -void Species::set_birth_state(const std::vector& times, +template +void Species::set_birth_state(const std::vector& times, const std::vector& patch_density, const std::vector& pr_patch_survival) { util::check_length(times.size(), size()); diff --git a/inst/include/plant/strategy.h b/inst/include/plant/strategy.h index 6fb74c44..bd1de00a 100644 --- a/inst/include/plant/strategy.h +++ b/inst/include/plant/strategy.h @@ -52,6 +52,15 @@ class Strategy { void compute_rates(const environment_type& environment, Internals& vars); + // Exact d(growth rate)/d(height) via AD, if this strategy provides one + // (#537 A1). The default signals "unavailable" (NA), so + // Node::growth_rate_gradient falls back to its finite difference; FF16 + // overrides this. Non-virtual: Node calls it on the concrete strategy type. + double growth_rate_gradient_height_ad(double /*height*/, + const environment_type& /*environment*/) { + return NA_REAL; + } + void update_dependent_aux(const int index, Internals& vars); // Seed strategy-specific initial ODE states for a newly introduced individual, diff --git a/man/birth_rate_gradient.Rd b/man/birth_rate_gradient.Rd new file mode 100644 index 00000000..594cdb76 --- /dev/null +++ b/man/birth_rate_gradient.Rd @@ -0,0 +1,71 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/emergent_gradient.R +\name{birth_rate_gradient} +\alias{birth_rate_gradient} +\title{Resident-coupled birth-rate gradient of emergent census metrics (FF16)} +\usage{ +birth_rate_gradient( + scm, + metrics = c("LAI", "biomass", "size_moment"), + species = 1L, + birth_rate = NULL +) +} +\arguments{ +\item{scm}{An \code{SCM} run with \code{save_RK45_cache = TRUE} (FF16 strategy).} + +\item{metrics}{Character vector, any of the census metrics \code{"LAI"}, +\code{"biomass"}, \code{"size_moment"} and -- single-species only -- +\code{"net_reproduction_ratio"} (the demographic-equilibrium / \eqn{R_0 = 1} axis: +\eqn{\partial R_0/\partial b}, with unit per-seed weights so \code{birth_rate} enters +only through the canopy -- the pure density feedback, i.e. the change in a +trait-identical mutant's fitness as the resident density moves; negative, which is +what makes an \eqn{R_0 = 1} equilibrium well-posed). (\code{offspring_production} is +excluded -- its birth-rate derivative is the trivial frozen identity above.)} + +\item{species}{Integer species index; on a multi-species stand this selects which +species' birth-rate driver to differentiate (the cross-species total). Default \code{1}.} + +\item{birth_rate}{The birth-rate driver to differentiate at; recovered from the run +by default (single-species only; the multi-species path recovers per species).} +} +\value{ +A list with \code{d_birth_rate} (named \code{d(metric)/d(birth_rate)}) and the + reconstructed \code{values} (plus \code{env_err}, the joint-canopy drift gauge, on a + multi-species stand). +} +\description{ +Reverse-mode derivative of an SCM's emergent census metrics w.r.t. the per-species +\emph{birth-rate driver} (#472 scope B, the birth-rate gradient AXIS). +} +\details{ +This differentiates a different axis than \code{\link{stand_gradient}}: not a trait, +but the constant birth-rate driver of a species. The use case is evolving a community +toward \strong{demographic equilibrium} -- birth_rate scales every cohort's +establishment density, so it re-shades the whole resident canopy, and its derivative +gives the Newton/gradient step that drives an emergent census target (e.g. a target +LAI / biomass, or self-replacement). + +Only the \strong{resident-coupled} reading needs a tape: with the canopy held frozen +every cohort's density is exactly linear in birth_rate, so the (rare-mutant) frozen +derivative is the trivial identity \eqn{\partial(\mathrm{metric})/\partial b = +\mathrm{metric}/b} (and \code{offspring_production} keeps that frozen reading even +under resident feedback -- the canopy a rare mutant invades is the resident's). So +this returns the resident TOTAL \eqn{\partial(\mathrm{census\ metric})/\partial b}: +birth_rate is registered as the sole input on the same coupled whole-stand replay the +resident trait gradient uses (one reverse sweep per metric). + +On a \strong{multi-species} stand this returns the CROSS-SPECIES total +\eqn{\partial(\mathrm{total\text{-}stand\ metric})/\partial b_s} for the chosen +\code{species} \eqn{s}: species \eqn{s}'s recruitment density re-shades the joint +canopy every species reads (the cross term the frozen reading sets to zero). It is +well-conditioned on a fixed node schedule and gated (with a clear error) for the stiff +refined-schedule case, like the cross-species trait gradient. + +FF16 only so far (its closed-form net keeps the coupled replay robust to long +horizons); TF24f's coupled replay is stiffness-gated (see +\code{stand_gradient(feedback = "resident")}). +} +\seealso{ +\code{\link{stand_gradient}}, \code{\link{offspring_production_gradient}}. +} diff --git a/man/grow_individual_to_size_gradient.Rd b/man/grow_individual_to_size_gradient.Rd new file mode 100644 index 00000000..8f848f81 --- /dev/null +++ b/man/grow_individual_to_size_gradient.Rd @@ -0,0 +1,79 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/emergent_gradient.R +\name{grow_individual_to_size_gradient} +\alias{grow_individual_to_size_gradient} +\title{Reverse-mode gradient of grow_individual_to_size (FF16 / TF24f)} +\usage{ +grow_individual_to_size_gradient( + individual, + sizes, + size_name, + env, + traits = NULL, + time_max = Inf, + warn = TRUE +) +} +\arguments{ +\item{individual}{An \code{Individual} object (FF16 or TF24f strategy), as passed to +\code{\link{grow_individual_to_size}}. For TF24f the trajectory carries the tracked +collar (\code{opt_root_psi_state}) as a 6th state and the gradient is dispatched to +the TF24f AD tape (same frozen-schedule replay + stopping-time IFT, with the collar +curvature-linearised); all other behaviour is identical.} + +\item{sizes}{A vector of target sizes to grow the plant to (increasing).} + +\item{size_name}{The size variable the targets refer to (one of the ODE state +names, e.g. \code{"height"}; FF16's monotonic size).} + +\item{env}{An \code{Environment} object (the fixed environment).} + +\item{traits}{Character vector of FF16 trait names; \code{NULL} (default) uses all +28 production-relevant parameters.} + +\item{time_max, warn}{Passed through to \code{\link{grow_individual_to_size}} for the +schedule-discovery pass.} +} +\value{ +A list with \code{$time} (the reconstructed \eqn{t^*} per size), + \code{$state} (a sizes x component matrix of the ODE state at \eqn{t^*}), + \code{$d_time} (a sizes x trait matrix \eqn{dt^*/d\theta}) and \code{$d_state} (a + sizes x component x trait array of TOTAL \eqn{d\,\mathrm{state}/d\theta}). The + reconstructed \code{$time}/\code{$state} match \code{grow_individual_to_size} to + the live \code{uniroot} tolerance. +} +\description{ +Reverse-mode trait gradient of \code{\link{grow_individual_to_size}} (#472 scope +B, FF16): a single plant grown in a FIXED environment to target size(s), +differentiated w.r.t. traits. +} +\details{ +The per-plant, fixed-environment counterpart of \code{\link{stand_gradient}}: there +is no resident feedback (the environment is given), so the gradient is the exact +derivative of the grow-to-size solve. For each target size it returns the derivative +of the stopping TIME \eqn{t^*} (the time the plant reaches the target) and of every +ODE STATE component at \eqn{t^*}, w.r.t. all FF16 traits, as an optional output of +the same call you already make. This is the gradient a growth-rate or +time-to-size calibration / optimisation consumes (e.g. +\code{\link{optimise_individual_rate_at_size_by_trait}}). + +It is a two-pass method mirroring the SCM gradients. Pass 1 runs the ordinary +\code{\link{grow_individual_to_size}} (its adaptive Cash-Karp step schedule and +per-node trajectory are harvested). Pass 2 replays the demographic ODE over that +FROZEN schedule with the trait active, reading the fixed environment with the +default deep-crown assimilation, to a single partial final step landing on +\eqn{t^*}; one reverse sweep per state component gives \eqn{\partial +\mathrm{state}/\partial\theta} at fixed \eqn{t^*}. The stopping time itself responds +to the trait through the implicit function theorem on the stopping condition +\eqn{\mathrm{size}(t^*,\theta) = \mathrm{target}}: +\deqn{dt^*/d\theta = -\,(\partial\,\mathrm{size}/\partial\theta\,|_{t^*}) / + \dot{\mathrm{size}}(t^*),} +so the TOTAL derivative of each returned component \eqn{y_c} is \eqn{dy_c/d\theta = +\partial y_c/\partial\theta|_{t^*} + \dot y_c(t^*)\,dt^*/d\theta} (for the size +component itself the two terms cancel, as it is pinned to the target). The seedling +size \eqn{h_0} (which solves \eqn{\mathrm{mass}(h_0) = \mathrm{seed\ mass}}) carries +its own \eqn{dh_0/d\theta} by the same implicit-function step. +} +\seealso{ +\code{\link{grow_individual_to_size}}, \code{\link{stand_gradient}}. +} diff --git a/man/offspring_production_gradient.Rd b/man/offspring_production_gradient.Rd new file mode 100644 index 00000000..68badf7f --- /dev/null +++ b/man/offspring_production_gradient.Rd @@ -0,0 +1,50 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/emergent_gradient.R +\name{offspring_production_gradient} +\alias{offspring_production_gradient} +\title{Reverse-mode gradient of emergent offspring_production} +\usage{ +offspring_production_gradient( + scm, + traits = NULL, + species = 1L, + birth_rate = NULL +) +} +\arguments{ +\item{scm}{An \code{SCM} run with \code{control(save_RK45_cache = TRUE)} (FF16, +TF24, or TF24f). The cached schedule + resident light are read from its patch +natively; the SCM is not re-run.} + +\item{traits}{Character vector of trait (parameter) names to differentiate. +\code{NULL} (default) uses the strategy's full production-trait set.} + +\item{species}{Integer index of the species (cohort family) to differentiate. +Default \code{1}.} + +\item{birth_rate}{The (constant) birth-rate driver used in the run. Recovered as +\code{offspring_production / net_reproduction_ratio} by default; pass it for a +time-varying driver.} +} +\value{ +A named numeric vector of trait derivatives, with attribute + \code{"offspring_production"} (the value reconstructed by the replay, which + should match \code{scm$offspring_production[[species]]}). +} +\description{ +Reverse-mode trait gradient of an SCM's emergent \code{offspring_production} +(#472 scope B). +} +\details{ +A convenience wrapper over \code{\link{stand_gradient}} for the single +\code{offspring_production} metric, returning a named vector instead of a +one-row Jacobian. It works for \strong{any} strategy that \code{stand_gradient} +supports (FF16, TF24, TF24f) -- the dispatch and the reverse-mode replay live in +\code{stand_gradient}; this just selects the metric and reshapes. It is the +rare-mutant / invasion-fitness (selection) gradient: the resident light is held +frozen, so for a multi-species stand it is species \code{s}'s gradient against +the shared frozen canopy of all species. +} +\seealso{ +\code{\link{stand_gradient}}, the canonical metrics x traits engine. +} diff --git a/man/stand_gradient.Rd b/man/stand_gradient.Rd new file mode 100644 index 00000000..0e3f01df --- /dev/null +++ b/man/stand_gradient.Rd @@ -0,0 +1,99 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/emergent_gradient.R +\name{stand_gradient} +\alias{stand_gradient} +\title{Reverse-mode Jacobian of emergent stand metrics (FF16)} +\usage{ +stand_gradient( + scm, + metrics = "offspring_production", + traits = NULL, + species = 1L, + birth_rate = NULL, + feedback = c("frozen", "resident") +) +} +\arguments{ +\item{scm}{An \code{SCM} run with \code{save_RK45_cache = TRUE} (FF16 strategy).} + +\item{metrics}{Character vector of stand-metric names, any of +\code{"offspring_production"} (the seed-rain integral), \code{"LAI"} (leaf-area +index = the SCM's \code{compute_competition(0)}), \code{"biomass"} (the +size-distribution integral of live + heartwood mass) and \code{"size_moment"} +(the first moment of the size distribution, \eqn{\int n(h)\,h\,dh}). All are +symmetric reductions over the replayed cohorts; none is privileged.} + +\item{traits}{Character vector of FF16 trait names. \code{NULL} (default) uses +all 28 production-relevant parameters.} + +\item{species}{Integer index of the species (cohort family); see +\code{\link{offspring_production_gradient}}. Default \code{1}.} + +\item{birth_rate}{The (constant) birth-rate driver; recovered from the run by +default.} + +\item{feedback}{How the resident light responds to the trait. \code{"frozen"} +(default) holds the canopy fixed -- the rare-mutant / invasion-fitness gradient, +correct for \code{offspring_production}. \code{"resident"} is the resident TOTAL +gradient via the COUPLED whole-stand replay: all cohorts are re-evolved together +over the frozen schedule and the canopy light is reconstructed each RK stage from +the active stand (cohort heights AND densities respond to the trait, odelia #32 +active-knot spline). So EVERY trait feeds back -- a trait that changes growth or +mortality moves the canopy everyone reads -- and the feedback routinely dominates +and flips the sign of the census metrics (LAI / biomass / size-moment) relative +to the frozen reading. Applies to the census metrics; \code{offspring_production} +stays \code{"frozen"} (the invasion gradient) even under \code{"resident"}. FF16 +only so far (#472 scope B). \strong{Recommended use} (scope decision): for the +census metrics (LAI / biomass / size_moment) \code{"resident"} is the correct +stand-level total gradient and \code{"frozen"} the rare-mutant invasion gradient; +for \code{offspring_production} the two coincide (it is always the invasion +gradient). The default is left at \code{"frozen"} so the function is a safe, +backward-compatible no-feedback derivative unless feedback is asked for. +\strong{Multi-species:} on a stand with more than one species \code{"resident"} +returns the CROSS-SPECIES total gradient \eqn{d(\mathrm{total\ stand\ metric})/ +d(\theta_{\mathrm{species}})} -- the joint canopy is rebuilt from every species' +re-evolved cohorts and the census metrics are summed over species, so the +differentiated species' traits feed back through the canopy that every species +reads (the cross term \code{"frozen"} sets to zero). The joint re-evolution is +well-conditioned on a \emph{fixed} node schedule; on a strongly clustered +adaptively-refined schedule it can go stiff, and the function then raises a clear +error (gated by a cheap baseline pass) asking for a fixed/uniform schedule rather +than returning a diverged result.} +} +\value{ +A list with \code{$jacobian} (a metrics x traits matrix) and + \code{$values} (the reconstructed metric values, which should match the SCM's + emergent outputs). +} +\description{ +Reverse-mode trait gradient of an SCM's emergent stand metrics (#472 scope B, +the calibration-facing generic engine, FF16). +} +\details{ +The generic counterpart of \code{\link{offspring_production_gradient}}: given a +resident \code{SCM} run with \code{control(save_RK45_cache = TRUE)}, it returns a +\strong{metrics x traits Jacobian} \eqn{d(\mathrm{metric}_m)/d(\theta_k)} for a +set of emergent stand metrics, computed from ONE resident baseline. Every metric +is a weighted reduction over the replayed cohorts, \eqn{\mathrm{metric} = \sum_i +w_i\, f(\mathrm{state}_i)}; \code{offspring_production} is just one such entry -- +none is privileged. The engine records one forward replay onto a single adjoint +tape and takes one cheap reverse sweep \emph{per metric}, so M metrics cost a +replay plus M sweeps, not M replays. This is the calibration core: \code{plant} +returns the Jacobian; a downstream package composes likelihoods (data never +enters here), which is what lets many likelihood terms share one stand baseline. + +Works for FF16, TF24 and TF24f residents (dispatched on the strategy). A + TF24/TF24f resident must have been run with \code{shading_model = "crown-centre"}. + FF16 supports all metrics (frozen + resident). TF24 currently supports + \code{"offspring_production"} (its census metrics need a leaf-optimisation + cross-sensitivity that is a follow-up). TF24f -- the fast-acclimation variant + whose tracked-collar leaf eval is analytic -- supports the CENSUS metrics + (\code{"LAI"} / \code{"biomass"} / \code{"size_moment"}) under \code{feedback = + "frozen"} via the reverse-mode AD census tape (the collar is carried as a taped + state with a curvature-linearised gradient-ascent rate); its resident (coupled) + census gradient and an offspring tape are follow-ups. \code{\link{stand_state_jacobian}} + works for FF16 and TF24. +} +\seealso{ +\code{\link{offspring_production_gradient}}, \code{\link{birth_rate_gradient}}. +} diff --git a/man/stand_state_jacobian.Rd b/man/stand_state_jacobian.Rd new file mode 100644 index 00000000..56e3425c --- /dev/null +++ b/man/stand_state_jacobian.Rd @@ -0,0 +1,43 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/emergent_gradient.R +\name{stand_state_jacobian} +\alias{stand_state_jacobian} +\title{Per-cohort state x trait Jacobian (FF16)} +\usage{ +stand_state_jacobian(scm, traits = NULL, species = 1L, birth_rate = NULL) +} +\arguments{ +\item{scm}{An \code{SCM} run with \code{save_RK45_cache = TRUE} (FF16 strategy).} + +\item{traits}{Character vector of FF16 trait names; \code{NULL} uses all 28.} + +\item{species}{Integer species index (see \code{\link{stand_gradient}}).} + +\item{birth_rate}{The (constant) birth-rate driver; recovered by default.} +} +\value{ +A list with \code{$states} (a cohort x component matrix of final-state + values) and \code{$jacobian} (a cohort x component x trait array). +} +\description{ +Per-cohort state x trait Jacobian of a resident SCM (#472 scope B, the +calibration-facing engine's escape hatch, FF16). +} +\details{ +The escape hatch for emergent metrics that are NOT a simple weighted reduction +(quantiles, ratios, bespoke statistics): rather than register a \code{(w, f)} +metric, this exposes \eqn{d(\mathrm{state}_{i,c})/d(\theta_k)} -- the derivative +of each replayed cohort's final demographic state component +\code{c} \eqn{\in} \{height, mortality, fecundity, area_heartwood, +mass_heartwood, offspring\} w.r.t. each trait. ANY smooth downstream metric over +the cohort states then has its gradient by the chain rule, with \code{plant} +never needing to know the metric -- the same boundary as "likelihoods live +downstream". Each cohort's final state is independent, so this tapes one cohort +at a time (one reverse sweep per state component). + +Works for both FF16 and TF24 residents (dispatched on the strategy). A + TF24 resident must have been run with \code{shading_model = "crown-centre"}. +} +\seealso{ +\code{\link{stand_gradient}}. +} diff --git a/notes/ad-refactor-optimize-roadmap.md b/notes/ad-refactor-optimize-roadmap.md new file mode 100644 index 00000000..9e065f7b --- /dev/null +++ b/notes/ad-refactor-optimize-roadmap.md @@ -0,0 +1,590 @@ +# AD gradient machinery: refactor, consolidate, optimize + +## Context + +Reverse-mode automatic differentiation of emergent trait gradients is **built and +validated end-to-end** for FF16, TF24, and TF24f on branch `spike-ff16-scm-emergent` +(101 commits, ~21,700 net lines over `develop`). The prototype phase is done; this plan +is the **refactor + optimize phase**. It does **not** add modelling capability — it +makes the validated machinery faster, smaller, more correct, and CI-covered. + +Three problems motivate it: + +1. **The R-side `Rcpp::as<>` env round-trip is the correctness ceiling.** The gradient + is two-pass: pass 1 runs a resident `run_scm(save_RK45_cache=TRUE)` and harvests a + frozen schedule + per-RK-stage environment; pass 2 is an AD replay that never re-runs + the solver. Pass-1 harvest happens **in R** (`ff16_harvest`, `tf24f_harvest`, …), + pulling `patch$environment_history` across the boundary and rebuilding each env via + `Rcpp::as<>`. That reconstruction is **not faithful** for the crown-sampled light + above cohort heights — it is the root cause of the TF24f census-recon fidelity floor + at patch lifetime > 4 and of a validation trap (a native-correct AD `g'` read as "18% + wrong" only through a round-tripped env). It also rebuilds the whole RcppR6 patch on + every access (O(stand size), ~1600× slower per `R/emergent_gradient.R:70-73`). + +2. **Duplication.** Three near-identical replay engines (`src/ff16_emergent.cpp` 1851, + `src/tf24_emergent.cpp` 479, `src/tf24f_emergent.cpp` 2266) plus three R harvest + seams. The Cash-Karp stepper and state containers are *already* shared; what is + triplicated is the **orchestration glue** (the `Frozen` struct, `build_frozen`/`as<>` + marshalling, per-strategy deriv/axpy/seed, `census_reduce`, coupled-canopy + reconstruction). + +3. **No benchmarks or CI coverage for the AD path.** There is no timing harness for the + gradient calls, and 7 FF16 AD test files `skip()` with "AD tape symbols unavailable + in this load_all session" — so the FF16 AD surface is not exercised in CI. ~31 + intermediate `ad_*`/`coupled_*`/`recon_*` validation scripts accumulated during the + prototype. + +**Intended outcome:** the AD path drives native environments (round-trip gone, the +long-horizon fidelity floor dissolved), one strategy-parameterized engine instead of +three, a benchmark + regression net that proves every step preserves the validated +Jacobians, and the spike merged to `develop` as a clean reviewable PR stack. + +### Decision on PR #541 / merge sequencing (resolved with Dan) + +**Do NOT gate this work on landing #541 first; close #541 without merging.** #541 +(`spike-ff16-hierarchy`, the `` templating) is already the ancestor base of this +spike — nothing depends on it being in `develop`, and its commits remain in the spike +history whether or not the PR is open. The only value of keeping it open was reviewer +ergonomics, which is better served by re-cutting fresh PRs from the *final* tree. +**Action: close PR #541 (do not merge)** — it isn't serving a purpose as a live review +object. We do all the work in place on the spike (no rebase, no merge conflict), then +decompose the final tree into a stacked PR sequence for `develop` (Phase 5), where a +foundation slice (the templating) can be cut fresh as PR #1. *(Closing the PR is a +mutating GitHub action — done on plan approval, not during planning.)* + +--- + +## Phase 1 — Safety net: benchmarks + correctness-regression fixture (do first) + +Establish the baseline *before* touching anything, so every later step is provably +non-regressing in both timing and value. Templates the existing +`scripts/bench_ab.R` / `scripts/bench_tf24.R` conventions. + +- **`scripts/bench_gradient.R`** — same-machine A/B harness (mirrors `bench_ab.R`: + `path`+`label` args, `pkgload::load_all(compile=FALSE)`, warm-up, one `RESULT|…` line + per build, interleave this-branch vs a worktree). The novel requirement: **time the + three sub-costs separately** by calling internals directly — + - `run_ms` = resident `run_scm(save_RK45_cache=TRUE)` (context), + - `harvest_ms` = the R-side `*_harvest()` call alone (**the refactor's target**), + - `impl_ms` = the C++ `*_impl` replay+sweep alone, + - headline `harvest_frac = harvest_ms/(harvest_ms+impl_ms)` and a `val=` + bit-identity check (à la `bench_tf24.R:31`). + Canonical cases lifted from test fixtures (fixed schedules): `ff16_frozen` (28 traits, + 4 metrics), `ff16_resident_coupled`, `ff16_resident_ms`, `ff16_offspring`, + `tf24f_census`, `tf24f_resident`, `grow_individual`. `BENCH_REPS`/`BENCH_CASES` env + vars; optional `BENCH_SCALING=1` to factor `impl_ms ≈ replay(traits) + M·sweep`. + +- **Correctness-regression fixture (AD-vs-AD baseline).** Snapshot the current validated + Jacobians (values + `d(metric)/d(trait)` matrices) to + `tests/testthat/fixtures/gradient-baseline.rds` for every distinct engine + (frozen / coupled-single / coupled-ms / state-jacobian / offspring / grow, across + FF16/TF24/TF24f). A new `tests/testthat/test-gradient-regression.R` reloads it and + asserts two-tier tolerance: **bit-identical** (`1e-12` rel) for pure-relocation steps, + **noise-floor** (`~5e-6` rel) for the coupled/ms paths that legitimately reorder FP + sums. This is complementary to the existing AD-vs-FD tests (those pin AD to physics at + ~1%; this pins AD to its own validated self at machine precision). Use fixed node + schedules so snapshot/check are reproducible. + +- **Extend `scripts/profile-benchmarks.R`** with gradient cases (`grad_ff16_frozen`, + `grad_ff16_resident`, `grad_tf24f_census`) reusing its Rprof + native `/usr/bin/sample` + machinery. Per the profile-plant skill: sample **ON** for hotspot localisation, **OFF** + for A/B ratios; only same-session ratios are trustworthy. + +- **Timing history:** recorded alongside the `bench_gradient.R` RESULT lines, columns + `date | step+sha | case | run_ms | harvest_ms | impl_ms | public_ms | harvest_frac | + cum_speedup | incr_speedup | sample | fixture(PASS/FAIL + worst rel_dev) | notes`. + No speedup row is recorded without its correctness verdict attached. + +**Validates:** the fixture itself is the validation infrastructure; confirm it passes on +the current HEAD and that `harvest_ms + impl_ms ≈ public_ms`. + +--- + +## Phase 2 — Native-env harvest in C++ (kills the `Rcpp::as<>` round-trip) + +> **UPDATE 2026-06-30 (measured, supersedes the fidelity rationale below).** The premise +> that the `Rcpp::as<>` env round-trip causes the TF24f census floor is **NOT supported by +> measurement**. A lossy-vs-native A/B of the census recon gives bit-identical LAI at every +> horizon (H=4 rel 6.3e-7; H=5 3.98e-3; H=8 1.19e-2, no abort) — `native == lossy` exactly. +> The TF24_Environment round-trip is faithful for the census path (knots fully determine the +> spline; `compute_competition(0)` is read at z=0, inside the domain). **The lifetime>4 floor +> is the g' backward-FD near-cancellation + frozen-schedule replay, not the env.** Phases 2a +> (FF16) and 2b (TF24f census) are DONE and committed as **bit-identical refactors** that +> remove a real serialization round-trip + the documented env-validation hazard and unlock the +> `harvest_frac` perf headroom — **not** a fidelity fix. Lifting the census floor is a separate +> item (exact-AD g' / refined-schedule replay). No false lifetime gate added. See memory +> `ad-env-roundtrip-faithful-for-census`. + +The highest-value change: it has a **correctness payoff** (removes the lifetime-floor and +the validation trap) and shrinks the per-engine surface Phase 3 then unifies. The +round-trip is structurally avoidable — confirmed: the replay reads the env only as a +*double* value + analytic derivative (`src/ff16_emergent.cpp:125-128`), lifting just those +scalars into the AD type; the AD type never touches the env. A faithful native env +**pointer** is sufficient and exact. + +**Approach (engine stays OUTSIDE the SCM object).** Do **not** add a `collect_gradient` +run mode to `SCM::run()` — the harvest is +already captured during the normal `save_RK45_cache=TRUE` run. Instead: + +- Add `inst/include/plant/gradient/resident_harvest.h` with a `ResidentHarvest` struct + and a free `harvest_resident(const Patch&, int species)` builder that gathers, + **by const-ref/pointer into the Patch's own storage**, the pieces R rebuilds today: + `step_history`, `&environment_history[n][s]` (native pointers, not copies), + `stand_*_stage_history`, `stand_newnode_*`, and per-species cohort metadata. Move the + two remaining R computations — trapezoid node weights (`emergent_gradient.R:94-97`) and + the per-RK-stage `pr_survival` matrix (`:99-103`) — into this C++ builder. +- Change each engine's `Frozen.eh` from owning `std::vector<…Environment>` to holding + `const Environment*` borrowed from the Patch. **Delete the `Rcpp::as<>` loop** + (`ff16_emergent.cpp:450`, and the TF24/TF24f equivalents). This line is the round-trip. +- The C++ entry takes the **live `Patch&`** from the RcppR6 external pointer (env pointers + are only valid while the SCM C++ object lives; RcppR6 keeps it alive in R). The R API + (`stand_gradient()` etc.) becomes a thin forwarder passing the live SCM xptr down. + +**Sub-steps, each bit-checkable against Phase 1's fixture:** +1. **2a** — `ResidentHarvest` + C++ gather + the two moved arithmetic computations, still + building `Frozen` from copied envs. **Bit-identical** to current `ff16_harvest` + outputs (assert `tw`, `ppsurv`, sampled env values equal). Isolates the gather move. +2. **2b** — switch `Frozen.eh` to borrowed native pointers; delete `as<>`. Bit-identical + at lifetime ≤ 4; at **lifetime 5–8** the recon LAI must now match + `compute_competition(0)` to ~1e-6 and TF24f must no longer abort at the spline edge — + add new long-horizon gates to `test-tf24f-census-gradient.R` that were previously + impossible. **This is the headline correctness win.** +3. **2c** — repeat for TF24 / TF24f env entries. + +**Files:** new `inst/include/plant/gradient/resident_harvest.h`; `inst/include/plant/patch.h` +(public accessors already exist, `:153`, `:176-218`); `src/ff16_emergent.cpp`, +`src/tf24_emergent.cpp`, `src/tf24f_emergent.cpp` (`Frozen`, `build_frozen`, `attach_*`); +`R/emergent_gradient.R`, `R/tf24_emergent_gradient.R`, `R/tf24f_emergent_gradient.R` +(harvests shrink to forwarders). `make rebuild` if an export signature changes. + +**Validates:** all `test-{ff16,tf24,tf24f}-*gradient*.R` stay green; the new lifetime-5/8 +TF24f gates pass; the Phase-1 fixture confirms ≤4-horizon Jacobians unchanged; record the +`harvest_frac` drop in the timing-history file. + +--- + +## Phase 3 — Unify the three replay engines (maximal dedup, per Dan) + +> **OUTCOME 2026-06-30 (safe dedup done; finding on the limit of further unification).** +> Executed the safe, high-value, bit-identical shares and measured the result. What was +> genuinely duplicated and is now shared in `inst/include/plant/gradient/`: +> - `scm_harvest.h`: `recover_birth_rate` (was inlined 5×), `birth_steps` (2×), +> `census_trapezium` (the descending-height census reduction, FF16 + TF24f). +> - `coupled_canopy.h`: `canopy_comp_at` (the Yokozawa light trapezium, was hand-copied +> as ff16 `coupled_comp_at` + tf24f `tf24f_comp_at`; 4 call sites). +> - The Cash-Karp stepper (`ff16_cashkarp_replay`) + state containers were **already** +> shared before this work. +> +> All bit-identical (fixture max_rel 0.0; one 9e-16 FP-reassociation on frozen census). +> **Key finding:** the engines shrank only ~107 lines for ~132 lines of shared header — +> i.e. the duplicated *orchestration* was small relative to the **irreducible +> strategy-specific physiology** (FF16's light-response hyperbola + deep-crown GK vs +> TF24f's tracked-collar leaf solve + curvature harvest; the coupled `deep_net_coupled` +> / deriv kernels). So the "4900 → 2000" estimate was optimistic: after sharing the +> stepper/harvest/census/canopy, the bulk of each engine is **distinct biology, not +> redundant copy**. A generic `replay_cohort_final` trait-class would add +> machinery to abstract over genuinely-different deriv kernels for modest line savings, +> and the multi-species coupled unification is high-risk (the stiffness/conditioning +> guards) for low dedup value (the MS engines are largely distinct). **Recommendation: +> stop the engine unification here** — the redundant duplication is removed; what remains +> is either irreducible or net-negative to force. The remaining `build_interp` +> single-vs-MS intra-file overlap is a small (~25-line) optional same-file cleanup. +> Original (pre-execution) plan kept below for context. + +Collapse the triplicated orchestration into one strategy-parameterized engine. The +stepper (`ff16_cashkarp_replay`) and state containers are already shared across all three +(`tf24_emergent.cpp` and `tf24f_emergent.cpp` already `#include ff16_production_kernel.h`). +The engines differ in `replay_cohort_final` at exactly three points: the env→net call, the +state-vector width (5/6/7), and the rate fill. + +**Minimal trait-class interface** each strategy supplies (mostly *binding* things already +templated — `FF16ProdPars`/`TF24ProdPars` and `*_compute_rates_from_net` already +exist with matching shapes): +``` +struct StrategyTraits { + using ProdPars = …; // existing + using LifeState = …; // strategy-supplied width (TF24f adds collar + log_density) + static S net_at(ProdPars, Env*, LifeState, /*harvest@n,stage*/); + static Rates rates_from_net(ProdPars, height, area_leaf, net); + static LifeState seed(ProdPars, Env* birth_env, h0, …); // establishment / IFT collar birth +}; +``` +TF24f's collar-as-6th-state + curvature harvest fits as a **hook**, not a fork: the collar +is an extra `LifeState` component with rate `k_acclim·dprofit_dpsi`; the curvature harvest +is a *harvest producer* orthogonal to the replay template. **Do not** merge TF24's +envelope-zeroed collar with TF24f's taped collar — they are different math; keep them as +two `seed`/`net_at` implementations behind the hook. + +**Ordered dedup (each gated by Phase-1 fixture); Dan chose to push through all ranks:** +1. **Shared `census_reduce` + `(w,f)` metric dispatch** — pure arithmetic, currently + hand-copied identically (`ff16:587`, `tf24f:226`). **Bit-identical.** (~300 lines) +2. **Generic `replay_cohort_final`** via the trait class (frozen/invasion path). + FF16/TF24 offspring are bit-exact references. **Bit-identical.** (~400-600 lines) +3. **Generic coupled-canopy reconstruction** — the Yokozawa trapezium + `Q=(1-(z/h)^eta)^2` is identical across strategies (`coupled_comp_at`/`build_interp` + vs `tf24f_comp_at`/`build_canopy`); only the per-cohort rate differs. Validated to the + ~1e-4 FD-noise floor + sign-flip invariants (not bit-identical — the coupled tape is + numerically delicate). (~500-800 lines) +4. **Multi-species engine** (`assemble_metrics_coupled_ms`/`FrozenMS`). **Highest risk** — + the MS stiffness/conditioning guards (joint-env-drift > 1e-3 gate, finiteness guard, + tape scoping) must be carried through **verbatim**; a refactor that drops them silently + re-opens divergence. Validate against the MS-R0/MS-cross-species-R1 gates on fast fixed + schedules. + +**Files:** `inst/include/plant/models/ff16_production_kernel.h` (shared `census_reduce`, +generic `replay_cohort_final`, trait base), new +`inst/include/plant/gradient/replay_engine.h` (strategy-agnostic orchestration); the three +`src/*_emergent.cpp` shrink to: include the generic engine + supply the strategy trait +class + the `[[Rcpp::export]]` entry points. + +**Validates:** ranks 1-2 bit-for-bit; rank 3 at ~1e-4 + sign-flips; rank 4 against the MS +gates. Run the **whole** test suite each step (generic tests loop over +`helper-plant.R` strategy lists; remember TF24f is the non-5-state variant excluded +there). Record line-count reduction + timing per step. + +--- + +## Phase 4 — Cleanup: fold validation scripts into tests, cover the AD surface in CI + +Per Dan: **fold the still-valuable validation scripts into proper tests**, then delete the +scripts (rather than deleting outright or leaving them as loose scripts). + +- **Migrate** the meaningful checks from `scripts/coupled_r0_check.R`, + `coupled_r1_check.R`, `recon_static_check.R`, `coupled_vs_full_scm_fd.R`, + `coupled_multispp_r0.R`, `coupled_multispp_r1.R`, `coupled_gap_decompose.R`, + `coupled_birthenv_channel.R`, etc. into `testthat` tests (or fold into the existing + `test-{ff16,tf24f}-*gradient*.R`). Delete the migrated scripts. Keep a small curated set + of `ad_*` demos as runnable examples; the `overstorey-staging` guide `.qmd` remains the + user-facing narrative. +- **Make the AD tests actually run in CI (enabler — otherwise folded tests just skip).** + The 7 FF16 AD test files skip with "AD tape symbols unavailable in this load_all + session"; TF24f tests already run live. Diagnose the asymmetry and make the FF16 AD + symbols available under the test load path — tie to the **odelia DLL load-ordering** + fix (plant needs `importFrom(odelia,…)` so odelia's `.onLoad` runs before plant's + `useDynLib`) and the build-optimized-then-`load_all(compile=FALSE)` path. The aim: + folded + existing AD tests execute in CI rather than skip. (Aligns with roadmap Phase C: + "make the headline result CI-testable in plain R".) + +**Validates:** the migrated tests pass under `make test` / `devtools::test()` *without* +skipping; the AD surface is genuinely covered. Confirm no orphaned `source(scripts/…)` +references remain. + +--- + +## Phase 5 — Decompose for merge to `develop` + +PR #541 is **closed (not merged)** at the start of this work — its commits stay as spike +ancestors. Re-cut the final spike tree into a stacked, reviewable PR sequence (this is +*merge* decomposition, independent of the in-place writing order above): + +1. **Foundation** (`` templating + exact AD growth-rate gradient) — additive, + bit-identical to current package; cut fresh from the final tree, the clean foundation + reviewers vet first (replaces the role the old #541 would have played). +2. **Native-env harvest** (Phase 2) — the correctness fix; land even if later phases slip. +3. **Engine unification** (Phase 3 ranks 1-2 bit-identical; rank 3 noise-gated; rank 4 MS). +4. **Cleanup + CI coverage** (Phase 4). + +Each PR carries its slice of the timing-history table and its fixture verdict. + +--- + +## Verification (end-to-end) + +- **Build discipline:** C++-only change → `make compile`; new export / RcppR6 field → + `make rebuild`; then `pkgload::load_all(".", compile=FALSE)` (optimized -O2, **never** + `devtools::load_all()` which is -O0). Run the **whole** plant suite before each PR + (`make test`), not just `test-strategy-*` — generic tests loop over `helper-plant.R`. +- **Correctness gate (every step):** `Rscript scripts/gradient_fixture.R check` (or the + `test-gradient-regression.R` testthat run) — bit-identical for relocation/rank-1-2 + steps, noise-floor for coupled/ms. Plus the existing AD-vs-FD tests stay green. +- **The Phase-2 headline proof:** new lifetime-5 and lifetime-8 TF24f census gates that + were impossible under the round-trip now pass to ~1e-6 vs `compute_competition(0)`. +- **Timing (same-session A/B):** build this branch + a worktree of the pre-step ref, + interleave `Rscript scripts/bench_gradient.R