library(shiny)
library(shinydashboard)
library(shinyjs)
library(shinyAce)
library(shinydashboard)
library(shinycssloaders)
library(shinyBS)
library(purrr)
library(gslnls)
library(tidyverse)
library(ggplot2)
library(reshape2)
library(openxlsx)
library(DT)
library(ggpubr)
library(gridExtra)
library(drc)
library(twopartm)
library(car)
library(dplyr)


#### reactive values ----
Dat <- reactiveValues()
REP <- reactiveValues()

#### functions ----
dilFUN2 <- function(cs_,dils,Faktor) {
  av <- cs_
  dils_av <- dils_av
  dils_avsc <- dils_av*Faktor
  dils2 <- dils_avsc+av
  dilfactors <- 1/exp(dils2-lag(dils2))
  return(dilfactors)
}

plot_f <- function(dat, sigmoid,det_sig, TransFlag=F) {
  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()  
  if(is.null(det_sig)) {
    
      startlist <- list(a=sigmoid[1], b=sigmoid[5],cs=sigmoid[7],
                        d=sigmoid[3],r=sigmoid[8])
   
  } else {
    startlist <- list(a=det_sig[5], b=det_sig[1],cs=det_sig[7],
                      d=det_sig[3],r=det_sig[7] - det_sig[8])
  }
  
  mr <- gsl_nls(fn = readout ~ a+(d-a)/(1+exp(b*((cs-r*isSample)-log_dose))),
                data=all_l,
                start=startlist,
                control=gsl_nls_control(xtol=1e-6,ftol=1e-6, gtol=1e-6))
  s_mr <- tryCatch({
    s_mr <- summary(mr)
  },
  error = function(err) {
    s_mr <- NULL
  })
  
  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)))
  
  if (is.null(det_sig)) {
    SAMPLEtrue <- sigmoid[2] + (sigmoid[4] -sigmoid[2])/(1+exp(sigmoid[6]*((sigmoid[7]-sigmoid[8]-seq_x))))
    REFtrue <- sigmoid[1] + (sigmoid[3] -sigmoid[1])/(1+exp(sigmoid[5]*((sigmoid[7]-seq_x))))
  } else {
    SAMPLEtrue <- det_sig[4] + (det_sig[6] -det_sig[4])/(1+exp(det_sig[2]*(det_sig[8]-seq_x)))
    REFtrue <- det_sig[3] + (det_sig[5] -det_sig[3])/(1+exp(det_sig[1]*(det_sig[7]-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.31696/b)
  Xbendu3 <- cs+(1.31696/b)
  XbendlT <- cs-r-(1.31696/b)
  XbenduT <- cs-r+(1.31696/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
  
  p <- ggplot(all_l2, aes(x=log_dose, y=readout, color=factor(isRef))) +
    geom_point(shape=factor(isRef), alpha=0.8) +
    labs(title = paste("restricted 4pl"),
         color="product") +
    scale_color_manual(labels=c("test","reference"), values=c("#C2173F","#4545BA")) +
    scale_shape_manual(labels=c("test","reference")) +
    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="#4545BA55",linetype=2) +
    geom_vline(xintercept=c(XasymplT, XasympuT), col="#C2173F55",linetype=2) +
    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), alpha=0.8) +
    labs(title = paste("restricted transformed 4pl"), color="product") +
    scale_color_manual(labels=c("test","reference"), values=c("#C2173F","#4545BA")) +
    theme_bw()
  
  mrt <- gsl_nls(fn = readouttrans ~ 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_mrt <- summary(mrt)
  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.31696/b_trans)
  XbenduTrans <- cs_trans+(1.31696/b_trans)
  XbendlTransT <- cs_trans-r_trans-(1.31696/b_trans)
  XbenduTransT <- cs_trans-r_trans+(1.31696/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()
    if(is.null(det_sig)) {
      
      startlistmu <- list(as=sigmoid[1], bs=sigmoid[5],cs=sigmoid[7],
                        ds=sigmoid[3], at=sigmoid[2], bt=sigmoid[6],
                        dt=sigmoid[4],r=sigmoid[8])
      
    } else {
      startlistmu <- list(as=det_sig[5], bs=det_sig[1],cs=det_sig[7],
                        ds=det_sig[3],at=det_sig[6], bt=det_sig[2],
                        dt=det_sig[4],r=det_sig[7] - det_sig[8])
    }
    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) })
    #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() +
    labs(title="unrestricted 4_pl-Model", color="product") +
    scale_color_manual(labels = c("test","reference"), values=c("#C2173F","#4545BA")) +
    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() +
    labs(title="unrestricted transformed 4_pl-Model", color="product") +
    scale_color_manual(labels = c("test","reference"), values=c("#C2173F","#4545BA")) +
    theme_bw()
  
  unrestr_trans <- drm(readouttrans ~ exp(log_dose), isSample, data=all_l2, fct=LL.4(),
                       pmodels=data.frame(isSample, isSample,isSample,isSample))
  Sum_ut <- summary(unrestr_trans)
  ast_t <- Sum_ut$coefficients[3,1]
  ate_t <- Sum_ut$coefficients[4,1]
  bst_t <- Sum_ut$coefficients[1,1]
  bte_t <- Sum_ut$coefficients[2,1]
  cst_t <- log(Sum_ut$coefficients[7,1])
  cte_t <- log(Sum_ut$coefficients[8,1])
  dst_t <- Sum_ut$coefficients[5,1]
  dte_t <- Sum_ut$coefficients[6,1]
  
  REFu_trans <- ast_t + (dst_t-ast_t)/(1+exp(bst_t*(seq_x-cst_t)))
  SAMPLEu_trans <- ate_t + (dte_t-ate_t)/(1+exp(bte_t*(seq_x-cte_t)))
  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)
}

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)
}

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

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
  
  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)
}

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)
  
  CORdat <- cor(ro_new[,1],ro_new[,ncol(ro_new)])
  if (CORdat<0) SLOPE <- -1 else SLOPE <- 1
  
  startlist <- list(a=min(ro_new[,2]), b=SLOPE, d=max(ro_new[,2]), cs=mean(all_l$log_dose),r=0)
  tryCatch({
    mr <- gsl_nls(fn = readout ~ a+(d-a)/(1+exp(b*(cs-r*isSample-log_dose))),
                  data=all_l,
                  start=startlist,
                  control=gsl_nls_control(xtol=1e-6,ftol=1e-6, gtol=1e-6))
  },
  error = function(err) {
    err$message
  })
  
  startlistmu <- list(as=max(ro_new[,2]), bs=SLOPE, ds=min(ro_new[,2]), cs=mean(all_l$log_dose),
                      at=max(ro_new[,2]), bt=SLOPE, dt=min(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(err) {
    err$message
  })
  
  if (!PureErrFlag) {
    pot_est <- exp(confintd(mr, "r", method="asymptotic"))
    potU_est <- exp(confintd(mu, "r", method="asymptotic"))
  } else {
    FitAnova <- anova(lm(readout ~ factor(Conc)*isSample, all_l))
    meanPureErr <- FitAnova[4,3]
    SU_mr <- tryCatch({
      SU_mr <- summary(mr)
    }, error = function(msg) {
      return()
    })
    
    browser() 
    if (length(SU_mr)>1) {
      s_mr <- SU_mr$coefficients
    } else { SU_mr <- rep(NA,5) }
    
    VCOV <- vcov(mr)
    V_V <- VCOV/SU_mr$sigma^2
    SEsPure <- sqrt(diag(V_V)*meanPureErr)
    pot_est <- c(exp(s_mr['r',1]),exp(s_mr['r',1]-qt(0.975,nrow(all_l)-5)*SEsPure['r']),
                 exp(s_mr['r',1]+qt(0.975,nrow(all_l)-5)*SEsPure['r']))
    # unrestricted
    s_mu <- summary(mu)$coefficients
    SU_mu <- summary(mu)
    VCOVu <- vcov(mu)
    V_Vu <- VCOVu/SU_mu$sigma^2
    SEsPureU <- sqrt(diag(V_Vu)*meanPureErr)
    potU_est <- c(exp(s_mu['r',1]),exp(s_mu['r',1]-qt(0.975,nrow(all_l)-8)*SEsPureU['r']),
                  +   exp(s_mu['r',1]+qt(0.975,nrow(all_l)-8)*SEsPureU['r']))
  } # PureErrFlag
  
  startlistmr_log <- list(a=max(all_l$readouttrans), b=SLOPE, d=min(all_l$readouttrans), cs=mean(all_l$log_dose),r=0)
  
  tryCatch({
    mr_log <- gsl_nls(fn = readouttrans ~ a+(d-a)/(1+exp(b*(cs-r*isSample-log_dose))),
                      data=all_l,
                      start=startlistmr_log,
                      control=gsl_nls_control(xtol=1e-6,ftol=1e-6, gtol=1e-6))
  },
  error = function(err) {
    err$message
  })
  
  
  startlistmu_log <- list(as=max(ro_new[,2]), bs=SLOPE, ds=min(ro_new[,2]), cs=mean(all_l$log_dose),
                          at=max(ro_new[,2]), bt=SLOPE, dt=min(ro_new[,2]), r=0)
  
  tryCatch({
    mu_log <- gsl_nls(fn = readouttrans ~ 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_log,
                      control=gsl_nls_control(xtol=1e-6,ftol=1e-6, gtol=1e-6))
  },
  error = function(err) {
    err$message
  })
  
  pot_est_log <- exp(confintd(mr_log, "r", method="asymptotic"))
  potU_est_log <- exp(confintd(mu_log, "r", method="asymptotic"))
  colnames(pot_est_log) <- c("estimate","lowerCI2","upperCI")
  colnames(potU_est_log) <- c("estimate","lowerCI2","upperCI")
  #browser()   
  su_mr_log <- summary(mr_log)
  Dat$RMSE_Rlog <- su_mr_log$sigma
  su_mu_log <- 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(pot_est, potU_est, pot_est_log, potU_est_log)
  
  potALL2 <- cbind(c("restricted","unrestricted","ln-transformed restr","ln-transformed unrestr"), potALL)
  return(potALL2)
}

ParamCI_F <- function(xt,xs,se_xt, se_xs, CoVarlog,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(bs)
  var_log_xt <- (se_xt/xt)^2
  se_log_ratio <- sqrt(var_log_xs + var_log_xt) #-2*CoVarlog)
  
  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)
}

tests_FUNC <- function(ro_new, Lim, PureErrFlag) {
 #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_l$isRef <- isRef
  all_l$isSample <- isSample
  all_l$Conc <- exp(all_l$log_dose)
  all_l$readout[all_l$readout < 0] <- 0.01
  tryCatch({
  pot <- drm(readout ~ Conc, isSample, data=all_l, fct=LL.4(names=c("b","d","a","c")),
             pmodels=data.frame(1,1,1,isSample))
  },
  error = function(msg){
    return(0) })
  compParm(pot, "c",display=T)
  ED50 <- ED(pot,c(50), interval="delta")
  PotEst <- ED50[1,1]/ED50[2,1]
  potAll <- EDcomp(pot, percVec=c(50,50), interval="delta", display=FALSE)
  potAll2 <- potAll[1:3]
  
  CORro <- cor(ro_new[,1], ro_new[,ncol(ro_new)])
  
  if (CORro<0) SLOPE <- -1 else SLOPE <- 1
  startlist <- list(a=max(ro_new[,2]), b=SLOPE, d=min(ro_new[,2]), cs=mean(all_l$log_dose),r=0)
  tryCatch({
    mr <- gsl_nls(fn = readout ~ a+(d-a)/(1+exp(b*(cs-r*isSample-log_dose))),
                  data=all_l,
                  start=startlist,
                  control=gsl_nls_control(xtol=1e-6,ftol=1e-6, gtol=1e-6))
  },
  error = function(msg){
    return(0) })
  
  startlistmu <- list(as=max(ro_new[,2]), bs=SLOPE, ds=min(ro_new[,2]), cs=mean(all_l$log_dose),
                      at=max(ro_new[,2]), bt=SLOPE, dt=min(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) })
  
  smu <- tryCatch({ summary(mu) },
                  error=function(msg){
                    return(0) })
  
  POTr_CI <- potAll2[2:3]
  
  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 <- vcovMU/smu$sigma^2
  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
  
  potU <- drm(readout ~ Conc, isSample, data=all_l, fct=LL.4(names=c("b","d","a","c")),
              pmodels=data.frame(isSample, isSample,isSample,isSample))
  potAllU <- EDcomp(potU, percVec=c(50,50), interval="delta", display=FALSE)
  potAllU2 <- potAllU[1:3]
  sum_potU <- summary(potU)
  coeffs <- potU$coefficients
  coeffs[1] <- ifelse(CORro<0, -coeffs[1], coeffs[1])
  coeffs[2] <- ifelse(CORro<0, -coeffs[2], coeffs[2])
  names(coeffs) <- c("bs","bt","ds","dt","as","at","cs","ct")
  
  e_c_ref <- coeffs[7]
  e_c_test <- coeffs[8]
  coeffs[7:8] <- log(coeffs[7:8])
  test_c <- logical()
  if((potAllU2[2] >Lim[[9]]/100 & potAllU2[3] <Lim[[10]]/100)) test_c <- 0 else test_c <- 1
  
  #### 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((predict(potU)-mean(all_l$readout))^2),5)
  SSregr <- round(sum((predict(pot)-mean(all_l$readout))^2),5)
  # non-parallelism
  SSnonparall <- round(sum(resid(pot)^2)-sum(resid(potU)^2),5)
  SSprep <- round(sum((predict(lm(readout ~ isSample, all_l))-mean(all_l$readout))^2),5)
  
  RSS <- round(sum(potU$predres[,2]^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))-predict(potU))^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(pot$predres[,2]^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)
  #browser()  
  ## EQ test on lower As diff
  ds <- coeffs["ds"]
  dt <- coeffs["dt"]
  lAs_diff <- (dt-ds)
  uCI_laDiff <- lAs_diff+qt(0.975,df.residual(mu))*sqrt(sum_potU$coefficients[3,2]^2+sum_potU$coefficients[4,2]^2)
  lCI_laDiff <- lAs_diff-qt(0.975,df.residual(mu))*sqrt(sum_potU$coefficients[3,2]^2+sum_potU$coefficients[4,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)-noConc
  uAsCI2 <- ParamCI_F(dt,ds,se_dt, se_ds,CoVarlog_d, DFs, Conf=0.9975)
  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 <- slopeCI
  
  #### 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 <- (at-dt)/(as-ds)
  
  at_dt <- (at-dt)
  as_ds <- (as-ds)
  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(at_dt,as_ds,se_dt_at, se_ds_as,CoVarlog=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(at_dt/as_ds, 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)
}

ANOVA4plUnresfunc <- function(ro_new, sigmoid) {
  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
  
  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((potU$predres[,1] - mean(all_l$readout))^2),5)
  SStreat_df <- length(unique(all_l$log_dose))-1
  SSregr <- round(sum((predict(pot)-mean(all_l$readout))^2),5)
  ## Non-parallel
  SSnonparallel <- round(sum(resid(pot)^2) - sum(resid(potU)^2),5)
  ## Preparation
  SSprep <- round(sum((predict(lm(readout ~ isSample, all_l)) - mean(all_l$readout))^2),5)
  ## Resid Err
  RSS <- round(sum(potU$predres[,2]^2),5)
  RSS_df <- nrow(all_l)-SStreat_df-1
  FitAnova <- anova(lm(readout ~ factor(Conc)*isSample, all_l))
  # PureErr
  SSE <- FitAnova[4,3]
  SSE_df <- FitAnova[4,1]
  # Non-Linearity
  SSnonlin <- round(sum((predict(lm(readout ~ factor(Conc)*isSample, all_l)) - predict(potU))^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)
}

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)
  
}

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)
}



#### ui ----

ui <- dashboardPage(
  dashboardHeader(title = "Plateflow"),
  dashboardSidebar(
    sidebarMenu(
    img(src="logo.png", width=230),
    menuItem("Home", tabName="home", icon=icon("home")),
    menuItem("Data template", tabName = "template", icon=icon("table"),
             menuSubItem( tags$li(a("EXCEL File",  target="self",href="TestFile.xlsx")))
    ),
    # menuItem("User Manual /Validation", tabName = "manual", icon=icon("book"), # tabName here and in dashboard body need to be identical
    #          menuSubItem(icon = NULL, tags$li(a("Document",  target="self",href="UserManual.pdf")))
    # ),
    menuItem("EXCEL upload", tabName="Dataupload", icon=icon("magnet", lib="glyphicon")),
    menuItem("4PL + report", tabName="fourPL", icon=icon("chart-line", lib="font-awesome")),
    #menuItem("XLSX diagnostics", tabName="XLdiagn", icon=icon("chart-bar", lib="font-awesome")),
    menuItem("Linear regression + report", tabName="pla", icon=icon("pencil", lib="glyphicon")),
    menuItem("Wizard", tabName="wizard", icon=icon("chart-column", lib="font-awesome")),
    menuItem("Documentation", tabName="documentation", icon=icon("chart-area", lib="font-awesome"))
  ),
  tags$footer(HTML(paste(tags$strong(tags$u("InnerAnalytics")), paste(rep("&nbsp",9), collapse=""),
                         "Developer:", paste(rep("&nbsp",9), collapse=""),
                         "Host on:", paste(rep("&nbsp",9), collapse=""),
                         "Version:")),
              align="left", style=
                "position:fixed; bottom:0;width=100%; background: #FFC337BB; font-family: Times New Roman; font-size:100%;
                  padding: 5px; color:#4545BA; box-sizing: border-box; z-index: 1000;")),
  dashboardBody(
    fluidPage(
      tabItems(
        tabItem(tabName = "home", htmlOutput("homePage")),
        tabItem(tabName = "Dataupload", uiOutput("Dataupload")),
        tabItem(tabName = "fourPL", uiOutput("fourPL")),
        #tabItem(tabName = "XLdiagn", uiOutput("XLdiagn")),
        tabItem(tabName = "pla", uiOutput("pla")),
        tabItem(tabName = "wizard", uiOutput("wizard")),
        tabItem(tabName = "documentation", uiOutput("docu"))
      )
      
    )
  ), skin="blue"
)


#### server ----
server <- function(input, output, session) {
  ReportParS <- reactiveValues()
  IPReportParS <- reactiveValues()
  
 #### renderUIs ---- 
  output$homePage <- renderUI({
    navbarPage("Home",
               tabPanel("Introduction",
                        tags$img(src="logo.png", class="adv_logo"),
                        h4("Introduction to the bioassay software"),
                        tags$mark("linear regression"), br(),
                        column(4,
                               tags$table(id="dose-table",
                                          numericInput("lEACdiffla","lower EAC for diff. of LA", -0.175, step=0.001),
                                          numericInput("uEACdiffla","upper EAC for diff. of LA", 0.189, step=0.001),
                                          numericInput("lEACratiola","lower EAC ratio of LAs", 0.005, step=0.001),
                                          numericInput("uEACratiola","upper EAC for ratio of LAs", 100, step=1),
                                          
                                          numericInput("lEACratioSlope","lower EAC for ratio of slopes", 0.55, step=0.01),
                                          numericInput("uEACratioSlope","upper EAC for ratio of slopes", 1.84, step=0.1),
                                          numericInput("lEACratioua","lower EAC for ratio of UAs", 0.75, step=0.1),
                                          numericInput("uEACratioua","upper EAC for ratio of UAs", 1.33, step=0.1),
                                          numericInput("lowerPot","lower EAC for potency", 75, step=1),
                                          numericInput("upperPot","upper EAC for potency", 133, step=1),
                                          numericInput("lEACratioAdiff","lower EAC of ratio of asymptote differences", 0.75, step=0.01),
                                          numericInput("uEACratioAdiff","upper EAC of ratio of asymptote differences", 1.33, step=0.01)
                               ))
                        
               ),
               tabPanel("Documentation",
                        h4("Introduction "),
                        h4("Information on dilution setting"),
                        "(for details see:", a(href="ADONIS.pdf","Whitepaper", download=NA, target="_blank"),")",tags$br(),
                        "Bend points are calculated according to following formula:",
                        withMathJax(" $$bp_{1/2} = \\pm\\frac{1.31696}{Hill's slope}$$")),
               tabPanel("Configuration",
                        verbatimTextOutput("sessioninfo"))
    )
  })
  
  output$Dataupload <- renderUI({
    navbarPage(title="Information",
               tabPanel(title = "Real data",
                        tabsetPanel(
                          tabPanel("Data input",
                            
                        
                        column(3,
                               #img(src="Screenshot.png", width=200),
                              box(title = "Upload", status="warning",solidHeader = T, width=12, "Please upload your EXCEL file here",
                                  fileInput("iFile",'',accept=".xlsx")),
                              uiOutput(outputId = "sheetName"),
                              "For data format in the EXCEL file see Data template",
                              "If no data are uploaded, the settings to the right are used for calculations.",
                              
                              tags$head(tags$style(HTML("label {font-size:80%;margin-bottom: 3px;margin-top: 3px;}"))),
                               
                               div(checkboxInput("PureErr", "Should pure error be used for calculation of CIs?", FALSE),
                                   style = "font-size: 24px !important;color: #C2173F"),
                               checkboxGroupInput("selectedSSTs", "Which suitability tests to be used?", choices= c("F-test on Regr."="1",
                                                   "EQ-test on lower asymptote difference"= "2",
                                                   "EQ-test on ratio of lower asymptote"= "3","EQ-test on ratio of Hill slopes"= "4",
                                                   "EQ-test on ratio of upper asymptote"= "5", "F-test on non-linearity"="6",
                                                  "EQ-test on ratio of asymptote differences"= "7"),
                                                  selected= c("1","2","3","4","5","6","7")),
                              #actionLink("selectall","SelectAll"),
                               h5("\n\n\n Author: Franz Innerbichler, InnerAnalytics")),
                        
                          ),
                        tabPanel("4pl-Analysis",
                                 tags$style(HTML("pre { color: black; background-color: #FFE1FF; 
                                             font-family: 'Helvetica Neue', Helvetica, Arial, sans-serif;font-size: 12px;} ")),
                                
                                   wellPanel(
                                     fluidRow(
                                       column(4, 
                                              #h5("Diagnostics only shown if EXCEL is uploaded"),
                                              htmlOutput("PureErrW2"), 
                                              tags$head(tags$style("#PureErrW2{color: red;
                                                             font-size: 16px;
                                                             font_style: italic;}")),
                                              tableOutput("Filesampl"),
                                              tableOutput("relpotTestTab"),
                                              plotOutput("relpotTestPlot", width="300px", height="150px"), # Pot CI plot
                                              h4("Akaike Information Criterion"),
                                              tableOutput("AIC"),
                                              h5("First row: restricted model; 2nd row: unrestricted model"),
                                              h5("Smaller values of AIC indicate better fit to the data"),
                                              tableOutput("VarDiagn")
                                       ),
                                       column(8,
                                              plotOutput("XLplot"),
                                              
                                              plotOutput("diagnplot"),
                                              
                                              tableOutput("ANOVAXLS"),
                                              DTOutput("EQtests"))
                                       
                                     ))),
                        tabPanel("linear Analysis",
                                 sidebarLayout(
                                   sidebarPanel(
                                     width=3,
                                     fluidRow(
                                       column(6,
                                              numericInput("Limits",p("limit to be >", bsButton("q4",label="", icon=icon("info"), style="primary", size="extra-small")),
                                                           bsPopover(id="q4", title="", content="The calculated limits ...")))
                                     )),
                                   mainPanel(
                                     tabsetPanel(id="tabs",
                                                 tabPanel("linear PLA",
                                                          box(title="ANOVA table", status="primary",solidHeader = T, width=12,
                                                              tableOutput("Anovatab")),
                                                          column(6,
                                                                 htmlOutput("PureErrW3"),
                                                                 tags$head(tags$style("#PureErrW3{color: red;
                                                             font-size: 16px;
                                                             font_style: italic;}")),
                                                                 plotOutput("plotLin"),
                                                                 "Delta method is used for potency CIs",
                                                                 DT::dataTableOutput("pottab"),
                                                                 h4("Unrestricted linear model (SSSI):"),
                                                                 tableOutput("SummaryModABu"),
                                                                 h4("Restricted linear model (CSSI):"),
                                                                 tableOutput("SummaryModAB")),
                                                          column(3,
                                                                 h3("Tests for linear PLA):"),
                                                                 DT::dataTableOutput("TESTSlin"),
                                                                 h5("The estimate is the p-value of the test"),
                                                                 h5("F-tests on regression, significance of slopes, and preparation need to have a p-value <0.05 to pass"),
                                                                 h5("All other tests pass if p-value > 0.05"),
                                                                 "SST CI for difference of slopes:",
                                                                 tableOutput("SlopeDiffCI")),
                                                          column(3,
                                                                 h3("ANOVA for parallel line assay"),
                                                                 DT::dataTableOutput("ANOVAlin"))),
                                                 tabPanel("Report",
                                                          h4("Settings for report")
                                                 ))
                                   )
                                 )
                                 ),
                                   tabPanel("parameter estimates",
                                            column(3,style = "background: #4FCBD922",
                                                   br(),
                                                   h4("Regression results restricted"),
                                                   tableOutput("coeffs_r"),
                                                   "Bend points restricted",
                                                   tableOutput("bends_r2")),
                                            column(3,style = "background: #B5C74022",       
                                                   br(),
                                                   h4("Regression results unrestricted"),
                                                   tableOutput("coeffs_unr")),
                                            column(3,style = "background: #F9545422", 
                                                   h4("Regression results (ln-transformed)"),
                                                   tableOutput("logcoeffs_r"),
                                                   tableOutput("bends_unr2"),
                                                   
                                                   tableOutput("logcoeffs_unr"))
                                   ),
                                   tabPanel("Report",
                                            h4("Settings for report"),
                                            downloadButton("downloadXLReport", label="Download PDF report", class="butt"),
                                            tags$style(type="text/css","#downloadXLReport {background-color: orange; color: black;font-family: COurier New}"),
                                   ) 
                        )
                        )
               )
  })
  
  
  
  
  output$fourPL <- renderUI({
    navbarPage(title="4PL",
               tabPanel("Analysis and Plots",
                        #sidebarLayout(
                          # sidebarPanel(
                          #   width=4,
                          #   fluidRow(
                          #            )
                          #   ),
                          mainPanel(width=12,
                            tabsetPanel(id="tabs",
                                        tabPanel("Settings",
                                                 h4("Settings of 4PL regression"),
                                                 div(checkboxInput("PureErr4pl", "Should pure error be used for calculation of CIs?", FALSE),
                                                     style = "font-size: 24px !important;color: #C2173F"),
                                                 
                                                 h4("User help"),
                                                 h5("If new dilutions are entered, please adjust EC50 to avoid calculation errors"),
                                                 # numericInput("Limits",p("limit to be >", bsButton("q4",label="", icon=icon("info"), style="primary", size="extra-small")),
                                                 #              bsPopover(id="q4", title="", content="The calculated limits ...")),
                                                 #h5("Diagnostics only shown if EXCEL is uploaded"),
                                                 column(2,style = "background: #7FAEFF",
                                                        #actionButton("StartCalc", "Click, when calculations to be started"),
                                                        h4("curve settings"),
                                                        numericInput("lowAsymptREF", "lower asymptote REF",10,step=1,min=0),
                                                        numericInput("lowAsymptTEST", "lower asymptote TEST",10,step=1,min=0),
                                                        numericInput("uppAsymptREF", "upper asymptote REF",110,step=1,min=0),
                                                        numericInput("uppAsymptTEST", "upper asymptote TEST",110,step=1,min=0)
                                                 ),
                                                 column(2,style = "background: #7FAEFF",
                                                        numericInput("slopeREF", "slope REF",1,step=0.1,min=-10),
                                                        numericInput("slopeTEST", "slope TEST",1,step=0.1,min=-10),
                                                        numericInput("EC50", "EC50 REF",-3.5),
                                                        numericInput("potDiff", "potency difference",0)
                                                 ),
                                                 column(2,style = "background: #627ADD",
                                                        h4("dilutions"),
                                                        numericInput("CONC1", "highest concentration",0.3, min=-3.5),
                                                        numericInput("CONC2", "2nd concentration",0.15),
                                                        numericInput("CONC3", "3rd concentration",0.075),
                                                        numericInput("CONC4", "4th concentration",0.0375),
                                                        numericInput("CONC5", "5th concentration",0.01875),
                                                        
                                                        numericInput("CONC6", "6th concentration",0.00938)
                                                 ),
                                                 column(2,style = "background: #627ADD",
                                                        
                                                        numericInput("CONC7", "7th concentration",0.00469),
                                                        numericInput("CONC8", "8thd concentration",0.00235),
                                                        numericInput("CONC9", "9thd concentration",value=NA),
                                                        numericInput("CONC10", "10th concentration",value=NA),
                                                        numericInput("CONC11", "11th concentration",value=NA),
                                                        
                                                        numericInput("CONC12", "lowest concentration",NA)
                                                 ),
                                                 column(2,style = "background: #4FCBD9",
                                                        h4("geometric dilution scheme"),
                                                        numericInput("ConcStart", "starting concentration",value=NA, min=0),
                                                        numericInput("dilutionFac", "dilution factor",value=NA, min=0, max=10),
                                                        numericInput("NoDil", "no. of dilutions",value=NA, min=8),
                                                        numericInput("NoDilSer", "no. of dil. series",value=NA),
                                                        verbatimTextOutput("dilutions")
                                                 ),

                                                     column(2,
                                                            h3("Settings"),
                                                            helpText("Vary the slider to see the effect of special cause"),
                                                            sliderInput("sdfac","Variability as % of lower to upper asymptote", max=10, value=3, min=0.1, step=0.1),
                                                            checkboxInput("heterosked","heteroskedastic noise", FALSE),
                                                            sliderInput("potencydiff","potency of test (%)", max=200, min=50, value=100, step=1),
                                                            sliderInput("outlL","outlier in lower asymptote",  min=0, max=1.5, value=0, step=0.1),
                                                            sliderInput("outlM","outlier in mid part",   min=0, max=1.5,value=0, step=0.1),
                                                            sliderInput("outlU","outlier in upper asymptote",   min=0, max=1.5,value=0, step=0.1)
                                                     ),
                                                 h4("log-dilutions from settings above"),
                                                 verbatimTextOutput("logdil")
                                                     
                                                 #)         
                                        ),
                                        #### 4pl fits ----
                                        tabPanel("4pl-fit",
                                                 tags$head(tags$style("#PureErrW2{color: red;
                                                             font-size: 10px;
                                                             font_style: italic;}")),
                                                 wellPanel(
                                                   fluidRow(
                                                     column(10,
                                                            tags$style(span(htmlOutput("PureErrW3"), style="color: red")),
                                                            htmlOutput("PureErrW3"),
                                                            
                                                            plotOutput("plot", width = "80%"),
                                                            DT::dataTableOutput("pottab4pl"),
                                                            "Footnote: test performed on relative CIs.",
                                                     
                                                            DTOutput("EQtests4pl"), # SSTs
                                                            h5("*...The estimate for F-test on regression and on non-linearity is the p-value"),
                                                            h5("F-test on regression passes if F-value > F-crit and thus p < 0.05"),
                                                            h5("F-test on non-linearity passes if F-value < F-crit and thus p > 0.05"),
                                                            h5("Test results outcome: 0 ... test passed (for EQ tests: CI within limits);
                                           1 ... test failed (for EQ tests CI not within limits);
                                           -1 ... calculations unbound/denominator too close to 0"),
                                                            #plotOutput("CIplot, height=50%")
                                                     ),
                                                     column(8,
                                                            "4 PL ANOVA unrestricted",
                                                            box(title = "ANOVA unrestricted", status="warning",solidHeader = T, width=12, "",
                                                                DT::dataTableOutput("ANOVA")),
                                                            h5("Please note that for the CIs of rel. potency the RSS of the restricted model is used"),
                                                            h5("RSS ... 'Residual error' in the SumSquares column"),
                                                            h5("MSE ... 'Residual error' in the MeanSquaress column"),
                                                            h5("SSE ... 'Pure error' in the SumSquares column"),
                                                            h5("RMSE ... Square root of the 'Residual Error' in the MeanSquares column"),
                                                            verbatimTextOutput("RMSE")
                                                     ),
                                                     column(8,
                                                            box(title = "Simulated data per log-concentration", status="warning",solidHeader = T, width=12, "incl. mean, sd and CV%",
                                                                DT::dataTableOutput("Conctab")))                  
                                                   ))
                                        ),
                                        
                                        tabPanel("ln-transformed y",
                                                 h4("ln-transformed y-axis plots"),
                                                 plotOutput("plot4plTrans", width = "80%"),
                                                 DT::dataTableOutput("pottab4plTrans"),
                                        ),
                                        tabPanel("Report",
                                                 h4("Settings for report"),
                                                 downloadButton("downloadXLReport", label="Download PDF report", class="butt"),
                                                 tags$style(type="text/css","#downloadXLReport {background-color: orange; color: black;font-family: COurier New}"),
                                        ) 
                            )
                          )
                        #)
               )
    )
  })
  
  
  output$pla <- renderUI({
    navbarPage(title="pla",
               tabPanel("Analysis and Plots",
                        
                        )
               )
  })
  
  output$wizard <- renderUI({
    navbarPage(title="Dilution setting",
               tabPanel("Plots",
                        sidebarLayout(
                          sidebarPanel(
                            width=3,
                            fluidRow(
                              column(6,
                                     numericInput("Limits",p("limit to be >", bsButton("q4",label="", icon=icon("info"), style="primary", size="extra-small")),
                                                  bsPopover(id="q4", title="", content="The calculated limits ...")))
                            )),
                          mainPanel(
                            tabsetPanel(id="tabs",
                                        tabPanel("4pl",
                                                 box(title="ANOVA table", status="primary",solidHeader = T, width=12,
                                                     tableOutput("Anovatab")),
                                                 column(4,
                                                        h3("Confidence intervals"),
                                                        tableOutput("CIs"),
                                                        "The confidence intrval table is interaactive for changes in: variability slider (%SD), potency of test-slider,
                                       and 'Adjust the dilutions'-slider",
                                                        tableOutput("optimalDils"),
                                                        selectInput(inputId="scenario", label= "Select an 'optimal' scenario:", choices = c("scenario 2","scenario 3","scenario 6","steep slope"))),
                                                 column(5,
                                                        plotOutput("plotfordilutions"),
                                                        h4("in grey: most extreme bend point lines of theoretical samples with 50% and 200% potency"),
                                                        sliderInput("dilslider", "Adjust the dilutions(+-change in %)", min = -100,max=100, value=0, step=1, round=0),
                                                        checkboxInput("fixupper","Fix highest concentration (if unticked, the center is fixed)",FALSE),
                                                        h5("Dilution factors"),
                                                        tableOutput("adjlogdil"),
                                                        "Short guidance: wider dilution ranges increase the CIs of rel. potency, and decrease the CIs of upper and lower asymptote ratios, as well as Hill's slope ratios", br(),
                                                        "Narrower dilution ranges decrease the CIs of rel. potency, and increase the CIs of upper and lower asymptote ratios, ands Hill's slope ratios",),
                                                 column(3,
                                                        h3("Bend points"),
                                                        tableOutput("bps"),
                                                        tableOutput("extremebps"),
                                                        h4("Explanation of the plots")
                                                 )),
                                        tabPanel("Report",
                                                 h4("Settings for report")
                                        ))
                          )
                        )))
  })
  

  
  v <- reactiveValues(num_dose=0, next.dose.t=0)
  
  sigmoid <- reactive({
    sig <- c(input$lowAsymptREF, input$lowAsymptTEST,input$uppAsymptREF,input$uppAsymptREF,
             input$slopeREF,input$slopeTEST,input$EC50,input$potDiff)
    sig
  })
  
  CONC <- reactive({
    Konz_ <- c(input$CONC1,input$CONC2,input$CONC3,input$CONC4,
               input$CONC5,input$CONC6,input$CONC7,input$CONC8,
               input$CONC9,input$CONC10,input$CONC11,input$CONC12)
    if (any(na.omit(Konz_)==0)) Konz_[Konz_ ==0] <- 0.0000001
    Konz <- na.omit(Konz_)
  })
  
  Dils <- reactive({
    Dilutions <- c(input$ConcStart,input$dilutionFac,input$NoDil,input$NoDilSer)
  })
  
  #### input EXCEL file ----
  
  observe({
    if (!is.null(input$iFile)) {
      inFile <- input$iFile
      ext <- tools::file_ext(inFile$name)
      file.rename(inFile$datapath, paste(inFile$datapath, ".xlsx",sep=""))
      t.filelocation <- gsub('\\\\','/',paste(inFile$datapath, ext,sep="."))
      sheets <- openxlsx::getSheetNames(t.filelocation)
      dat <- lapply(sheets, openxlsx::read.xlsx, xlsxFile = t.filelocation)
      names(dat) <- sheets
      Dat$wb <- dat
      names(Dat$wb) <- sheets
      Dat$sheets <- sheets
      Dat$FileName <- input$iFile[["name"]]
      
    }
  })
  output$sheetName <- renderUI({
    if (!is.null(Dat$wb)) {
      #browser()
      cnSheets <- Dat$sheets
      cnSheets2 <- c("please choose", cnSheets)
      selectInput(inputId = "sheet", label="Select a sheet:",choices = cnSheets)
    }
  })
  observeEvent(input$sign_out, {
    unlink(input$iFile$datapath)
    reset(id = "") # from shinyjs package
  })
  
  
  #### process XLSX file ----
  observe({
    if (!is.null(input$iFile)) {
      if (!is.null(input$sheet)) {
        if (input$sheet != "please choose") {
          #browser()       
          
          XLdat <- Dat$wb[input$sheet][[1]]
          if (is.null(XLdat)) XLdat <- Dat$wb[Dat$sheets[1]][[1]]
          cn <- colnames(XLdat)
          logI <- grep("log", cn)
          logDoseI <- grep("log_dose", cn)
          if (length(logI)>0 & length(logDoseI)==0) {
            XLdat$log_dose <- XLdat[,logI]
            XLdat2 <- XLdat[,-logI]
            CORro <- cor(XLdat$log_dose, XLdat[,3])
          } else if (length(logI)==0 & length(logDoseI)==0) {
            Ind <- grep("dilu|dose|Dose|Conc|conc",cn)
            XLdat$log_dose <- log(XLdat[,Ind])
            CORro <- cor(XLdat[,Ind], XLdat[,3])
            XLdat2 <- XLdat[,-Ind]
          } else if (length(logI)>0 & length(logDoseI)>0) {
            XLdat2 <- XLdat
            CORro <- cor(XLdat[,logI], XLdat[,3])
          }
          Dat$EXCEL <- XLdat2
          PureErrFlag <- input$PureErr
          warning_text2 <- reactive({
            ifelse(PureErrFlag, 'Pure Error is selected', '')
          })
          output$PureErrW2 <- renderText(warning_text2())
          
          noDilSeries <-(ncol(XLdat2)-1)/2
          noDils <- nrow(XLdat2)
          Dat$noDilSeriesXL <- noDilSeries
          
          all_l <- melt(data.frame(XLdat2), id.vars="log_dose",variable.name = "replname", value.name = "readout")
          isRef <- rep(c(1,0),1,each=nrow(XLdat2)*noDilSeries) 
          isSample <- rep(c(0,1),1,each=nrow(XLdat2)*noDilSeries)
          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
          REP$all_l <- all_l
          
          #### XLSX eval ----
          if (CORro<0) SLOPE <- -1 else SLOPE <- 1
          ec50est <- (max(all_l$log_dose)+min(all_l$log_dose))/2
          startlist <- list(a=min(all_l$readout), b=SLOPE, d=max(all_l$readout), cs=ec50est,r=0)
          tryCatch({
            mr <- gsl_nls(fn = readout ~ a+(d-a)/(1+exp(b*(cs-r*isSample-log_dose))),
                          data=all_l,
                          start=startlist,
                          control=gsl_nls_control(xtol=1e-6,ftol=1e-6, gtol=1e-6))
          },
          error = function(err) {
            err$message
          })
          
          startlistmu <- list(as=min(all_l$readout), bs=SLOPE, ds=max(all_l$readout), cs=ec50est,
                              at=min(all_l$readout), bt=SLOPE, dt=max(all_l$readout), 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(err) {
            err$message
          })
          
          Smr <- summary(mr)
          Smu <- summary(mu)
          coeffsMR <- Smr$coefficients[,1]
          coeffsMU <- Smu$coefficients[,1]
          Dat$coeffsMRes <- coeffsMR
          names(coeffsMU) <- c("lowAsym REF", "slope REF","upperAsym REF","EC50 REF","lowAsym TEST","slope TEST","upperAsym TEST","r")
          
          if (!PureErrFlag) {
            pot_est <- exp(confintd(mr, "r", method="asymptotic"))
            potU_est <- exp(confintd(mu, "r", method="asymptotic"))
            colnames(pot_est) <- c("estimate","lowerCI","upperCI")
            colnames(potU_est) <- c("estimate","lowerCI","upperCI")
          } else {
            FitAnova <- anova(lm(readout ~ factor(log_dose)*isSample, all_l))
            meanPureErr <- FitAnova[4,3]
            DFsPure <- FitAnova[4,1]
            VCOV <- vcov(mr)
            V_V <- VCOV/Smr$sigma^2
            VCOVpure <- V_V*meanPureErr
            SEsPure <- sqrt(diag(V_V)*meanPureErr)
            pot_est <- data.frame(estimate=exp(coeffsMR[5]), lowerCI = exp(coeffsMR[5]-qt(0.975,DFsPure)*SEsPure[5]),
                                  upperCI = exp(coeffsMR[5]+qt(0.975,DFsPure)*SEsPure[5]))
            VCOVu <- vcov(mu)
            V_Vu <- VCOVu/Smu$sigma^2
            #VCOVpure <- V_Vu*meanPureErr
            SEsPureU <- sqrt(diag(V_Vu)*meanPureErr)
            potU_est <- data.frame(estimate=exp(coeffsMU[7]), lowerCI = exp(coeffsMU[7]-qt(0.975,DFsPure)*SEsPureU[7]),
                                   upperCI = exp(coeffsMU[7]+qt(0.975,DFsPure)*SEsPureU[7]))
            
          }
          
          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))
          SR <- summary(pot)
          SU <- summary(potU)
          coeffs_UN <- potU$coefficients
          coeffs_UN[1] <- ifelse(xor(CORro>0, coeffs_UN[1]>0), -coeffs_UN[1],coeffs_UN[1])
          coeffs_UN[2] <- ifelse(xor(CORro>0, coeffs_UN[2]>0), -coeffs_UN[2],coeffs_UN[2])
          coeffs_UN[7:8] <- log(coeffs_UN[7:8])
          POTU <- EDcomp(potU, percVec = c(50,50), interval="delta",display=F)
          Dat$potDiffXL <- POTU[1]*100
          RMSE_unr_diagn <- sqrt(SU$resVar)
          RMSE_res_diagn <- sqrt(SR$resVar)
          up_lowDiffDiagn <- SU$coefficients[5,1] - SU$coefficients[3,1]
          ProzSD_diagn <- RMSE_unr_diagn*100/up_lowDiffDiagn
          Dat$ProzSD_XL <- ProzSD_diagn
          
          observe({
            pot_est3 <- data.frame(pot_est*100)
            MaxPl <- max(input$upperPot, pot_est3$upperCI)
            MinPl <- min(input$lowerPot, pot_est3$lowerCI)
            MaxPl_ <- MaxPl*1.2
            MinPl_ <- MinPl*0.8
        #browser()    
            p_relCI <- ggplot(data=pot_est3, aes(xmin=lowerCI, xmax=upperCI, y=1)) +
              geom_linerange(size=4, col="darkseagreen",alpha=0.5) +
              geom_point(aes(x=estimate, y=1), col="grey15", shape=13, size=10) +
              geom_vline(xintercept = c(input$lowerPot, input$upperPot), col="indianred") +
              annotate("text", x=input$lowerPot-13, y=1.040, label=paste("lower EAC:", input$lowerPot), col="indianred") +
              annotate("text", x=input$upperPot+13, y=1.040, label=paste("upper EAC:", input$upperPot), col="indianred") +
              annotate("text", x=pot_est3$lowerCI-10, y=1.020, label=paste("lower CL:", round(pot_est3$lowerCI,1)), col="darkgreen") +
              annotate("text", x=pot_est3$upperCI+10, y=1.020, label=paste("upper CL:", round(pot_est3$upperCI,1)), col="darkgreen") +
              annotate("text", x=pot_est3$estimate, y=0.98, label=paste("rel. potency:", round(pot_est3$estimate,1)), col="black") +
              ylim(c(0.95, 1.05)) +
              xlim(c(MinPl_,MaxPl_)) +
              xlab("relative potency + confidence interval") +
              theme_bw() +
              theme(axis.title.y=element_blank(),
                    axis.text.y=element_blank(),
                    axis.ticks.y=element_blank())
            output$relpotTestPlot <- renderPlot({
              p_relCI
            })
            REP$relpotTestPlot <- p_relCI
            
            output$relpotTestTab <- renderTable({ pot_est3 })
          })
          
          SStreat <- round(sum((potU$predres[,1] - mean(all_l$readout))^2),5)
          SStreat_df <- length(unique(all_l$log_dose))-1
          SSregr <- round(sum((predict(pot)-mean(all_l$readout))^2),5)
          ## Non-parallel
          SSnonparallel <- round(sum(resid(pot)^2) - sum(resid(potU)^2),5)
          ## Preparation
          SSprep <- round(sum((predict(lm(readout ~ isSample, all_l)) - mean(all_l$readout))^2),5)
          ## Resid Err
          RSS <- round(sum(potU$predres[,2]^2),5)
          RSS_df <- nrow(all_l)-SStreat_df-1
          FitAnova <- anova(lm(readout ~ factor(Conc)*isSample, all_l))
          # PureErr
          SSE <- FitAnova[4,3]
          SSE_df <- FitAnova[4,1]
          # Non-Linearity
          SSnonlin <- round(sum((predict(lm(readout ~ factor(Conc)*isSample, all_l)) - predict(potU))^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)
          
          ANOVAtab2 <- data.frame(Source = c("Treatment","Preparation","Regression",
                                             "Non-Parallelism","Residual Error","Non-linearity",
                                             "Pure Error","Total"),
                                  DF = round(AnovaDFs,0),
                                  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,"","")
          )
          
          output$ANOVAXLS <- renderTable({ ANOVAtab2 })
          
          REP$ANOVAXLS <- ANOVAtab2
          
          #browser()   
          
          startlistlog <- list(a=min(log(all_l$readout)), b=SLOPE, d=max(log(all_l$readout)), cs=ec50est,r=0)
          tryCatch({
            mrlog <- gsl_nls(fn = log(readout) ~ a+(d-a)/(1+exp(b*(cs-r*isSample-log_dose))),
                             data=all_l,
                             start=startlistlog,
                             control=gsl_nls_control(xtol=1e-6,ftol=1e-6, gtol=1e-6))
          },
          error = function(err) {
            print("Error in mrlog gsl_nls")
          })
          
          startlistmulog <- list(as=min(log(all_l$readout)), bs=SLOPE, ds=max(log(all_l$readout)), cs=ec50est,
                                 at=min(log(all_l$readout)), bt=SLOPE, dt=max(log(all_l$readout)), r=0)
          tryCatch({
            mulog <- 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=startlistmulog,
                             control=gsl_nls_control(xtol=1e-6,ftol=1e-6, gtol=1e-6))
          },
          error = function(err) {
            print("Error in murlog gsl_nls")
          })
          logpot <- drm(log(readout) ~ Conc, isSample, data=all_l, fct=LL.4(names=c("b","d","a","c")),
                        pmodels=data.frame(1,1,1,isSample))
          logpotU <- drm(log(readout) ~ Conc, isSample, data=all_l, fct=LL.4(names=c("b","d","a","c")),
                         pmodels=data.frame(isSample, isSample,isSample,isSample))
          Smrlog <- summary(mrlog)$coefficients
          Smulog <- summary(mulog)$coefficients
          SUlog <- summary(logpotU)
          SRlog <- summary(logpot)
          RMSE_unrlog_diagn <- sqrt(SUlog$resVar)
          RMSE_reslog_diagn <- sqrt(SRlog$resVar)
          
          up_lowDifflogDiagn <- SUlog$coefficients[5, 1] - SUlog$coefficients[3, 1]
          ProzSDlog_diagn    <- RMSE_unrlog_diagn * 100 / up_lowDifflogDiagn
          
          #### Diagnostic RMSE table ####
          DiagnTable <- data.frame(parameter = c("RMSE unrestricted", "RMSE_restr.", "Diff_upper-lowerAsymp", "%SD (unrestricted)", 
                                                 "RMSE log_unrestricted", "RMSE log_restr", "diff_up-lowAsymp_log", "%SD (log unrestricted)"),
                                   result = c(round(RMSE_unr_diagn, 4), round(RMSE_res_diagn, 4), 
                                              round(up_lowDiffDiagn, 4), round(ProzSD_diagn, 4),
                                              round(RMSE_unrlog_diagn, 4), round(RMSE_reslog_diagn, 4),
                                              round(up_lowDifflogDiagn, 4), round(ProzSDlog_diagn, 4)))
          
          Dat$DiagnTable <- DiagnTable
          REP$DiagnTable <- DiagnTable
          
          logpotest  <- exp(confintd(mrlog, "r", method = "asymptotic")) # compParm(logpot, "c")
          logpotuest <- exp(confintd(mulog, "r", method = "asymptotic")) # compParm(logpotu, "c")
          
          # Berechnung der Konfidenzintervalle (CI)
          # logpotCI <- c(exp(Smrlog[5,1] - qt(0.975, nrow(all_1)-5) * Smrlog[5,2]), exp(Smrlog[5,1]), exp(Smrlog[5,1] + qt(0.975, nrow(all_1)-5) * Smrlog[5,2]))
          colnames(logpotest)  <- c("estimate", "lowerCI", "upperCI")
          
          colnames(logpotuest) <- c("estimate", "lowerCI", "upperCI")
      #browser()    
          cnXL <- colnames(XLdat2)
          Filesample <- data.frame(Test = c("File name", "samples"), Test2=c(Dat$FileName, paste(cnXL[1], " vs ", cnXL[4])))
          colnames(Filesample) <- c("", "")
          output$Filesampl <- renderTable({ Filesample }, rownames = F)
          
          UnRPLAausw <- data.frame(Information = c("model", "lower asymptote Ref", "Hill's slope Ref", "upper asymptote Ref","EC50 Ref", 
                                                   "lower asymptote Test", "Hill's slope Test", 
                                                   "upper asymptote Test","EC50 Difference",
                                                   "relative potency", "lower CI", "upper CI"),
                                   Results = unlist(c("UNRESTRICTED", round(coeffsMU, 3), round(potU_est*100, 3)))) # von psl_nls
          
          # "log relative potency", "log lower CI", "log upper CI", round(logpotest, 3), round(compParm(potu, "c", display = F), 3)
          
          output$coeffs_unr <- renderTable({
            UnRPLAausw
          })
          
      #browser()
          UnRPLAausw2 <- data.frame(Dat$bendpointsTRANS)
          if (length(UnRPLAausw2) > 0) {
            colnames(UnRPLAausw2) <- c("bendpoints log")
            UnrBendLog <- data.frame(Bendpoint = c("REF_lower","REF_upper",
                                                   "TEST_lower","REF_lower"),
                                     bendpoints_logscale = UnRPLAausw2)
            
            output$bends_unr2 <- renderTable({
              UnrBendLog
            })
          }
          REP$UnRPLAausw  <- UnRPLAausw
          REP$UnRPLAausw2 <- UnRPLAausw2
          
          # browser()
          coeffs_R  <- coeffsMR # pot$coefficients
          coeffs_R[5] <- coeffs_R[4] - coeffs_R[5]
          names(coeffs_R) <- c("lower A", "slope",  "upper A", "EC50 REF", "EC50 TEST")
          # coeffs_R[4] <- log(coeffs_R[4])
          # coeffs_R[5] <- log(coeffs_R[5])
          # --- Ergebnistabelle: RESTRICTED (Eingeschränktes Modell) ---
          PLAAusw <- data.frame(
            Information = c("model", "lower asymptote", "Hill's slope", "upper asymptote","EC50 Ref", 
                            "EC50 Test", "relative potency", 
                            "lower CI", "upper CI"),
            Results = unlist(c("RESTRICTED", round(coeffs_R, 3), 
                               round(pot_est[1, ] * 100, 3)))) # von gs1_nls
          output$coeffs_r <- renderTable({ PLAAusw  })
          
          PLAAusw2 <- data.frame(Dat$bendpoints)
          output$bends_r2 <- renderTable({ PLAAusw2 }, digits = 3, rownames = T)
          
          REP$PLAausw <- PLAAusw
          REP$PLBend <- PLAAusw2
          
          # --- Koeffizienten-Extraktion ---
          logcoeffs_R <- Smrlog[, 1] # logpot$coefficients
          names(logcoeffs_R) <- c("lower A", "Hill's slope",  "upper A", "EC50 REF","EC50 DIFF")
          
          # --- Ergebnistabelle: LOG RESTRICTED ---
          
          LogPLAAusw <- data.frame(
            Information = c("model", "lower asymptote", "Hill's slope", "upper asymptote","EC50 Ref", 
                            "EC50 difference", "log relative potency", 
                            "log lower CI", "log upper CI"),
            Results = unlist(c("LOG RESTRICTED", round(logcoeffs_R, 3), 
                               round(logpotest * 100, 3)))) # von gs1_nls
          
          output$logcoeffs_r <- renderTable({ LogPLAAusw  })
          REP$LogPLAausw <- LogPLAAusw
          
          logcoeffs_UNR <- Smulog[,1]
          names(logcoeffs_UNR) <- c("lower asymptote Ref", "Hill's slope Ref",  "upper asymptote Ref","EC50 Ref",
                                    "lower asymptote Test", "Hill's slope Test", "upper asymptote Test","EC50 Diff" 
          )
          
          # --- Ergebnistabelle: LOG UNRESTRICTED ---
          
          LogUnrPLAAusw <- data.frame(
            Information = c("model", "lower asymptote Ref", "Hill's slope Ref",  "upper asymptote Ref","EC50 Ref",
                            "lower asymptote Test", "Hill's slope Test", "upper asymptote Test","EC50 Diff" ,
                            "relative potency", "lower CI", "upper CI"),
            
            Results = unlist(c("LOG UNRESTRICTED", round(logcoeffs_UNR, 3), 
                               round(logpotest * 100, 3)))) # von gs1_nls
          
          output$logcoeffs_unr <- renderTable({
            LogUnrPLAAusw
          })
          REP$LogUnrPLAausw <- LogUnrPLAAusw
          #browser()        
          Dat$coeffs_UN <- coeffs_UN
          
          
          if (exists("Ind")) {
            Dat$dilution <- XLdat[,Ind]
          } else Dat$dilution <- exp(XLdat[,logI])
          # --- Plot-Ausgabe ---
          output$XLplot <- renderPlot({
            plot_f(XLdat2, sigmoid = NULL, det_sig = coeffs_UN, TransFlag=F)
          })
          
          REP$XLdat2 <- XLdat2
          
          # --- Diagnose-Plots (Residualanalyse) ---
          output$diagnplot <- renderPlot({
            op <- par(mfrow = c(2, 2), mar = c(3.2, 3.2, 2, .5), mgp = c(2, .7, 0)) 
            
            # 1. Residuals vs Fitted
            plot(residuals(pot) ~ fitted(pot), main = "Residuals restricted")
            abline(h = 0)
            qqnorm(residuals(pot))
            qqline(residuals(pot))
            plot(residuals(potU) ~ fitted(potU), main = "Residuals unrestricted")
            abline(h = 0)
            
            qqnorm(residuals(potU))
            qqline(residuals(potU))
            
            par(op) # Parameter zurücksetzen
          })
          
          output$AIC <- renderTable({
            AIC <- AIC(pot, potU)
          })
          
          output$VarDiagn <- renderTable({
            DiagnTable
          }, digits=4)
          
          output$relpotplot <- renderPlot({
            relpot(potU, intervall="fieller", bty="l",
                   main="Quality of rel. potency over response")
          })
          
          
          
        } # !please choose
      } # input$sheet
    } # input$iFile
  })
  
  #### make geomDils reactive ----
  observe({
    #browser()
    if (is.null(input$ConcStart)) return(NULL)
    if (!is.na(input$ConcStart) & !is.na(input$dilutionFac) &!is.na(input$NoDil) &!is.na(input$NoDilSer)) {
      upR <- input$ConcStart
      noDil <- input$NoDil
      noDilSer <- input$NoDilSer
      Div <- input$dilutionFac
      res <- c()
      N_ <- 1
      Conc <- c(upR, divFUN(upR,Div,N=N_,res,noDil))
      
      Dat$MetaConc <- Conc
    } 
  })
  
  #### updateSlider on XLSX ----
  observe({
    if (!is.null(Dat$potDiffXL)) {
      updateSliderInput(session, "potencydiff",
                        value=round(as.numeric(Dat$potDiffXL[[1]]),5))
    }
  })
  observeEvent(input$potencydiff, {
    if (!is.null(Dat$potDiffXL)) {
      updateSliderInput(session, "potencydiff",
                        value=round(as.numeric(input$potencydiff),5))
    }
  })
  observe({
    if (!is.null(Dat$ProzSD_XL)) {
      updateSliderInput(session, "sdfacf",
                        value=round(as.numeric(Dat$ProzSD_XL[[1]]),5))
    }
  })
  observeEvent(input$sdfac, {
    if (!is.null(Dat$ProzSD_XL)) {
      updateSliderInput(session, "sdfac",
                        value=round(as.numeric(Dat$ProzSD_XL[[1]]),5))
    }
  })
  
  #### updaterNumeric Input ----
  observe({
    if(!is.null(Dat$coeffs_UN)) {
      updateNumericInput(session, "lowAsymptREF",
                         value=round(as.numeric(Dat$coeffs_UN[3]),5), min=0)
      updateNumericInput(session, "lowAsymptTEST",
                         value=round(as.numeric(Dat$coeffs_UN[4]),5), min=0)
      updateNumericInput(session, "uppAsymptREF",
                         value=round(as.numeric(Dat$coeffs_UN[5]),5), min=0)
      updateNumericInput(session, "uppAsymptTEST",
                         value=round(as.numeric(Dat$coeffs_UN[6]),5), min=0)
      updateNumericInput(session, "slopeREF",
                         value=round(as.numeric(Dat$coeffs_UN[1]),5))
      updateNumericInput(session, "slopeTEST",
                         value=round(as.numeric(Dat$coeffs_UN[2]),5))
      updateNumericInput(session, "EC50",
                         value=round(as.numeric(Dat$coeffs_UN[7]),5))
      updateNumericInput(session, "potDiff",
                         value=round(as.numeric(Dat$coeffs_UN[7])- as.numeric(Dat$coeffs_UN[8]),5))
    }
  })
  
  observe({
    if(!is.null(Dat$dilution)) {
      updateNumericInput(session, "CONC1",
                         value=as.numeric(Dat$dilution[1]))
      updateNumericInput(session, "CONC2",
                         value=as.numeric(Dat$dilution[2]))
      updateNumericInput(session, "CONC3",
                         value=as.numeric(Dat$dilution[3]))
      updateNumericInput(session, "CONC4",
                         value=as.numeric(Dat$dilution[4]))
      updateNumericInput(session, "CONC5",
                         value=as.numeric(Dat$dilution[5]))
      updateNumericInput(session, "CONC6",
                         value=as.numeric(Dat$dilution[6]))
      updateNumericInput(session, "CONC7",
                         value=as.numeric(Dat$dilution[7]))
      updateNumericInput(session, "CONC8",
                         value=as.numeric(Dat$dilution[8]))
      updateNumericInput(session, "CONC9",
                         value=as.numeric(Dat$dilution[9]))
      updateNumericInput(session, "CONC10",
                         value=as.numeric(Dat$dilution[10]))
      updateNumericInput(session, "CONC11",
                         value=as.numeric(Dat$dilution[11]))
      updateNumericInput(session, "CONC12",
                         value=as.numeric(Dat$dilution[12]))
      
    }
  })
  
  observe({
    if(!is.null(Dat$MetaConc)) {
      updateNumericInput(session, "CONC1",
                         value=as.numeric(Dat$MetaConc[1]))
      updateNumericInput(session, "CONC2",
                         value=as.numeric(Dat$MetaConc[2]))
      updateNumericInput(session, "CONC3",
                         value=as.numeric(Dat$MetaConc[3]))
      updateNumericInput(session, "CONC4",
                         value=as.numeric(Dat$MetaConc[4]))
      updateNumericInput(session, "CONC5",
                         value=as.numeric(Dat$MetaConc[5]))
      updateNumericInput(session, "CONC6",
                         value=as.numeric(Dat$MetaConc[6]))
      updateNumericInput(session, "CONC7",
                         value=as.numeric(Dat$MetaConc[7]))
      updateNumericInput(session, "CONC8",
                         value=as.numeric(Dat$MetaConc[8]))
      updateNumericInput(session, "CONC9",
                         value=as.numeric(Dat$MetaConc[9]))
      updateNumericInput(session, "CONC10",
                         value=as.numeric(Dat$MetaConc[10]))
      updateNumericInput(session, "CONC11",
                         value=as.numeric(Dat$MetaConc[11]))
      updateNumericInput(session, "CONC12",
                         value=as.numeric(Dat$MetaConc[12]))
      
    }
  })
  
  #### render logDilsText ----
  output$logdil <- renderText({
    if (!is.null(Dat$MetaConc)) {
      Conc <- Dat$MetaConc
    } else Conc <- CONC()
    logdilu <-log(Conc)
    logdilu
  })
  
  
  #### reactive dataset sim ----
  
  sim <- reactive({
    #browser()
    if(is.null(sigmoid())) return(NULL)
    sd_fac_ <- as.numeric(input$sdfac)
    r_ <- log(as.numeric(input$potencydiff)/100)
    as = sigmoid()[1]; bs = sigmoid()[5];cs = sigmoid()[7];ds = sigmoid()[3];at = sigmoid()[2];
    bt = sigmoid()[6];r = sigmoid()[8]; ct = cs-r_; dt = sigmoid()[4];
    if (!is.null(Dat$MetaConc)) Conc <- Dat$MetaConc else Conc <- CONC()
    log_conc <- log(Conc)
    av_test <- as + (ds-as)/(1+exp(bs*(cs - log_conc)))
    av_ref <- at + (dt-at)/(1+exp(bt*(ct - log_conc)))
    #browser()  
    if (!is.na(input$NoDilSer)) {
      noDilSer <- input$NoDilSer
    } else if (!is.null(Dat$noDilSeriesXL)) noDilSer <- Dat$noDilSeriesXL else noDilSer <- 3
    if (!is.na(input$NoDil)) noDil <- input$NoDil else noDil <- length(log_conc)
    isRef <- rep(c(1,0), 1,each=noDilSer*noDil)
    isSample <- rep(c(0,1), 1,each=noDilSer*noDil)
    
    #if (is.null(Dat$EXCEL)) {
      ro_new <- Calc_DilRes(as=as,at=at,ds=ds,dt=dt,cs=cs,ct=ct,r=r_,bt=bt,bs=bs, log_conc = log_conc,
                            sd_fac=sd_fac_,
                            # auslenkU=outlierU,
                            # auslenkM=outlierM,
                            # auslenkL=outlierL,
                            heteroNoise = input$heterosked, noDilSeries = noDilSer, noDils = noDil)
    #} else ro_new <- Dat$EXCEL
  })
  # })
  
  ####sim2 ----
  sim2 <- reactive({
    tab <- sim()
    if (is.null(Dat$EXCEL)) return(tab) else return(Dat$EXCEL)
  })
  
  
  #### Plot 4pl ----
  output$plot <- renderPlot({
    #browser()
    sigmoid <- sigmoid()
    det_sig=NULL
    plot_f(sim2(),sigmoid, det_sig, TransFlag = F)
  })
  
  #### Plot 4pl Transformed ----
  output$plot4plTrans <- renderPlot({
    #browser()
    sigmoid <- sigmoid()
    det_sig=NULL
    plot_f(sim2(),sigmoid, det_sig, TransFlag = T)
  })
  
  
  #### Testergebnisse für 4PL ----
  observe({
    if (is.null(sim2())) return(NULL)
    if (is.null(input$PureErr4pl)) return(NULL)
    #observeEvent(input$StartCalc,{  
    PureErrFlag <- input$PureErr4pl
    warning_text3 <- reactive({
      ifelse(PureErrFlag, 'Pure error selected','')
    })
    
   #browser() 
    output$PureErrW3 <- renderText(warning_text3())
    
    Limite <- list(as.numeric(input$lEACdiffla), as.numeric(input$uEACdiffla), 
                   as.numeric(input$lEACratiola), as.numeric(input$uEACratiola), 
                   as.numeric(input$lEACratioSlope), as.numeric(input$uEACratioSlope), 
                   as.numeric(input$lEACratioua), as.numeric(input$uEACratioua), 
                   as.numeric(input$lowerPot), as.numeric(input$upperPot), 
                   as.numeric(input$lEACratioAdiff), as.numeric(input$uEACratioAdiff))
    Dat$limite <- Limite
    #browser()
    
    tab <- tests_FUNC(sim2(), Limite, PureErrFlag = PureErrFlag)
    if (length(tab)>1) {
      tab[1,6:7] <- c("-","-")
      Dat$tests_FUNC <- tab
      REP$testsTab <- tab
      tab2 <- tab[1:7,]
      
      dat <- datatable(tab2,options = list(
        paging=TRUE,
        dom="t",
        rownames=FALSE
      )) %>% formatStyle("test_results",
                         target='row',
                         backgroundColor = styleEqual(c(-1,0,1),
                                                      c("pink",'lightgreen','lightgrey')))
    } else { dat <- datatable(data.frame(test_results = "Convergeance failed for the uploaded dataset")) }
    #browser()
    output$EQtests4pl <- renderDT({ dat})
    
  }) # observe
  
  #### Testergebnisse für XLSX ----
  observe({
    if (is.null(Dat$EXCEL)) return(NULL)
    if (is.null(input$PureErr)) return(NULL)
    #observeEvent(input$StartCalc,{  
    PureErrFlag <- input$PureErr
    warning_text3 <- reactive({
      ifelse(PureErrFlag, 'Pure error selected','')
    })
    
    
    output$PureErrW3 <- renderText(warning_text3())
    
    Limite <- list(as.numeric(input$lEACdiffla), as.numeric(input$uEACdiffla), 
                   as.numeric(input$lEACratiola), as.numeric(input$uEACratiola), 
                   as.numeric(input$lEACratioSlope), as.numeric(input$uEACratioSlope), 
                   as.numeric(input$lEACratioua), as.numeric(input$uEACratioua), 
                   as.numeric(input$lowerPot), as.numeric(input$upperPot), 
                   as.numeric(input$lEACratioAdiff), as.numeric(input$uEACratioAdiff))
    Dat$limite <- Limite
    #browser()
    SelTests <- as.numeric(input$selectedSSTs)
    tab <- tests_FUNC(Dat$EXCEL, Limite, PureErrFlag = PureErrFlag)
    
    tab[1,6:7] <- c("-","-")
    Dat$tests_FUNC <- tab
    REP$testsTab <- tab
    tab2 <- tab[SelTests,]
    
    dat <- datatable(tab2,options = list(
      paging=TRUE,
      dom="t",
      rownames=FALSE
    )) %>% formatStyle("test_results",
                       target='row',
                       backgroundColor = styleEqual(c(-1,0,1),
                                                    c("pink",'lightgreen','lightgrey')))
    
    output$EQtests <- renderDT({ dat })
    
  }) # observe
  
  ####plot CIs ----
  observe({
    tab <- Dat$tests_FUNC
    if (is.null(tab)) return(NULL)
    
    tab2 <- tab[-c(1,2,3,6),]
    tab2[,3:7] <- apply(tab2[,3:7],2,as.numeric)
    tab2[4:5,3:7] <- tab2[4:5,3:7]/100
    
    p_CIs <- ggplot(tab2,aes(x=test,y=estimate, color=test,group=test)) +
      geom_point() +
      geom_errorbar(aes(ymin=lower_CI, ymax=upper_CI), width=0.4) +
      geom_crossbar(aes(ymin=lower_limit, ymax=upper_limit), size=0.8) + 
      coord_flip() + 
      theme_bw() +
      theme(legend.position = "none",text = element_text(size=20))
    
    output$CIplot <- renderPlot({ p_CIs}, height=200)
    
    REP$CIplot <- p_CIs
  })
  
  
  output$simdat <- DT::renderDataTable({
    tab <- sim2()
    if (is.character(tab)) stop(tab)
    
    tab2 <- round(tab, 5)
    colnames(tab2) <- c(paste("T", seq(1,(ncol(tab2)-1)/2)),
                        paste("R", seq(1,(ncol(tab2)-1)/2)), "log_conc" )
    dat <- datatable(tab2, options=list(
      paging=T,
      pageLength=20,
      dom="t"
    ))                    
  })
  
  output$Conctab <- DT::renderDataTable({
    if (!is.na(Dils()[1]) & is.na(Dils()[4])) return(NULL)
    tab <- sim2()
    if (is.character(tab)) stop(tab)
    if (!is.na(Dils()[4])) {
      noDilSer <- Dils()[4]
    } else if (!is.null(Dat$noDilSeriesXL)) {
      noDilSer <- Dat$noDilSeriesXL
    } else { noDilSer <- 3 }
    
    Conc <- CONC()
    Conctab <- perConcTab(tab, noDilSeries = noDilSer)
    Dat$Conctab <- Conctab
    
    dat <- datatable(Conctab, options=list(
      paging=T,
      pageLength=12,
      dom="t"
    ))    %>% formatStyle(0,
                          target='row',
                          backgroundColor = styleEqual(c("avs","sds","cv", "avs_test","sds_test","cv_test"),
                                                       c('lightgrey','lightgreen','pink','lightgrey','lightgreen','pink'))
    )    %>% formatRound(columns=colnames(Conctab), digits=3)
  })
  
  #### linear Plot output ----
  
  output$plotLin <- renderPlot({
    tab <- Dat$EXCEL
    
    # tab <- sim2()
    # if (is.character(tab)) stop(tab)
    # #browser()  
    # if (!is.na(Dils()[4])) noDilSer <- Dils()[4] else noDilSer = (ncol(tab)-1)/2
    # Conc <- CONC()
    # Conctab <- Dat$Conctab
    # if (!is.na(Dils()[3])) noDil <- Dils()[3] else noDil = length(Conc)
    # 
    # slopeSt <- slopeTe <- matrix(NA, nrow=noDil-2,ncol=2)
    # for (i in 1:(noDil-2)) {
    #   avs <- Conctab[noDilSer+1,]
    #   threes <- data.frame(lnC=log(Conc[i:(i+2)]), resp=avs[i:(i+2)])
    #   lm3St <- lm(resp ~ lnC, data=threes)
    #   slopeSt[i,] <- lm3St$coefficients
    #   avt <- Conctab[noDilSer*2+4,]
    #   threet <- data.frame(lnC=log(Conc[i:(i+2)]), resp=avt[i:(i+2)])
    #   lm3Te <- lm(resp ~ lnC, data=threet)
    #   slopeTe[i,] <- lm3Te$coefficients
    # }
    # 
    # indS <- which(abs(slopeSt[,2]) == max(abs(slopeSt[,2])))
    # indT <- which(abs(slopeTe[,2]) == max(abs(slopeTe[,2])))
    # 
    # pl_ <- slopeSt[indS,1]+slopeSt[indS,2]*log(Conc)
    # pl_T <- slopeTe[indT,1]+slopeTe[indT,2]*log(Conc)
    # pl_df <- data.frame(lnC=log(Conc), plotS=pl_, plotT=pl_T)
    # 
    # all_l <- melt(data.frame(tab), 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)
    # all_l2S <- all_l2[all_l2$isRef == 1,]
    # all_l2T <- all_l2[all_l2$isRef == 0,]
    # all_mS <- all_l2S[order(all_l2S$log_dose, decreasing=TRUE),]
    # all_mT <- all_l2T[order(all_l2T$log_dose, decreasing=TRUE),]
    # 
    # circleS <- all_mS[(indS*noDilSer-(noDilSer-1)):((indS+2)*noDilSer),]
    # circleT <- all_mT[(indT*noDilSer-(noDilSer-1)):((indT+2)*noDilSer),]
    # circle <- rbind(circleS,circleT)
    # 
    # Dat$circles <- circle
    # sigmoid <- sigmoid()
    # log_dose <- unique(all_l$log_dose)
    # seq_x <- seq(min(log_dose), max(log_dose),0.1)
    # SAMPLEtrue <- sigmoid[2] + (sigmoid[4]-sigmoid[2])/(1+exp(sigmoid[6]*((sigmoid[7]-log(input$potencydiff/100)-seq_x))))
    # REFtrue <- sigmoid[1] + (sigmoid[3]-sigmoid[1])/(1+exp(sigmoid[5]*((sigmoid[7]-seq_x))))
    # 
    # truePL_df <- cbind(seq_x,SAMPLEtrue, REFtrue)
    # 
    # p <- ggplot(all_l2,aes(x=log_dose,y=readout, color=factor(isRef))) +
    #   geom_point() +
    #   labs(title=paste("linear regression model", indS,indT), color="product") +
    #   scale_colour_manual(labels = c("test","reference"), values=c("red","blue")) +
    #   ylim(min(all_l2$readout),max(all_l2$readout)) +
    #   theme_bw()
    # p2 <- p + geom_line(data=pl_df,aes(x=lnC,y=plotS),color="blue",
    #                     inherit.aes = F) +
    #   geom_line(data=pl_df,aes(x=lnC,y=plotT),color="red",
    #             inherit.aes = F) +
    #   geom_line(data=data.frame(truePL_df),aes(x=seq_x,y=SAMPLEtrue),color="red", linetype=2,alpha=0.4,
    #             inherit.aes = F) +
    #   geom_line(data=data.frame(truePL_df),aes(x=seq_x,y=REFtrue),color="blue", linetype=2,alpha=0.4,
    #             inherit.aes = F) +
    #   labs(title = paste("unrestricted linear regression model",indS,indT), color="product") +
    #   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), inherit.aes = FALSE) +
    #   scale_shape_manual(labels=c("test","reference"), values=c(21,21))
    # 
    # 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)
    # 
    # pl_restS <- sum_mLin$coefficients[1,1] + sum_mLin$coefficients[2,1]*log(Conc)
    # pl_restT <- sum_mLin$coefficients[1,1] + sum_mLin$coefficients[3,1] + sum_mLin$coefficients[2,1]*log(Conc)
    # pl_rest <- data.frame(lnC=log(Conc), plotS=pl_restS, plotT=pl_restT)
    # 
    # pr2 <- p + geom_line(data=pl_rest,aes(x=lnC,y=plotS),color="blue",
    #                      inherit.aes = F) +
    #   geom_line(data=pl_rest,aes(x=lnC,y=plotT),color="red",
    #             inherit.aes = F) +
    #   geom_line(data=data.frame(truePL_df),aes(x=seq_x,y=SAMPLEtrue),color="red", linetype=2,alpha=0.4,
    #             inherit.aes = F) +
    #   geom_line(data=data.frame(truePL_df),aes(x=seq_x,y=REFtrue),color="blue", linetype=2,alpha=0.4,
    #             inherit.aes = F) +
    #   labs(title = paste("restricted linear regression model",indS,indT), color="product") +
    #   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), inherit.aes = FALSE) +
    #   scale_shape_manual(labels=c("test","reference"), values=c(21,21))
    # grid.arrange(p3,pr3,nrow=1)
  })
  
  
  output$plotLin2 <- renderPlot({
    tab <- sim2()
    if (is.character(tab)) stop(tab)
  #browser()  
    if (!is.na(Dils()[4])) noDilSer <- Dils()[4] else noDilSer = (ncol(tab)-1)/2
    Conc <- CONC()
    Conctab <- Dat$Conctab
    if (!is.na(Dils()[3])) noDil <- Dils()[3] else noDil = length(Conc)
    
    slopeSt <- slopeTe <- matrix(NA, nrow=noDil-2,ncol=2)
    for (i in 1:(noDil-2)) {
      avs <- Conctab[noDilSer+1,]
      threes <- data.frame(lnC=log(Conc[i:(i+2)]), resp=avs[i:(i+2)])
      lm3St <- lm(resp ~ lnC, data=threes)
      slopeSt[i,] <- lm3St$coefficients
      avt <- Conctab[noDilSer*2+4,]
      threet <- data.frame(lnC=log(Conc[i:(i+2)]), resp=avt[i:(i+2)])
      lm3Te <- lm(resp ~ lnC, data=threet)
      slopeTe[i,] <- lm3Te$coefficients
    }
    
    indS <- which(abs(slopeSt[,2]) == max(abs(slopeSt[,2])))
    indT <- which(abs(slopeTe[,2]) == max(abs(slopeTe[,2])))
    
    pl_ <- slopeSt[indS,1]+slopeSt[indS,2]*log(Conc)
    pl_T <- slopeTe[indT,1]+slopeTe[indT,2]*log(Conc)
    pl_df <- data.frame(lnC=log(Conc), plotS=pl_, plotT=pl_T)
    
    all_l <- melt(data.frame(tab), 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)
    all_l2S <- all_l2[all_l2$isRef == 1,]
    all_l2T <- all_l2[all_l2$isRef == 0,]
    all_mS <- all_l2S[order(all_l2S$log_dose, decreasing=TRUE),]
    all_mT <- all_l2T[order(all_l2T$log_dose, decreasing=TRUE),]
    
    circleS <- all_mS[(indS*noDilSer-(noDilSer-1)):((indS+2)*noDilSer),]
    circleT <- all_mT[(indT*noDilSer-(noDilSer-1)):((indT+2)*noDilSer),]
    circle <- rbind(circleS,circleT)
    
    Dat$circles <- circle
    sigmoid <- sigmoid()
    log_dose <- unique(all_l$log_dose)
    seq_x <- seq(min(log_dose), max(log_dose),0.1)
    SAMPLEtrue <- sigmoid[2] + (sigmoid[4]-sigmoid[2])/(1+exp(sigmoid[6]*((sigmoid[7]-log(input$potencydiff/100)-seq_x))))
    REFtrue <- sigmoid[1] + (sigmoid[3]-sigmoid[1])/(1+exp(sigmoid[5]*((sigmoid[7]-seq_x))))
    
    truePL_df <- cbind(seq_x,SAMPLEtrue, REFtrue)
    
    p <- ggplot(all_l2,aes(x=log_dose,y=readout, color=factor(isRef))) +
      geom_point() +
      labs(title=paste("linear regression model", indS,indT), color="product") +
      scale_colour_manual(labels = c("test","reference"), values=c("red","blue")) +
      ylim(min(all_l2$readout),max(all_l2$readout)) +
      theme_bw()
    p2 <- p + geom_line(data=pl_df,aes(x=lnC,y=plotS),color="blue",
                        inherit.aes = F) +
      geom_line(data=pl_df,aes(x=lnC,y=plotT),color="red",
                inherit.aes = F) +
      geom_line(data=data.frame(truePL_df),aes(x=seq_x,y=SAMPLEtrue),color="red", linetype=2,alpha=0.4,
                inherit.aes = F) +
      geom_line(data=data.frame(truePL_df),aes(x=seq_x,y=REFtrue),color="blue", linetype=2,alpha=0.4,
                inherit.aes = F) +
      labs(title = paste("unrestricted linear regression model",indS,indT), color="product") +
      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), inherit.aes = FALSE) +
      scale_shape_manual(labels=c("test","reference"), values=c(21,21))
    
    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)
    
    pl_restS <- sum_mLin$coefficients[1,1] + sum_mLin$coefficients[2,1]*log(Conc)
    pl_restT <- sum_mLin$coefficients[1,1] + sum_mLin$coefficients[3,1] + sum_mLin$coefficients[2,1]*log(Conc)
    pl_rest <- data.frame(lnC=log(Conc), plotS=pl_restS, plotT=pl_restT)
    
    pr2 <- p + geom_line(data=pl_rest,aes(x=lnC,y=plotS),color="blue",
                         inherit.aes = F) +
      geom_line(data=pl_rest,aes(x=lnC,y=plotT),color="red",
                inherit.aes = F) +
      geom_line(data=data.frame(truePL_df),aes(x=seq_x,y=SAMPLEtrue),color="red", linetype=2,alpha=0.4,
                inherit.aes = F) +
      geom_line(data=data.frame(truePL_df),aes(x=seq_x,y=REFtrue),color="blue", linetype=2,alpha=0.4,
                inherit.aes = F) +
      labs(title = paste("restricted linear regression model",indS,indT), color="product") +
      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), inherit.aes = FALSE) +
      scale_shape_manual(labels=c("test","reference"), values=c(21,21))
    grid.arrange(p3,pr3,nrow=1)
  })
  
  #### linear PLA tests ----
  output$TESTSlin <- DT::renderDataTable({
    tab <- sim2()
    if (is.character(tab)) stop(tab)
    Conc <- CONC()
    Limite <- list(as.numeric(input$lEACdiffla), as.numeric(input$uEACdiffla), 
                   as.numeric(input$lEACratiola), as.numeric(input$uEACratiola), 
                   as.numeric(input$lEACratioSlope), as.numeric(input$uEACratioSlope), 
                   as.numeric(input$lEACratioua), as.numeric(input$uEACratioua), 
                   as.numeric(input$lowerPot), as.numeric(input$upperPot), 
                   as.numeric(input$lEACratioAdiff), as.numeric(input$uEACratioAdiff))
    
    circles <- Dat$circles
    PureErrFlag <- input$PureErr
    warning_text <- reactive({
      ifelse(PureErrFlag, 'Pure error is selected','')
    })
    output$PureErrW <- renderText(warning_text())
    
    LIN <- ANOVAlintests(tab,circles,Limite,PureErrFlag=PureErrFlag)
    df <- LIN[[1]]
    su_modU <- LIN[[2]]
    su_mod2 <- LIN[[4]]
    
    output$SummaryModABu <- renderTable({ su_modU }, digits=5)
    output$SummaryModAB <- renderTable({ su_mod2 }, digits=5)
    
    slopeDiffCI <- t(data.frame(LIN[[3]]))
    colnames(slopeDiffCI) <- c("slope difference","lower CI","upper CI")
    output$SlopeDiffCI <- renderTable({ slopeDiffCI },digits=5)
    #browser()
    Dat$ANOVA <- df[,4:length(df)]
    dat <- datatable(df[,1:3],
                     options=list(
                       paging=T, dom="t",rownames=F
                     )) %>% formatStyle("test_results", target="row",backgroundColor = styleEqual(c(-1,0,1),
                                                                                                  c("pink","lightgreen","lightgrey")))
    
  })
  
  #### output 4PL ANOVA tests ---
  output$ANOVA <- DT::renderDataTable({
    sigmoid <- sigmoid()
    tab <- ANOVA4plUnresfunc(sim2(),sigmoid)
    dat <- datatable(tab,
                     options=list(
                       dom="t",rownames=F
                     )) %>% formatStyle("p_value", target="row",
                                        backgroundColor = styleEqual(c("p_value"),
                                                                     c("lightgrey")))
  })
  
  #### output RMSEs ----
  output$RMSE <- renderText({
    paste("RMSE (unrestricted model):", Dat$RMSE_unr, "(~ entered % upper-lower asymptote)\n",
          "RMSE restricted model:", Dat$RMSE_r, "\n",
          "Pure RMSE unrestricted model:", Dat$RMSE_pure, "\n",
          "%SD (unr model): ", Dat$RMSE_unr*100/Dat$up_lowAs, "(calculated as: RMSE/(upper-lower Asymptote)*100\n",
          "RMSE (log restr. model): ", Dat$RMSE_Rlog, "\n",
          "RMSE (log unrestr. model): ", Dat$RMSE_Ulog, "\n",
          "%SDlog (unr model): ", Dat$RMSE_Ulog*100/Dat$up_lowAslog )
  })
  
  output$ANOVAlin <- DT::renderDataTable({
    ANOVAlin <- Dat$ANOVA
    dat <- datatable(ANOVAlin,
                     options=list(
                       dom="t",rownames=F
                     )) %>% formatStyle("p.value", target='cell',
                                        backgroundColor = styleEqual(c("p.value"),
                                                                     c("lightgrey")))
  })
  ### output pot tab ----
  output$pottab <- DT::renderDataTable({
    
    Lim <- list(as.numeric(input$lEACdiffla), as.numeric(input$uEACdiffla), 
                as.numeric(input$lEACratiola), as.numeric(input$uEACratiola), 
                as.numeric(input$lEACratioSlope), as.numeric(input$uEACratioSlope), 
                as.numeric(input$lEACratioua), as.numeric(input$uEACratioua), 
                as.numeric(input$lowerPot), as.numeric(input$upperPot))
    
    circles <- Dat$circles
    PureErrFlag <- input$PureErr
    pottab <- LinPotTab(circles,Lim,PureErrFlag = PureErrFlag)
    #browser()
    dat <- datatable(pottab,
                     options=list(
                       dom="t",rownames=F
                     )) %>% formatStyle("test_result", target='row',
                                        backgroundColor = styleEqual(c(0,1), c("lightgrey")))
  }) 
  
  #### 4pl potency table ----
 observe({
    #browser()
   if (is.null(sim2()) | is.null(Dils())) return(NULL)
    ro_new <- sim2()
    Dils_ <- Dils()
    if (!is.na(Dils()[4])) noDilSer <- Dils()[4] else noDilSer <- 3
    PureErrFl <- input$PureErr4pl
    pottab4 <- pot4plFUNC(ro_new = ro_new, PureErrFlag = PureErrFl)
#browser()
    Lim <- list(as.numeric(input$lEACdiffla), as.numeric(input$uEACdiffla), 
                as.numeric(input$lEACratiola), as.numeric(input$uEACratiola), 
                as.numeric(input$lEACratioSlope), as.numeric(input$uEACratioSlope), 
                as.numeric(input$lEACratioua), as.numeric(input$uEACratioua), 
                as.numeric(input$lowerPot), as.numeric(input$upperPot))
    
    
    pottab4_ <- data.frame(pottab4)
    pottab4_$potency <- as.numeric(pottab4[,2])*100
    pottab4_$`lower95%CI` <- as.numeric(pottab4[,3])*100
    pottab4_$`upper95%CI` <- as.numeric(pottab4[,4])*100
    pottab4_$relative_lowerCL <- round(pottab4_[,6]/pottab4_[,5]*100,3)
    pottab4_$relative_upperCL <- round(pottab4_[,7]/pottab4_[,5]*100,3)
    
    if (as.numeric(pottab4_$relative_lowerCL[1]) > Lim[[9]] & as.numeric(pottab4_$relative_upperCL[1]) < Lim[[10]] ) {
      test_potCI <- 0
    } else {test_potCI <- 1 }
    if (as.numeric(pottab4_$relative_lowerCL[2]) > Lim[[9]] & as.numeric(pottab4_$relative_upperCL[2]) < Lim[[10]] ) {
      test_potUCI <- 0
    } else {test_potUCI <- 1 }
    if (as.numeric(pottab4_$relative_lowerCL[3]) > Lim[[9]] & as.numeric(pottab4_$relative_upperCL[3]) < Lim[[10]]  ) {
      test_potCI_t <- 0
    } else {test_potCI_t <- 1 }
    if (as.numeric(pottab4_$relative_lowerCL[4]) > Lim[[9]] & as.numeric(pottab4_$relative_upperCL[4]) < Lim[[10]] ) {
      test_potUCI_t <- 0
    } else {test_potUCI_t <- 1 }
    pottab4_ <- cbind(pottab4_[,-(2:4)], data.frame(tests=c(test_potCI, test_potUCI,test_potCI_t,test_potUCI_t)))
    colnames(pottab4_) <- c("model","potency","lower95%CI","upper95%CI","relative_lower95%CI","relative_upper95%CI","test_result")
    
    output$pottab4pl <- DT::renderDataTable({
    dat <- datatable(pottab4_[1:2,],
                     options=list(
                       paging=T, dom="t",rownames=F
                     )) %>% formatStyle("test_result", target="row",backgroundColor = styleEqual(c(0,1),
                                                                                                 c("lightgreen","pink")))
    })
    output$pottab4plTrans <- DT::renderDataTable({
      dat <- datatable(pottab4_[3:4,],
                       options=list(
                         paging=T, dom="t",rownames=F
                       )) %>% formatStyle("test_result", target="row",backgroundColor = styleEqual(c(0,1),
                                                                                                   c("lightgreen","pink")))
    })
})
  
  #### Dilutions Simulator ----
  output$plotfordilutions <- renderPlot({
    tab <- sim2()
 #browser()   
    tab <- as.data.frame(tab)
    dils <- tab$log_dose
    min_y <- min(tab[,1:3])
    max_y <- max(tab[,1:3])
    
    if (input$fixupper) {
      dils_av <- dils-max(dils)
      dils_av_ <- dils_av*(input$dilslider/100+1)
      dils2 <- round(dils_av_ + max(dils),4)
      dilfactors <- 1/exp(dils2-lag(dils2))
    } else {
      if (!is.null(Dat$cfordils)) {
        av <- Dat$cfordils
      } else { av <- (min(dils) + max(dils))/2 }
      dils_av <- dils-av
      dils_avsc <- dils_av*(input$dilslider/100+1)
      dils2 <- dils_avsc+av
      dilfactors <- 1/exp(dils2-lag(dils2))
    }
    
    Dat$newDils <- dils2
    
    sigmoid <- sigmoid()
    
    #browser()  
    BPs <- Dat$bendpoints
    EC50REF <- (BPs[2]+BPs[1])/2
    Einh <- abs((BPs[2]-BPs[1])/5)
    asyml <- EC50REF-2*(EC50REF-BPs[1])
    asymu <- EC50REF+2*(EC50REF-BPs[1])
    
    det_sig <- Dat$coeffs_UN
    
    if (is.null(Dat$coeffs_UN)) {
      SAMPLE50 <- sigmoid[1] + (sigmoid[3] - sigmoid[1])/(1+exp(sigmoid[5]*( (sigmoid[7]+0.693147)- dils2)))
      SAMPLE200 <- sigmoid[1] + (sigmoid[3] - sigmoid[1])/(1+exp(sigmoid[5]*( (sigmoid[7]-0.693147)-dils2)))
      Xbend50l <- sigmoid[7] + 0.693147-1.31696/sigmoid[5]
      Xbend200l <- sigmoid[7] - 0.693147-1.31696/sigmoid[5]
      Xbend50u <- sigmoid[7] + 0.693147+1.31696/sigmoid[5]
      Xbend200u <- sigmoid[7] - 0.693147+1.31696/sigmoid[5]
      Xbend50 <- max(Xbend50l, Xbend50u)
      Xbend200 <- min(Xbend200l, Xbend200u)
      dummy <- plot_f(tab,sigmoid,det_sig=NULL)
    } else {
      
      #browser()
      SAMPLE50 <- det_sig[3] + (det_sig[5] - det_sig[3])/(1+exp(det_sig[1]*(det_sig[7]+0.693147-dils2)))
      SAMPLE200 <- det_sig[3] + (det_sig[5] - det_sig[3])/(1+exp(det_sig[1]*(det_sig[7]-0.693147-dils2)))
      Xbend50l <- det_sig[7] + 0.693147-1.31696/det_sig[1]
      Xbend200l <- det_sig[7] - 0.693147-1.31696/det_sig[1]
      Xbend50u <- det_sig[7] + 0.693147+1.31696/det_sig[1]
      Xbend200u <- det_sig[7] - 0.693147+1.31696/det_sig[1]
      Xbend50 <- max(Xbend50l, Xbend50u)
      Xbend200 <- min(Xbend200l, Xbend200u)
      dummy <- plot_f(tab,sigmoid=NULL,det_sig=det_sig)
    }
    
    
    
    
    
    pl_df <- cbind(dils2, SAMPLE50, SAMPLE200)
    
    #browser()
    # scenario2
    eqSpac <- abs((BPs[1]-BPs[2])/5)
    optdils <- c((asyml+BPs[1])/2, BPs[1], BPs[1]+1*eqSpac, BPs[1]+2*eqSpac,BPs[1]+3*eqSpac,BPs[1]+4*eqSpac,BPs[2], (asymu+BPs[2])/2)
    # scenario 3
    eqSpac_3 <- abs((BPs[1]-BPs[2])/3)
    optdils_3 <- c(BPs[1]-2*eqSpac_3, BPs[1]-eqSpac_3, BPs[1], BPs[1]+1*eqSpac_3, BPs[1]+2*eqSpac_3,BPs[2], BPs[2]+eqSpac_3, BPs[2]+2*eqSpac_3)
    # scenario 6
    Einh2 <- abs(((BPs[2]-BPs[1])*0.7)/5)
    eqSpac2 <- (2*0.7/Einh)/3
    optdils2 <- c((asyml+BPs[1])/2, BPs[1], EC50REF-1.5*Einh2, EC50REF-0.5*Einh2,EC50REF+0.5*Einh2,EC50REF+1.5*Einh2, BPs[2], (asymu+BPs[2])/2)
    # steep slope
    eqSpac3 <- (abs(Xbend200-Xbend50))/5
    optdils3 <- c(Xbend200-eqSpac3,Xbend200, Xbend200+1*eqSpac3, Xbend200+2*eqSpac3,Xbend200+3*eqSpac3,Xbend200+4*eqSpac3,Xbend50, Xbend50+eqSpac3)
    
    output$extremebps <- renderTable({
      ExtremeBPs <- c(Xbend50,Xbend200)
      DF2 <- data.frame(sample=c("50% sample (right)", "200% sample (left)"), Extreme_BPs=ExtremeBPs)
      DF2
    })
    
    optD <- data.frame(cbind(optdils, optdils_3,optdils2, optdils3))
    colnames(optD) <- c("scenario2","scenario3","scenario6","steep slope")
    
    output$optimalDils <- renderTable({ optD })
    
    output$adjlogdil <- renderTable({ 
      adjlogdilfactors <- round(dilfactors,3)
      adjlogdils <- round(dils2,3)
      adjdils <- round(exp(dils2),3)
      DilsTable <- data.frame('adjusted ln(dilutions)' = adjlogdils,
                              'adjusted ln_dilution_factors' = adjlogdilfactors,
                              'adjusted dilutions' = adjdils)
      DilsTable
    })
    
    if (!is.null(Dat$p2)) {
      p2 <- Dat$p2
      p_dil <- p2 +
        annotate("pointrange",x=dils2,y=rep(min_y, length(dils2)), xmin=min(dils2), xmax=max(dils2)) +
        annotate("text", x=dils2,y=rep(min_y+(max_y-min_y)*0.05, length(dils2)), label=as.character(round(dils2,3))) +
        annotate("text", x=dils2[-1]+(max(dils2)-min(dils2))*0.05,
                 y=rep(min_y+(max_y-min_y)*0.1, length(dils2[-1])),
                 label=as.character(round(dilfactors[-1],3))) +
        geom_line(data=as.data.frame(pl_df),aes(x=dils2,y=SAMPLE50), color="grey15", linetype=2,
                  inherit.aes = F) +
        geom_line(data=as.data.frame(pl_df),aes(x=dils2,y=SAMPLE200), color="grey15", linetype=2,
                  inherit.aes = F) +
        geom_vline(xintercept=c(Xbend50,Xbend200), col="grey15", linetype=2) +
        {if (input$scenario =="scenario 6") annotate("pointrange",x=optdils2,y=rep(min_y+(max_y-min_y)*0.2, length(optdils2)),
                                                     xmin=min(optdils2), xmax=max(optdils2), color="seagreen")} +
        {if (input$scenario =="scenario 6") annotate("text",x=optdils2,y=rep(min_y+(max_y-min_y)*0.25, length(optdils2)),
                                                     label=as.character(round(optdils2,3)), color="seagreen")} +
        {if (input$scenario =="scenario 2") annotate("pointrange",x=optdils,y=rep(min_y+(max_y-min_y)*0.2, length(optdils)),
                                                     xmin=min(optdils), xmax=max(optdils), color="seagreen")} +
        {if (input$scenario =="scenario 2") annotate("text",x=optdils,y=rep(min_y+(max_y-min_y)*0.25, length(optdils)),
                                                     label=as.character(round(optdils,3)), color="seagreen")} +
        {if (input$scenario =="scenario 3") annotate("pointrange",x=optdils_3,y=rep(min_y+(max_y-min_y)*0.2, length(optdils_3)),
                                                     xmin=min(optdils_3), xmax=max(optdils_3), color="seagreen")} +
        {if (input$scenario =="scenario 3") annotate("text",x=optdils_3,y=rep(min_y+(max_y-min_y)*0.25, length(optdils_3)),
                                                     label=as.character(round(optdils_3,3)), color="seagreen")} +
        {if (input$scenario =="steep slope") annotate("pointrange",x=optdils3,y=rep(min_y+(max_y-min_y)*0.2, length(optdils3)),
                                                      xmin=min(optdils3), xmax=max(optdils3), color="seagreen")} +
        {if (input$scenario =="steep slope") annotate("text",x=optdils3,y=rep(min_y+(max_y-min_y)*0.25, length(optdils3)),
                                                      label=as.character(round(optdils3,3)), color="seagreen")} +
        annotate("text",x=optdils[1],y=(max_y+min_y)*0.5,
                 label=paste("in green: optimal \n dilutions acc. to Whitepaper\n", input$scenario), color="seagreen",
                 size=14/.pt,fontface="bold")
    }
    print(p_dil)
    
  })
  
  #### Dilutions CI table ----
  observe({
    if (is.null(input$potencydiff)) return(NULL)
    output$CIs <- renderTable({
      PureErrFlag <- input$PureErr
      if (is.null(Dat$coeffs_UN)) {
        # checks if an EXCEL was uploaded
        sigmoid <- sigmoid()
        det_sig=NULL
        
        ast = sigmoid()[1];bst = sigmoid()[5];cst = sigmoid()[7];dst = sigmoid()[3];ate = sigmoid()[2];
        bte = sigmoid()[6];r_ = sigmoid()[8];
        cte = cst-r_;dte = sigmoid()[4];
      } else {
        sigmoid <- NULL
        det_sig <- Dat$coeffs_UN
        ast <- det_sig[3]
        ate <- det_sig[4]
        bst <- det_sig[1]
        bte <- det_sig[2]
        cst <- det_sig[7]
        cte <- det_sig[7] -log(input$potencydiff/100)
        dst <- det_sig[5]
        dte <- det_sig[6]
        r_ <- log(input$potencydiff/100)
        
      }
      if (!is.na(input$NoDilSer)) {
        noDilSer <- input$NoDilSer
      } else if (!is.null(Dat$NoDilSeriesXL)) noDilSer <- Dat$noDilSeriesXL else noDilSer <- 3
      if (!is.na(input$NoDil)) noDil <- input$NoDil else noDil <- length(Dat$newDils)
      #browser()     
      tab <- Calc_DilRes(as=ast,at=ate,ds=dst,dt=dte,cs=cst,ct=cte,r=r_,bt=bte,bs=bst, 
                         sd_fac=input$sdfac,log_conc=Dat$newDils,
                         # auslenkU=outlierU,
                         # auslenkM=outlierM,
                         # auslenkL=outlierL,
                         heteroNoise = FALSE, noDilSeries = noDilSer, noDils = noDil)
      Limite <- list(as.numeric(input$lEACdiffla), as.numeric(input$uEACdiffla), 
                     as.numeric(input$lEACratiola), as.numeric(input$uEACratiola), 
                     as.numeric(input$lEACratioSlope), as.numeric(input$uEACratioSlope), 
                     as.numeric(input$lEACratioua), as.numeric(input$uEACratioua), 
                     as.numeric(input$lowerPot), as.numeric(input$upperPot), 
                     as.numeric(input$lEACratioAdiff), as.numeric(input$uEACratioAdiff))
      
      CItable <- tests_FUNC(tab,Limite,PureErrFlag=PureErrFlag)
      
      CItable_ <- CItable[-c(1,2,6,8,9),-c(2,4,5)]
      potAll <- pot4plFUNC(tab, input$PureErr)
      restrPot <- potAll[1,1:4]
      restrPot[2:4] <- round(as.numeric(restrPot[2:4]),5)
      potAll_ <- rbind(CItable_, restrPot)
      potAll_$CIwidth <- as.numeric(potAll_[,4])-as.numeric(potAll_[,3])
      potAll_[,1] <- c("ratio of lower asymptotes","ratio of slopes","ratio of upper asymptotes", "ratio of asympt. differences","restricted potency")
      
      output$bps <- renderTable({
        DF <- data.frame(sample=names(Dat$bendpoints),BPs=Dat$bendpoints)
        DF
      })
      return(potAll_)
    })
  })
  
  #### simulations ----
  observe({
    observeEvent(input$goSim,{
      sd_fac_ <- as.numeric(input$sdfac)
      
      r_ <- log(as.numeric(input$potencydiff)/100)
      Conc <- Dat$MetaConc
      as = sigmoid()[1]; bs = sigmoid()[5];cs = sigmoid()[7];ds = sigmoid()[3];at = sigmoid()[2];
      bt = sigmoid()[6];r = sigmoid()[8]; ct = cs-r_; dt = sigmoid()[4]
      
      if (!is.null(Dat$MetaConc)) {
        Conc <- Dat$MetaConc
      } else {
        Conc <- CONC()
      }
      log_dose <- log(Conc)
      yAxfac <- (ds-as)
      
      if (!is.na(input$NoDilSer)) {
        noDilSer <- input$NoDilSer
      } else if (!is.null(Dat$NoDilSeriesXL)) noDilSer <- Dat$noDilSeriesXL else noDilSer <- 3
      if (!is.na(input$NoDil)) noDil <- input$NoDil else noDil <- length(Conc)
      isRef <- rep(c(1,0),1,each=noDilSer*noDil) 
      isSample <- rep(c(0,1),1,each=noDilSer*noDil)
      N <- as.numeric(input$simN)
      
      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)))
      
      resHist <- matrix(NA,nrow=N, ncol=13)
      residualsList <- list()
      start.time2 <- Sys.time()
      withProgress(message = 'Making plot', value=0, {
        for (i in 1:N) {
          if (input$heterosked) {
            # heterosc noise
            ro_jit <- matrix(unlist(map(av, function(x) x+rnorm(1,0,x*sd_fac_/100))), nrow=noDil, ncol=noDilSer*2)
          } else {
            # homosc noise
            ro_jit <- matrix(unlist(map(av, function(x) x+rnorm(1,0,sd_fac_*yAxfac/100))), nrow=noDil, ncol=noDilSer*2)
          }
          # browser()
          ro_jit <- abs(ro_jit)
          ro_new <- cbind(ro_jit, log_dose)
          all_l <- melt(data.frame(ro_new), id.vars="log_dose", variable.name = "replname", value.name = "readout")
          all_l$isRef <- isRef
          all_l$isSample <- isSample
          all_l$Conc <- exp(all_l$log_dose)
          pot <- drm(readout ~ Conc, isSample, data=all_l, fct=LL.4(names=c("b","d","a","c")),
                     pmodels=data.frame(1,1,1,isSample))
          potAll <- EDcomp(pot, percVec=c(50,50), interval="delta", display=FALSE)
          potAll2 <- potAll[1:3]
          RSS <- sum(pot$predres[,2]^2)
          dfreed <- nrow(all_l)-5
          MSE <- RSS/dfreed
          
          potU <- drm(readout ~ Conc, isSample, data=all_l, fct=LL.4(names=c("b","d","a","c")),
                      pmodels=data.frame(isSample, isSample,isSample,isSample))
          DF_U <- nrow(all_l)-8
          
          uAsratio <- compParm(potU, "a",display=F)
          uCIuAs <- uAsratio[1]+qt(0.975,DF_U)*uAsratio[2]
          lCIuAs <- uAsratio[1]-qt(0.975,DF_U)*uAsratio[2]
          lAsratio <- compParm(potU, "d",display=F)
          uCIlAs <- lAsratio[1]+qt(0.975,DF_U)*lAsratio[2]
          lCIlAs <- lAsratio[1]-qt(0.975,DF_U)*lAsratio[2]
          Sloperatio <- compParm(potU, "b",display=F)
          uCISlo <- Sloperatio[1]+qt(0.975,DF_U)*Sloperatio[2]
          lCISlo <- Sloperatio[1]-qt(0.975,DF_U)*Sloperatio[2]
          su <- summary(potU)
          
          v <- vcov(potU)[c(5,6),c(5,6)]
          Vd <- vcov(potU)[c(3,4),c(3,4)]
          Va_d <- v+Vd
          
          A_DTEST <- su$coefficients[6,1]-su$coefficients[4,1]
          A_DREF <- su$coefficients[5,1]-su$coefficients[3,1]
          if (abs(at/(sqrt(Va_d[2,2]/3))) > qt(0.95,2)) {
            try(Fie_ad <- round(FiellerRatio(A_DREF,A_DTEST, Va_d),5))
          }
          if (!exists("Fie_ad")) Fie_ad <- NA
          
          resHist[i,] <-  c(potAll2, sqrt(MSE),Sloperatio[1],lCISlo, uCISlo,
                            uAsratio[1], lCIuAs, uCIuAs, Fie_ad[1],Fie_ad[2],Fie_ad[3])
          colnames(resHist) <- c("pot4pl","lCI4pl","uCI4pl","RMSE","estSlope_ratio",
                                 "lCISlope_ratio","uCISlope_ratio","estuAs_ratio",
                                 "lCIuAs_ratio","uCIuAs_ratio","estAsyDiff_ratio",
                                 "lCIAsyDiff_ratio", "uCIAsyDiff_ratio")
          
          incProgress(1/N, detail=paste("Doing simulations",i))
        } # withProgress    
        
      })
      end.time2 <- Sys.time()
      
      Dat$resHist <- resHist
    })
  })
  
  
  #### simulation Histograms output ----
  
  output$plotHistuAs <- renderPlot({
    if (!is.null(Dat$resHist)) {
      
      resHist <- Dat$resHist
      #browser()     
      resHistuAs <- as.data.frame(resHist[,8:10])
      resHistuAs_l <- melt(data.frame(resHistuAs), variable.name="ratio_CIs", value.name = "readout")
      #browser()  
      lowquant_uAs <- quantile(resHistuAs[,2], probs=as.numeric(input$lowQuant)/100)
      upquant_uAs <- quantile(resHistuAs[,3], probs=as.numeric(input$uppQuant)/100)
      
      p_uAs <- ggplot(resHistuAs_l) +
        geom_histogram(aes(readout, fill=ratio_CIs),alpha=0.5,position="identity") +
        labs(title = paste("upper asymptote ratio EACs:", round(lowquant_uAs,3), " to ", round(upquant_uAs,3))) +
        geom_vline(xintercept = c(lowquant_uAs, upquant_uAs), color="black", linetype="dashed", linewidth=1) +
        geom_vline(xintercept = c(input$lEACratioua , input$uEACratioua), color="red", linetype="dashed", linewidth=1) +
        theme_bw()
      
      # asympt diff ratio
      resHistAsDiff <- as.data.frame(resHist[,11:13])
      resHistAsDiff_l <- melt(data.frame(resHistAsDiff), variable.name="ratio_CIs", value.name = "readout")
      
      lowquant_AsDiff <- quantile(resHistAsDiff[,2], probs=as.numeric(input$lowQuant)/100)
      upquant_AsDiff <- quantile(resHistAsDiff[,3], probs=as.numeric(input$uppQuant)/100)
      
      p_AsDiff <- ggplot(resHistAsDiff_l, aes(readout, fill=ratio_CIs)) +
        geom_histogram(alpha=0.5,position="identity") +
        labs(title = paste("asymptote diff. ratio EACs:", round(lowquant_AsDiff,3), " to ", round(upquant_AsDiff,3))) +
        geom_vline(xintercept = c(lowquant_AsDiff, upquant_AsDiff), color="black", linetype="dashed", linewidth=1) +
        geom_vline(xintercept = c(input$lEACratioAdiff , input$uEACratioAdiff), color="red", linetype="dashed", linewidth=1) +
        theme_bw()
      
      # Slope ratio
      resHistSlo <- as.data.frame(resHist[,5:7])
      resHistSlo_l <- melt(data.frame(resHistSlo), variable.name="ratio_CIs", value.name = "readout")
      
      lowquant_Slo <- quantile(resHistSlo[,2], probs=as.numeric(input$lowQuant)/100)
      upquant_Slo <- quantile(resHistSlo[,3], probs=as.numeric(input$uppQuant)/100)
      
      p_Slo <- ggplot(resHistSlo_l, aes(readout, fill=ratio_CIs)) +
        geom_histogram(alpha=0.5,position="identity") +
        labs(title = paste("Slope ratio EACs:", round(lowquant_Slo,3), " to ", round(upquant_Slo,3))) +
        geom_vline(xintercept = c(lowquant_Slo, upquant_Slo), color="black", linetype="dashed", linewidth=1) +
        geom_vline(xintercept = c(input$lEACratioSlope , input$uEACratioSlope), color="red", linetype="dashed", linewidth=1) +
        theme_bw()
      
      # poency ratio
      resHistPot <- as.data.frame(resHist[,1:3])
      resHistPot_l <- melt(data.frame(resHistPot), variable.name="ratio_CIs", value.name = "readout")
      
      lowquant_Pot <- quantile(resHistPot[,2], probs=as.numeric(input$lowQuant)/100)
      upquant_Pot <- quantile(resHistPot[,3], probs=as.numeric(input$uppQuant)/100)
      #browser()    
      p_Pot <- ggplot(resHistPot_l, aes(readout, fill=ratio_CIs)) +
        geom_histogram(alpha=0.5,position="identity") +
        labs(title = paste("Poency ratio EACs:", round(lowquant_Pot,3), " to ", round(upquant_Pot,3))) +
        geom_vline(xintercept = c(lowquant_Pot, upquant_Pot), color="black", linetype="dashed", linewidth=1) +
        geom_vline(xintercept = c(input$lowerPot/100, input$upperPot/100), color="red", linetype="dashed", linewidth=1) +
        theme_bw()
      
      grid.arrange(p_Slo, p_AsDiff, p_uAs, p_Pot, nrow=1)
      
    }
  })
  
  #### download XL report----
  
  output$downloadXLReport <- downloadHandler(
    filename= paste0("Report_4PLEvaluation", Dat$FileName,".pdf"),
    
    content = function(file) {
      tpdr <- tempdir()
      tempReport <- file.path(tpdr,"Doc_BioassayReport.Rmd")
      file.copy("Doc_BioassayReport.Rmd", tempReport, overwrite = T)
      
      tempReportc <- file.path(tpdr,"logo.png")
      file.copy("logo.png", tempReportc, overwrite = T)
      
      rmarkdown::render(tempReport, output_file = file,
                        params = list(FileName = Dat$FileName,
                                      author = Dat$author,
                                      REP = REP,
                                      coeffs = Dat$coeffs_UN),
                        envir = new.env(parent = globalenv()))
      
    }
  )
  
  
}

shinyApp(ui, server)
