More tweaks

NEWS.md

@@ -34,6 +34,9 @@ be directly added to this file to describe the related changes.
 
 ## UNRELEASED
 
+- Added new options for adding penalties to the error function during Variable J
+  fits, that enable the user to selectively penalize negative or unreasonably
+  large values of mesophyll conductance.
 - Made a few improvements to C3 curve fitting functions (`fit_c3_aci` and
   `fit_c3_variable_j`):
   - Two more parameters can now be fit: `alpha` (related to TPU) and
@@ -64,11 +67,34 @@ be directly added to this file to describe the related changes.
       fixed values or to values from a column of an exdf object.
     - Unknown parameters are now provided in alphabetical order when applicable
       (such as the first input arguments to `calculate_c3_assimilation`).
+  - The functions are more tolerant to curves with severe problems (such as
+    negative Ci) that prevent a good fit from being found; rather than throwing
+    an error, the fit functions now silently return `NA` for all results, along
+    with a message explaining the issue.
 - Added a function for estimating `Rd` with the Laisk method:
   `calculate_rd_laisk`
 - A "unit dictionary" was added for internal use; this may be expanded and used
   more often in the future.
-- Tests for `calculate_c3_assimilation` were added.
+- Renamed several variables and input arguments:
+  - Licor files contain a column called `alpha`, and several different "alphas"
+    were used throughout `PhotoGEA`. To avoid confusion, the values in
+    `PhotoGEA` were renamed as follows:
+    - `alpha_g`: used in C3 assimilation calculations
+    - `alpha_pr`: used in Gamma_star calculations
+    - `alpha_psii`: used in C4 assimilation calculations
+  - The acronym "TPU" was used to refer to a process (triose phosphate
+    utilization) and the maximum rate of that process. To avoid confusion, the
+    rate parameter was renamed to `Tp`
+- Tests were added for several functions:
+  - `calculate_c3_assimilation`
+  - `fit_c3_aci`
+  - `fit_c3_variable_j`
+  - `fit_c4_aci`
+  - `calculate_c3_limitations_grassi`
+  - `calculate_c3_limitations_warren`
+- The `read_gasex_file` function now automatically includes the filename as a
+  column in the resulting `exdf` object; this helps with troubleshooting
+  problematic curves or files.
 - PRs related to creating this version:
   - https://github.com/eloch216/PhotoGEA/pull/85
   - https://github.com/eloch216/PhotoGEA/pull/86

R/calculate_c3_assimilation.R

@@ -1,10 +1,10 @@
 calculate_c3_assimilation <- function(
     exdf_obj,
-    alpha,        # dimensionless      (this value is sometimes being fitted)
+    alpha_g,      # dimensionless      (this value is sometimes being fitted)
     Gamma_star,   # micromol / mol     (this value is sometimes being fitted)
     J_at_25,      # micromol / m^2 / s (at 25 degrees C; typically this value is being fitted)
     Rd_at_25,     # micromol / m^2 / s (at 25 degrees C; typically this value is being fitted)
-    TPU,          # micromol / m^2 / s (typically this value is being fitted)
+    Tp,           # micromol / m^2 / s (typically this value is being fitted)
     Vcmax_at_25,  # micromol / m^2 / s (at 25 degrees C; typically this value is being fitted)
     POc = 210000, # microbar           (typically this value is known from the experimental setup)
     atp_use = 4.0,
@@ -38,11 +38,11 @@
         required_variables[[vcmax_norm_column_name]]     <- 'normalized to Vcmax at 25 degrees C'
 
         flexible_param <- list(
-            alpha = alpha,
+            alpha_g = alpha_g,
             Gamma_star = Gamma_star,
             J_at_25 = J_at_25,
             Rd_at_25 = Rd_at_25,
-            TPU = TPU,
+            Tp = Tp,
             Vcmax_at_25 = Vcmax_at_25
         )
 
@@ -70,11 +70,11 @@
     }
 
     # Retrieve values of flexible parameters as necessary
-    if (!value_set(alpha))       {alpha       <- exdf_obj[, 'alpha']}
+    if (!value_set(alpha_g))     {alpha_g     <- exdf_obj[, 'alpha_g']}
     if (!value_set(Gamma_star))  {Gamma_star  <- exdf_obj[, 'Gamma_star']}
     if (!value_set(J_at_25))     {J_at_25     <- exdf_obj[, 'J_at_25']}
     if (!value_set(Rd_at_25))    {Rd_at_25    <- exdf_obj[, 'Rd_at_25']}
-    if (!value_set(TPU))         {TPU         <- exdf_obj[, 'TPU']}
+    if (!value_set(Tp))          {Tp          <- exdf_obj[, 'Tp']}
     if (!value_set(Vcmax_at_25)) {Vcmax_at_25 <- exdf_obj[, 'Vcmax_at_25']}
 
     # Extract a few columns from the exdf object to make the equations easier to
@@ -92,23 +92,27 @@
     Rd_tl <- Rd_at_25 * exdf_obj[, rd_norm_column_name]          # micromol / m^2 / s
     J_tl <- J_at_25 * exdf_obj[, j_norm_column_name]             # micromol / m^2 / s
 
-    # Make sure key inputs have reasonable values; these checks cannot be
-    # bypassed
+    # Make sure key inputs have reasonable values
     msg <- character()
 
-    if (any(alpha < 0 | alpha > 1)) {msg <- append(msg, 'alpha must be >= 0 and <= 1')}
-    if (any(Cc < 0))                {msg <- append(msg, 'Cc must be >= 0')}
-    if (any(Gamma_star < 0))        {msg <- append(msg, 'Gamma_star must be >= 0')}
-    if (any(J_at_25 < 0))           {msg <- append(msg, 'J_at_25 must be >= 0')}
-    if (any(Kc < 0))                {msg <- append(msg, 'Kc must be >= 0')}
-    if (any(Ko < 0))                {msg <- append(msg, 'Ko must be >= 0')}
-    if (any(pressure < 0))          {msg <- append(msg, 'pressure must be >= 0')}
-    if (any(Rd_at_25 < 0))          {msg <- append(msg, 'Rd_at_25 must be >= 0')}
-    if (any(TPU < 0))               {msg <- append(msg, 'TPU must be >= 0')}
-    if (any(Vcmax_at_25 < 0))       {msg <- append(msg, 'Vcmax_at_25 must be >= 0')}
-
-    if (length(msg) > 0) {
-        stop(paste(msg, collapse = '. '))
+    if (any(alpha_g < 0 | alpha_g > 1, na.rm = TRUE)) {msg <- append(msg, 'alpha_g must be >= 0 and <= 1')}
+    if (any(Cc < 0, na.rm = TRUE))                    {msg <- append(msg, 'Cc must be >= 0')}
+    if (any(Gamma_star < 0, na.rm = TRUE))            {msg <- append(msg, 'Gamma_star must be >= 0')}
+    if (any(J_at_25 < 0, na.rm = TRUE))               {msg <- append(msg, 'J_at_25 must be >= 0')}
+    if (any(Kc < 0, na.rm = TRUE))                    {msg <- append(msg, 'Kc must be >= 0')}
+    if (any(Ko < 0, na.rm = TRUE))                    {msg <- append(msg, 'Ko must be >= 0')}
+    if (any(pressure < 0, na.rm = TRUE))              {msg <- append(msg, 'pressure must be >= 0')}
+    if (any(Rd_at_25 < 0, na.rm = TRUE))              {msg <- append(msg, 'Rd_at_25 must be >= 0')}
+    if (any(Tp < 0, na.rm = TRUE))                    {msg <- append(msg, 'Tp must be >= 0')}
+    if (any(Vcmax_at_25 < 0, na.rm = TRUE))           {msg <- append(msg, 'Vcmax_at_25 must be >= 0')}
+
+    msg <- paste(msg, collapse = '. ')
+
+    # We only bypass these checks if !perform_checks && return_exdf
+    if (perform_checks || !return_exdf) {
+        if (msg != '') {
+            stop(msg)
+        }
     }
 
     # Rubisco-limited carboxylation (micromol / m^2 / s)
@@ -118,8 +122,8 @@
     Wj <- PCc * J_tl / (atp_use * PCc + nadph_use * Gamma_star)
 
     # TPU-limited carboxylation (micromol / m^2 / s)
-    Wp <- PCc * 3 * TPU / (PCc - Gamma_star * (1 + 3 * alpha))
-    Wp[PCc <= Gamma_star * (1 + 3 * alpha)] <- Inf
+    Wp <- PCc * 3 * Tp / (PCc - Gamma_star * (1 + 3 * alpha_g))
+    Wp[PCc <= Gamma_star * (1 + 3 * alpha_g)] <- Inf
 
     # Co-limitation between Wc and Wj
     a_cj <- curvature_cj
@@ -158,11 +162,11 @@
         # Make a new exdf object from the calculated variables and make sure units
         # are included
         output <- exdf(data.frame(
-            alpha = alpha,
+            alpha_g = alpha_g,
             Gamma_star = Gamma_star,
             J_tl = J_tl,
             Rd_tl = Rd_tl,
-            TPU = TPU,
+            Tp = Tp,
             Vcmax_tl = Vcmax_tl,
             Ac = Ac,
             Aj = Aj,
@@ -171,25 +175,27 @@
             Wc = Wc,
             Wj = Wj,
             Wp = Wp,
-            Vc = Wcjp
+            Vc = Wcjp,
+            c3_assimilation_msg = msg
         ))
 
         document_variables(
             output,
-            c('calculate_c3_assimilation', 'alpha',      'dimensionless'),
-            c('calculate_c3_assimilation', 'Gamma_star', 'micromol mol^(-1)'),
-            c('calculate_c3_assimilation', 'J_tl',       'micromol m^(-2) s^(-1)'),
-            c('calculate_c3_assimilation', 'Rd_tl',      'micromol m^(-2) s^(-1)'),
-            c('calculate_c3_assimilation', 'TPU',        'micromol m^(-2) s^(-1)'),
-            c('calculate_c3_assimilation', 'Vcmax_tl',   'micromol m^(-2) s^(-1)'),
-            c('calculate_c3_assimilation', 'Ac',         'micromol m^(-2) s^(-1)'),
-            c('calculate_c3_assimilation', 'Aj',         'micromol m^(-2) s^(-1)'),
-            c('calculate_c3_assimilation', 'Ap',         'micromol m^(-2) s^(-1)'),
-            c('calculate_c3_assimilation', 'An',         'micromol m^(-2) s^(-1)'),
-            c('calculate_c3_assimilation', 'Wc',         'micromol m^(-2) s^(-1)'),
-            c('calculate_c3_assimilation', 'Wj',         'micromol m^(-2) s^(-1)'),
-            c('calculate_c3_assimilation', 'Wp',         'micromol m^(-2) s^(-1)'),
-            c('calculate_c3_assimilation', 'Vc',         'micromol m^(-2) s^(-1)')
+            c('calculate_c3_assimilation', 'alpha_g',             'dimensionless'),
+            c('calculate_c3_assimilation', 'Gamma_star',          'micromol mol^(-1)'),
+            c('calculate_c3_assimilation', 'J_tl',                'micromol m^(-2) s^(-1)'),
+            c('calculate_c3_assimilation', 'Rd_tl',               'micromol m^(-2) s^(-1)'),
+            c('calculate_c3_assimilation', 'Tp',                  'micromol m^(-2) s^(-1)'),
+            c('calculate_c3_assimilation', 'Vcmax_tl',            'micromol m^(-2) s^(-1)'),
+            c('calculate_c3_assimilation', 'Ac',                  'micromol m^(-2) s^(-1)'),
+            c('calculate_c3_assimilation', 'Aj',                  'micromol m^(-2) s^(-1)'),
+            c('calculate_c3_assimilation', 'Ap',                  'micromol m^(-2) s^(-1)'),
+            c('calculate_c3_assimilation', 'An',                  'micromol m^(-2) s^(-1)'),
+            c('calculate_c3_assimilation', 'Wc',                  'micromol m^(-2) s^(-1)'),
+            c('calculate_c3_assimilation', 'Wj',                  'micromol m^(-2) s^(-1)'),
+            c('calculate_c3_assimilation', 'Wp',                  'micromol m^(-2) s^(-1)'),
+            c('calculate_c3_assimilation', 'Vc',                  'micromol m^(-2) s^(-1)'),
+            c('calculate_c3_assimilation', 'c3_assimilation_msg', '')
         )
     } else {
         return(list(An = An, Ac = Ac, Aj = Aj, Ap = Ap, Wc = Wc, Wj = Wj))

R/calculate_c3_limitations_warren.R

@@ -7,21 +7,16 @@ calculate_c3_limitations_warren <- function(
     POc = 210000,
     atp_use = 4.0,
     nadph_use = 8.0,
-    alpha = 0.0,
     curvature_cj = 1.0,
     curvature_cjp = 1.0,
     ca_column_name = 'Ca',
     cc_column_name = 'Cc',
     ci_column_name = 'Ci',
-    j_column_name = 'J_at_25',
     j_norm_column_name = 'J_norm',
     kc_column_name = 'Kc',
     ko_column_name = 'Ko',
-    rd_column_name = 'Rd_at_25',
     rd_norm_column_name = 'Rd_norm',
     total_pressure_column_name = 'total_pressure',
-    tpu_column_name = 'TPU',
-    vcmax_column_name = 'Vcmax_at_25',
     vcmax_norm_column_name = 'Vcmax_norm'
 )
 {
@@ -31,29 +26,29 @@ calculate_c3_limitations_warren <- function(
     }
 
     # Make sure the required variables are defined and have the correct units.
-    # We don't need to check all of them, since calculate_c3_assimilation will
-    # perform checks.
     required_variables <- list()
-    required_variables[['Gamma_star']]      <- unit_dictionary[['Gamma_star']]
-    required_variables[[ca_column_name]]    <- 'micromol mol^(-1)'
-    required_variables[[cc_column_name]]    <- 'micromol mol^(-1)'
-    required_variables[[ci_column_name]]    <- 'micromol mol^(-1)'
-    required_variables[[j_column_name]]     <- 'micromol m^(-2) s^(-1)'
-    required_variables[[rd_column_name]]    <- 'micromol m^(-2) s^(-1)'
-    required_variables[[tpu_column_name]]   <- 'micromol m^(-2) s^(-1)'
-    required_variables[[vcmax_column_name]] <- 'micromol m^(-2) s^(-1)'
+    required_variables[['alpha_g']]                  <- unit_dictionary[['alpha_g']]
+    required_variables[['Gamma_star']]               <- unit_dictionary[['Gamma_star']]
+    required_variables[['J_at_25']]                  <- unit_dictionary[['J_at_25']]
+    required_variables[['Rd_at_25']]                 <- unit_dictionary[['Rd_at_25']]
+    required_variables[['Tp']]                       <- unit_dictionary[['Tp']]
+    required_variables[['Vcmax_at_25']]              <- unit_dictionary[['Vcmax_at_25']]
+    required_variables[[ca_column_name]]             <- unit_dictionary[['Ca']]
+    required_variables[[cc_column_name]]             <- unit_dictionary[['Cc']]
+    required_variables[[ci_column_name]]             <- unit_dictionary[['Ci']]
+    required_variables[[j_norm_column_name]]         <- unit_dictionary[['J_norm']]
+    required_variables[[kc_column_name]]             <- unit_dictionary[['Kc']]
+    required_variables[[ko_column_name]]             <- unit_dictionary[['Ko']]
+    required_variables[[rd_norm_column_name]]        <- unit_dictionary[['Rd_norm']]
+    required_variables[[total_pressure_column_name]] <- unit_dictionary[['total_pressure']]
+    required_variables[[vcmax_norm_column_name]]     <- unit_dictionary[['Vcmax_norm']]
 
     check_required_variables(exdf_obj, required_variables)
 
     # Extract key variables to make the following equations simpler
-    Ca         <- exdf_obj[, ca_column_name]    # micromol / mol
-    Cc         <- exdf_obj[, cc_column_name]    # micromol / mol
-    Ci         <- exdf_obj[, ci_column_name]    # micromol / mol
-    Gamma_star <- exdf_obj[, 'Gamma_star']      # micromol / mol
-    J          <- exdf_obj[, j_column_name]     # micromol / m^2 / s
-    Rd         <- exdf_obj[, rd_column_name]    # micromol / m^2 / s
-    TPU        <- exdf_obj[, tpu_column_name]   # micromol / m^2 / s
-    Vcmax      <- exdf_obj[, vcmax_column_name] # micromol / m^2 / s
+    Ca <- exdf_obj[, ca_column_name]    # micromol / mol
+    Cc <- exdf_obj[, cc_column_name]    # micromol / mol
+    Ci <- exdf_obj[, ci_column_name]    # micromol / mol
 
     # If gsc is as measured and gmc is infinite, we have the measured drawdown
     # across the stomata, but no drawdown across the mesophyll:
@@ -75,31 +70,29 @@ calculate_c3_limitations_warren <- function(
     # Make a helping function for calculating assimilation rates for different
     # Cc column names
     an_from_cc <- function(cc_name) {
-        sapply(seq_len(nrow(exdf_obj)), function(i) {
-            calculate_c3_assimilation(
-                exdf_obj[i, , TRUE],
-                J[i],
-                Gamma_star[i],
-                Rd[i],
-                TPU[i],
-                Vcmax[i],
-                POc = POc,
-                atp_use = atp_use,
-                nadph_use = nadph_use,
-                alpha = alpha,
-                curvature_cj = curvature_cj,
-                curvature_cjp = curvature_cjp,
-                cc_column_name = cc_name,
-                j_norm_column_name = j_norm_column_name,
-                kc_column_name = kc_column_name,
-                ko_column_name = ko_column_name,
-                rd_norm_column_name = rd_norm_column_name,
-                total_pressure_column_name = total_pressure_column_name,
-                vcmax_norm_column_name = vcmax_norm_column_name,
-                perform_checks = TRUE,
-                return_exdf = FALSE
-            )[['An']]
-        })
+        calculate_c3_assimilation(
+            exdf_obj,
+            '', # alpha_g
+            '', # Gamma_star
+            '', # J_at_25
+            '', # Rd_at_25
+            '', # Tp
+            '', # Vcmax_at_25
+            POc = POc,
+            atp_use = atp_use,
+            nadph_use = nadph_use,
+            curvature_cj = curvature_cj,
+            curvature_cjp = curvature_cjp,
+            cc_column_name = cc_name,
+            j_norm_column_name = j_norm_column_name,
+            kc_column_name = kc_column_name,
+            ko_column_name = ko_column_name,
+            rd_norm_column_name = rd_norm_column_name,
+            total_pressure_column_name = total_pressure_column_name,
+            vcmax_norm_column_name = vcmax_norm_column_name,
+            perform_checks = FALSE,
+            return_exdf = TRUE
+        )[, 'An']
     }
 
     # Calculate the net assimilation rate assuming gmc and gsc are as measured

R/calculate_c3_variable_j.R

@@ -59,20 +59,24 @@
 
     Rd_tl <- Rd_at_25 * exdf_obj[, rd_norm_column_name]         # micromol / m^2 / s
 
-    # Make sure key inputs have reasonable values; these checks cannot be
-    # bypassed
+    # Make sure key inputs have reasonable values
     msg <- character()
 
-    if (any(Ci < 0))         {msg <- append(msg, 'Ci must be >= 0')}
-    if (any(Gamma_star < 0)) {msg <- append(msg, 'Gamma_star must be >= 0')}
-    if (any(PhiPS2 < 0))     {msg <- append(msg, 'PhiPS2 must be >= 0')}
-    if (any(pressure < 0))   {msg <- append(msg, 'pressure must be >= 0')}
-    if (any(Qin < 0))        {msg <- append(msg, 'Qin must be >= 0')}
-    if (any(Rd_at_25 < 0))   {msg <- append(msg, 'Rd_at_25 must be >= 0')}
-    if (any(tau < 0))        {msg <- append(msg, 'tau must be >= 0')}
+    if (any(Ci < 0, na.rm = TRUE))         {msg <- append(msg, 'Ci must be >= 0')}
+    if (any(Gamma_star < 0, na.rm = TRUE)) {msg <- append(msg, 'Gamma_star must be >= 0')}
+    if (any(PhiPS2 < 0, na.rm = TRUE))     {msg <- append(msg, 'PhiPS2 must be >= 0')}
+    if (any(pressure < 0, na.rm = TRUE))   {msg <- append(msg, 'pressure must be >= 0')}
+    if (any(Qin < 0, na.rm = TRUE))        {msg <- append(msg, 'Qin must be >= 0')}
+    if (any(Rd_at_25 < 0, na.rm = TRUE))   {msg <- append(msg, 'Rd_at_25 must be >= 0')}
+    if (any(tau < 0, na.rm = TRUE))        {msg <- append(msg, 'tau must be >= 0')}
 
-    if (length(msg) > 0) {
-        stop(paste(msg, collapse = '. '))
+    msg <- paste(msg, collapse = '. ')
+
+    # We only bypass these checks if !perform_checks && return_exdf
+    if (perform_checks || !return_exdf) {
+        if (msg != '') {
+            stop(msg)
+        }
     }
 
     # Calculate J_F (actual RuBP regeneration rate as estimated from
@@ -101,17 +105,19 @@
             Rd_tl = Rd_tl,
             J_F = J_F,
             gmc = gmc,
-            Cc = Cc
+            Cc = Cc,
+            c3_variable_j_msg = msg
         ))
 
         document_variables(
             output,
-            c('calculate_c3_variable_j', 'Gamma_star', 'micromol mol^(-1)'),
-            c('calculate_c3_variable_j', 'tau',        'dimensionless'),
-            c('calculate_c3_variable_j', 'Rd_tl',      'micromol m^(-2) s^(-1)'),
-            c('calculate_c3_variable_j', 'J_F',        'micromol m^(-2) s^(-1)'),
-            c('calculate_c3_variable_j', 'gmc',        'mol m^(-2) s^(-1) bar^(-1)'),
-            c('calculate_c3_variable_j', 'Cc',         'micromol mol^(-1)')
+            c('calculate_c3_variable_j', 'Gamma_star',        'micromol mol^(-1)'),
+            c('calculate_c3_variable_j', 'tau',               'dimensionless'),
+            c('calculate_c3_variable_j', 'Rd_tl',             'micromol m^(-2) s^(-1)'),
+            c('calculate_c3_variable_j', 'J_F',               'micromol m^(-2) s^(-1)'),
+            c('calculate_c3_variable_j', 'gmc',               'mol m^(-2) s^(-1) bar^(-1)'),
+            c('calculate_c3_variable_j', 'Cc',                'micromol mol^(-1)'),
+            c('calculate_c3_variable_j', 'c3_variable_j_msg', '')
         )
     } else {
         return(list(gmc = gmc, Cc = Cc))