#' Installation Test
#' =========================================================================================
#' This example script is intended to 
#' to test installation of software and tools 
#' using the DDMoRe Standalone Execution Environment (SEE). 
#'   
#' The following steps are implemented in this workflow: 
#' * Initialisation of the R console
#' * Exploratory graphical analysis
#' * Parameter estimation with Monolix and model evaluation in Xpose
#' * Parameter estimation in NONMEM using FOCE and model evaluation in Xpose
#' * Parameter estimation in Winbugs and MCMC diagnostics using coda
#' * Updating parameter estimates in the MDL Parameter Object 
#'   using MLE values from NONMEM
#' * Performing a Visual Predictive Check (VPC) in PsN
#' * Using MLE values from NONMEM to simulate new observed values using the simulx function
#'   in the mlxR package.
#' * Updating parameter estimates with NONMEM and simulx output
#' * Performing evaluation of study design characteristics using PFIM and PopED

#' 
#' Initialisation
#' =========================
#' Clear workspace and set working directory 
#+ Initialisation
rm(list=ls(all=F))
getwd()
#' Some useful functions for working with the SEE and workspaces in the MDL-IDE
#' A system variable has been set which retrieves the directory of the workspaces
Sys.getenv("MDLIDE_WORKSPACE_HOME")
#' By typing "~" you can get to this directory
mydir <- file.path("~/InstallationTest/models")
setwd(mydir)
#' Another function from the ddmore R package retrieves the directory that the SEE is installed in.
#' You can then navigate relative to these directories.
ddmore:::DDMORE.checkConfiguration()

#' Set name of .mdl file 
mdlfile <- "InstallationTest_COMPARTMENT.mdl"
model <- tools::file_path_sans_ext(mdlfile)

#' Read the different MDL objects needed for upcoming tasks.
#' An MDL file may contain more than one object of any type 
#' (though typically only one model!) so these functions return
#' a list of all objects of that type found in the target MDL file.
#' To pick out the first, use the double square bracket notation in R
#' to retrieve the first item of the list. 
myModelObj <- getModelObjects(mdlfile)[[1]] 
myDataObj <- getDataObjects(mdlfile)[[1]]
myParameterObj <- getParameterObjects(mdlfile)[[1]]
myPriorObj <- getPriorObjects(mdlfile)[[1]]
myDesignObj <- getDesignObjects(mdlfile)[[1]]

#' You can also refer  to objects by their name.
myTaskPropertiesObj_saem <- getTaskPropertiesObjects(mdlfile)[["SAEM_task"]]
myTaskPropertiesObj_focei <- getTaskPropertiesObjects(mdlfile)[["NONMEM_task"]]
myTaskPropertiesObj_mcmc <- getTaskPropertiesObjects(mdlfile)[["BUGS_task"]]
myTaskPropertiesObj_Evaluate<- getTaskPropertiesObjects(mdlfile)[["EVALUATE_task"]]
#myTaskPropertiesObj_Optimise <- getTaskPropertiesObjects(mdlfile)[["warfarin_Optimise_Explicit"]]
 
#' Exploratory Data Analysis
#' =========================
#' getDataObjects reads and parses the MDL. It doesn't actually read the data file.
#' To do so, use the ddmore function readDataObj() to create an R object (data frame) 
#' based on the MDL data object. This uses data column headers defined in the Data Object.
myData <- readDataObj(myDataObj)

#' Let's look at the first 6 lines of the data set
head(myData)
myEDAData <- myData[myData$MDV==0,]

library(ggplot2)
#' Plot the data using ggplot2 library 
plot1 <- ggplot(aes(y=DV,x=TIME, group=ID), data=myEDAData) +
				geom_line()+
				geom_point() + 
				labs(y="Concentration" ,x="Time (hours)")
#+ ExploratoryAnalysis
print(plot1)


#' Model Development
#' =========================

#' ESTIMATE model parameters using Monolix
#' -------------------------
#' You can pass MDL files (or PharmML files) into ddmore functions
#' to perform tasks. But you can also define new Model Object Groups (MOGs)
#' and then write these to file. This allows you to combine different elements
#' from the MOG to do specific tasks.
#'  
#' To define the new MOG for estimation in Monolix:
MLX.MOG <- createMogObj(dataObj = myDataObj, 
		parObj = myParameterObj, 
		mdlObj = myModelObj, 
		taskObj = myTaskPropertiesObj_saem)

#' We can then write the MOG back out to an .mdl file.
mdlfile.MLX <- paste0(model,"_MLX.mdl")
writeMogObj(MLX.MOG,mdlfile.MLX)

#' The ddmore "estimate" function translates the contents of the .mdl file to 
#' MLXTRAN and then estimates parameters using Monolix. After estimation, the output 
#' from Monolix is converted to a Standardised Output object which is saved in a .SO.xml file.
#'   
#' Translated files and Monolix output will be returned in the ./Monolix subfolder.
#' The Standardised Output object (.SO.xml) is read and parsed into an R object called "mlx" 
#' of (S4) class "StandardOutputObject".
#+ Monolix, cache=TRUE
MLX <- estimate(mdlfile.MLX, target="MONOLIX", subfolder="Monolix")
slotNames(MLX)

#' Retrieve estimated parameters
#' ------------------------------
#' The ddmore function getPopulationParameters extracts the Population Parameter values from 
#' an R object of (S4) class "StandardOutputObject" and returns the MLE estimates. 
#' See documentation for getPopulationParameters to see other arguments and settings for this function.  
getPopulationParameters(MLX)$MLE

parameters_mlx <- getPopulationParameters(MLX, what="estimates")$MLE
parameters_mlx

#' Other functions exist to extract information from the SO, but the SO is also available as an R
#' object and so it is relatively straightforward to retrieve elements directly using R commands.
#'  
#' Perform basic model diagnostics 
#' --------------------------------
#' Use ddmore function as.xpdb() to create an Xpose database object from an R object of (S4) class "StandardOutputObject".
#' We can then call Xpose functions referencing this mlx.xpdb object as the input.
datafile <- myDataObj@SOURCE$srcfile$file
mlx.xpdb<-as.xpdb(MLX,datafile)

#+ Basic_Goodness_of_Fit_following_Monolix_Estimation
print(basic.gof(mlx.xpdb))
mlx.xpdb@Data$ID <- as.numeric(as.character(mlx.xpdb@Data$ID))
#+ Individual_Plots_following_Monolix_Estimation
print(ind.plots(mlx.xpdb))

#' ESTIMATE model parameters in NONMEM
#' -------------------------

#' Assemble the new Modelling Object Group (MOG) using Nonmem task properties.
#' Note that the only change is in the Task Properties Object.
#' All other MDL objects remain the same.
NM.MOG <- createMogObj(dataObj = myDataObj, 
		parObj = myParameterObj, 
		mdlObj = myModelObj, 
		taskObj = myTaskPropertiesObj_focei)

#' We can then write the MOG back out to an .mdl file as before.
mdlfile.NM <- paste0(model,"_NM.mdl")
writeMogObj(NM.MOG,mdlfile.NM)

#' The ddmore "estimate" function translates the contents of the .mdl file to 
#' NMTRAN and then estimates parameters using NONMEM. After estimation, the output 
#' from NONMEM is converted to a Standardised Output object which is saved in a .SO.xml file.
#'   
#' Translated files and NONMEM output will be returned in the ./NONMEM subfolder.
#' The Standardised Output object (.SO.xml) is read and parsed into an R object called "nm" 
#' of (S4) class "StandardOutputObject".
#+ NONMEM, cache=TRUE
NM <- estimate(mdlfile.NM, target="NONMEM", subfolder="NONMEM")

#' Once the model has run and the Standard Output (SO) object has returned
#' it is possible to read that object in to R at any time.
# NM <-LoadSOObject("NONMEM/Simeoni_PAGE_Estimation_NM.SO.xml")

#' Results from NONMEM should be comparable to previous results
parameters_nm <- getPopulationParameters(NM,  what="estimates")$MLE
parameters_nm
            
#' ESTIMATE model parameters in WINBUGS
#' -------------------------
#' Assembling the new MOG for Winbugs. Note that we reuse the data and model. 
#' Two main changes are made: selecting the appropriate Task Properties Object for BUGS
#' and exchanging the Prior Object instead of the Parameter Object.
BUGS.MOG <- createMogObj(dataObj = myDataObj, 
		priorObj = myPriorObj, 
		mdlObj = myModelObj, 
		taskObj = myTaskPropertiesObj_mcmc)

#' We can then write the MOG back out to an .mdl file.
mdlfile.BUGS <- paste0(model,"_BUGS.mdl")
writeMogObj(BUGS.MOG,mdlfile.BUGS)
#mdlfile.BUGS <- "Simeoni_PAGE_Estimation_BUGS.mdl"

#' Winbugs converter and function call has not been integrated into the SEE yet,
#' so we cannot run:
#+ BUGSestimate, eval=FALSE
# estimate(mdlfile.BUGS, target="BUGS")
#' in the meantime, an R script has been prepared to call the converter, generate Winbugs input. 
#' This involves conversion of the NONMEM format dataset to BUGS (R list) format and preparation
#' of the model file, data file, initial values for the MCMC and creation of the script file
#' for running BUGS in batch mode.
#+ BUGS, cache=TRUE

source(file.path(ddmore:::DDMORE.checkConfiguration(),"BlackBoxWinBUGS","ddmorewinbugs.R"))
BUGS<-run.winbugs(mdlfile.BUGS,subfolder="WinBUGSNEW")

getPopulationParameters(BUGS)$Bayesian
getPopulationParameters(BUGS, what="estimates")$Bayesian # o other options,  what= precisions, intervals, all
# it is possible also to omitt $Bayesian
getIndividualParameters(BUGS)

#' Here we use the coda package to examine and assess MCMC convergence.
library(coda)

#' By default the parameters monitored by BUGS are all parameters specified with
#' random variable distributions in the PRIOR_VARIABLE_DEFINITION. In this example
#' these are specified on the log scale. The following code transforms them back to
#' the natural scale. Note that in the public release, MCMC nodes for monitoring 
#' will be specified via the Task Properties Object.

parameters_BUGS<- getPopulationParameters(BUGS, what="estimates")$Bayesian

#' coda MCMC trace and density plots
#+ MCMC_Diagnostics
#to read only one chain
BUGSoutputPath <- file.path(getwd(),"WinBUGS")
coda_out <- read.coda(output.file=file.path(BUGSoutputPath,"output1.txt"), index.file=file.path(BUGSoutputPath,"outputIndex.txt"))

#to read more chains (3 in this case) use read.coda.interactive()
#coda_out <- read.coda.interactive()
# Enter in the order:
  #outputIndex.txt
  #output1.txt  
  #output2.txt  
  #output2.txt
## summary statistics
summary(coda_out)
## trace and density plots
plot(coda_out, ask=F)
## gelman-rubin diagnostics (if at least 2 chains)
#gelman.diag(coda_out)
#gelman.plot(coda_out)

#' Compare estimated parameters in the 3 tools
#' ----------------------------------

estPars <- intersect(rownames(parameters_BUGS), names(parameters_nm))
cbind(BUGS=parameters_BUGS[estPars,"Median"],
		NM=parameters_nm[estPars],
		MLX=parameters_mlx[estPars]) 

#' VPC of model 
#' --------------------------------
#' Update the (initial) values in the MDL Parameter Object
#' Extract the structural and variability parameters from the SO object
#' following NONMEM estimation.
structuralPar <- getPopulationParameters(NM, what="estimates",block='structural')$MLE
variabilityPar <- getPopulationParameters(NM, what="estimates",block='variability')$MLE

#' Update the parameter object using the ddmore "updateParObj" function.
#' This function updates an R object of (S4) class "parObj".
#' The user chooses which block to update, what items within that block,
#' and what to replace those items with. 
#' NOTE: that updateParObj can only update attributes which ALREADY EXIST in the MDL Parameter Object for that item.
#' This ensures that valid MDL is preserved.   

myParObjUpdated <- updateParObj(myParameterObj,block="STRUCTURAL",
		item=names(structuralPar),
		with=list(value=structuralPar))

myParObjUpdated <- updateParObj(myParObjUpdated,block="VARIABILITY",
		item=names(variabilityPar),
		with=list(value=variabilityPar))

#' Check that the appropriate initial values have been updated to the MLE values from the previous fit.
#print(myParObjUpdated@STRUCTURAL)
#print(myParObjUpdated@VARIABILITY)

#' Assembling the new MOG. Note that we reuse the data, model and tasks from the previous NONMEM run. 
VPC.MOG <- createMogObj(dataObj = myDataObj, 
		parObj = myParObjUpdated, 
		mdlObj = myModelObj, 
		taskObj = myTaskPropertiesObj_focei)

#' We can then write the MOG back out to an .mdl file.
mdlfile.VPC <- paste0(model,"_VPC.mdl")
writeMogObj(VPC.MOG,mdlfile.VPC)

#' Similarly as above, ddmore "VPC.PsN" function can be used to run a VPC using PsN as target tool
#+ VPC, cache=TRUE
vpcFiles <- VPC.PsN(mdlfile.VPC,samples=200, seed=12345,
		vpcOptions ="-n_simulation=10 -auto_bin=5, -min_points_in_bin=2",
		subfolder="VPC", plot=FALSE) 

#' To replay the visualisation using information from the VPC SO file
#+ VPC_using_NONMEM_estimates
plot(xpose.VPC(vpc.info= file.path("./VPC","vpc_results.csv"),
				vpctab= file.path("./VPC","vpctab"),
				main="VPC TGI model"))

#' Simulation using simulx
#' --------------------------------
#' The mlxR package has been developed to visualize and explore models that are encoded in
#' MLXTRAN or PharmML.
#' The ddmore function as.PharmML translates an MDL file (extension .mdl) to its PharmML
#' representation. The output file (extension .xml) is saved in the working directory.

myPharmML <- as.PharmML(mdlfile)
model.mlxtran <- pharmml2mlxtran(myPharmML, "InstallationTest_Model.txt")

library(mlxR)
#' mlxR contains a function pharmml2mlxtran(pharmMLFile) that performs the appropriate translation
#' if needed.

#model.mlxtran <- file.path(getwd(),"Monolix","InstallationTest_MLX","InstallationTest_MLX_model.txt")

#' In future, we will be able to read Design Object information from MDL directly
#' and create inputs to simulx by specifying as.simulx(mdlfile). For now though,
#' we create the relevant inputs for simulx by hand.
p<-parameters_nm
print(p)

#' adm : dose administration times
adm <- list( time  = 0, 
              amount = 100, 
              target = 'INPUT_CENTRAL')

#' f : what to sample ("observe") and at what times?
f <- list( name='C', time=seq(0,24,by=0.1))
y <- list( name='Y', time=c(0.001, 0.5, 1, 4, 8 , 12, 24))

#' g : define group sizes and assign treatments.
g <- list(size=12, level="longitudinal", treatment=adm)

#' Call simulx.
#' simulx can take PharmML as an input - so the model does NOT need to be translated to MLXTRAN.

res <- simulx( model     = myPharmML,
               parameter = p,
               group = g,
               output    = list(f,y))

#' Plot simulx output
#+ Plot_of_simulated_data_using_NONMEM_estimates
ggplot() + 
		geom_line(data=res$WTOT, aes(x=time, y=WTOT, colour=id))+
		geom_point(data=res$Y, aes(x=time, y=Y, colour=id))+
		scale_colour_discrete(name="Experimental\n condition",
				breaks=c("1", "2"),
				labels=c("Control", "Treatment"))


#' Design evaluation using PFIM
#' --------------------------------
#' Update the (initial) values in the MDL Parameter Object for PFIM evaluation
#' Note that due to current limitations in PFIM, evaluation and optimisation of
#' multiple doses (as in the Simeoni estimation model) is not straightforward to encode.
#' For this demonstration, we will examine one single dose at
#' time 8 (days) so that the tumour has time to grow before treatment, so that we 
#' can see the effect of treatment.
#'  
#' Although the multiple dose design is not supported by PFIM, the model stays the same
#' and we can update the parameter values used for evaluating the study design using the
#' results of estimation using NONMEM and use simulation results to provide the prediction
#' of the total tumour weight (WTOT) at Day 8.
#'  

##' Assembling the new MOG.  
#PFIM.MOG <- createMogObj(designObj = myDesignObj, 
#		parObj = myParObjUpdated, 
#		mdlObj = myModelObj, 
#		taskObj = myTaskPropertiesObj_Evaluate)
##
##' We can then write the MOG back out to an .mdl file.
mdlfile.PFIM <- "InstallationTest_Design.mdl"
#writeMogObj(PFIM.MOG,mdlfile.PFIM)

#' Again, as for BUGS, the PFIM converter and connector have not yet fully been integrated
#' into the SEE environment. Here we call a small script which calls the converter which
#' creates the required PFIM input R scripts.
#'  
#' The Design Object settings and Task Properties settings dictate whether evaluation or 
#' optimisation of the trial design is performed.
source("callPFIMconverter_MKS_v2.R")

#' The runPFIM function takes an MDL file, converts to PharmML, runs the converter and then
#' runs a created batch script to run PFIM.
pharmmlfile <- as.PharmML(mdlfile.PFIM)

#+ PFIM, cache=TRUE
runPFIM(pharmmlfile=pharmmlfile.PFIM, jarLocation=file.path(ddmore:::DDMORE.checkConfiguration(),"PFIM","converter"))
#+ readPFIM
readLines(file.path(getwd(),"PFIM","stdout.out"))

#' Design evaluation using PopED
#' --------------------------------
#' PopED function as.poped takes a PharmML file and creates 
#' a poped.db database object ready for use with PopED.
#' We can then use PopED functions directly (natively) in R.
library(PopED)

as.poped(pharmMLFile)

#'  create plot of model without variability 
#+ Plot_model_predictions_from_PopED
plot_model_prediction(poped.db)

#' get predictions from model
model_prediction(poped.db)

#' evaluate initial design
FIM <- evaluate.fim(poped.db) 
FIM
det(FIM)
get_rse(FIM,poped.db)
