Integrate over height or leaf area of an individual?

[dfalster]

384 seeks to simplify methods used to create spline of light_availability and integrate over it.

In thinking about this, I came to think of an older method we used to calculate assimilation and wondering if we should return to this. I'm not sure when it was removed, but it's present here https://github.com/traitecoevo/plant/blob/34de32261b18089b24b7676bd1aa1f6e1a3649d9/src/ff16_strategy.cpp#L265

The choice is between

Method 1 (current):

  // Define an anonymous function to integrate
  // For given height in crown z, take photosynthesis at depth z multiplied by 
  //   amount of leaf at that depth
  std::function<double(double)> f1 = [&](double z) -> double {
    return assimilation_leaf(environment.get_environment_at_height(z)) * q(z, height);
  };

  // Integrate over crown depth
  A1 = function_integrator.integrate(f1, 0.0, height);

or

Method 2 (proposed):

  // For given proportion of leaf area, calculate the height z where that occurs in plant height =height (via function Qp), then light at that depth and photosynthesis at that point
  //   amount of leaf at that depth
  std::function<double(double)> f2 = [&](double p) -> double {
    return assimilation_leaf(environment.get_environment_at_height(Qp(p, height)));
  };

  double A2 = function_integrator.integrate(f2, 0.0, 1.0);

the key difference is that

1. method 1 integrates over the height of the plant, , i.e. from 0-height. The points used in the integration are distributed in proportion to height. So for a top weighted canopy, more effort is being used at lower heights, where there isn't much leaf area.
2. method 1 integrates over the proportion of leaf area in the plant, i.e. from 0-1. The points used in the integration area distributed in proportion to leaf area. So no matter what the leaf area distribution in the plant canopy, effort is deployed at the heights where the leaf area actually occurs.

I confirm that the two methods give almost identical results and run with similar speed, but they do differ subtly in results when in competition, as we would expect them to. So values in tests need to be adjusted for tests to pass. E.g.

Failure (test-scm-support.R:100:3): collect_auxiliary_variables
as.numeric(state[8:9, 142, 141]) not equal to c(0.0007211209, 0.0001707292).
2/2 mismatches (average diff: 2.34e-07)
[1] 0.000721 - 0.000721 == -3.22e-07
[2] 0.000171 - 0.000171 == -1.45e-07

Failure (test-strategy-ff16-reference-comparison.R:115:3): Reference comparison
p$rate("height") not equal to `cmp_height_dt`.
1/1 mismatches
[1] 0.732 - 0.732 == 4.02e-05

Any thoughts @itowers1 @aornugent ? I prefer method 2 and propose we switch to that, in a seperate commit after #384 is merged.

Alternatively, we could implement this for new growth model, that will incorporate a variety of changes, including water usage. See #387

[aornugent]

I support the the change but don't totally understand the motivation.

Your point about expending compute time in proportion to leaf area makes sense, but that doesn't seem relevant if the timing is similar.

Conceptual clarity is a perfectly reasonable reason if that's what you're after.

On Wed, 31 Jan 2024, 8:54 am Daniel Falster, @.***> wrote:

384 https://github.com/traitecoevo/plant/pull/384 seeks to simplify

methods used to create spline of light_availability and integrate over it.

In thinking about this, I came to think of an older method we used to calculate assimilation and wondering if we should return to this.

The choice is between

Method 1 (current):

// Define an anonymous function to integrate // For given height in crown z, take photosynthesis at depth z multiplied by // amount of leaf at that depth std::function<double(double)> f1 = [&](double z) -> double { return assimilation_leaf(environment.get_environment_at_height(z)) * q(z, height); };

// Integrate over crown depth A1 = function_integrator.integrate(f1, 0.0, height);

or

Method 2 (proposed):

// For given proportion of leaf area, calculate the height z where that occurs in plant height =height (via function Qp), then light at that depth and photosynthesis at that point // amount of leaf at that depth std::function<double(double)> f2 = [&](double p) -> double { return assimilation_leaf(environment.get_environment_at_height(Qp(p, height))); };

double A2 = function_integrator.integrate(f2, 0.0, 1.0);

the key difference is that

1. method 1 integrates over the height of the plant, , i.e. from 0-height. The points used in the integration are distributed in proportion to height. So for a top weighted canopy, more effort is being used at lower heights, where there isn't much leaf area.
2. method 1 integrates over the proportion of leaf area in the plant, i.e. from 0-1. The points used in the integration area distributed in proportion to leaf area. So no matter what the leaf area distribution in the plant canopy, effort is deployed at the heights where the leaf area actually occurs.

I confirm that the two methods give almost identical results and run with similar speed, but they do differ subtly in results when in competition, as we would expect them to. So values in tests need to be adjusted for tests to pass. E.g.

Failure (test-scm-support.R:100:3): collect_auxiliary_variables as.numeric(state[8:9, 142, 141]) not equal to c(0.0007211209, 0.0001707292). 2/2 mismatches (average diff: 2.34e-07) [1] 0.000721 - 0.000721 == -3.22e-07 [2] 0.000171 - 0.000171 == -1.45e-07

Failure (test-strategy-ff16-reference-comparison.R:115:3): Reference comparison p$rate("height") not equal to cmp_height_dt. 1/1 mismatches [1] 0.732 - 0.732 == 4.02e-05

Any thoughts @itowers1 https://github.com/itowers1 @aornugent https://github.com/aornugent ? I prefer method 2 and propose we switch to that, in a seperate commit after #384 https://github.com/traitecoevo/plant/pull/384 is merged.

— Reply to this email directly, view it on GitHub https://github.com/traitecoevo/plant/issues/386, or unsubscribe https://github.com/notifications/unsubscribe-auth/AE3RPMK5Y5X3IEUXJIZTGQDYRFT2PAVCNFSM6AAAAABCR7ZXWKVHI2DSMVQWIX3LMV43ASLTON2WKOZSGEYDQOBVGE4DEMA . You are receiving this because you were mentioned.Message ID: @.***>

[itowers1]

I would support moving to option 2, although results are similar in this case, I would suspect you might get greater divergence in outcomes for more extreme canopy shapes (such as an almost flat-topped tree). Presumably, this is also what we want to do for roots, so good to be consistent as well.

[dfalster]

Thanks. I'll try to explain the difference better. Both are using the same number of points, hence the same speed. The difference is where the points are situated within the canopy of the plant. In the plot below, the line shows the density of leaf area at height (z) vs height (z). This is the default crown profile with eta -12. The red points are option 1. There's a lot of points where there is very little leaf area. The green is option 2, points are distributed with respect to leaf area. The last green point is at 3m, but 99.7% of leaf area is above that.

With option 1 we only have 10 points from 3-5m, so arguably this is more like a 10 point integration, with 11 wasted points. In option 2, all 21 points are inlacing the end integral equally.

Note, in both cases, the spacing is what is actually used in the integration, determined by the Gauss-Kronrod quadrature routine (with rule 21). which is not evenly spaced.

# abscissa of QK21 routine
ab <- c(
  0.995657163025808080735527280689003,
  0.973906528517171720077964012084452,
  0.930157491355708226001207180059508,
  0.865063366688984510732096688423493,
  0.780817726586416897063717578345042,
  0.679409568299024406234327365114874,
  0.562757134668604683339000099272694,
  0.433395394129247190799265943165784,
  0.294392862701460198131126603103866,
  0.148874338981631210884826001129720,
  0.000000000000000000000000000000000
)

# This is how spacing is calculated 
x_GK21 <- 0.5 * c(1 - ab, 1 + ab) |> sort() |> unique()

eta <- 12

Qp <- function(p, height) { 
  return ((1 - sqrt(p))^(1/eta) * height);
}

Q <- function(z, height) {
  if (z > height) {
    return(0.0);
  }
  tmp <- 1.0-pow(z / height, eta);
  return(tmp * tmp);
}

q <- function(z, height) {
  tmp <- (z / height)^eta
  q <- 2 * eta * (1 - tmp) * tmp / z
  q[z==0] <- 0
  return(q)
}

height <- 5
z <- seq(0, height, length.out = 100)

# Leaf area distribution
plot(q(z, height), z, type = "l", ylim = c(-1, 6), col= "grey")

# options 1 - Intergrate wrt z, points evenly spaced along z
z_s <- x_GK21 * height
points(q(z_s, height), z_s, col = "red", pch = 16)

# options 2 - Integrate wrt p, proportion of leaf area
p <- x_GK21
z_s <- Qp(p, height)
points(q(z_s, height), z_s, col = "green", pch = 16)