################################################################################
# All functions for plateflow
################################################################################


#' Levenberg Marquard fit of 4 pl
#'
#' Returns list of following: summary of restricted and unrestricted model fits of the 4pl function,
#' potency estimates and respective CIs of restricted and unrestricted models, and the predictions thereof.
#'
#' @param ro_new The dilution series readouts for REF and TEST + 1 column log(loncentration) or concentration.
#' @param TransFlag A boolean if the plot should be with natural log as readout (TRUE) or not (FALSE).
#' @returns Returns list of following: summary of restricted and unrestricted model fits of the 4pl function,
#' potency estimates and respective CIs of restricted and unrestricted models, and the predictions thereof.
#' @export
#' @examples
#' dat <- data.frame(
#'   REF1 = c(1547, 1620, 1644, 2504, 3426, 3512, 3401, 3787), REF2 = c(1492, 1536, 1384, 2286, 3046, 3479, 3516, 3497),
#'   REF3 = c(1468, 1827, 1558, 2252, 3002, 3349, 2945, 3665),
#'   TEST1 = c(1405, 1523, 1502, 1474, 2383, 3221, 3589, 3445), TEST2 = c(1420, 1516, 1544, 1512, 2226, 3219, 3327, 3591),
#'   TEST3 = c(1399, 1376, 1588, 1475, 2148, 3083, 2942, 3466), log_dose = c(5.01, 3.401, 2.708, 2.015, 1.32176, 0.62861, -0.0645385, -1.6739764)
#' )
#' TransF <- FALSE
#' Dat <- list()
#' te <- Fitting_FUNC(dat, TransF)
#' print(te)
Fitting_FUNC <- function(ro_new, TransFlag = F) {
  CORro <- cor(ro_new[, 1], ro_new[, ncol(ro_new)])
  # browser()
  all_l <- melt(data.frame(ro_new), id.vars = "log_dose", variable.name = "replname", value.name = "readout")
  isRef <- rep(c(1, 0), 1, each = nrow(all_l) / 2)
  isSample <- rep(c(0, 1), 1, each = nrow(all_l) / 2)
  all_l2 <- cbind(all_l, isRef, isSample)
  # browser()
  if (CORro < 0) SLOPE <- -1 else SLOPE <- 1
  if (!TransFlag) {
    startlist <- list(a = min(ro_new[, 2]), b = SLOPE, d = max(ro_new[, 2]), cs = mean(all_l$log_dose), r = 0)

    mr <- tryCatch(
      {
        gsl_nls(
          fn = readout ~ a + (d - a) / (1 + exp(b * ((cs - r * isSample) - log_dose))),
          data = all_l2,
          start = startlist, # race=T,
          control = gsl_nls_control(xtol = 1e-6, ftol = 1e-6, gtol = 1e-6)
        )
      },
      warning = function(e) {
        mr <<- "In nlsModel singular gradient matrix"
      }
    )
    # Stop if singular gradient matrix
    if (is.character(mr)) {
      return(mr)
    }

    s_mr <- tryCatch(
      {
        s_mr <- summary(mr)
      },
      error = function(err) {
        s_mr <- NULL
      }
    )
  } else {
    startlist <- list(a = log(min(ro_new[, 2])), b = SLOPE, d = log(max(ro_new[, 2])), cs = mean(all_l$log_dose), r = 0)
    mrT <- gsl_nls(
      fn = log(readout) ~ a + (d - a) / (1 + exp(b * ((cs - r * isSample) - log_dose))),
      data = all_l2,
      start = startlist,
      control = gsl_nls_control(xtol = 1e-6, ftol = 1e-6, gtol = 1e-6)
    )
    s_mr <- summary(mrT)
  }

  if (!TransFlag) {
    startlistmu <- list(
      as = min(ro_new[, 2]), bs = SLOPE, ds = max(ro_new[, 2]), cs = mean(all_l$log_dose),
      at = min(ro_new[, 2]), bt = SLOPE, dt = max(ro_new[, 2]), r = 0
    )
    tryCatch(
      {
        mu <- gsl_nls(
          fn = readout ~ as * isRef + at * isSample + (ds * isRef + dt * isSample - as * isRef - at * isSample) /
            (1 + isRef * exp(bs * (cs - log_dose)) + isSample * exp(bt * (cs - r * isSample - log_dose))),
          data = all_l,
          start = startlistmu,
          control = gsl_nls_control(xtol = 1e-6, ftol = 1e-6, gtol = 1e-6)
        )
      },
      error = function(msg) {
        return(0)
      }
    )

    Sum_u <- tryCatch(
      {
        summary(mu)
      },
      error = function(msg) {
        return(0)
      }
    )
  } else {
    startlistmu <- list(
      as = log(min(ro_new[, 2])), bs = SLOPE, ds = log(max(ro_new[, 2])), cs = mean(all_l$log_dose),
      at = log(min(ro_new[, 2])), bt = SLOPE, dt = log(max(ro_new[, 2])), r = 0
    )
    tryCatch(
      {
        muT <- gsl_nls(
          fn = log(readout) ~ as * isRef + at * isSample + (ds * isRef + dt * isSample - as * isRef - at * isSample) /
            (1 + isRef * exp(bs * (cs - log_dose)) + isSample * exp(bt * (cs - r * isSample - log_dose))),
          data = all_l,
          start = startlistmu,
          control = gsl_nls_control(xtol = 1e-6, ftol = 1e-6, gtol = 1e-6)
        )
      },
      error = function(msg) {
        return(0)
      }
    )

    Sum_u <- tryCatch(
      {
        summary(muT)
      },
      error = function(msg) {
        return(0)
      }
    )
  }
  if (!TransFlag) {
    pot_est <- exp(confintd(mr, "r", method = "asymptotic"))
    potU_est <- exp(confintd(mu, "r", method = "asymptotic"))
    PRED <- predict(mr)
    PREDu <- predict(mu)
  } else {
    pot_est <- exp(confintd(mrT, "r", method = "asymptotic"))
    potU_est <- exp(confintd(muT, "r", method = "asymptotic"))
    PRED <- predict(mrT)
    PREDu <- predict(muT)
  }
  return(list(s_mr, Sum_u, pot_est, potU_est, PRED, PREDu))
}


#' Plot when the fitting has a singularity
#'
#' Returns the scatter plot of the data, with REF in blue and TEST in red.
#' One plots are returned.
#'
#' @param dat The dilution series readouts for REF and TEST + 1 column log(loncentration) or concentration.
#' @returns A grid object with 2 linearity plots, restricted and unrestricted model.
#' @export
#' @examples
#' data.frame(
#'   R_dil1 = c(
#'     10.0651024695491, 10.9844983291817, 10.7635586089293, 10.4597656321327, 10.3898668457823, 10.8171761349909,
#'     10.319758021908, 10.1304854046653
#'   ),
#'   R_dil2 = c(
#'     10.9649145494504, 10.0202868589385, 10.8424145955735, 10.9311360356894, 10.3284659026404,
#'     10.6890147558796, 10.3014450252305, 10.9594838595181
#'   ),
#'   R_dil3 = c(
#'     10.4630510824383, 10.4566715089363, 10.2350765290036, 10.3300581874798, 10.9648088137065,
#'     10.286893755805, 10.4856643841389, 10.5275521552307
#'   ),
#'   T_dil1 = c(
#'     12.732175566336, 12.7756403995095, 12.1672539684741, 12.7060603907892, 12.8000685682832,
#'     12.8800092157515, 12.7160581291873, 12.6996878912416
#'   ),
#'   T_dil2 = c(
#'     12.3923194313831, 12.0943488144175, 12.7955302154828, 12.4825917078735, 12.6856540203788,
#'     12.7348548498556, 12.9222470610476, 12.1186618671252
#'   ),
#'   T_dil3 = c(
#'     12.7899182255274, 12.9722600411128, 12.7078445380891, 12.4913523531941, 12.1718281909609,
#'     12.5313873615133, 12.952802332772, 12.5960321394342
#'   ),
#'   log_dose = c(
#'     0, -1.09861228866811, -2.19722457733622, -3.29583686600433, -4.39444915467244,
#'     -5.49306144334055, -6.59167373200866, -7.69028602067677
#'   )
#' )
#'
#'
#' p <- plotSingularity(dat)
#' print(p)
plotSingularity <- function(dat) { # sigmoid,det_sig,
  CORdat <- cor(dat[, 1], dat[, ncol(dat)])
  # browser()
  all_l <- melt(data.frame(dat), id.vars = "log_dose", variable.name = "replname", value.name = "readout")
  isRef <- rep(c(1, 0), 1, each = nrow(all_l) / 2)
  isSample <- rep(c(0, 1), 1, each = nrow(all_l) / 2)
  all_l2 <- cbind(all_l, isRef, isSample)
  # browser()

  log_dose <- unique(all_l$log_dose)
  seq_x <- seq(min(log_dose), max(log_dose), 0.1)


  # browser()
  # all_l2$readout[all_l2$readout < 0] <- 0.01
  all_l2$readouttrans <- log(all_l2$readout)

  # browser()
  pSing <- ggplot(all_l2, aes(x = log_dose, y = readout, color = factor(isRef))) +
    geom_point(shape = factor(isRef), size = 3, alpha = 0.8) +
    labs(
      title = paste("No 4pl fit possible"),
      color = "product"
    ) +
    scale_color_manual(labels = c("test", "reference"), values = c("#C2173F", "#4545BA")) +
    scale_shape_manual(labels = c("test", "reference")) +
    scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
    scale_y_continuous(breaks = scales::pretty_breaks(n = 10)) +
    # theme_bw() +
    theme(axis.text = element_text(size = 14))

  return(pSing)
}


#' Plot sigmoidal curve
#'
#' Returns the final plots of the 4pl function as sigmoidal lines, and the single readouts as scatter, with REF in blue and TEST in red.
#' Two plots are returned in a grid object: the unrestricted model and the restricted model.
#'
#' @param dat The dilution series readouts for REF and TEST + 1 column log(loncentration) or concentration.
#' @param sigmoid The parameters of the 4pl fit of unrestricted model from EXCEL input.
#' @param det_sig The parameters of the 4pl fit of unrestricted model from the metadata input.
#' @param TransFlag A boolean if the plot should be with natural log as readout (TRUE) or not (FALSE).
#' @returns A grid object either of the original scale or the natural log of the readouts.
#' @export
#' @examples
#' dat <- data.frame(
#'   REF1 = c(1547, 1620, 1644, 2504, 3426, 3512, 3401, 3787), REF2 = c(1492, 1536, 1384, 2286, 3046, 3479, 3516, 3497),
#'   REF3 = c(1468, 1827, 1558, 2252, 3002, 3349, 2945, 3665),
#'   TEST1 = c(1405, 1523, 1502, 1474, 2383, 3221, 3589, 3445), TEST2 = c(1420, 1516, 1544, 1512, 2226, 3219, 3327, 3591),
#'   TEST3 = c(1399, 1376, 1588, 1475, 2148, 3083, 2942, 3466), log_dose = c(5.01, 3.401, 2.708, 2.015, 1.32176, 0.62861, -0.0645385, -1.6739764)
#' )
#' sigmoid <- c(0.7163324, 0.5636804, 10.6156340, 9.9784160, -0.7504673, -0.7108692, -3.5788141, -0.6662962)
#' det_sig <- FALSE
#' TransF <- FALSE
#' Dat <- list()
#' p <- plot_f(dat, sigmoid, det_sig, TransF)
#' print(p)
plot_f <- function(dat, TransFlag = F) { # sigmoid,det_sig,
  CORdat <- cor(dat[, 1], dat[, ncol(dat)])
  # browser()
  all_l <- melt(data.frame(dat), id.vars = "log_dose", variable.name = "replname", value.name = "readout")
  isRef <- rep(c(1, 0), 1, each = nrow(all_l) / 2)
  isSample <- rep(c(0, 1), 1, each = nrow(all_l) / 2)
  all_l2 <- cbind(all_l, isRef, isSample)
  # browser()
  MODLS <- Fitting_FUNC(dat, TransFlag = F)
  s_mr <- MODLS[[1]]
  a <- s_mr$coefficients["a", 1]
  b <- s_mr$coefficients["b", 1]
  cs <- s_mr$coefficients["cs", 1]
  d <- s_mr$coefficients["d", 1]
  r <- s_mr$coefficients["r", 1]

  log_dose <- unique(all_l$log_dose)
  seq_x <- seq(min(log_dose), max(log_dose), 0.1)
  SAMPLE <- a + (d - a) / (1 + exp(b * ((cs - r) - seq_x)))
  REF <- a + (d - a) / (1 + exp(b * ((cs) - seq_x)))

  s_mu <- MODLS[[2]]

  as <- s_mu$coefficients["as", 1]
  bs <- s_mu$coefficients["bs", 1]
  cs <- s_mu$coefficients["cs", 1]
  ds <- s_mu$coefficients["ds", 1]
  at <- s_mu$coefficients["at", 1]
  bt <- s_mu$coefficients["bt", 1]
  ct <- s_mu$coefficients["cs", 1] - s_mu$coefficients["r", 1]
  dt <- s_mu$coefficients["dt", 1]
  r <- s_mu$coefficients["r", 1]


  SAMPLEtrue <- at + (dt - at) / (1 + exp(bt * ((cs - r - seq_x))))
  REFtrue <- as + (ds - as) / (1 + exp(bs * ((cs - seq_x))))

  # browser()
  pl_df <- cbind(seq_x, SAMPLE, REF, SAMPLEtrue, REFtrue)
  all_l2$readout[all_l2$readout < 0] <- 0.01
  all_l2$readouttrans <- log(all_l2$readout)
  slopeEC50 <- b * (d - a) / 4
  Intercept <- a + (d - a) / 2 - b * (d - a) / 4 * cs
  # browser()
  Xbendl3 <- cs - (1.5434 / b)
  Xbendu3 <- cs + (1.5434 / b)
  XbendlT <- cs - r - (1.5434 / b)
  XbenduT <- cs - r + (1.5434 / b)
  XasymplS <- cs - (3 / b)
  XasympuS <- cs + (3 / b)
  XasymplT <- cs - r - (3 / b)
  XasympuT <- cs - r + (3 / b)
  bendpoints <- c(
    bendREF_lower = round(Xbendl3, 3), bendREF_upper = round(Xbendu3, 3),
    bendSAMPLE_lower = round(XbendlT, 3), bendSAMPLE_upper = round(XbenduT, 3),
    asympREF_lower = round(XasymplS, 3), asympREF_upper = round(XasympuS, 3),
    asympSAMPLE_lower = round(XasymplT, 3), asympSAMPLE_upper = round(XasympuT, 3)
  )
  Dat$bendpoints <- bendpoints
  Dat$cfordils <- cs
  # browser()
  p <- ggplot(all_l2, aes(x = log_dose, y = readout, color = factor(isRef))) +
    geom_point(shape = factor(isRef), size = 3, alpha = 0.8) +
    labs(
      title = paste("restricted 4pl model"),
      color = "product"
    ) +
    scale_color_manual(labels = c("test", "reference"), values = c("#C2173F", "#4545BA")) +
    scale_shape_manual(labels = c("test", "reference")) +
    scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
    scale_y_continuous(breaks = scales::pretty_breaks(n = 10)) +
    # theme_bw() +
    theme(axis.text = element_text(size = 14))

  p2 <- p + geom_line(
    data = as.data.frame(pl_df), aes(x = seq_x, y = SAMPLE), color = "#C2173F",
    inherit.aes = F
  ) +
    geom_line(
      data = as.data.frame(pl_df), aes(x = seq_x, y = REF), color = "#4545BA",
      inherit.aes = F
    ) +
    geom_line(
      data = as.data.frame(pl_df), aes(x = seq_x, y = SAMPLEtrue), color = "#C2173F", linetype = 2, alpha = 0.4,
      inherit.aes = F
    ) +
    geom_line(
      data = as.data.frame(pl_df), aes(x = seq_x, y = REFtrue), color = "#4545BA", linetype = 2, alpha = 0.4,
      inherit.aes = F
    ) +
    geom_vline(xintercept = c(Xbendl3, Xbendu3), col = "#4545BA", linetype = 2) +
    geom_vline(xintercept = c(XbendlT, XbenduT), col = "#C2173F", linetype = 2) +
    geom_vline(xintercept = c(XasymplS, XasympuS), col = "#4545BABB", linetype = 3) +
    geom_vline(xintercept = c(XasymplT, XasympuT), col = "#C2173FBB", linetype = 3) +
    annotate("text", x = cs, y = a + (d - a) / 2, label = "0", size = 5) +
    geom_abline(slope = slopeEC50, intercept = Intercept) +
    theme(legend.position = "none")
  Dat$p2 <- p2

  # transformed plots
  p_rt <- ggplot(all_l2, aes(x = log_dose, y = readouttrans, color = factor(isRef))) +
    geom_point(shape = factor(isRef), size = 3, alpha = 0.8) +
    labs(title = paste("restricted transformed 4pl"), color = "product") +
    scale_color_manual(labels = c("test", "reference"), values = c("#C2173F", "#4545BA")) +
    theme_bw()

  MODLStrans <- Fitting_FUNC(dat, TransFlag = TRUE)
  s_mrt <- MODLStrans[[1]]
  a_trans <- s_mrt$coefficients["a", 1]
  b_trans <- s_mrt$coefficients["b", 1]
  cs_trans <- s_mrt$coefficients["cs", 1]
  d_trans <- s_mrt$coefficients["d", 1]
  r_trans <- s_mrt$coefficients["r", 1]

  XbendlTrans <- cs_trans - (1.5434 / b_trans)
  XbenduTrans <- cs_trans + (1.5434 / b_trans)
  XbendlTransT <- cs_trans - r_trans - (1.5434 / b_trans)
  XbenduTransT <- cs_trans - r_trans + (1.5434 / b_trans)
  bendpointsTRANS <- c(
    bendREF_lower = round(XbendlTrans, 3), bendREF_upper = round(XbenduTrans, 3),
    bendSAMPLE_lower = round(XbendlTransT, 3), bendSAMPLE_upper = round(XbenduTransT, 3)
  )
  Dat$bendpointsTRANS <- bendpointsTRANS
  SAMPLEtrans <- a_trans + (d_trans - a_trans) / (1 + exp(b_trans * ((cs_trans - r_trans) - seq_x)))
  REFtrans <- a_trans + (d_trans - a_trans) / (1 + exp(b_trans * ((cs_trans) - seq_x)))

  pl_df_trans <- cbind(seq_x, SAMPLEtrans, REFtrans)
  p_rt2 <- p_rt + geom_line(
    data = as.data.frame(pl_df_trans), aes(x = seq_x, y = SAMPLEtrans), color = "#C2173F",
    inherit.aes = F
  ) +
    geom_line(
      data = as.data.frame(pl_df_trans), aes(x = seq_x, y = REFtrans), color = "#4545BA",
      inherit.aes = F
    ) +
    geom_vline(xintercept = c(XbendlTrans, XbenduTrans), col = "#4545BA", linetype = 2) +
    geom_vline(xintercept = c(XbendlTransT, XbenduTransT), col = "#C2173F", linetype = 2) +
    theme(legend.position = "none", axis.text = element_text(size = 14))

  # browser()
  Sum_u <- MODLS[[2]]
  # browser()
  # if (is.null(det_sig)) {
  ast <- Sum_u$coefficients["as", 1]
  ate <- Sum_u$coefficients["at", 1]
  bst <- Sum_u$coefficients["bs", 1]
  bte <- Sum_u$coefficients["bt", 1]
  cst <- Sum_u$coefficients["cs", 1]
  cte <- Sum_u$coefficients["cs", 1] - Sum_u$coefficients["r", 1]
  dst <- Sum_u$coefficients["ds", 1]
  dte <- Sum_u$coefficients["dt", 1]
  # } else {
  #   ast <- det_sig[5]
  #   ate <- det_sig[6]
  #   bst <- det_sig[1]
  #   bte <- det_sig[2]
  #   cst <- det_sig[7]
  #   cte <- det_sig[8]
  #   dst <- det_sig[3]
  #   dte <- det_sig[4]
  # }
  REFu <- ast + (dst - ast) / (1 + exp(bst * (cst - seq_x)))
  SAMPLEu <- ate + (dte - ate) / (1 + exp(bte * (cte - seq_x)))
  pl_df2 <- cbind(seq_x, SAMPLEu, REFu)
  # browser()
  pu <- ggplot(all_l2, aes(x = log_dose, y = readout, color = factor(isRef))) +
    geom_point(size = 3) +
    labs(title = "unrestricted 4pl model", color = "product") +
    scale_color_manual(labels = c("test", "reference"), values = c("#C2173F88", "#4545BA88")) +
    theme_bw()
  pu2 <- pu + geom_line(
    data = as.data.frame(pl_df2), aes(x = seq_x, y = SAMPLEu),
    color = "#C2173F", inherit.aes = F
  ) +
    geom_line(
      data = as.data.frame(pl_df2), aes(x = seq_x, y = REFu),
      color = "#4545BA", inherit.aes = F,
      show.legend = F
    )
  pu2_ <- pu2 +
    theme(legend.position = "none", axis.text = element_text(size = 14))

  putrans <- ggplot(all_l2, aes(x = log_dose, y = readouttrans, color = factor(isRef))) +
    geom_point(size = 3) +
    labs(title = "unrestricted transformed 4_pl-Model", color = "product") +
    scale_color_manual(labels = c("test", "reference"), values = c("#C2173FAA", "#4545BAAA")) +
    theme_bw()

  Sum_ut <- MODLStrans[[2]]

  ast_t <- Sum_ut$coefficients["as", 1]
  ate_t <- Sum_ut$coefficients["at", 1]
  bst_t <- Sum_ut$coefficients["bs", 1]
  bte_t <- Sum_ut$coefficients["bt", 1]
  cst_t <- Sum_ut$coefficients["cs", 1]
  cte_t <- Sum_ut$coefficients["cs", 1] - Sum_ut$coefficients["r", 1]
  dst_t <- Sum_ut$coefficients["ds", 1]
  dte_t <- Sum_ut$coefficients["dt", 1]

  REFu_trans <- ast_t + (dst_t - ast_t) / (1 + exp(bst_t * (cst_t - seq_x)))
  SAMPLEu_trans <- ate_t + (dte_t - ate_t) / (1 + exp(bte_t * (cte_t - seq_x)))
  pl_df2u_t <- cbind(seq_x, SAMPLEu_trans, REFu_trans)

  pu2_t <- putrans + geom_line(
    data = as.data.frame(pl_df2u_t), aes(x = seq_x, y = SAMPLEu_trans),
    color = "#C2173F", inherit.aes = F
  ) +
    geom_line(
      data = as.data.frame(pl_df2u_t), aes(x = seq_x, y = REFu_trans),
      color = "#4545BA", inherit.aes = F,
      show.legend = F
    )
  pu3_t <- pu2_t
  if (TransFlag) grid.arrange(p_rt2, pu3_t, nrow = 1) else grid.arrange(p2, pu2_, nrow = 1)
}

#' Simulate a well-plate readout
#'
#' Returns a possible well-plate result based on the metadata provided and a variability.
#' Homo- and heteroscedasticity can be simulated.
#'
#' @param as, bs, cs, ds, at, bt,  dt,r,ct, gt, gs are the parameter of the 5pl logistic regression.
#' @param sd_fac The standard deviation with units in % of upper-lower asymptote.
#' @param log_conc The natural log of the concentrations.
#' @param heteroNoise Boolean if heteroscedstic noise should be added.
#' @param noDilSeries A number >1 indicating how many dilution series per product should be simulated.
#' @param noDils A number >7 indicating how many dilutions steps per product should be simulated.
#' @returns A data-frame with readouts and natural log of concentrations.
#' @export
#' @examples
#' as <- 3
#' bs <- 1
#' cs <- -4
#' ds <- 10
#' at <- 3
#' bt <- 1
#' dt <- 10
#' r <- 0.0001
#' ct <- cs - r
#' sd_fac <- 0.1
#' gt <- 1
#' gs <- 1
#' lnConc <- c(1, 0, -1, -2, -3, -4, -5, -6)
#' heteroNoise <- FALSE
#' noDilS <- 3
#' noD <- 8
#' Calc_DilRes(as, bs, cs, ds, at, bt, dt, r, ct, sd_fac, gt, gs, log_conc = lnConc, heteroNoise, noDilS, noD)
Calc_DilRes <- function(as = 3, bs = 1, cs = -4, ds = 10, at = 3, bt = 1, dt = 10, r = 0.0001, ct = cs - r,
                        sd_fac = 0.1, gt = 1, gs = 1, log_conc,
                        heteroNoise = FALSE, noDilSeries, noDils) {
  yAxfac <- (ds - as)
  log_dose <- log_conc
  isRef <- rep(c(1, 0), 1, each = length(log_conc) * noDilSeries)
  isSample <- rep(c(0, 1), 1, each = length(log_conc) * noDilSeries)
  # browser()
  av <- as * isRef + at * isSample + (ds * isRef + dt * isSample - as * isRef - at * isSample) /
    (1 + isRef * exp(bs * (cs - log_dose)) + isSample * exp(bt * (ct - log_dose)))
  if (heteroNoise) {
    # heterosc noise
    ro_jit <- matrix(unlist(map(av, function(x) x + rnorm(1, 0, x * sd_fac / 100))), nrow = noDils, ncol = noDilSeries * 2)
  } else {
    # homosc noise
    ro_jit <- matrix(unlist(map(av, function(x) x + rnorm(1, 0, sd_fac * yAxfac / 100))), nrow = noDils, ncol = noDilSeries * 2)
  }
  ro_jit <- abs(ro_jit)

  ro_jit2 <- cbind(ro_jit, log_dose)
  if (noDilSeries == 3) {
    colnames(ro_jit2) <- c("R_dil1", "R_dil2", "R_dil3", "T_dil1", "T_dil2", "T_dil3", "log_dose")
  } else {
    colnames(ro_jit2) <- c("R_dil1", "R_dil2", "T_dil1", "T_dil2", "log_dose")
  }
  return(ro_jit2)
}

#' Calculate the potency of linear PLA
#'
#' The delta Method is used for calculating the potency using a restricted linear model.
#' Absolute and relative CIs are calculated. The model RMSE or the Pure Error can be used.
#'
#' @param circles Data frame with the 3 consecutive concentrations and the respective readouts.
#' @param Lim The limits to check if the relative CI is within the limits.
#' @param PureErrFlag Boolean if Pure Error should be used.
#' @returns A data-frame with potency estimate, absolute CIs, test result, relative CIs.
#' @export
#' @examples
#' CIRC <- data.frame(
#'   log_dose = c(
#'     -2.5, -2.5, -2.5, -3.2, -3.2, -3.2, -3.9, -3.9, -3.9,
#'     -3.2, -3.2, -3.2, -3.9, -3.9, -3.9, -4.7, -4.7, -4.7
#'   ),
#'   replname = c(
#'     "R_dil1", "R_dil1", "R_dil1", "R_dil2", "R_dil2", "R_dil2", "R_dil3", "R_dil3", "R_dil3",
#'     "T_dil1", "T_dil1", "T_dil1", "T_dil2", "T_dil2", "T_dil2", "T_dil3", "T_dil3", "T_dil3"
#'   ),
#'   readout = c(72.1, 75.8, 76.04, 59.8, 61, 62.7, 43.6, 45, 41.5, 53.5, 62.2, 65.9, 48.3, 43.8, 43.14, 28.17, 29.2, 31.2),
#'   isRef = c(rep(1, 9), rep(0, 9)),
#'   isSample = c(rep(0, 9), rep(1, 9))
#' )
#' Lim <- c(rep(0, 8), 70, 130) # only Lim 9 and 10 relevant
#' PureErrF <- TRUE
#'
#'
#' LinPotTab(circles = CIRC, Lim, PureErrF)
LinPotTab <- function(circles, Lim, PureErrFlag) {
  circ_ABl <- circles
  circ_Al <- circ_ABl[circ_ABl$isSample == 1, ]
  circ_Al <- circ_ABl[circ_ABl$isSample == 0, ]
  # restr CSSI model
  modAB <- lm(readout ~ log_dose + isSample, circ_ABl)
  coeffs <- modAB$coefficients
  SU_modAB <- tryCatch(
    {
      SU_modAB <- summary(modAB)
    },
    error = function(msg) {
      return(NA)
    }
  )
  # Intercept diff/slope modAB
  linPot <- exp(modAB$coefficients[3] / modAB$coefficients[2])

  if (PureErrFlag) {
    FitAnova <- anova(lm(readout ~ factor(log_dose) * isSample, circ_ABl))
    meanPureErr <- FitAnova[4, 3]
    DFsPure <- FitAnova[4, 1]
    VCOV <- vcov(modAB)
    V_V <- VCOV / SU_modAB$sigma^2
    VCOVpure <- V_V * meanPureErr
    SEsPure <- sqrt(diag(V_V) * meanPureErr)
  }

  log_pot_delta <- deltaMethod(modAB, "isSample/log_dose")
  if (PureErrFlag) {
    V_ <- log_pot_delta$SE^2 / SU_modAB$sigma^2
    V_p <- V_ * meanPureErr
    potDeltaPureSE <- sqrt(V_p)
    CI_log_low <- log_pot_delta$Estimate - qt(0.975, DFsPure) * potDeltaPureSE
    CI_log_up <- log_pot_delta$Estimate + qt(0.975, DFsPure) * potDeltaPureSE
  } else {
    CI_log_low <- log_pot_delta$Estimate - qt(0.975, df.residual(modAB)) * log_pot_delta$SE
    CI_log_up <- log_pot_delta$Estimate + qt(0.975, df.residual(modAB)) * log_pot_delta$SE
  }
  # browser()
  ExpLinPot <- exp(c(log_pot_delta$Estimate, CI_log_low, CI_log_up))

  # Rel pot CI
  relLinpotCI <- ExpLinPot / linPot * 100
  if (relLinpotCI[2] > Lim[[9]] & relLinpotCI[3] < Lim[[10]]) test_potCI <- 0 else test_potCI <- 1

  pottab <- cbind(
    round(linPot * 100, 3), round(ExpLinPot[2] * 100, 3), round(ExpLinPot[3] * 100, 3),
    round(test_potCI, 3), round(relLinpotCI[2], 3), round(relLinpotCI[3], 3)
  )
  colnames(pottab) <- c("Potency", "lower 95%CI", "upper 95%CI", "test_result", "lowerRel95%CI", "upperRel95%CI")
  return(pottab)
}


#' Calculate the ANOVA of linear PLA and restricted and unrestr. model summaries
#'
#' The delta Method is used for calculating the potency using a restricted linear model.
#' Absolute and relative CIs are calculated. The model RMSE or the Pure Error can be used.
#'
#' @param circles Data frame with the 3 consecutive concentrations and the respective readouts.
#' @param Lim The limits to check if the relative CI is within the limits.
#' @param PureErrFlag Boolean if Pure Error should be used.
#' @returns A data-frame with 1) data.frame with F-tests ANOVA and test results
#' 2) summary of unrestricted linear model 3) Slope difference +CIs
#' 4) summary of restricted linear model.
#' @export
#' @examples
#' dat <- data.frame(
#'   R_dil1 = c(10221, 18258, 31993, 49336, 68332, 83527, 95584, 102229),
#'   R_dil2 = c(10136, 19078, 31925, 49003, 68034, 83776, 95495, 101608),
#'   T_dil1 = c(10830, 19891, 33915, 52131, 70617, 85784, 95937, 102791),
#'   T_dil2 = c(11169, 20153, 34007, 52179, 69962, 85543, 96439, 102655),
#'   log_dose = c(-1.2029, -1.89712, -2.590267, -3.2834, -3.97656, -4.66917, -5.362323, -6.05334)
#' )
#' CIRC <- data.frame(
#'   log_dose = c(-2.590267, -2.590267, -3.2834, -3.2834, -3.97656, -3.97656, -2.590267, -2.590267, -3.2834, -3.2834, -3.97656, -3.97656),
#'   replname = c("R_dil1", "R_dil2", "R_dil1", "R_dil2", "R_dil1", "R_dil2", "T_dil1", "T_dil2", "T_dil1", "T_dil2", "T_dil1", "T_dil2"),
#'   readout = c(31993, 31925, 49336, 49003, 68332, 68034, 33915, 34007, 52131, 52179, 70617, 69962),
#'   isRef = c(rep(1, 6), rep(0, 6)),
#'   isSample = c(rep(0, 6), rep(1, 6))
#' )
#' Lim <- c(rep(0, 8), 70, 130) # only Lim 9 and 10 relevant
#' PureErrF <- TRUE
#'
#'
#' ANOVAlintests(ro_new, circles, Lim, PureErrF)
ANOVAlintests <- function(ro_new, circles, Lim, PureErrFlag) {
  all_l <- melt(data.frame(ro_new), id.vars = "log_dose", variable.name = "replname", value.name = "readout")
  isRef <- rep(c(1, 0), 1, each = nrow(all_l) / 2)
  isSample <- rep(c(0, 1), 1, each = nrow(all_l) / 2)
  all_l$isRef <- isRef
  all_l$isSample <- isSample
  all_l$Conc <- exp(all_l$log_dose)
  all_lA <- all_l[all_l$isSample == 1, ] # TEST
  all_lB <- all_l[all_l$isSample == 0, ] # REF
  # browser()
  circ_ABl <- circles
  circ_Al <- circ_ABl[circ_ABl$isSample == 1, ]
  circ_Bl <- circ_ABl[circ_ABl$isSample == 0, ]
  # restr CSSI model
  modAB <- lm(readout ~ log_dose + isSample, circ_ABl)
  # unrestr with interact term SSSI
  modABu <- lm(readout ~ log_dose + isSample + log_dose * isSample, circ_ABl)
  modA <- lm(readout ~ log_dose + isSample, circ_Al)
  modB <- lm(readout ~ log_dose + isSample, circ_Bl)

  log_pot_delta <- deltaMethod(modAB, "isSample/log_dose")
  CI_log_low <- log_pot_delta$Estimate - qt(0.975, df.residual(modAB)) * log_pot_delta$SE
  CI_log_up <- log_pot_delta$Estimate + qt(0.975, df.residual(modAB)) * log_pot_delta$SE
  ExpLinPot <- exp(c(log_pot_delta$Estimate, CI_log_low, CI_log_up))
  if (ExpLinPot[2] * 100 > Lim[9] & ExpLinPot[3] * 100 > Lim[10]) test_potCI <- 0 else test_potCI <- 1

  su_mod <- summary(modAB)$coefficients
  su_mod2 <- cbind(data.frame(parameter = c("intercept REF", "slope REF", "intercepts diff.")), su_mod)
  su_modU <- summary(modABu)$coefficients
  su_modU2 <- cbind(data.frame(parameter = c("intercept REF", "slope REF", "intercepts diff.", "slope difference")), su_modU)

  uCI_SloDiff <- su_modU[4, 1] + qt(0.975, 8) * su_modU[4, 2]
  lCI_SloDiff <- su_modU[4, 1] - qt(0.975, 8) * su_modU[4, 2]
  SlopeDiffCI <- c(su_modU[4, 1], lCI_SloDiff, uCI_SloDiff)

  lenCirc <- nrow(circ_ABl)
  dfTreat <- 3
  dfPrep <- 1
  dfReg <- 1
  dfnonP <- 1
  dfRMSE <- c(lenCirc - 3 - 1)
  dfTotal <- lenCirc - 1
  dfPureE <- lenCirc - 6
  dfNonLin <- dfRMSE - dfPureE

  RSS <- sum(resid(lm(readout ~ log_dose * isSample, circ_ABl))^2)
  MSE <- RSS / dfRMSE
  SSE <- sum(resid(lm(readout ~ factor(log_dose) * isSample, circ_ABl))^2)
  MSpure <- SSE / dfPureE

  if (PureErrFlag) {
    FitAnova <- anova(lm(readout ~ factor(log_dose) * isSample, circ_ABl))
    meanPureErr <- FitAnova[4, 3]
    SU_modAB <- tryCatch(
      {
        SU_modAB <- summary(modAB)
      },
      error = function(msg) {
        return(NA)
      }
    )
    if (length(SU_modAB) > 1) s_modABcoeffs <- summary(modAB)$coefficients

    DFsPure <- FitAnova[4, 1]
    VCOV <- vcov(modAB)
    V_V <- VCOV / SU_modAB$sigma^2
    VCOVpure <- V_V * meanPureErr
    SEsPure <- sqrt(diag(V_V) * meanPureErr)
    su_mod2[, 3] <- SEsPure
    su_mod2[, 4] <- su_mod2[, 2] / su_mod2[, 3]
    su_mod2[, 5] <- 2 * (1 - pt(abs(su_mod[, 4]), FitAnova[4, 1]))

    s_mu <- summary(modABu)$coefficients
    SU_modABu <- summary(modABu)
    VCOVu <- vcov(modABu)
    V_Vu <- VCOVu / SU_modABu$sigma^2
    SEsPureU <- sqrt(diag(V_Vu) * meanPureErr)

    su_modU2[, 3] <- SEsPureU
    su_modU2[, 4] <- su_modU2[, 2] / su_modU2[, 3]
    su_modU2[, 5] <- 2 * (1 - pt(abs(su_modU2[, 4]), FitAnova[4, 1]))

    uCI_SloDiffP <- su_modU[4, 1] + qt(0.975, 8) * SEsPureU[4]
    lCI_SloDiffP <- su_modU[4, 1] - qt(0.975, 8) * SEsPureU[4]
    SlopeDiffCI <- c(su_modU[4, 1], lCI_SloDiffP, uCI_SloDiffP)

    SSRes <- SSE
    dfRes <- dfPureE
  } else {
    SSRes <- RSS
    dfRes <- dfRMSE
  }

  # treatment
  SStreat <- print(sum((predict(lm(readout ~ factor(log_dose) * isSample, circ_ABl)) - mean(circ_ABl$readout))^2))
  F_treat <- (SStreat / dfTreat) / (SSRes / dfRes)
  # Preparation
  SSprep <- print(sum((predict(lm(readout ~ isSample, circ_ABl)) - mean(circ_ABl$readout))^2))
  F_prep <- (SSprep / dfTreat) / (SSRes / dfRes)
  # Regression
  # ANOVA tape II SS of regression
  SSreg <- Anova(lm(readout ~ log_dose + isSample, circ_ABl))[1, 1]
  # Non-parallelism
  # diff of RSS of restricted and unrestricted model
  SSnonpar <- sum(resid(modAB)^2) - sum(resid(modABu)^2)
  F_nonpar <- SSnonpar / (sum(resid(lm(readout ~ factor(log_dose) * isSample, circ_ABl))^2) / (lenCirc - 4))

  # non-linearity
  SSnonlin <- sum((predict(modABu) - predict(lm(readout ~ as.factor(log_dose) * isSample, circ_ABl)))^2)
  # = RSS-SSE
  # Total SS
  SStot <- sum((circ_ABl$readout - mean(circ_ABl$readout))^2)
  # Significance of R^2 F-ratio
  # MSR/MSE
  # sample A
  F_R2_A <- sum((predict(lm(readout ~ log_dose + I(log_dose^2), circ_Al)) - mean(predict(modA)))^2 - (predict(modA) - mean(circ_Al$readout))^2) /
    (sum((predict(lm(readout ~ log_dose + I(log_dose^2), circ_Al)) - circ_Al$readout)^2) / (nrow(circ_Al) - 3))
  pFR2_A <- round(pf(F_R2_A, 1, 6), 4)
  # sample B
  F_R2_B <- sum((predict(lm(readout ~ log_dose + I(log_dose^2), circ_Bl)) - mean(predict(modB)))^2 - (predict(modB) - mean(circ_Bl$readout))^2) /
    (sum((predict(lm(readout ~ log_dose + I(log_dose^2), circ_Bl)) - circ_Bl$readout)^2) / (nrow(circ_Bl) - 3))
  pFR2_B <- round(pf(F_R2_B, 1, 6), 4)
  # sign of non-lin with pure error: MSSnonlin/MSSE
  F_nonlin <- (SSnonlin / 2) / (SSE / dfPureE)

  # sign of slope
  F_slope_B <- sum((predict(modB) - mean(circ_Bl$readout))^2) / (sum((circ_Bl$readout - predict(modB))^2) / (nrow(circ_Bl) - 2))
  F_slope_A <- sum((predict(modA) - mean(circ_Al$readout))^2) / (sum((circ_Al$readout - predict(modA))^2) / (nrow(circ_Al) - 2))
  # F-test on regression: MSSreg/MSSE
  if (is.na(F_nonlin)) F_nonlin <- 0
  if (F_nonlin > 0) {
    p_F_nonlin <- round(pf(F_nonlin, 2, dfPureE, lower.tail = F), 5)
  } else {
    p_F_nonlin <- "SSnonlin neg or 0"
  }

  # significances
  F_regr <- (SSreg / 1) / (SSRes / dfRes)
  p_F_regr <- round(pf(F_regr, 1, dfRes, lower.tail = F), 3)
  p_F_treat <- round(pf(F_treat, 3, dfRes, lower.tail = F), 3)
  p_F_prep <- round(pf(F_prep, 1, dfRes, lower.tail = F), 3)
  p_F_slope_A <- round(pf(F_slope_A, 1, (nrow(circ_Al) - 2), lower.tail = F), 3)
  p_F_slope_B <- round(pf(F_slope_B, 1, (nrow(circ_Bl) - 2), lower.tail = F), 3)
  p_F_nonp <- round(pf(F_nonpar, 1, dfRes, lower.tail = F), 3)
  p_F_LoF <- p_F_nonlin

  res_tab_lin <- data.frame(
    test = c(
      "F-test on sign. of regression", "F_test on non-lin",
      "F-test on R^2 A", "F_test on R^2 B",
      "F-test on slope A", "F-test on slope B",
      "F-test on non-parallelism", "F-test on preparation"
    ),
    test_results = c(
      ifelse(p_F_regr < 0.05, 0, 1), ifelse(p_F_nonlin < 0.05, 1, 0),
      ifelse(pFR2_A < 0.05, 1, 0), ifelse(pFR2_B < 0.05, 1, 0),
      ifelse(p_F_slope_A < 0.05, 0, 1), ifelse(p_F_slope_B < 0.05, 0, 1),
      ifelse(p_F_nonp < 0.05, 1, 0), ifelse(p_F_prep < 0.05, 0, 1)
    ),
    estimate = c(
      p_F_regr, p_F_nonlin, pFR2_A, pFR2_B, p_F_slope_A,
      p_F_slope_B, p_F_nonp, p_F_prep
    ),
    Source = c(
      "Treatment", "Preparation", "Regression", "Non-parallelism",
      "Resid Error", "Non-linearity", "Pure error", "Total"
    ),
    df = c(dfTreat, 1, 1, 1, dfRMSE, 2, dfPureE, lenCirc - 1),
    SumSquares = c(
      round(SStreat, 5), round(SSprep, 5), round(SSreg, 5),
      round(SSnonpar, 5), round(RSS, 5), round(SSnonlin, 5),
      round(SSE, 5), round(SStot, 5)
    ),
    MS = c(
      round(SStreat / dfTreat, 5), round(SSprep, 5), round(SSreg, 5),
      round(SSnonpar, 5), round(RSS / dfRMSE, 5), round(SSnonlin / 2, 5),
      round(SSE / dfPureE, 5), round(SStot / dfTotal, 5)
    ),
    "F-value" = c(
      round(F_treat, 5), round(F_prep, 5), round(F_regr, 5),
      round(F_nonpar, 5), "", round(F_nonlin, 5), "", ""
    ),
    "p-value" = c(p_F_treat, p_F_prep, p_F_regr, p_F_nonp, "", p_F_LoF, "", "")
  )
  RET <- list(res_tab_lin, su_modU2, SlopeDiffCI, su_mod2)
  return(RET)
}


#' Plots linear regression on the scatterplot.
#'
#' Restricted and unrestricted model are plotted. The number of dilution steps used for the regression
#' are printed in the title. Reference is blue, red is the sample. The points used for regression
#' are highlighted with a circle.
#'
#' @param circle Data frame with the 3 consecutive concentrations and the respective readouts.
#' @param sigmoid Parameters of the unrestricted fit of the full dataset for drawing the curves.
#' @param all_l2 Long format data.frame of the readout data incl log(concentration) and isRef (0s and 1s).
#' @param pl_df  Data.frame for plotting the regression lines.
#' @param indS Index of dilution, where regression starts for reference sample.
#' @param indT Index of dilution, where regression starts for test sample.
#' @returns A grid object of 2 plots of unrestricted and restricted linear PLA models.
#' @export
#' @examples
#' circle <- read.csv("~/plateflow/circle.csv")
#' all_l2 <- read.csv("~/plateflow/all_l2.csv")
#' sigmoid <- c(10.0, 10.0, 110.0, 110.0, 1.0, 1.0, -3.5, 0.0)
#' indS <- 3
#' indT <- 3
#' pl_df <- data.frame(
#'   lnC = c(-1.203973, -1.897120, -2.590267, -3.283414, -3.976562, -4.669176, -5.362323, -6.053340),
#'   plotS = c(113.772511, 97.668371, 81.564231, 65.460091, 49.355952, 33.264200, 17.160060, 1.105405),
#'   plotT = c(114.213375, 97.588663, 80.963951, 64.339239, 47.714527, 31.102604, 14.477892, -2.095735)
#' )
#'
#'
#' PlotLinPLA_FUNC(circle, sigmoid, all_l2, pl_df, indS, indT)
PlotLinPLA_FUNC <- function(circle, sigmoid, all_l2, pl_df, indS, indT) {
  # browser()
  mLin <- gsl_nls(readout ~ (intS + r) * isSample + intS * isRef + k * log_dose,
    data = circle,
    start = list(intS = 0, k = 1, r = 0),
    control = gsl_nls_control(xtol = 1e-10, ftol = 1e-10, gtol = 1e-10)
  )
  # alternativ: modAB <- lm(readout ~ log_dose+isSample, circle)
  sum_mLin <- summary(mLin)

  log_dose <- unique(all_l2$log_dose)
  seq_x <- seq(min(log_dose), max(log_dose), 0.1)
  if (!is.null(sigmoid)) {
    SAMPLEtrue <- sigmoid[5] + (sigmoid[7] - sigmoid[5]) / (1 + exp(sigmoid[6] * ((sigmoid[4] - sigmoid[8] - seq_x))))
    REFtrue <- sigmoid[1] + (sigmoid[3] - sigmoid[1]) / (1 + exp(sigmoid[2] * ((sigmoid[4] - seq_x))))
    truePL_df <- cbind(seq_x, SAMPLEtrue, REFtrue)
  } else {
    SAMPLEtrue <- NULL
    REFtrue <- NULL
    truePL_df <- NULL
  }


  p <- ggplot(all_l2, aes(x = log_dose, y = readout, color = factor(isRef))) +
    geom_point(size = 2) +
    # labs(title=paste("linear regression model", indS,indT), color="product") +
    scale_colour_manual(labels = c("test", "reference"), values = c("#C2173F", "#4545BA")) +
    ylim(min(all_l2$readout), max(all_l2$readout)) +
    scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
    scale_y_continuous(breaks = scales::pretty_breaks(n = 10)) +
    theme_bw()
  p2 <- p + geom_line(
    data = pl_df, aes(x = lnC, y = plotS), color = "#4545BA",
    inherit.aes = F
  ) +
    geom_line(
      data = pl_df, aes(x = lnC, y = plotT), color = "#C2173F",
      inherit.aes = F
    ) +
    {
      if (!is.null(truePL_df)) {
        geom_line(
          data = data.frame(truePL_df), aes(x = seq_x, y = SAMPLEtrue), color = "#C2173F", linetype = 2, alpha = 0.4,
          inherit.aes = F
        )
      }
    } +
    {
      if (!is.null(truePL_df)) {
        geom_line(
          data = data.frame(truePL_df), aes(x = seq_x, y = REFtrue), color = "#4545BA", linetype = 2, alpha = 0.4,
          inherit.aes = F
        )
      }
    } +
    labs(title = paste("unrestricted PLA model"), subtitle = paste("Regression starts for reference sample:", indS, "for test sample:", indT)) +
    theme(legend.position = "none", axis.text = element_text(size = 14))
  p3 <- p2 + geom_point(circle, mapping = aes(
    x = log_dose, y = readout, shape = factor(isRef),
    size = 5, alpha = 0.2
  ), col = c("black"), inherit.aes = FALSE) +
    scale_shape_manual(labels = c("test", "reference"), values = c(21, 21))
  # fit intercept for test and ref and common slope


  pl_restS <- sum_mLin$coefficients[1, 1] + sum_mLin$coefficients[2, 1] * log_dose
  pl_restT <- sum_mLin$coefficients[1, 1] + sum_mLin$coefficients[3, 1] + sum_mLin$coefficients[2, 1] * log_dose
  pl_rest <- data.frame(lnC = log_dose, plotS = pl_restS, plotT = pl_restT)

  pr2 <- p + geom_line(
    data = pl_rest, aes(x = lnC, y = plotS), color = "#4545BA",
    inherit.aes = F
  ) +
    geom_line(
      data = pl_rest, aes(x = lnC, y = plotT), color = "#C2173F",
      inherit.aes = F
    ) +
    {
      if (!is.null(truePL_df)) {
        geom_line(
          data = data.frame(truePL_df), aes(x = seq_x, y = SAMPLEtrue), color = "#C2173F", linetype = 2, alpha = 0.4,
          inherit.aes = F
        )
      }
    } +
    {
      if (!is.null(truePL_df)) {
        geom_line(
          data = data.frame(truePL_df), aes(x = seq_x, y = REFtrue), color = "#4545BA", linetype = 2, alpha = 0.4,
          inherit.aes = F
        )
      }
    } +
    labs(
      title = paste("restricted linear regression model"),
      subtitle = paste("Regression on highlighted points")
    ) +
    theme(legend.position = "none", axis.text = element_text(size = 14))
  pr3 <- pr2 + geom_point(circle, mapping = aes(
    x = log_dose, y = readout, shape = factor(isRef),
    size = 5, alpha = 0.2
  ), col = c("black"), inherit.aes = FALSE) +
    scale_shape_manual(labels = c("test", "reference"), values = c(21, 21))
  return(grid.arrange(p3, pr3, nrow = 1))
}


#' Calculates the potency of 4PL PLA of all models model
#'
#' The gradient method is used for calculating the potency for a restricted model, an unrestricteed model,
#' and the log-transformations thereof.
#' Absolute and relative CIs are calculated. The model RMSE or the Pure Error can be used.
#'
#' @param ro_new The dilution series readouts for REF and TEST + 1 column log(concentration).
#' @param PureErrFlag Boolean if Pure Error should be used.
#' @returns A data-frame with 1) data.frame with F-tests ANOVA and test results
#' 2) summary of unrestricted linear model 3) Slope difference +CIs
#' 4) summary of restricted linear model.
#' @export
#' @examples
#' ro_new <- data.frame(
#'   R_dil1 = c(10221, 18258, 31993, 49336, 68332, 83527, 95584, 102229),
#'   R_dil2 = c(10136, 19078, 31925, 49003, 68034, 83776, 95495, 101608),
#'   T_dil1 = c(10830, 19891, 33915, 52131, 70617, 85784, 95937, 102791),
#'   T_dil2 = c(11169, 20153, 34007, 52179, 69962, 85543, 96439, 102655),
#'   log_dose = c(-1.2029, -1.89712, -2.590267, -3.2834, -3.97656, -4.66917, -5.362323, -6.05334)
#' )
#' PureErrF <- TRUE
#'
#'
#' pot4plFUNC(ro_new, PureErrF)
pot4plFUNC <- function(ro_new, PureErrFlag) {
  all_l <- melt(data.frame(ro_new), id.vars = "log_dose", variable.name = "replname", value.name = "readout")
  isRef <- rep(c(1, 0), 1, each = nrow(all_l) / 2)
  isSample <- rep(c(0, 1), 1, each = nrow(all_l) / 2)
  all_l$isRef <- isRef
  all_l$isSample <- isSample
  all_l$Conc <- exp(all_l$log_dose)
  all_l$readout[all_l$readout < 0] <- 0.01
  all_l$readouttrans <- log(all_l$readout)
  # browser()
  CORdat <- cor(ro_new[, 1], ro_new[, ncol(ro_new)])
  if (CORdat < 0) SLOPE <- -1 else SLOPE <- 1
  #
  FITs <- Fitting_FUNC(ro_new, TransFlag = F)
  if (!PureErrFlag) {
    pot_est <- FITs[[3]]
    potU_est <- FITs[[4]]
  } else {
    FitAnova <- anova(lm(readout ~ factor(Conc) * isSample, all_l))
    meanPureErr <- FitAnova[4, 3]
    DFsPure <- FitAnova[4, 1]
    # SU_mr <- tryCatch({
    #   SU_mr <- summary(mr)
    # }, error = function(msg) {
    #   return()
    # })
    SU_mr <- FITs[[1]]
    SU_mrCoeff <- SU_mr$coefficients
    # if (length(SU_mr)>1) {
    #   SU_mrCoeff <- SU_mr$coefficients
    # } else { SU_mrCoeff <- rep(NA,5) }
    # te2$cov.unscaled[1,1]*te2$sigma^2
    # VCOV <- vcov(mr)
    # V_V <- VCOV/SU_mr$sigma^2
    V_V <- SU_mr$cov.unscaled
    SEsPure <- sqrt(diag(V_V) * meanPureErr)
    pot_est <- c(
      exp(SU_mrCoeff["r", 1]), exp(SU_mrCoeff["r", 1] - qt(0.975, DFsPure) * SEsPure["r"]),
      exp(SU_mrCoeff["r", 1] + qt(0.975, DFsPure) * SEsPure["r"])
    )
    # unrestricted
    SU_mu <- FITs[[2]]
    s_muCoeff <- SU_mu$coefficients
    # SU_mu <- summary(mu)
    # VCOVu <- vcov(mu)
    V_Vu <- SU_mu$cov.unscaled
    SEsPureU <- sqrt(diag(V_Vu) * meanPureErr)
    potU_est <- c(
      exp(s_muCoeff["r", 1]), exp(s_muCoeff["r", 1] - qt(0.975, DFsPure) * SEsPureU["r"]),
      +exp(s_muCoeff["r", 1] + qt(0.975, DFsPure) * SEsPureU["r"])
    )
  } # PureErrFlag

  FITstrans <- Fitting_FUNC(ro_new, TransFlag = TRUE)


  # pot_est_log <- exp(confintd(mr_log, "r", method="asymptotic"))
  # potU_est_log <- exp(confintd(mu_log, "r", method="asymptotic"))
  pot_est_log <- FITstrans[[3]]
  potU_est_log <- FITstrans[[4]]
  colnames(pot_est_log) <- c("estimate", "lowerCI2", "upperCI")
  colnames(potU_est_log) <- c("estimate", "lowerCI2", "upperCI")
  # browser()
  su_mr_log <- FITstrans[[1]] # summary(mr_log)
  Dat$RMSE_Rlog <- su_mr_log$sigma
  su_mu_log <- FITstrans[[2]] # summary(mu_log)
  Dat$RMSE_Ulog <- su_mu_log$sigma
  Dat$up_lowAslog <- su_mu_log$coefficients[1, 1] - su_mu_log$coefficients[4, 1]
  potALL <- rbind(round(pot_est, 6), round(potU_est, 6), round(pot_est_log, 6), round(potU_est_log, 6))

  potALL2 <- cbind(c("restricted", "unrestricted", "ln-transformed restr", "ln-transformed unrestr"), potALL)
  return(potALL2)
}

#' Calculates logarithmized CIs for ratios
#'
#' Ratio CIs can be unbound, and not calulatable. Therefore the ratios are logarithmized.
#' The Covariance is not used, as it is 0 for comparisons between products.
#'
#' @param xs The estimate of the reference sample.
#' @param xt The estimate of the test sample.
#' @param se_xs The standard error of the estimate of the reference sample.
#' @param se_xt The standard error of the estimate of the test sample.
#' @param CoVar The covariance between test and reference estimate (set to 0).
#' @param DFs Degrees of freedom, either from pure error term (no of readouts - no of concentrations*2), or RMSE term (no of readouts - 8 parameters).
#' @param Conf The confidence of the Confidence Interval (default 95% which requires 0.975).
#' @returns A data-frame with the lower and upper CI in anti-log form.
#' @export
#' @examples
#' xs <- 2
#' xt <- 3.2
#' se_xt <- 0.34
#' se_xs <- 0.23
#' DFs <- 32 - 16
#'
#'
#' ParamCI_F(xt, xs, se_xt, se_xs, CoVar, DFs, Conf = 0.975)
ParamCI_F <- function(xt, xs, se_xt, se_xs, CoVar, DFs, Conf = 0.975) {
  log_xs <- log(abs(xs))
  log_xt <- log(abs(xt))
  var_log_xs <- (se_xs / xs)^2 # approximate variance of log(xs)
  var_log_xt <- (se_xt / xt)^2
  se_log_ratio <- sqrt(var_log_xs + var_log_xt) #-2*CoVar/(xs*xt)

  lower_log_ratio <- log_xt - log_xs - qt(Conf, DFs) * se_log_ratio
  upper_log_ratio <- log_xt - log_xs + qt(Conf, DFs) * se_log_ratio
  ci_ratio <- exp(c(lower_log_ratio, upper_log_ratio))
  return(ci_ratio)
}


#' Calculates EQ tests and F-Tests for 4P
#'
#' CIs of asmptote ratios,slope ratios, asympt-difference ratio, and lower asympt-diff. calculated.
#' In addition the F-test on significance of regression and non-linearity are calculated.
#'
#' @param ro_new Data frame with the concentrations and the respective readouts.
#' @param Lim The limits to check if the relative CI is within the limits.
#' @param PureErrFlag Boolean if Pure Error should be used.
#' @returns A data-frame with the lower and upper CI in anti-log form.
#' @export
#' @examples
#'
#' dat <- data.frame(
#'   REF1 = c(1547, 1620, 1644, 2504, 3426, 3512, 3401, 3787), REF2 = c(1492, 1536, 1384, 2286, 3046, 3479, 3516, 3497),
#'   REF3 = c(1468, 1827, 1558, 2252, 3002, 3349, 2945, 3665),
#'   TEST1 = c(1405, 1523, 1502, 1474, 2383, 3221, 3589, 3445), TEST2 = c(1420, 1516, 1544, 1512, 2226, 3219, 3327, 3591),
#'   TEST3 = c(1399, 1376, 1588, 1475, 2148, 3083, 2942, 3466), log_dose = c(5.01, 3.401, 2.708, 2.015, 1.32176, 0.62861, -0.0645385, -1.6739764)
#' )
#' Lim <- c(-1, 1, 0.005, 2, 0.5, 2, 0.5, 2, 75, 133, 75, 133)
#' PureErrF <- FALSE
#'
#' tests_FUNC(ro_new, Lim, PureErrF)
tests_FUNC <- function(ro_new, Lim, PureErrFlag) {
  all_l <- melt(data.frame(ro_new), id.vars = "log_dose", variable.name = "replname", value.name = "readout")
  isRef <- rep(c(1, 0), 1, each = nrow(all_l) / 2)
  isSample <- rep(c(0, 1), 1, each = nrow(all_l) / 2)
  all_l$isRef <- isRef
  all_l$isSample <- isSample
  all_l$Conc <- exp(all_l$log_dose)
  all_l$readout[all_l$readout < 0] <- 0.01
  # browser()
  FITs <- Fitting_FUNC(ro_new = ro_new, TransFlag = FALSE)
  if (is.character(FITs)) {
    return(FITs)
  } # if singularity

  POTr_CI <- FITs[[3]][2:3]
  potAll2 <- FITs[[3]]

  smr <- FITs[[1]]
  smu <- FITs[[2]]
  FitAnova <- anova(lm(readout ~ factor(Conc) * isSample, all_l))
  # pure error
  pureSS <- FitAnova[4, 2]
  pureSS_df <- FitAnova[4, 1]
  meanPureErr <- FitAnova[4, 3]
  # vcovMU <- vcov(mu)
  V_V <- smu$cov.unscaled
  SEsPure <- sqrt(diag(V_V) * meanPureErr)
  VCOVpure <- V_V * meanPureErr
  DFsPure <- FitAnova[4, 1]


  testPOTr <- logical()
  if (POTr_CI[1] * 100 > Lim[[9]] & POTr_CI[2] * 100 < Lim[[10]]) testPOTr <- 0 else testPOTr <- 1

  potAllU2 <- FITs[[4]]
  test_c <- logical()
  if ((potAllU2[2] > Lim[[9]] / 100 & potAllU2[3] < Lim[[10]] / 100)) test_c <- 0 else test_c <- 1
  predPot <- FITs[[5]]
  predPotU <- FITs[[6]]

  vcovMU <- V_V * smu$sigma^2

  #### ANOVA ----
  noConc <- length(unique(all_l$Conc))
  nofitted <- noConc
  AnovaDFs <- c(nofitted - 1, 1, 3, nofitted - 4 - 1, nrow(all_l) - nofitted, nofitted, nrow(all_l) - 2 * nofitted, nrow(all_l) - 1)
  SStreat <- round(sum((predPotU - mean(all_l$readout))^2), 5)
  SSregr <- round(sum((predPot - mean(all_l$readout))^2), 5)
  # non-parallelism
  SSnonparall <- round(sum(smr$residuals^2) - sum(smu$residuals^2), 5)
  SSprep <- round(sum((predict(lm(readout ~ isSample, all_l)) - mean(all_l$readout))^2), 5)

  RSS <- round(sum(smu$residuals^2), 5)
  RSS_df <- AnovaDFs[5]
  MSEunr <- RSS / RSS_df
  RMSEunr <- sqrt(RSS / RSS_df)
  # Pure Err
  FitAnova <- anova(lm(readout ~ factor(Conc) * isSample, all_l))
  SSE <- sum(resid(lm(readout ~ factor(Conc) * isSample, all_l))^2) # =FitAnova[4,2]
  SSE_df <- FitAnova[4, 1]
  PureMSE <- SSE / SSE_df
  RMSE_pure <- sqrt(PureMSE)
  ## non-lin = LoF
  if (PureErrFlag) {
    ERR <- PureMSE
    ERR_df <- SSE_df
  } else {
    ERR <- MSEunr
    ERR_df <- RSS_df
  }
  SSnonlin <- sum((predict(lm(readout ~ factor(Conc) * isSample, all_l)) - predPotU)^2)
  LoF_df <- FitAnova[1, 1] + FitAnova[2, 1]
  F_regr <- (SSregr / AnovaDFs[3]) / ERR
  p_F_regr <- round(pf(F_regr, AnovaDFs[3], ERR_df, lower.tail = F), 5)
  if (ncol(ro_new) < 4) F_nonlin <- 0 else F_nonlin <- (SSnonlin / AnovaDFs[6]) / ERR
  if (F_nonlin > 0) {
    p_F_nonlin <- round(pf(F_nonlin, AnovaDFs[6], ERR_df, lower.tail = F), 5)
  } else {
    p_F_nonlin <- "SSnonlin neg or single dilutions"
  }

  test_a <- test_b <- test_d <- test_ad <- logical()

  RSS_r <- round(sum(smr$residuals^2), 5)
  MSE_r <- RSS_r / (nrow(all_l) - 5)
  RMSE_r <- round(sqrt(MSE_r), 6)
  Dat$RMSE_r <- RMSE_r
  Dat$RMSE_pure <- RMSE_pure
  Dat$RMSE_unr <- round(RMSEunr, 6)

  coeffs <- smu$coefficients[, 1]
  # browser()
  ## EQ test on lower As diff
  as <- coeffs["as"]
  at <- coeffs["at"]
  lAs_diff <- (at - as)
  uCI_laDiff <- lAs_diff + qt(0.975, smu$df[2]) * sqrt(smu$coefficients["ds", 2]^2 + smu$coefficients["dt", 2]^2)
  lCI_laDiff <- lAs_diff - qt(0.975, smu$df[2]) * sqrt(smu$coefficients["ds", 2]^2 + smu$coefficients["dt", 2]^2)
  if (uCI_laDiff < Lim[[2]] & lCI_laDiff > Lim[[1]]) test_la_diff <- 0 else test_la_diff <- 1

  #### EQ test on upper asymptote ratio ----
  # as <- coeffs["as"]
  # at <- coeffs["at"]
  # uAsRatio <- compParm(potU, "a","/",display=F)
  # uAsCI <- c(uAsRatio[1]-qt(0.975,RSS_df)*uAsRatio[2], uAsRatio[1]+qt(0.975,RSS_df)*uAsRatio[2])
  # browser()
  ds <- smu$coefficients["ds", 1]
  dt <- smu$coefficients["dt", 1]
  if (PureErrFlag) se_ds <- sqrt(VCOVpure["ds", "ds"]) else se_ds <- smu$coefficients["ds", 2]
  if (PureErrFlag) se_dt <- sqrt(VCOVpure["dt", "dt"]) else se_dt <- smu$coefficients["dt", 2]
  if (PureErrFlag) CoVarlog_d <- VCOVpure["dt", "ds"] else CoVarlog_d <- vcovMU["dt", "ds"]
  if (PureErrFlag) DFs <- DFsPure else DFs <- nrow(all_l) - 8
  uAsCI2 <- ParamCI_F(dt, ds, se_dt, se_ds, CoVarlog_d, DFs, Conf = 0.975)
  if (uAsCI2[1] > Lim[[7]] & uAsCI2[2] < Lim[[8]]) test_a <- 0 else test_a <- 1
  estUppA <- round(at / as, 5)

  Dat$uAsCI <- uAsCI2

  #### EQ test on slope ratio ----
  # bs <- coeffs["bs"]
  # bt <- coeffs["bt"]
  # slopeRatio <- compParm(potU, "b","/",display=F)
  # slopeCI <- c(slopeRatio[1,1]-qt(0.975,RSS_df)*slopeRatio[1,2], slopeRatio[1,1]+qt(0.975,RSS_df)*slopeRatio[1,2])

  bs <- smu$coefficients["bs", 1]
  bt <- smu$coefficients["bt", 1]
  if (PureErrFlag) se_bs <- sqrt(VCOVpure["bs", "bs"]) else se_bs <- smu$coefficients["bs", 2]
  if (PureErrFlag) se_bt <- sqrt(VCOVpure["bt", "bt"]) else se_bt <- smu$coefficients["bt", 2]
  if (PureErrFlag) CoVarlog_b <- VCOVpure["bt", "bs"] else CoVarlog_b <- vcovMU["bt", "bs"]
  slopeCI2 <- ParamCI_F(bt, bs, se_bt, se_bs, CoVarlog_b, DFs, Conf = 0.975)
  if (slopeCI2[1] > Lim[[5]] & slopeCI2[2] < Lim[[6]]) test_b <- 0 else test_b <- 1
  estUppA <- round(at / as, 5)

  Dat$slopeRatioCI <- slopeCI2

  #### EQ test on lower As ratio ----

  # lAsRatio <- compParm(potU, "d","/",display=F)
  # slopeCI <- c(lAsRatio[1,1]-qt(0.975,RSS_df)*lAsRatio[1,2], lAsRatio[1,1]+qt(0.975,RSS_df)*lAsRatio[1,2])

  as <- smu$coefficients["as", 1]
  at <- smu$coefficients["at", 1]
  if (PureErrFlag) se_as <- sqrt(VCOVpure["as", "as"]) else se_as <- smu$coefficients["as", 2]
  if (PureErrFlag) se_at <- sqrt(VCOVpure["at", "at"]) else se_at <- smu$coefficients["at", 2]
  if (PureErrFlag) CoVarlog_a <- VCOVpure["at", "as"] else CoVarlog_a <- vcovMU["at", "as"]
  lAsCI2 <- ParamCI_F(at, as, se_at, se_as, CoVarlog_a, DFs, Conf = 0.975)
  if (lAsCI2[1] > Lim[[3]] & lAsCI2[2] < Lim[[4]]) test_d <- 0 else test_d <- 1
  estLowA <- round(at / as, 5)

  Dat$lAsCI <- lAsCI2

  #### EQtest on ratio of As difference ----
  AsDiffRatio <- (dt - at) / (ds - as)

  dt_at <- (dt - at)
  ds_as <- (ds - as)
  se_ds_asPure <- sqrt(VCOVpure["as", "as"] + VCOVpure["ds", "ds"] - 2 * VCOVpure["as", "ds"])
  se_dt_atPure <- sqrt(VCOVpure["at", "at"] + VCOVpure["dt", "dt"] - 2 * VCOVpure["at", "dt"])
  se_ds_asRMSE <- sqrt(vcovMU["as", "as"] + vcovMU["ds", "ds"] - 2 * vcovMU["as", "ds"])
  se_dt_atRMSE <- sqrt(vcovMU["at", "at"] + vcovMU["dt", "dt"] - 2 * vcovMU["at", "dt"])
  if (PureErrFlag) se_ds_as <- se_ds_asPure else se_ds_as <- se_ds_asRMSE
  if (PureErrFlag) se_dt_at <- se_dt_atPure else se_dt_at <- se_dt_atRMSE

  AsDiffCI2 <- ParamCI_F(dt_at, ds_as, se_dt_at, se_ds_as, CoVar = 0, DFs, Conf = 0.975)
  if (AsDiffCI2[1] > Lim[[11]] & AsDiffCI2[2] < Lim[[12]]) test_ad <- 0 else test_ad <- 1
  estLowA <- round(at / as, 5)

  Dat$up_lowAs <- abs(ds - as)

  lowerCIlowerA <- lAsCI2[1]
  lowerCIupperA <- uAsCI2[1]
  upperCIlowerA <- lAsCI2[2]
  upperCIupperA <- uAsCI2[2]
  test_lowA <- test_d
  test_uppA <- test_a
  # browser()
  res_tab <- data.frame(
    test = c(
      "F-test on sign. of regression*",
      "EQ test on lower asymptotes difference",
      "EQ test ratio of lower asymptotes",
      "EQ test ratio of Hill slopes",
      "EQ test ratio of upper asymptotes",
      "F-test on non-linearity*",
      "EQ test ratio of asymptote difference",
      "geom. rel. CI restr. model",
      "geom. rel. CI unrestr. model"
    ),
    test_results = c(
      ifelse(p_F_regr < 0.05, 0, 1), test_la_diff, test_lowA, test_b, test_uppA,
      ifelse(p_F_nonlin > 1, 1, ifelse(p_F_nonlin < 0.05, 1, 0)), test_ad,
      testPOTr, test_c
    ),
    estimate = c(
      round(p_F_regr, 3), round(lAs_diff, 5),
      estLowA, round(bs / bt, 5), estUppA, p_F_nonlin,
      round(dt_at / ds_as, 5), round(potAll2[1] * 100, 2), round(potAllU2[1] * 100, 2)
    ),
    lower_limit = c("-", Lim[[1]], Lim[[3]], Lim[[5]], Lim[[7]], "-", Lim[[11]], Lim[[9]], Lim[[9]]),
    upper_limit = c("-", Lim[[2]], Lim[[4]], Lim[[6]], Lim[[8]], "-", Lim[[12]], Lim[[10]], Lim[[10]]),
    lower_CI = c(
      RMSE_r, round(lCI_laDiff, 3), round(lAsCI2[1], 5), round(slopeCI2[1], 5),
      round(uAsCI2[1], 5), "-", round(AsDiffCI2[1], 5), round(potAll2[2], 2), round(potAllU2[2], 2)
    ),
    upper_CI = c(
      RMSE_pure, round(uCI_laDiff, 3), round(lAsCI2[2], 5), round(slopeCI2[2], 5),
      round(uAsCI2[2], 5), "-", round(AsDiffCI2[2], 5), round(potAll2[3], 2), round(potAllU2[3], 2)
    )
  )
  return(res_tab)
}

#' ANOVA unrestricted model
#'
#' ANOVA with unrestricted model.
#'
#'
#' @param ro_new Data frame with the concentrations and the respective readouts.
#' @returns A data-frame with the analysis of variance.
#' @export
#' @examples
#'
#' ro_new <- data.frame(
#'   REF1 = c(1547, 1620, 1644, 2504, 3426, 3512, 3401, 3787), REF2 = c(1492, 1536, 1384, 2286, 3046, 3479, 3516, 3497),
#'   REF3 = c(1468, 1827, 1558, 2252, 3002, 3349, 2945, 3665),
#'   TEST1 = c(1405, 1523, 1502, 1474, 2383, 3221, 3589, 3445), TEST2 = c(1420, 1516, 1544, 1512, 2226, 3219, 3327, 3591),
#'   TEST3 = c(1399, 1376, 1588, 1475, 2148, 3083, 2942, 3466), log_dose = c(5.01, 3.401, 2.708, 2.015, 1.32176, 0.62861, -0.0645385, -1.6739764)
#' )
#'
#'
#' ANOVA4plUnresfunc(ro_new)
ANOVA4plUnresfunc <- function(ro_new) {
  all_l <- melt(data.frame(ro_new), id.vars = "log_dose", variable.name = "replname", value.name = "readout")
  all_len <- nrow(all_l)
  isRef <- rep(c(1, 0), 1, each = all_len / 2)
  isSample <- rep(c(0, 1), 1, each = all_len / 2)
  all_l$isRef <- isRef
  all_l$isSample <- isSample
  all_l$Conc <- exp(all_l$log_dose)
  all_l$readout[all_l$readout < 0] <- 0.01

  FITs <- Fitting_FUNC(ro_new = ro_new, TransFlag = FALSE)
  smr <- FITs[[1]]
  smu <- FITs[[2]]
  smrPREDs <- FITs[[5]]
  smuPREDs <- FITs[[6]]
  # pot <- drm(readout ~ Conc, isSample, data=all_l, fct=LL.4(names=c("b","d","a","c")),
  #            pmodels=data.frame(1,1,1,isSample))
  # potU <- drm(readout ~ Conc, isSample, data=all_l, fct=LL.4(names=c("b","d","a","c")),
  #             pmodels=data.frame(isSample, isSample,isSample,isSample))
  SStreat <- round(sum((smuPREDs - mean(all_l$readout))^2), 5)
  SStreat_df <- length(unique(all_l$log_dose)) - 1
  SSregr <- round(sum((smrPREDs - mean(all_l$readout))^2), 5)
  ## Non-parallel
  SSnonparallel <- round(sum(smr$residuals^2) - sum(smu$residuals^2), 5)
  ## Preparation
  SSprep <- round(sum((predict(lm(readout ~ isSample, all_l)) - mean(all_l$readout))^2), 5)
  ## Resid Err
  RSS <- round(sum(smu$residuals^2), 5)
  RSS_df <- nrow(all_l) - SStreat_df - 1
  FitAnova <- anova(lm(readout ~ factor(Conc) * isSample, all_l))
  # PureErr
  SSE <- round(FitAnova[4, 3], 5)
  SSE_df <- FitAnova[4, 1]
  # Non-Linearity
  SSnonlin <- round(sum((predict(lm(readout ~ factor(Conc) * isSample, all_l)) - smuPREDs)^2), 4)
  LoF_df <- FitAnova[1, 1] + FitAnova[2, 1]
  ## Total
  SStot <- round(sum((all_l$readout - mean(all_l$readout))^2), 5)
  MSE <- RSS / RSS_df
  noConc <- length(unique(all_l$Conc))
  AnovaDFs <- c(noConc - 1, 1, 3, noConc - 4 - 1, nrow(all_l) - noConc, noConc, nrow(all_l) - noConc - noConc, nrow(all_l) - 1)
  p_SStreat <- round(pf((SStreat / AnovaDFs[1]) / MSE, AnovaDFs[1], RSS_df, lower.tail = F), 3)
  p_SSprep <- round(pf((SSprep / AnovaDFs[2]) / MSE, AnovaDFs[2], RSS_df, lower.tail = F), 3)
  p_SSregr <- round(pf((SSregr / AnovaDFs[3]) / MSE, AnovaDFs[3], RSS_df, lower.tail = F), 3)
  p_SSnonp <- round(pf((SSnonparallel / AnovaDFs[4]) / MSE, AnovaDFs[3], RSS_df, lower.tail = F), 3)
  p_SSLoF <- round(pf((SSnonlin / LoF_df) / (SSE / SSE_df), LoF_df, SSE_df, lower.tail = F), 5)

  ANOVAtab <- data.frame(
    Source = c(
      "Treatment", "Preparation", "Regression",
      "Non-Parallelism", "Residual Error", "Non-linearity",
      "Pure Error", "Total"
    ),
    DF = AnovaDFs,
    SumSquares = c(
      SStreat, SSprep, SSregr, SSnonparallel,
      RSS, SSnonlin, SSE, SStot
    ),
    MeanSquares = c(
      round(SStreat / AnovaDFs[1], 3), SSprep, round(SStreat / AnovaDFs[3], 3), round(SSnonparallel / AnovaDFs[4], 3),
      round(MSE, 5), round(SSnonlin / LoF_df, 5), round(SSE / SSE_df, 5), ""
    ),
    "F-value" = c(
      round((SStreat / AnovaDFs[1]) / MSE, 5), round((SSprep / AnovaDFs[2]) / MSE, 5),
      round((SSregr / AnovaDFs[3]) / MSE, 5), round((SSnonparallel / AnovaDFs[4]) / MSE, 5),
      "", round((SSnonlin / LoF_df) / (SSE / SSE_df), 5), "", ""
    ),
    "p_value" = c(round(p_SStreat, 3), p_SSprep, round(p_SSregr, 3), p_SSnonp, "", p_SSLoF, "", "")
  )

  return(ANOVAtab)
}

#' Calculates descriptive statistics per concentration
#'
#' Mean, SD, and CV% per concentration and product..
#'
#'
#' @param ro_new Data frame with the concentrations and the respective readouts.
#' @returns A data-frame with the concentrations and metadata.
#' @export
#' @examples
#'
#' ro_new <- data.frame(
#'   REF1 = c(1547, 1620, 1644, 2504, 3426, 3512, 3401, 3787), REF2 = c(1492, 1536, 1384, 2286, 3046, 3479, 3516, 3497),
#'   REF3 = c(1468, 1827, 1558, 2252, 3002, 3349, 2945, 3665),
#'   TEST1 = c(1405, 1523, 1502, 1474, 2383, 3221, 3589, 3445), TEST2 = c(1420, 1516, 1544, 1512, 2226, 3219, 3327, 3591),
#'   TEST3 = c(1399, 1376, 1588, 1475, 2148, 3083, 2942, 3466), log_dose = c(5.01, 3.401, 2.708, 2.015, 1.32176, 0.62861, -0.0645385, -1.6739764)
#' )
#'
#'
#' perConcTab(ro_new, noDilSeries = 3)
perConcTab <- function(ro_new, noDilSeries) {
  Reftab <- ro_new[, c(1:noDilSeries)]
  Testtab <- ro_new[, c((noDilSeries + 1):(2 * noDilSeries))]
  tReftab <- t(Reftab)
  colnames(tReftab) <- round(ro_new[, ncol(ro_new)], 5)

  avs <- apply(tReftab, 2, mean)
  sds <- apply(tReftab, 2, sd)
  cv <- sds / avs * 100
  tReftab2 <- rbind(tReftab, avs, sds, cv)

  tTesttab <- t(Testtab)
  colnames(tTesttab) <- round(ro_new[, ncol(ro_new)], 5)

  avs_test <- apply(tTesttab, 2, mean)
  sds_test <- apply(tTesttab, 2, sd)
  cv_test <- sds_test / avs_test * 100
  tTesttab2 <- rbind(tTesttab, avs_test, sds_test, cv_test)
  concTab <- rbind(tReftab2, tTesttab2)
  return(concTab)
}

#' Calculates dilution series.
#'
#' Recursive function to calcuate a geometric dilution series.
#'
#' @param x Starting concentration.
#' @param Div Divisor to get next concentration.
#' @param N Number of dilution step, increases.
#' @param res Vector of results.
#' @param noDil Final number of concentrations -1 (the initial one).
#' @returns A data-frame with the concentrations and metadata.
#' @export
#' @examples
#'
#' x <- 1
#' Div <- 3
#' N <- 0
#' res <- c()
#' noDil <- 7
#'
#' divFUN(x, Div, N, res, noDil)
divFUN <- function(x, Div, N, res, noDil) {
  N <- N + 1
  y <- x / Div
  res <- c(res, y)
  if (N == noDil) {
    return(res)
  }
  divFUN(y, Div, N, res, noDil)
}