R/1TM_Qsun_regression_functions.R

In paulnorthrop/walls: Statistical Inference for the Thermophysical Properties of Walls from In-Situ Surveys

#' @export
one_tm_sun_reg_fn <- function(x1, x2, x3, x4, x5, Qp1, R1, R2, R3, C1, in_out,
                              tau) {
  # Total number of Q values (Qin and Qout)
  n <- length(Qp1)
  m <- n / 2
  # Indices for inside and outside
#  ind_i <- which(in_out == "inside")
#  ind_o <- which(in_out == "outside")
  ind_i <- 1:m
  ind_o <- (m + 1): n
  # Quantities involved in the regression function
  r1 <- 1 / R1
  r23 <- 1 / (R2 + R3)
  r3 <- 1 / R3
  C1tau <- C1 / tau[1]
  dsplus <- 2 * C1tau + r1 + r23
  dsminus <- 2 * C1tau - r1 - r23
  bs <- dsminus / dsplus
  # Extract covariates for inside and for outside
  Qinp1 <- Qp1[ind_i]
  Qoutp1 <- Qp1[ind_o]
  x1i <- x1[ind_i]
  x2i <- x2[ind_i]
  x3i <- x3[ind_i]
  x4i <- x4[ind_i]
  x5i <- x5[ind_i]
  x1o <- x1[ind_o]
  x2o <- x2[ind_o]
  x3o <- x3[ind_o]
  x4o <- x4[ind_o]
  x5o <- x5[ind_o]
  # Qin regression function
  Qin_reg <- 2 * r1 *(r23 * x1i + C1tau * x2i - r23 * x4i / r3) / dsplus
  # Qout regression function
  Qout_reg <- 2 * r23 * (r1 * x1o - C1tau * x3o - r1 * x4o / r3
                         - C1tau * x5o / r3) / dsplus
  return(c(Qin_reg, Qout_reg))
}

#' @export
one_tm_sun_grad_fn <- function(x1, x2, x3, x4, x5, Qp1, R1, R2, R3, C1, in_out,
                               tau) {
  # R1, R2, R3, C1 and C2 are vectors when gnls() calls this, but not when
  # nls() calls.
  # Use R1[1], R2[1], R3[1], C1[1]
  R1 <- R1[1]
  R2 <- R2[1]
  R3 <- R3[1]
  C1 <- C1[1]
  func <- function(x) {
    R1 <- x[1]
    R2 <- x[2]
    R3 <- x[3]
    C1 <- x[4]
    val <- one_tm_sun_reg_fn(x1, x2, x3, x4, x5, Qp1, R1, R2, R3, C1, in_out,
                             tau)
    return(val)
  }
  reg_fn <- one_tm_sun_reg_fn(x1, x2, x3, x4, x5, Qp1, R1, R2, R3, C1, in_out,
                              tau)
  pars <- c(R1, R2, R3, C1)
  for_jacobian <- list(func = func, x = pars)
  grad_mat <- do.call(numDeriv::jacobian, for_jacobian)
  dimnames(grad_mat) <- list(NULL, c("R1", "R2", "R3", "C1"))
  attr(reg_fn, "gradient") <- grad_mat
  return(reg_fn)
}

# Regression function using the parameterization
# (R1 + R2 + R3, R1 / (R1 + R2 + R3), R2 / (R1 + R2 + R3), C1)

#' @export
R_one_tm_sun_reg_fn <- function(x1, x2, x3, x4, x5, Qp1, sumR, fracR1, fracR2,
                                C1, in_out, tau) {
  # Transform back to the parameterization (R1, R2, R3, C1)
  R1 <- fracR1 * sumR
  R2 <- fracR2 * sumR
  R3 <- sumR - R1 - R2
  return(one_tm_sun_reg_fn(x1, x2, x3, x4, x5, Qp1, R1, R2, R3, C1, in_out,
                           tau))
}

#' @export
R_one_tm_sun_grad_fn <- function(x1, x2, x3, x4, x5, Qp1, sumR, fracR1, fracR2,
                                 C1, in_out, tau) {
  # sumR, fracR1, fracR2, C1 and C2 are vectors when gnls() calls this,
  # but not when nls() calls.
  # Use sumR[1], fracR1[1], fracR2[1], C1[1], C2[1]
  sumR <- sumR[1]
  fracR1 <- fracR1[1]
  fracR2 <- fracR2[1]
  C1 <- C1[1]
  func <- function(x) {
    sumR <- x[1]
    fracR1 <- x[2]
    fracR2 <- x[3]
    C1 <- x[4]
    val <- R_one_tm_sun_reg_fn(x1, x2, x3, x4, x5, Qp1, sumR, fracR1, fracR2,
                               C1, in_out, tau)
    return(val)
  }
  reg_fn <- R_one_tm_sun_reg_fn(x1, x2, x3, x4, x5, Qp1, sumR, fracR1, fracR2,
                                C1, in_out, tau)
  pars <- c(sumR, fracR1, fracR2, C1)
  for_jacobian <- list(func = func, x = pars)
  grad_mat <- do.call(numDeriv::jacobian, for_jacobian)
  dimnames(grad_mat) <- list(NULL, c("sumR", "fracR1", "fracR2", "C1"))
  attr(reg_fn, "gradient") <- grad_mat
  return(reg_fn)
}

# Regression function using the parameterization
# (U = 1 / (R1 + R2 + R3), R1 / (R1 + R2 + R3), R2 / (R1 + R2 + R3), C1, C2)

#' @export
U_one_tm_sun_reg_fn <- function(x1, x2, x3, x4, x5, Qp1, Uvalue, fracR1,
                                fracR2, C1, in_out, tau) {
  # Transform back to the parameterization (R1, R2, R3, C1)
  sumR <- 1 / Uvalue
  R1 <- fracR1 * sumR
  R2 <- fracR2 * sumR
  R3 <- sumR - R1 - R2
  return(one_tm_sun_reg_fn(x1, x2, x3, x4, x5, Qp1, R1, R2, R3, C1, in_out,
                           tau))
}

#' @export
U_one_tm_sun_grad_fn <- function(x1, x2, x3, x4, x5, Qp1, Uvalue, fracR1,
                                 fracR2, C1, in_out, tau) {
  # sumR, fracR1, fracR2, C1 and C2 are vectors when gnls() calls this,
  # but not when nls() calls.
  # Use sumR[1], fracR1[1], fracR2[1], C1[1], C2[1]
  Uvalue <- Uvalue[1]
  fracR1 <- fracR1[1]
  fracR2 <- fracR2[1]
  C1 <- C1[1]
  func <- function(x) {
    Uvalue <- x[1]
    fracR1 <- x[2]
    fracR2 <- x[3]
    C1 <- x[4]
    val <- U_one_tm_sun_reg_fn(x1, x2, x3, x4, x5, Qp1, Uvalue, fracR1, fracR2,
                               C1, in_out, tau)
    return(val)
  }
  reg_fn <- U_one_tm_sun_reg_fn(x1, x2, x3, x4, x5, Qp1, Uvalue, fracR1,
                                fracR2, C1, in_out, tau)
  pars <- c(Uvalue, fracR1, fracR2, C1)
  for_jacobian <- list(func = func, x = pars)
  grad_mat <- do.call(numDeriv::jacobian, for_jacobian)
  dimnames(grad_mat) <- list(NULL, c("Uvalue", "fracR1", "fracR2", "C1"))
  attr(reg_fn, "gradient") <- grad_mat
  return(reg_fn)
}