## ----setup, include = FALSE----------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE
  , comment = "#>"
  #, fig.width=6
  #, cache = TRUE
)

## ----eval=FALSE----------------------------------------------------------
#  system.file("examples", package="PopED")

## ----eval=FALSE----------------------------------------------------------
#  ex_dir <- system.file("examples", package="PopED")
#  list.files(ex_dir)

## ----eval=FALSE----------------------------------------------------------
#  file_name <- "ex.1.a.PK.1.comp.oral.md.intro.R"
#  ex_file <- system.file("examples",file_name,package="PopED")
#  file.copy(ex_file,tempdir(),overwrite = T)
#  file.edit(file.path(tempdir(),file_name))

## ------------------------------------------------------------------------
library(PopED)
f_pkpdmodel <- function(model_switch,xt,parameters,poped.db){
  with(as.list(parameters),{
    y=xt
    MS <- model_switch
    
    # PK model
    CONC = DOSE/V*exp(-CL/V*xt) 
    
    # PD model
    EFF = E0 + CONC*EMAX/(EC50 + CONC)
    
    y[MS==1] = CONC[MS==1]
    y[MS==2] = EFF[MS==2]
    
    return(list( y= y,poped.db=poped.db))
  })
}

## ------------------------------------------------------------------------
## -- Residual Error function
## -- Proportional PK + additive PD
f_Err <- function(model_switch,xt,parameters,epsi,poped.db){
  returnArgs <- do.call(poped.db$model$ff_pointer,list(model_switch,xt,parameters,poped.db)) 
  y <- returnArgs[[1]]
  poped.db <- returnArgs[[2]]
  
  MS <- model_switch
  
  prop.err <- y*(1+epsi[,1])
  add.err <- y+epsi[,2]
  
  y[MS==1] = prop.err[MS==1]
  y[MS==2] = add.err[MS==2]
  
  return(list( y= y,poped.db =poped.db )) 
}

## ---- echo=FALSE---------------------------------------------------------
f_etaToParam <- function(x,a,bpop,b,bocc){
  parameters=c( 
    CL=bpop[1]*exp(b[1]),
    V=bpop[2]*exp(b[2]),
    E0=bpop[3]*exp(b[3]),
    EMAX=bpop[4],
    EC50=bpop[5]*exp(b[4]),
    DOSE=a[1]
  )
  return( parameters ) 
}

## ------------------------------------------------------------------------
poped.db <- create.poped.database(ff_fun=f_pkpdmodel,
                                  fError_fun=f_Err,
                                  fg_fun=f_etaToParam,
                                  sigma=diag(c(0.15,0.015)),
                                  bpop=c(CL=0.5,V=0.2,E0=1,EMAX=1,EC50=1),  
                                  d=c(CL=0.09,V=0.09,E0=0.04,EC50=0.09), 
                                  groupsize=20,
                                  m=3,
                                  xt=c( 0.33,0.66,0.9,5,0.1,1,2,5),
                                  minxt=0,
                                  maxxt=5,
                                  bUseGrouped_xt=1,
                                  model_switch=c(1,1,1,1,2,2,2,2),
                                  a=list(c(DOSE=0),c(DOSE=1),c(DOSE=2)),
                                  maxa=c(DOSE=10),
                                  mina=c(DOSE=0))



## ----simulate_multi-response_model---------------------------------------
plot_model_prediction(poped.db,IPRED=TRUE,DV=TRUE,facet_scales="free",
                      separate.groups=TRUE,
                      model.names=c("PK","PD")) 

## ----eval=TRUE-----------------------------------------------------------
library(Rcpp)
library(deSolve)

## ------------------------------------------------------------------------
PK.2.comp.oral.ode <- function(Time, State, Pars){
  with(as.list(c(State, Pars)), {    
    dA1 <- -KA*A1 
    dA2 <- KA*A1 + A3* Q/V2 -A2*(CL/V1+Q/V1)
    dA3 <- A2* Q/V1-A3* Q/V2
    return(list(c(dA1, dA2, dA3)))
  })
}

## ------------------------------------------------------------------------
ff.PK.2.comp.oral.md.ode <- function(model_switch, xt, parameters, poped.db){
  with(as.list(parameters),{
    A_ini <- c(A1=0, A2=0, A3=0)
    times_xt <- drop(xt)
    dose_times = seq(from=0,to=max(times_xt),by=TAU)
    eventdat <- data.frame(var = c("A1"), 
                           time = dose_times,
                           value = c(DOSE), method = c("add"))
    times <- sort(c(times_xt,dose_times))
    
    # Here "PK.2.comp.oral.ode" is equivalent to the pre-compiled version "two_comp_oral_ode_Rcpp".
    out <- ode(A_ini, times, PK.2.comp.oral.ode, parameters,
               events = list(data = eventdat))#atol=1e-13,rtol=1e-13)
    y = out[, "A2"]/(V1/Favail)
    y=y[match(times_xt,out[,"time"])]
    y=cbind(y)
    return(list(y=y,poped.db=poped.db))
  })
}


## ---- echo=FALSE---------------------------------------------------------
fg <- function(x,a,bpop,b,bocc){
  parameters=c( CL=bpop[1]*exp(b[1]),
                V1=bpop[2],
                KA=bpop[3]*exp(b[2]),
                Q=bpop[4],
                V2=bpop[5],
                Favail=bpop[6],
                DOSE=a[1],
                TAU=a[2])
  return( parameters ) 
}

## ------------------------------------------------------------------------
discrete_a <- cell(2,1)

poped.db <- create.poped.database(ff_fun="ff.PK.2.comp.oral.md.ode",
                                  fError_fun="feps.add.prop",
                                  fg_fun="fg",
                                  groupsize=20,
                                  m=1,      #number of groups
                                  sigma=c(prop=0.1^2,add=0.05^2),
                                  bpop=c(CL=10,V1=100,KA=1,Q= 3.0, V2= 40.0, Favail=1),
                                  d=c(CL=0.15^2,KA=0.25^2),
                                  notfixed_bpop=c(1,1,1,1,1,0),
                                  xt=c( 48,50,55,65,70,85,90,120), 
                                  minxt=0,
                                  maxxt=144,
                                  discrete_xt = list(0:144),
                                  a=c(DOSE=100,TAU=24),
                                  maxa=c(DOSE=1000,TAU=24),
                                  mina=c(DOSE=0,TAU=8),
                                  discrete_a = list(DOSE=seq(0,1000,by=100),TAU=8:24))

## ----simulate_ODE_model--------------------------------------------------
plot_model_prediction(poped.db)

## ------------------------------------------------------------------------
cppFunction('List two_comp_oral_ode_Rcpp(double Time, NumericVector A, NumericVector Pars) {
            int n = A.size();
            NumericVector dA(n);
            
            double CL = Pars[0];
            double V1 = Pars[1];
            double KA = Pars[2];
            double Q  = Pars[3];
            double V2 = Pars[4];
            
            dA[0] = -KA*A[0];
            dA[1] = KA*A[0] - (CL/V1)*A[1] - Q/V1*A[1] + Q/V2*A[2];
            dA[2] = Q/V1*A[1] - Q/V2*A[2];
            return List::create(dA);
            }')

## ---- echo=FALSE---------------------------------------------------------
ff.PK.2.comp.oral.md.ode.Rcpp <- function(model_switch, xt, parameters, poped.db){
  with(as.list(parameters),{
    A_ini <- c(A1=0, A2=0, A3=0)
    times_xt <- drop(xt)
    dose_times = seq(from=0,to=max(times_xt),by=TAU)
    eventdat <- data.frame(var = c("A1"), 
                           time = dose_times,
                           value = c(DOSE), method = c("add"))
    times <- sort(c(times_xt,dose_times))
    
    # Here "PK.2.comp.oral.ode" is equivalent to the pre-compiled version "two_comp_oral_ode_Rcpp".
    out <- ode(A_ini, times, two_comp_oral_ode_Rcpp, parameters,
               events = list(data = eventdat))#atol=1e-13,rtol=1e-13)
    y = out[, "A2"]/(V1/Favail)
    y=y[match(times_xt,out[,"time"])]
    y=cbind(y)
    return(list(y=y,poped.db=poped.db))
  })
}

## ------------------------------------------------------------------------
poped.db.Rcpp <- create.poped.database(poped.db,
                                       ff_fun="ff.PK.2.comp.oral.md.ode.Rcpp")

## ------------------------------------------------------------------------
tic(); eval <- evaluate_design(poped.db); toc()
tic(); eval <- evaluate_design(poped.db.Rcpp); toc()

## ---- echo=FALSE, results="hide"-----------------------------------------
cppFunction('List tmdd_qss_one_target_model_ode
            (double Time, NumericVector A, NumericVector Pars) {
            int n = A.size();
            NumericVector dA(n);
            
            double CL = Pars[0];
            double V1 = Pars[1];
            double Q  = Pars[2];
            double V2 = Pars[3];
            double FAVAIL = Pars[4];
            double KA = Pars[5];
            // double VMAX = Pars[6];
            // double KMSS = Pars[7];
            double R0 = Pars[8];
            double KSSS = Pars[9];
            double KDEG = Pars[10];
            double KINT = Pars[11];
            // double DOSE = Pars[12];
            // double SC_FLAG = Pars[13];
            
            double RTOT, CTOT, CFREE;
            RTOT = A[3];
            CTOT= A[1]/V1;
            CFREE = 0.5*((CTOT-RTOT-KSSS)+sqrt((CTOT-RTOT-KSSS)*(CTOT-RTOT-KSSS)+4*KSSS*CTOT));
            
            dA[0] = -KA*A[0];
            dA[1] = FAVAIL*KA*A[0]+(Q/V2)*A[2]-
            (CL/V1+Q/V1)*CFREE*V1-RTOT*KINT*CFREE*V1/(KSSS+CFREE);
            dA[2] = (Q/V1)*CFREE*V1 - (Q/V2)*A[2];
            dA[3] = R0*KDEG - KDEG*RTOT - (KINT-KDEG)*(RTOT*CFREE/(KSSS+CFREE));
            return List::create(dA);
            }')

## ----echo=FALSE, results="hide"------------------------------------------
sfg <- function(x,a,bpop,b,bocc){
  parameters=c( CL=bpop[1]*exp(b[1])  ,
                V1=bpop[2]*exp(b[2])	,
                Q=bpop[3]*exp(b[3])	,
                V2=bpop[4]*exp(b[4])	,
                FAVAIL=bpop[5]*exp(b[5])	,
                KA=bpop[6]*exp(b[6])	,                       
                VMAX=bpop[7]*exp(b[7])	,
                KMSS=bpop[8]*exp(b[8])	,
                R0=bpop[9]*exp(b[9])	,
                KSSS=bpop[10]*exp(b[10])	,
                KDEG=bpop[11]*exp(b[11])	,
                KINT=bpop[12]*exp(b[12])	,
                DOSE=a[1]	,
                SC_FLAG=a[2])   
  return(parameters) 
}

## ------------------------------------------------------------------------
tmdd_qss_one_target_model_compiled <- function(model_switch,xt,parameters,poped.db){
  with(as.list(parameters),{
    y=xt
    
    #The initialization vector for the compartment
    A_ini <- c(A1=DOSE*SC_FLAG,
               A2=DOSE*(1-SC_FLAG),
               A3=0,
               A4=R0)
    
    #Set up time points for the ODE
    times_xt <- drop(xt)
    times <- sort(times_xt)
    times <- c(0,times) ## add extra time for start of integration
    
    # solve the ODE
    out <- ode(A_ini, times, tmdd_qss_one_target_model_ode, parameters)#,atol=1e-13,rtol=1e-13)
    
    
    # extract the time points of the observations
    out = out[match(times_xt,out[,"time"]),]
    
    # Match ODE output to measurements
    RTOT = out[,"A4"]
    CTOT = out[,"A2"]/V1
    CFREE = 0.5*((CTOT-RTOT-KSSS)+sqrt((CTOT-RTOT-KSSS)^2+4*KSSS*CTOT))
    COMPLEX=((RTOT*CFREE)/(KSSS+CFREE))
    RFREE= RTOT-COMPLEX
    
    y[model_switch==1]= RTOT[model_switch==1]
    y[model_switch==2] =CFREE[model_switch==2]
    #y[model_switch==3]=RFREE[model_switch==3]
    
    return(list( y=y,poped.db=poped.db))
  })
}


## ---- echo=FALSE---------------------------------------------------------
tmdd_qss_one_target_model_ruv <- function(model_switch,xt,parameters,epsi,poped.db){
  returnArgs <- do.call(poped.db$model$ff_pointer,list(model_switch,xt,parameters,poped.db))
  y <- returnArgs[[1]]
  poped.db <- returnArgs[[2]]
  
  y[model_switch==1] = log(y[model_switch==1])+epsi[,1]
  y[model_switch==2] = log(y[model_switch==2])+epsi[,2]
  
  return(list(y=y,poped.db=poped.db))
}


## ------------------------------------------------------------------------
xt <- zeros(6,30)
study_1_xt <- matrix(rep(c(0.0417,0.25,0.5,1,3,7,14,21,28,35,42,49,56),8),nrow=4,byrow=TRUE)
study_2_xt <- matrix(rep(c(0.0417,1,1,7,14,21,28,56,63,70,77,84,91,98,105),4),nrow=2,byrow=TRUE)
xt[1:4,1:26] <- study_1_xt
xt[5:6,] <- study_2_xt

model_switch <- zeros(6,30)
model_switch[1:4,1:13] <- 1
model_switch[1:4,14:26] <- 2
model_switch[5:6,1:15] <- 1
model_switch[5:6,16:30] <- 2

G_xt <- zeros(6,30)
study_1_G_xt <- matrix(rep(c(1:13),8),nrow=4,byrow=TRUE)
study_2_G_xt <- matrix(rep(c(14:28),4),nrow=2,byrow=TRUE)
G_xt[1:4,1:26] <- study_1_G_xt
G_xt[5:6,] <- study_2_G_xt

## ------------------------------------------------------------------------
poped.db.2 <- create.poped.database(ff_fun=tmdd_qss_one_target_model_compiled,
                                    fError_fun=tmdd_qss_one_target_model_ruv,
                                    fg_fun=sfg,
                                    groupsize=rbind(6,6,6,6,100,100),
                                    m=6,      #number of groups
                                    sigma=c(0.04,0.0225),
                                    bpop=c(CL=0.3,V1=3,Q=0.2,V2=3,FAVAIL=0.7,KA=0.5,VMAX=0,
                                           KMSS=0,R0=0.1,KSSS=0.015,KDEG=10,KINT=0.05),
                                    d=c(CL=0.09,V1=0.09,Q=0.04,V2=0.04,FAVAIL=0.04,
                                        KA=0.16,VMAX=0,KMSS=0,R0=0.09,KSSS=0.09,KDEG=0.04,
                                        KINT=0.04),
                                    notfixed_bpop=c( 1,1,1,1,1,1,0,0,1,1,1,1),
                                    notfixed_d=c( 1,1,1,1,1,1,0,0,1,1,1,1),
                                    xt=xt,
                                    model_switch=model_switch,
                                    ni=rbind(26,26,26,26,30,30),
                                    bUseGrouped_xt=1,
                                    G_xt=G_xt,
                                    a=list(c(DOSE=100, SC_FLAG=0),
                                           c(DOSE=300, SC_FLAG=0),
                                           c(DOSE=600, SC_FLAG=0),
                                           c(DOSE=1000, SC_FLAG=1),
                                           c(DOSE=600, SC_FLAG=0),
                                           c(DOSE=1000, SC_FLAG=1)),
                                    discrete_a = list(DOSE=seq(100,1000,by=100),
                                                      SC_FLAG=c(0,1)))


## ----simulate_different_dose_regimen-------------------------------------
plot_model_prediction(poped.db.2,facet_scales="free")

## ------------------------------------------------------------------------
eval_2 <- evaluate_design(poped.db.2);
round(eval_2$rse)

## ----model---------------------------------------------------------------
mod_1 <- function(model_switch,xt,parameters,poped.db){
  with(as.list(parameters),{
    y=xt
    
    CL=CL*(WT/70)^(WT_CL)
    V=V*(WT/70)^(WT_V)
    DOSE=1000*(WT/70)
    y = DOSE/V*exp(-CL/V*xt) 
    
    return(list( y= y,poped.db=poped.db))
  })
}

par_1 <- function(x,a,bpop,b,bocc){
  parameters=c( CL=bpop[1]*exp(b[1]),
                V=bpop[2]*exp(b[2]),
                WT_CL=bpop[3],
                WT_V=bpop[4],
                WT=a[1])
  return( parameters ) 
}

## ----design--------------------------------------------------------------
poped_db <- 
  create.poped.database(
    ff_fun=mod_1,
    fg_fun=par_1,
    fError_fun=feps.add.prop,
    groupsize=50,
    m=1,
    sigma=c(0.015,0.0015),
    notfixed_sigma = c(1,0),
    bpop=c(CL=3.8,V=20,WT_CL=0.75,WT_V=1), 
    d=c(CL=0.05,V=0.05), 
    xt=c( 1,2,4,6,8,24),
    minxt=0,
    maxxt=24,
    bUseGrouped_xt=1,
    a=c(WT=70)
  )

## ------------------------------------------------------------------------
plot_model_prediction(poped_db)

## ------------------------------------------------------------------------
evaluate_design(poped_db)

## ------------------------------------------------------------------------
poped_db_2 <- 
  create.poped.database(
    ff_fun=mod_1,
    fg_fun=par_1,
    fError_fun=feps.add.prop,
    groupsize=1,
    m=50,
    sigma=c(0.015,0.0015),
    notfixed_sigma = c(1,0),
    bpop=c(CL=3.8,V=20,WT_CL=0.75,WT_V=1), 
    d=c(CL=0.05,V=0.05), 
    xt=c(1,2,4,6,8,24),
    minxt=0,
    maxxt=24,
    bUseGrouped_xt=1,
    a=as.list(rnorm(50, mean = 70, sd = 10))
  )


## ------------------------------------------------------------------------
ev <- evaluate_design(poped_db_2) 
round(ev$ofv,1)
round(ev$rse)

## ------------------------------------------------------------------------
sfg <- function(x,a,bpop,b,bocc){
  parameters=c( CL_OCC_1=bpop[1]*exp(b[1]+bocc[1,1]),
                CL_OCC_2=bpop[1]*exp(b[1]+bocc[1,2]),
                V=bpop[2]*exp(b[2]),
                KA=bpop[3]*exp(b[3]),
                DOSE=a[1],
                TAU=a[2])
  return( parameters ) 
}

## ------------------------------------------------------------------------
cppFunction(
  'List one_comp_oral_ode(double Time, NumericVector A, NumericVector Pars) {
   int n = A.size();
   NumericVector dA(n);
            
   double CL_OCC_1 = Pars[0];
   double CL_OCC_2 = Pars[1];
   double V = Pars[2];
   double KA = Pars[3];
   double TAU = Pars[4];
   double N,CL;
            
   N = floor(Time/TAU)+1;
   CL = CL_OCC_1;
   if(N>6) CL = CL_OCC_2;
   
   dA[0] = -KA*A[0];
   dA[1] = KA*A[0] - (CL/V)*A[1];
   return List::create(dA);
   }'
)

ff.ode.rcpp <- function(model_switch, xt, parameters, poped.db){
  with(as.list(parameters),{
    A_ini <- c(A1=0, A2=0)
    times_xt <- drop(xt) #xt[,,drop=T] 
    dose_times = seq(from=0,to=max(times_xt),by=TAU)
    eventdat <- data.frame(var = c("A1"), 
                           time = dose_times,
                           value = c(DOSE), method = c("add"))
    times <- sort(c(times_xt,dose_times))
    out <- ode(A_ini, times, one_comp_oral_ode, c(CL_OCC_1,CL_OCC_2,V,KA,TAU), 
               events = list(data = eventdat))#atol=1e-13,rtol=1e-13)
    y = out[, "A2"]/(V)
    y=y[match(times_xt,out[,"time"])]
    y=cbind(y)
    return(list(y=y,poped.db=poped.db))
  })
}


## ------------------------------------------------------------------------
poped.db <- 
  create.poped.database(
    ff_fun=ff.ode.rcpp,
    fError_fun=feps.add.prop,
    fg_fun=sfg,
    bpop=c(CL=3.75,V=72.8,KA=0.25), 
    d=c(CL=0.25^2,V=0.09,KA=0.09), 
    sigma=c(0.04,5e-6),
    notfixed_sigma=c(0,0),
    docc = matrix(c(0,0.09,0),nrow = 1),
    m=2,
    groupsize=20,
    xt=c( 1,2,8,240,245),
    minxt=c(0,0,0,240,240),
    maxxt=c(10,10,10,248,248),
    bUseGrouped_xt=1,
    a=list(c(DOSE=20,TAU=24),c(DOSE=40, TAU=24)),
    maxa=c(DOSE=200,TAU=24),
    mina=c(DOSE=0,TAU=24)
  )

## ----simulate_IOV_without_IIV--------------------------------------------
plot_model_prediction(poped.db, model_num_points = 300)

## ----simulate_IOV_with_IIV-----------------------------------------------
library(ggplot2)
set.seed(123)
plot_model_prediction(poped.db, PRED=F,IPRED=F, 
                      separate.groups=T, model_num_points = 300, 
                      groupsize_sim = 1,
                      IPRED.lines = T, alpha.IPRED.lines=0.6,
                      sample.times = F
                      ) + geom_vline(xintercept = 24*6,color="red")

## ------------------------------------------------------------------------
ev <- evaluate_design(poped.db)
round(ev$rse)

## ---- echo=FALSE, results="hide"-----------------------------------------
ff <- function(model_switch,xt,parameters,poped.db){
  ##-- Model: One comp first order absorption
  with(as.list(parameters),{
    y=xt
    y=(DOSE*Favail*KA/(V*(KA-CL/V)))*(exp(-CL/V*xt)-exp(-KA*xt))
    return(list(y=y,poped.db=poped.db))
  })
}

sfg <- function(x,a,bpop,b,bocc){
  ## -- parameter definition function 
  parameters=c(CL=bpop[1]*exp(b[1]),
               V=bpop[2]*exp(b[2]),
               KA=bpop[3]*exp(b[3]),
               Favail=bpop[4],
               DOSE=a[1])
  return(parameters) 
}

feps <- function(model_switch,xt,parameters,epsi,poped.db){
  ## -- Residual Error function
  ## -- Proportional 
  returnArgs <- ff(model_switch,xt,parameters,poped.db) 
  y <- returnArgs[[1]]
  poped.db <- returnArgs[[2]]
  y = y*(1+epsi[,1])
  
  return(list(y=y,poped.db=poped.db)) 
}

## -- Define initial design  and design space
poped.db <- 
  create.poped.database(
    ff_file="ff",
    fg_file="sfg",
    fError_file="feps",
    bpop=c(CL=0.15, V=8, KA=1.0, Favail=1), 
    notfixed_bpop=c(1,1,1,0),
    d=c(CL=0.07, V=0.02, KA=0.6), 
    sigma=0.01,
    groupsize=32,
    xt=c( 0.5,1,2,6,24,36,72,120),
    minxt=0,
    maxxt=120,
    a=70
  )

## ------------------------------------------------------------------------
poped.db_with <- 
  create.poped.database(
    ff_file="ff",
    fg_file="sfg",
    fError_file="feps",
    bpop=c(CL=0.15, V=8, KA=1.0, Favail=1), 
    notfixed_bpop=c(1,1,1,0),
    d=c(CL=0.07, V=0.02, KA=0.6), 
    covd = c(.03,.1,.09),
    sigma=0.01,
    groupsize=32,
    xt=c( 0.5,1,2,6,24,36,72,120),
    minxt=0,
    maxxt=120,
    a=70
  )

## ------------------------------------------------------------------------
(IIV <- poped.db_with$parameters$param.pt.val$d)
cov2cor(IIV)

## ----simulate_with_cov_matrix--------------------------------------------
library(ggplot2)
p1 <- plot_model_prediction(poped.db,IPRED=TRUE)+ylim(0,12)
p2 <- plot_model_prediction(poped.db_with,IPRED=TRUE) +ylim(0,12)
gridExtra::grid.arrange(p1, p2, nrow = 1)

## ------------------------------------------------------------------------
ev1 <- evaluate_design(poped.db)
ev2 <- evaluate_design(poped.db_with)

## ------------------------------------------------------------------------
round(ev1$ofv,1)
round(ev1$rse,1)

round(ev2$ofv,1)
round(ev2$rse,1)

## ------------------------------------------------------------------------
ev1 <- evaluate_design(poped.db, fim.calc.type=0)
ev2 <-evaluate_design(poped.db_with, fim.calc.type=0)

round(ev1$ofv,1)
round(ev1$rse,1)

round(ev2$ofv,1)
round(ev2$rse,1)

## ---- echo = FALSE-------------------------------------------------------
##-- Model: One comp first order absorption
## -- Analytic solution for both mutiple and single dosing
ff <- function(model_switch,xt,parameters,poped.db){
  with(as.list(parameters),{
    y=xt
    N = floor(xt/TAU)+1
    y=(DOSE*Favail/V)*(KA/(KA - CL/V)) * 
      (exp(-CL/V * (xt - (N - 1) * TAU)) * (1 - exp(-N * CL/V * TAU))/(1 - exp(-CL/V * TAU)) - 
         exp(-KA * (xt - (N - 1) * TAU)) * (1 - exp(-N * KA * TAU))/(1 - exp(-KA * TAU)))  
    return(list( y=y,poped.db=poped.db))
  })
}

## ------------------------------------------------------------------------
## -- parameter definition function 
## -- names match parameters in function ff
## -- note, covariate on clearance for pediatrics (using isPediatric x[1])
sfg <- function(x,a,bpop,b,bocc){
  parameters=c( 
    V=bpop[1]*exp(b[1]),
    KA=bpop[2]*exp(b[2]),
    CL=bpop[3]*exp(b[3])*bpop[5]^x[1], # add covariate for pediatrics
    Favail=bpop[4],
    isPediatric = x[1],
    DOSE=a[1],
    TAU=a[2])
  return( parameters ) 
}

## ---- echo=FALSE, results="hide"-----------------------------------------
## -- Residual unexplained variablity (RUV) function
## -- Additive + Proportional  
feps <- function(model_switch,xt,parameters,epsi,poped.db){
  returnArgs <- do.call(poped.db$model$ff_pointer,list(model_switch,xt,parameters,poped.db)) 
  y <- returnArgs[[1]]
  poped.db <- returnArgs[[2]]
  
  y = y*(1+epsi[,1])+epsi[,2]
  
  return(list( y= y,poped.db =poped.db )) 
}

## ------------------------------------------------------------------------
poped.db <- 
  create.poped.database(
    ff_fun="ff",
    fg_fun="sfg",
    fError_fun="feps",
    bpop=c(V=72.8,KA=0.25,CL=3.75,Favail=0.9,pedCL=0.8), 
    notfixed_bpop=c(1,1,1,0,1),
    d=c(V=0.09,KA=0.09,CL=0.25^2), 
    sigma=c(0.04,5e-6),
    notfixed_sigma=c(0,0),
    m=2,
    groupsize=20,
    xt=c( 1,8,10,240,245),
    bUseGrouped_xt=1,
    a=list(c(DOSE=20,TAU=24),c(DOSE=40, TAU=24)),
    x=list(isPediatric = 0)
  )

## ----simulate_adult------------------------------------------------------
plot_model_prediction(poped.db, model_num_points = 300)


## ------------------------------------------------------------------------
(outAdult = evaluate_design(poped.db))

## ------------------------------------------------------------------------
evaluate_design(create.poped.database(poped.db, notfixed_bpop=c(1,1,1,0,0)))

## ------------------------------------------------------------------------
poped.db.ped <- 
  create.poped.database(
    ff_fun="ff",
    fg_fun="sfg",
    fError_fun="feps",
    bpop=c(V=72.8,KA=0.25,CL=3.75,Favail=0.9,pedCL=0.8), 
    notfixed_bpop=c(1,1,1,0,1),
    d=c(V=0.09,KA=0.09,CL=0.25^2), 
    sigma=c(0.04,5e-6),
    notfixed_sigma=c(0,0),
    m=1,
    groupsize=6,
    xt=c( 1,2,6,240),
    bUseGrouped_xt=1,
    a=list(c(DOSE=40,TAU=24)),
    x=list(isPediatric = 1)
  )

## ----simulate_pediatrics-------------------------------------------------
plot_model_prediction(poped.db.ped, model_num_points = 300)

## ------------------------------------------------------------------------
evaluate_design(poped.db.ped)


## ------------------------------------------------------------------------
poped.db.all <- create.poped.database(
  poped.db.ped,
  prior_fim = outAdult$fim
)

(out.all <- evaluate_design(poped.db.all))

## ------------------------------------------------------------------------
evaluate_power(poped.db.all, bpop_idx=5, h0=1, out=out.all)

## ---- echo=FALSE, results="hide"-----------------------------------------
sfg <- function(x,a,bpop,b,bocc){
  ## -- parameter definition function 
  parameters=c(
    CL=bpop[1]*exp(b[1]),
    V=bpop[2]*exp(b[2]),
    KA=bpop[3]*exp(b[3]),
    Favail=bpop[4],
    DOSE=a[1])
  return(parameters) 
}

ff <- function(model_switch,xt,parameters,poped.db){
  ##-- Model: One comp first order absorption
  with(as.list(parameters),{
    y=xt
    y=(DOSE*Favail*KA/(V*(KA-CL/V)))*(exp(-CL/V*xt)-exp(-KA*xt))
    return(list(y=y,poped.db=poped.db))
  })
}

## ------------------------------------------------------------------------
bpop_vals <- c(CL=0.15, V=8, KA=1.0, Favail=1)
bpop_vals_ed <- 
  cbind(ones(length(bpop_vals),1)*4, # log-normal distribution
        bpop_vals,
        ones(length(bpop_vals),1)*(bpop_vals*0.1)^2) # 10% of bpop value
bpop_vals_ed["Favail",] <- c(0,1,0)
bpop_vals_ed

## ------------------------------------------------------------------------
poped.db <- 
  create.poped.database(
    ff_fun=ff,
    fg_fun=sfg,
    fError_fun=feps.add.prop,
    bpop=bpop_vals_ed, 
    notfixed_bpop=c(1,1,1,0),
    d=c(CL=0.07, V=0.02, KA=0.6), 
    sigma=c(0.01,0.25),
    groupsize=32,
    xt=c( 0.5,1,2,6,24,36,72,120),
    minxt=0,
    maxxt=120,
    a=70,
    mina=0,
    maxa=100,
    ED_samp_size=20)

## ------------------------------------------------------------------------
tic();evaluate_design(poped.db,d_switch=FALSE,ED_samp_size=20); toc()

## ------------------------------------------------------------------------
tic();evaluate_design(poped.db,d_switch=FALSE,ED_samp_size=20); toc()

## ---- echo=FALSE---------------------------------------------------------
sfg <- function(x,a,bpop,b,bocc){
  ## -- parameter definition function 
  parameters=c(CL=bpop[1]*exp(b[1]),
               V=bpop[2]*exp(b[2]),
               KA=bpop[3]*exp(b[3]),
               Favail=bpop[4],
               DOSE=a[1])
  return(parameters) 
}

ff <- function(model_switch,xt,parameters,poped.db){
  ##-- Model: One comp first order absorption
  with(as.list(parameters),{
    y=xt
    y=(DOSE*Favail*KA/(V*(KA-CL/V)))*(exp(-CL/V*xt)-exp(-KA*xt))
    return(list(y=y,poped.db=poped.db))
  })
}

## ------------------------------------------------------------------------
poped.db <- 
  create.poped.database(
    ff_fun=ff,
    fg_fun=sfg,
    fError_fun=feps.add.prop,
    bpop=c(CL=0.15, V=8, KA=1.0, Favail=1), 
    notfixed_bpop=c(1,1,1,0),
    d=c(CL=0.07, V=0.02, KA=0.6), 
    sigma=c(0.01,0.25),
    groupsize=32,
    xt=c( 0.5,1,2,6,24,36,72,120),
    minxt=0,
    maxxt=120,
    a=70,
    mina=0,
    maxa=100,
    ds_index=c(0,0,0,1,1,1,1,1), # size is number_of_non_fixed_parameters
    ofv_calc_type=6) # Ds OFV calculation

## ------------------------------------------------------------------------
evaluate_design(poped.db)

## ---- echo=FALSE, results="hide"-----------------------------------------

sfg <- function(x,a,bpop,b,bocc){
  ## -- parameter definition function 
  parameters=c(KA=bpop[1]*exp(b[1]),
               K=bpop[2]*exp(b[2]),
               V=bpop[3]*exp(b[3]),
               DOSE=a[1])
  return(parameters) 
}

ff <- function(model_switch,xt,parameters,poped.db){
  ##-- Model: One comp first order absorption
  with(as.list(parameters),{
    y<-(DOSE/V*KA/(KA-K)*(exp(-K*xt)-exp(-KA*xt)))
    return(list(y=y,poped.db=poped.db))
  })
}

## -- Residual unexplained variablity (RUV) function
## -- Additive + Proportional  
feps <- function(model_switch,xt,parameters,epsi,poped.db){
  returnArgs <- do.call(poped.db$model$ff_pointer,
                        list(model_switch,xt,parameters,poped.db)) 
  f <- returnArgs[[1]]
  y = f + (0.5 + 0.15*f)*epsi[,1]
  
  return(list( y= y,poped.db =poped.db )) 
}
## -- Define initial design  and design space
poped.db <- create.poped.database(
  ff_fun=ff,
  fg_fun=sfg,
  fError_fun =feps,
  bpop=c(KA=2, K=0.25, V=15), 
  d=c(KA=1, V=0.25,0.1), 
  sigma=c(1),
  notfixed_sigma = c(0),
  groupsize=1,
  xt=c( 1,3,8),
  minxt=0,
  maxxt=10,
  a=100)

plot_model_prediction(poped.db)


## ------------------------------------------------------------------------
shrinkage(poped.db)

