library(tidyverse)
library(gslnls)
library(ggplot2)
library(reshape2)
library(drc)
library(ggridges)

library(car)



# Pseudoreplication
# simulation of replication strategies

Calc_DilRes <- function(as = 10, bs = 1, cs = -6, ds = 110, at = 10, bt = 1, r=0.01,
                        ct = cs, dt = 110, sd_fac=0.1, gt=1, gs=1, log_concREF, log_concTEST,  # -r
                        THEOconc ,
                        heteroNoise=FALSE, noDilSeries, noDils) {   
  
  # browser()
  yAxfac <- (as-ds)
  log_cREF_l <- melt(log_concREF)
  log_doseREF <- log_cREF_l$value
  log_cTEST_l <- melt(log_concTEST)
  log_doseTEST <- log_cTEST_l$value
  isRef <- rep(c(1,0),1,each=length(log_concREF))
  isSample <- rep(c(0,1),1,each=length(log_concTEST))
  #um with equal sample without 
  ro_jit <- as*isRef + at*isSample + (ds*isRef + dt*isSample - as*isRef - at*isSample)/
    (1 + isRef*exp(bs*(cs - log_doseREF)) + isSample*exp(bt*(ct -log_doseTEST)))
  
  ro_jit_ <- abs(ro_jit)
  all_var <- rbind(log_cREF_l, log_cTEST_l)
  
  ro_jit2 <- cbind(all_var, ro_jit_, rep(THEOconc, 6), isRef, isSample)
  # if (noDilSeries==3) {
  #   colnames(ro_jit2) <- c("R_dil1","R_dil2","R_dil3","T_dil1","T_dil2","T_dil3", "log_dose")
  # } else if (noDilSeries==2) {
  #   colnames(ro_jit2) <- c("R_dil1","R_dil2","T_dil1","T_dil2", "log_dose")
  # }
  
  colnames(ro_jit2) <- c("dil_no","repl","trueLogConc" ,"readout","log_dose","isRef","isSample")
  return(ro_jit2)
}


RestrMF <- function(data, M) {
  
  all_l <- data
  #plot(all_l$log_dose, all_l$readout, col=all_l$sample)
  CORro <- cor(all_l$log_dose[1:8], all_l$readout[1:8])
  if (CORro<0) B <- -1 else B <- 1
  
  # browser()
  startlist <- list(a = min(all_l$readout), b = B, cs=mean(all_l$log_dose),
                    r=0.0001, d = max(all_l$readout))
  mr <- tryCatch({
    gsl_nls(fn = readout ~ a + (d - a)/(1 + exp(b*((cs-r*isSample)-log_dose))),
            data = all_l,
            start = startlist, #trace=T,
            control = gsl_nls_control(xtol = 1e-6, ftol = 1e-6, gtol = 1e-6))
  }, error = function(msg){
    return(0)
  } )
  
  if (length(mr)>1) {
    s_mr_coef <- summary(mr)$coefficients
    s_mr <- summary(mr)
  } else { s_mr_coef <- rep(0,4); s_mr <- 0 }
  
  if (M == "R") { 
    return(s_mr_coef) 
  } else if (M =="CI") {
    if (length(mr)==1) return(rep(0,4))
    (RelPot_mr <- exp(s_mr_coef[4,1])) # r exponentiated
    
    POTr_CI <- exp(confint(mr, "r", method = "asymptotic", level=0.95))
    # Pure error
    #browser()
    suppressWarnings(FitAnova <- anova(lm(readout ~ factor(log_dose)*isSample, all_l)))
    meanPureErr <- FitAnova[4,3]
    SEsPure <- sqrt(diag(vcov(mr)/s_mr$sigma^2)*meanPureErr)
    
    # PLA software: 1.0252 - 1.4242
    lCIpure <- exp(s_mr_coef[4,1]-qt(0.975,nrow(all_l)-16)*SEsPure[4])
    uCIpure <- exp(s_mr_coef[4,1]+qt(0.975,nrow(all_l)-16)*SEsPure[4])
    if (is.infinite(uCIpure)) uCIpure <- NA
    if (is.infinite(POTr_CI[2])) POTr_CI[2] <- NA
    relCIpure <- c(lCIpure/RelPot_mr, uCIpure/RelPot_mr)
    CIsDF <- c(RelPot_mr,POTr_CI, lCIpure, uCIpure)
    names(CIsDF) <- c("potency","lCI","uCI","lpureCI","upureCI")
    return(CIsDF)
  } else if (M=="U") {
    if (length(mr)==1) return(rep(0,4))
    
    startlstmu <- list(as = min(all_l$readout), bs = B, cs = mean(all_l$log_dose), ds = max(all_l$readout), 
                       at = min(all_l$readout), bt = B, ct = mean(all_l$log_dose), dt = max(all_l$readout))
    mu <- tryCatch({
      gsl_nls(readout ~ 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))),
              data = all_l,
              start = startlstmu, trace=T,
              control = gsl_nls_control(xtol = 1e-10, ftol = 1e-10, gtol = 1e-10))
    }, error = function(msg){
      return(0)
    } )
    if (length(mu)>1) {
      smu_coef <- tryCatch({
        summary(mu)$coefficients
      }, error = function(msg){
        return(rep(NA, 8))
      } )
      
      #s_mu <- summary(mu)
      return(smu_coef)
    } else return(rep(0,8))
    
  } else return()
}




Dil8F <- function(Cn1) {
  CSer <- matrix(NA, nrow=8, ncol=3)
  for (i in 1:8) {
    Cn1 <- Cn1/dilfactor+rnorm(3,0,SD*Cn1)
    CSer[i,] <- Cn1
  }
  return(CSer)
}

DilDirF <- function(OrigMenge) {
  CDir <- matrix(NA, nrow=10, ncol=3)
  for (DilExp in 1:10) {
    C_d <- OrigMenge/dilfactor^DilExp+rnorm(3,0,SD*OrigMenge/dilfactor^DilExp)
    CDir[DilExp,] <- C_d
  }
  return(CDir)
}

homhetCIF <- function(AssD_) {
  # Confidence level
  alpha <- 0.05
  z <- qnorm(1 - alpha/2)  # 1.96 for 95% CI
  
  ### 1. Homogeneously Weighted CI ###
  mean_logRP <- mean(AssD_$logRP)
  var_mean <- mean(AssD_$SE^2) / nrow(AssD_)
  SE_mean <- sqrt(var_mean)
  
  CI_homo <- c(mean_logRP - z * SE_mean, mean_logRP + z * SE_mean)
  CI_homo_exp <- exp(CI_homo)
  # cat("Homogeneous weighting:\n")
  # cat(round(exp(mean_logRP),3), "95% CI =", round(CI_homo_exp[1], 3), "-", round(CI_homo_exp[2], 3), "\n\n")
  
  ### 2. Heterogeneously Weighted CI ###
  weights <- 1 / (AssD_$SE^2)
  weighted_mean <- sum(weights * AssD_$logRP) / sum(weights)
  var_weighted <- 1 / sum(weights)
  SE_weighted <- sqrt(var_weighted)
  
  CI_hetero <- c(weighted_mean - z * SE_weighted, weighted_mean + z * SE_weighted)
  CI_hetero_exp <- exp(CI_hetero)
  res <- c(exp(mean_logRP) , CI_homo_exp, exp(weighted_mean), CI_hetero_exp )
  return(res)
  
  # cat("Heterogeneous weighting:\n")
  # cat("Combined RP =", round(exp(weighted_mean), 3), "\n")
  # cat("95% CI =", round(CI_hetero_exp[1], 6), "-", round(CI_hetero_exp[2], 6), "\n")
}

Conc <- c(1)
for (x in 1:10) Conc <- c(Conc, (1/3^x)) 


N <- 10000
OrigMenge <- 1
dilfactor <- 3
SD <- 0.02

rm(CIsM_all, DevTargC3_)
for (x in 1:3) {
  if (x==1) { Repl1 <- x; Repl2 <- x } # Anzahl Replikate für C2 und C3
  if (x==2) { Repl1 <- 1; Repl2 <- 3 }
  if (x==3) { Repl1 <- x; Repl2 <- x }
  
  Flags <- vector()
  SDs1 <- SDs2 <- POT <- DevTargC3 <- hhCIV <- c()
  for (nTrials in 1:N) {
    C2 <- OrigMenge/dilfactor+rnorm(Repl1,0,SD*OrigMenge) 
    C3 <- C2/dilfactor+rnorm(Repl1,0,C2*SD) 
    CSer <- Dil8F(C3)
    
    # "true" replication
    C2_1 <- OrigMenge/dilfactor+rnorm(Repl2,0,SD*OrigMenge) 
    C3_1 <- C2_1/dilfactor+rnorm(Repl2,0,SD*C2_1) 
    CSer_3 <-Dil8F(C3_1)
    
    
    # Flag <- sd(CSer)<sd(CSer_3)
    # SDs1 <- c(SDs1, sd(CSer))
    # SDs2 <- c(SDs2, sd(CSer_3))
    #Flags <- c(Flags, Flag)
    
    # all from stem solution directly
    CDir <- DilDirF(OrigMenge)
    CDir2 <- DilDirF(OrigMenge)
    if (x==1) { 
      DevTargC3 <- c(DevTargC3, mean(C3), mean(CDir[2,])) 
      log_concREF <- log(CSer)
      log_concTEST <- log(CDir[3:10,])
      testC <- "pseudo"; refC <- "direct"
    } # Anzahl Replikate für C2 und C3
    if (x==2) { 
      DevTargC3 <- c(DevTargC3, mean(C3_1), mean(CDir[2,])) 
      log_concREF <- log(CSer_3)
      log_concTEST <- log(CDir[3:10,])
      testC <- "true"; refC <- "direct"
    }
    if (x==3) { 
      DevTargC3 <- c(DevTargC3, mean(CDir[2,]), mean(CDir2[2,]))
      log_concREF <- log(CDir[3:10,])
      log_concTEST <- log(CDir2[3:10,])
      testC <- "true"; refC <- "true"
    }
    # DevTargC3 <- c(DevTargC3, mean(C3), mean(CDir[2,])) # c(DevTargC3, mean(CSer[4,]), mean(CSer_3[4,]))
    
    ro_new <- Calc_DilRes( log_concREF=log_concREF, 
                           log_concTEST=log_concTEST,THEOconc = log(Conc)[4:11])  
    CIs <- RestrMF(ro_new, M="CI")
    RM <- RestrMF(ro_new, M="R")
    
    POT <- c(POT, CIs)
    
    CIsV <- c()
    for (z in 0:2) {
      CL <- RestrMF(ro_new[c((z*8+1):(8*(z+1)),(z*8+25):(8*(z+1)+24)),], M="CI")
      #print(RestrMF(ro_new[c((z*8+1):(8*(z+1)),(z*8+25):(8*(z+1)+24)),], M="R"))
      CIsV <- c(CIsV, CL[1:3])
    }
    AssD <- t(matrix(CIsV, nrow=3))
    logAssD <- log(AssD)
    AssD_ <- as.data.frame(AssD)
    colnames(AssD_) <- c("RP", "lCI", "uCI")
    AssD_$SE <- (logAssD[,3]-logAssD[,1])/qt(0.975,11)
    AssD_$logRP <- logAssD[,1]
    
    hhCIs <-  homhetCIF(AssD_)
    hhCIV <- c(hhCIV, hhCIs)
    # p <- ggplot(ro_new, aes(x=trueLogConc, y=readout, col=as.factor(isRef))) +
    #   geom_point() +
    #   geom_line() +
    #   ggtitle(paste(nTrials, x)) +
    #   theme_bw()
    # print(p)
    # 
    # p2 <- ggplot(ro_new, aes(x=log_dose, y=readout, col=as.factor(isRef))) +
    #   geom_point(size=4) +
    #   ggtitle(paste(nTrials, x)) +
    #   geom_vline(xintercept = c(RM[3,1],(RM[3,1]-RM[4,1])), col=c("turquoise","red")) +
    #   theme_bw()
    # print(p2)
  }
  CIsM <- t(matrix(POT, nrow=5))
  CIsM <- cbind(CIsM, rep(x,nrow(CIsM)))
  print(paste( testC, refC,sd(CIsM[,1])))
  if (!exists("CIsM_all"))  CIsM_all <- CIsM else CIsM_all <- rbind(CIsM_all, CIsM)
  
  DevTargC3DF <- as.data.frame(t(matrix(DevTargC3, nrow=2)))
  if (!exists("DevTargC3_"))  DevTargC3_ <- DevTargC3DF  else DevTargC3_ <- rbind(DevTargC3_, DevTargC3DF )
  
  hhCIDF <- as.data.frame(t(matrix(hhCIV, nrow=6)))
  if (!exists("hhCIDF_"))  hhCIDF_ <- hhCIDF  else hhCIDF_ <- rbind(hhCIDF_, hhCIDF )
}

CIsM_allDF <- as.data.frame(CIsM_all)
colnames(CIsM_allDF) <- c("potency","lCI","uCI","lpureCI","upureCI", "scheme")
CIsM_allDF$widthCI <- CIsM_allDF$uCI-CIsM_allDF$lCI

CIsM_allDF_ <- cbind(CIsM_allDF, DevTargC3_)
cn <- colnames(CIsM_allDF_)
cn[8:9] <- c("C3_REF","C3_TEST")
colnames(CIsM_allDF_) <- cn
CIsM_allDF_$C3_DIFF <- CIsM_allDF_$C3_TEST-CIsM_allDF_$C3_REF

cor(CIsM_allDF_$potency, CIsM_allDF_$C3_DIFF)
#cor(abs(log(CIsM_allDF_$potency)), (CIsM_allDF_[,9]-CIsM_allDF_[,8]))
plot( CIsM_allDF_$C3_DIFF,CIsM_allDF_$potency)



for (n in 0:2) {
  te <- CIsM_allDF_[(n*N+1):((n+1)*N),]
  print(paste("normal CI coverage", sum(log(te$lCI)<0 & log(te$uCI)>0) /N*100))
  print(paste("pure CI coverage", sum(log(te$lpureCI)<0 & log(te$upureCI)>0) /N*100))
  print(cor(te$potency, (te$C3_DIFF)))
  te$logpot <- log(te$potency)
  # plot(te$C3_DIFF,te$potency, pch=19, xlab="Difference at C3")
  # abline(te$potency ~ te$C3_DIFF, col = "blue")
  # lines(lowess(te$potency ~ te$C3_DIFF), col = "green")
  scatterplot(logpot ~ C3_DIFF, data=te)
}

plot(te$C3_DIFF,te$logpot, pch=19, xlab="Difference at C3")
abline(te$logpot ~ te$C3_DIFF, col = "blue")
lines(lowess(te$logpot ~ te$C3_DIFF), col = "green")
scatterplot(logpot ~ C3_DIFF, data=te)

colnames(hhCIDF_) <- c("homomean", "lhomoCI", "uhomoCI","heteromean","lheteroCI","uheteroCI")
te <- cbind(CIsM_allDF, hhCIDF_)


p_histCI <- ggplot(CIsM_allDF, aes(x=widthCI, y=as.factor(scheme), fill=as.factor(scheme))) + 
  geom_density_ridges() +
  theme_ridges() 
p_histCI

p_histCI2 <- ggplot(CIsM_allDF, aes(x=widthCI, y=as.factor(scheme), fill=as.factor(scheme))) + 
  geom_density_ridges() +
  theme_ridges() 
p_histCI2  


sum(Flags)/N*100
mean(SDs1)
mean(SDs2)

ro_new <- Calc_DilRes(sd_fac=0.1, log_concREF=log_concREF, 
                      log_concTEST=log_concTEST,THEOconc = log(Conc)[4:11])