R/gradients.R

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

# Simplify, so that functions call both_reg_fn ?

#' @export
both_reg_fn <- function(x1, x23, Qsp1, R1, R2, C1, in_out, tau) {
  # R1, R2 and C1 are vectors when gnls() calls this, but not when nls() calls
  # Use R1[1], R2[1], C1[1]
  R1 <- R1[1]
  R2 <- R2[1]
  C1 <- C1[1]
  #
  r1 <- 1 / R1
  r2 <- 1 / R2
  dplus <- 2 * C1 / tau[1] + r1 + r2
  dminus <- 2 * C1 / tau[1] - r1 - r2
  b <- dminus / dplus
  term1 <- 2 * r1 * r2 * x1
  term2a <- 2 * r1 * C1 * x23[in_out == "inside"] / tau[1]
  term2b <- 2 * r2 * C1 * x23[in_out == "outside"] / tau[1]
  term2 <- c(term2a, term2b)
  term3 <- b * Qsp1
  reg_fn <- (term1 + term2) / dplus + term3
  #
  x2 <- x23[in_out == "inside"]
  x3 <- x23[in_out == "outside"]
  C1t <- C1 / tau[1]
  den <- (2 * C1 / tau[1] + r1 + r2) ^ 2
  dterm1dr1 <- 2 * r2 * (2 * C1t + r2) * x1 / den
  dterm1dr2 <- 2 * r1 * (2 * C1t + r1) * x1 / den
  dterm1dC1 <- -4 * r1 * r2 / tau[1] * x1 / den
  #
  dterm2adr1 <- 2 * C1t * (2 * C1t + r2) * x2 / den
  dterm2bdr1 <- -2 * r2 * C1t * x3 / den
  dterm2dr1 <- c(dterm2adr1, dterm2bdr1)
  dterm2adr2 <- -2 * r1 * C1t * x2 / den
  dterm2bdr2 <- 2 * C1t * (2 * C1t + r1) * x3 / den
  dterm2dr2 <- c(dterm2adr2, dterm2bdr2)
  dterm2adC1 <- 2 * r1 * (r1 + r2) * x2 / tau[1] / den
  dterm2bdC1 <- 2 * r2 * (r1 + r2) * x3 / tau[1] / den
  dterm2dC1 <- c(dterm2adC1, dterm2bdC1)
  #
  dterm3dr1 <- -4 * C1t * Qsp1 / den
  dterm3dr2 <- -4 * C1t * Qsp1 / den
  dterm3dC1 <- 4 * (r1 + r2) * Qsp1 / tau[1] / den
  #
  gradr1 <- dterm1dr1 + dterm2dr1 + dterm3dr1
  gradr2 <- dterm1dr2 + dterm2dr2 + dterm3dr2
  gradC1 <- dterm1dC1 + dterm2dC1 + dterm3dC1
  #
#  grad_mat <- cbind(-gradr1 / R1 ^ 2, -gradr2 / R2 ^ 2, gradC1)
#  dimnames(grad_mat) <- list(NULL, c("R1", "R2", "C1"))
  grad_mat <- array(0, c(length(x1), 3), list(NULL, c("R1", "R2", "C1")))
  grad_mat[, "R1"] <- -gradr1 / R1 ^ 2
  grad_mat[, "R2"] <- -gradr2 / R2 ^ 2
  grad_mat[, "C1"] <- gradC1
  attr(reg_fn, "gradient") <- grad_mat
  return(reg_fn)
}

#' @export
sum_both_reg_fn <- function(x1, x23, Qsp1, sum_R, frac_R, C1, in_out, tau) {
  sum_R <- sum_R[1]
  frac_R <- frac_R[1]
  C1 <- C1[1]
  R1 <- sum_R * frac_R
  R2 <- sum_R * (1 - frac_R)
  #
  r1 <- 1 / R1
  r2 <- 1 / R2
  dplus <- 2 * C1 / tau[1] + r1 + r2
  dminus <- 2 * C1 / tau[1] - r1 - r2
  b <- dminus / dplus
  term1 <- 2 * r1 * r2 * x1
  term2a <- 2 * r1 * C1 * x23[in_out == "inside"] / tau[1]
  term2b <- 2 * r2 * C1 * x23[in_out == "outside"] / tau[1]
  term2 <- c(term2a, term2b)
  term3 <- b * Qsp1
  reg_fn <- (term1 + term2) / dplus + term3
  #
  x2 <- x23[in_out == "inside"]
  x3 <- x23[in_out == "outside"]
  C1t <- C1 / tau[1]
  den <- (2 * C1 / tau[1] + r1 + r2) ^ 2
  dterm1dr1 <- 2 * r2 * (2 * C1t + r2) * x1 / den
  dterm1dr2 <- 2 * r1 * (2 * C1t + r1) * x1 / den
  dterm1dC1 <- -4 * r1 * r2 / tau[1] * x1 / den
  #
  dterm2adr1 <- 2 * C1t * (2 * C1t + r2) * x2 / den
  dterm2bdr1 <- -2 * r2 * C1t * x3 / den
  dterm2dr1 <- c(dterm2adr1, dterm2bdr1)
  dterm2adr2 <- -2 * r1 * C1t * x2 / den
  dterm2bdr2 <- 2 * C1t * (2 * C1t + r1) * x3 / den
  dterm2dr2 <- c(dterm2adr2, dterm2bdr2)
  dterm2adC1 <- 2 * r1 * (r1 + r2) * x2 / tau[1] / den
  dterm2bdC1 <- 2 * r2 * (r1 + r2) * x3 / tau[1] / den
  dterm2dC1 <- c(dterm2adC1, dterm2bdC1)
  #
  dterm3dr1 <- -4 * C1t * Qsp1 / den
  dterm3dr2 <- -4 * C1t * Qsp1 / den
  dterm3dC1 <- 4 * (r1 + r2) * Qsp1 / tau[1] / den
  gradr1 <- dterm1dr1 + dterm2dr1 + dterm3dr1
  gradr2 <- dterm1dr2 + dterm2dr2 + dterm3dr2
  gradC1 <- dterm1dC1 + dterm2dC1 + dterm3dC1
  #
#  grad_mat <- cbind(-gradr1 / R1 ^ 2, -gradr2 / R2 ^ 2, gradC1)
#  dimnames(grad_mat) <- list(NULL, c("R1", "R2", "C1"))
  grad_mat <- array(0, c(length(x1), 3), list(NULL, c("sum_R", "frac_R", "C1")))
  detadR1 <- -gradr1 / R1 ^ 2
  detadR2 <- -gradr2 / R2 ^ 2
#  grad_mat[, "sum_R"] <- detadR1 + detadR2  WRONG!!!!
  grad_mat[, "sum_R"] <- frac_R * detadR1 + (1 - frac_R) * detadR2
  grad_mat[, "frac_R"] <- sum_R * (detadR1 - detadR2)
  grad_mat[, "C1"] <- gradC1
  attr(reg_fn, "gradient") <- grad_mat
  return(reg_fn)
}

#' @export
log_sum_both_reg_fn <- function(x1, x23, Qsp1, log_sum_R, log_frac_R, log_C1,
                                in_out, tau) {
  log_sum_R <- log_sum_R[1]
  log_frac_R <- log_frac_R[1]
  log_C1 <- log_C1[1]
  sum_R <- exp(log_sum_R)
  frac_R <- exp(log_frac_R)
  C1 <- exp(log_C1)
  R1 <- sum_R * frac_R
  R2 <- sum_R * (1 - frac_R)
  #
  r1 <- 1 / R1
  r2 <- 1 / R2
  dplus <- 2 * C1 / tau[1] + r1 + r2
  dminus <- 2 * C1 / tau[1] - r1 - r2
  b <- dminus / dplus
  term1 <- 2 * r1 * r2 * x1
  term2a <- 2 * r1 * C1 * x23[in_out == "inside"] / tau[1]
  term2b <- 2 * r2 * C1 * x23[in_out == "outside"] / tau[1]
  term2 <- c(term2a, term2b)
  term3 <- b * Qsp1
  reg_fn <- (term1 + term2) / dplus + term3
  #
  x2 <- x23[in_out == "inside"]
  x3 <- x23[in_out == "outside"]
  C1t <- C1 / tau[1]
  den <- (2 * C1 / tau[1] + r1 + r2) ^ 2
  dterm1dr1 <- 2 * r2 * (2 * C1t + r2) * x1 / den
  dterm1dr2 <- 2 * r1 * (2 * C1t + r1) * x1 / den
  dterm1dC1 <- -4 * r1 * r2 / tau[1] * x1 / den
  #
  dterm2adr1 <- 2 * C1t * (2 * C1t + r2) * x2 / den
  dterm2bdr1 <- -2 * r2 * C1t * x3 / den
  dterm2dr1 <- c(dterm2adr1, dterm2bdr1)
  dterm2adr2 <- -2 * r1 * C1t * x2 / den
  dterm2bdr2 <- 2 * C1t * (2 * C1t + r1) * x3 / den
  dterm2dr2 <- c(dterm2adr2, dterm2bdr2)
  dterm2adC1 <- 2 * r1 * (r1 + r2) * x2 / tau[1] / den
  dterm2bdC1 <- 2 * r2 * (r1 + r2) * x3 / tau[1] / den
  dterm2dC1 <- c(dterm2adC1, dterm2bdC1)
  #
  dterm3dr1 <- -4 * C1t * Qsp1 / den
  dterm3dr2 <- -4 * C1t * Qsp1 / den
  dterm3dC1 <- 4 * (r1 + r2) * Qsp1 / tau[1] / den
  gradr1 <- dterm1dr1 + dterm2dr1 + dterm3dr1
  gradr2 <- dterm1dr2 + dterm2dr2 + dterm3dr2
  gradC1 <- dterm1dC1 + dterm2dC1 + dterm3dC1
  #
  #  grad_mat <- cbind(-gradr1 / R1 ^ 2, -gradr2 / R2 ^ 2, gradC1)
  #  dimnames(grad_mat) <- list(NULL, c("R1", "R2", "C1"))
  grad_mat <- array(0, c(length(x1), 3), list(NULL, c("log_sum_R", "log_frac_R",
                                                      "log_C1")))
  detadR1 <- -gradr1 / R1 ^ 2
  detadR2 <- -gradr2 / R2 ^ 2
#  grad_mat[, "log_sum_R"] <- sum_R * (detadR1 + detadR2) WRONG !!!
  grad_mat[, "log_sum_R"] <- sum_R * (frac_R * detadR1 + (1 - frac_R) * detadR2)
  grad_mat[, "log_frac_R"] <- frac_R * sum_R * (detadR1 - detadR2)
  grad_mat[, "log_C1"] <- C1 * gradC1
  attr(reg_fn, "gradient") <- grad_mat
  return(reg_fn)
}

#' @export
log_both_reg_fn <- function(x1, x23, Qsp1, logR1, logR2, logC1, in_out, tau) {
  R1 <- exp(logR1)[1]
  R2 <- exp(logR2)[1]
  C1 <- exp(logC1)[1]
  #
  r1 <- 1 / R1
  r2 <- 1 / R2
  dplus <- 2 * C1 / tau[1] + r1 + r2
  dminus <- 2 * C1 / tau[1] - r1 - r2
  b <- dminus / dplus
  term1 <- 2 * r1 * r2 * x1
  term2a <- 2 * r1 * C1 * x23[in_out == "inside"] / tau[1]
  term2b <- 2 * r2 * C1 * x23[in_out == "outside"] / tau[1]
  term2 <- c(term2a, term2b)
  term3 <- b * Qsp1
  reg_fn <- (term1 + term2) / dplus + term3
  #
  x2 <- x23[in_out == "inside"]
  x3 <- x23[in_out == "outside"]
  C1t <- C1 / tau[1]
  den <- (2 * C1 / tau[1] + r1 + r2) ^ 2
  dterm1dr1 <- 2 * r2 * (2 * C1t + r2) * x1 / den
  dterm1dr2 <- 2 * r1 * (2 * C1t + r1) * x1 / den
  dterm1dC1 <- -4 * r1 * r2 / tau[1] * x1 / den
  #
  dterm2adr1 <- 2 * C1t * (2 * C1t + r2) * x2 / den
  dterm2bdr1 <- -2 * r2 * C1t * x3 / den
  dterm2dr1 <- c(dterm2adr1, dterm2bdr1)
  dterm2adr2 <- -2 * r1 * C1t * x2 / den
  dterm2bdr2 <- 2 * C1t * (2 * C1t + r1) * x3 / den
  dterm2dr2 <- c(dterm2adr2, dterm2bdr2)
  dterm2adC1 <- 2 * r1 * (r1 + r2) * x2 / tau[1] / den
  dterm2bdC1 <- 2 * r2 * (r1 + r2) * x3 / tau[1] / den
  dterm2dC1 <- c(dterm2adC1, dterm2bdC1)
  #
  dterm3dr1 <- -4 * C1t * Qsp1 / den
  dterm3dr2 <- -4 * C1t * Qsp1 / den
  dterm3dC1 <- 4 * (r1 + r2) * Qsp1 / tau[1] / den
  #
  gradr1 <- dterm1dr1 + dterm2dr1 + dterm3dr1
  gradr2 <- dterm1dr2 + dterm2dr2 + dterm3dr2
  gradC1 <- dterm1dC1 + dterm2dC1 + dterm3dC1
  #
  grad_mat <- cbind(-gradr1 / R1, -gradr2 / R2, C1 * gradC1)
  dimnames(grad_mat) <- list(NULL, c("logR1", "logR2", "logC1"))
  attr(reg_fn, "gradient") <- grad_mat
  return(reg_fn)
}

#' @export
trans_sum_both_reg_fn <- function(x1, x23, Qsp1, log_sum_R, logit_frac_R,
                                  log_C1, in_out, tau) {
  log_sum_R <- log_sum_R[1]
  logit_frac_R <- logit_frac_R[1]
  log_C1 <- log_C1[1]
  sum_R <- exp(log_sum_R)
  frac_R <- exp(logit_frac_R) / (1 + exp(logit_frac_R))
  C1 <- exp(log_C1)
  R1 <- sum_R * frac_R
  R2 <- sum_R * (1 - frac_R)
  #
  r1 <- 1 / R1
  r2 <- 1 / R2
  dplus <- 2 * C1 / tau[1] + r1 + r2
  dminus <- 2 * C1 / tau[1] - r1 - r2
  b <- dminus / dplus
  term1 <- 2 * r1 * r2 * x1
  term2a <- 2 * r1 * C1 * x23[in_out == "inside"] / tau[1]
  term2b <- 2 * r2 * C1 * x23[in_out == "outside"] / tau[1]
  term2 <- c(term2a, term2b)
  term3 <- b * Qsp1
  reg_fn <- (term1 + term2) / dplus + term3
  #
  x2 <- x23[in_out == "inside"]
  x3 <- x23[in_out == "outside"]
  C1t <- C1 / tau[1]
  den <- (2 * C1 / tau[1] + r1 + r2) ^ 2
  dterm1dr1 <- 2 * r2 * (2 * C1t + r2) * x1 / den
  dterm1dr2 <- 2 * r1 * (2 * C1t + r1) * x1 / den
  dterm1dC1 <- -4 * r1 * r2 / tau[1] * x1 / den
  #
  dterm2adr1 <- 2 * C1t * (2 * C1t + r2) * x2 / den
  dterm2bdr1 <- -2 * r2 * C1t * x3 / den
  dterm2dr1 <- c(dterm2adr1, dterm2bdr1)
  dterm2adr2 <- -2 * r1 * C1t * x2 / den
  dterm2bdr2 <- 2 * C1t * (2 * C1t + r1) * x3 / den
  dterm2dr2 <- c(dterm2adr2, dterm2bdr2)
  dterm2adC1 <- 2 * r1 * (r1 + r2) * x2 / tau[1] / den
  dterm2bdC1 <- 2 * r2 * (r1 + r2) * x3 / tau[1] / den
  dterm2dC1 <- c(dterm2adC1, dterm2bdC1)
  #
  dterm3dr1 <- -4 * C1t * Qsp1 / den
  dterm3dr2 <- -4 * C1t * Qsp1 / den
  dterm3dC1 <- 4 * (r1 + r2) * Qsp1 / tau[1] / den
  gradr1 <- dterm1dr1 + dterm2dr1 + dterm3dr1
  gradr2 <- dterm1dr2 + dterm2dr2 + dterm3dr2
  gradC1 <- dterm1dC1 + dterm2dC1 + dterm3dC1
  #
  #  grad_mat <- cbind(-gradr1 / R1 ^ 2, -gradr2 / R2 ^ 2, gradC1)
  #  dimnames(grad_mat) <- list(NULL, c("R1", "R2", "C1"))
  grad_mat <- array(0, c(length(x1), 3), list(NULL, c("log_sum_R",
                                                      "logit_frac_R",
                                                      "log_C1")))
  detadR1 <- -gradr1 / R1 ^ 2
  detadR2 <- -gradr2 / R2 ^ 2
#  detadphi1 <- frac_R * detadR1 + (1 - frac_R) * detadR2
  grad_mat[, "log_sum_R"] <- sum_R * (frac_R * detadR1 + (1 - frac_R) * detadR2)
#  detadphi2 <- sum_R * (detadR1 - detadR2)
  grad_mat[, "logit_frac_R"] <- frac_R * (1 - frac_R) * sum_R * (detadR1 - detadR2)
  grad_mat[, "log_C1"] <- C1 * gradC1
  attr(reg_fn, "gradient") <- grad_mat
  return(reg_fn)
}