#' DDMoRe: Demonstration of interoperability
#' =========================================================================================
#' Authors: Letizia Carrara, Celine Sarr, Mike K Smith
#' Date: 18 July 2016
#'  
#' This example script is intended to demonstate interoperability
#' 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

#' 
#' Initialisation
#' =========================
#' Clear workspace and set working directory 

rm(list=ls(all=F))
#' 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")

#' 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()

#' The "~" directory is the default MDL-IDE workspace directory
projectDir <- file.path("~","NovartisDemo","models")
exampleDir <- file.path(projectDir, "simeoni_pop")
setwd(exampleDir)
getwd()

#' Nominate the MDL file for analysis. The MDL file below
#' describes the tumour growth of individual xenograft mice
#' following innoculation with tumour cells.
mdlfile <- "Simeoni_PAGE_Estimation_POP.mdl"
model <- tools::file_path_sans_ext(mdlfile)

#' The ddmore R package has been written to allow users to read, write and work
#' with MDL files and objects. getMDLObjects(<filename>) returns all MDL objects
#' from that file as an R list.
myMDLObj <- getMDLObjects(mdlfile) 
length(myMDLObj)
names(myMDLObj)
#str(myMDLObj) 

#' Next, 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. 
myDataObj <- getDataObjects(mdlfile)[[1]]
datafile <- myDataObj@SOURCE$scrfile$file

myParObj <- getParameterObjects(mdlfile)[[1]]
myModObj <- getModelObjects(mdlfile)[[1]]

#' You can also refer  to objects by their name. 
myTaskPropertiesObj_NM <- getTaskPropertiesObjects(mdlfile)[["simeoni2004_NONMEM_task"]]
myTaskPropertiesObj_MLX <- getTaskPropertiesObjects(mdlfile)[["simeoni2004_Monolix_task"]]

#' 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)
head(myData)

#' Extract only observation records
myEDAData<-myData[myData$MDV==0,]

library(ggplot2)
#' Plot the data using ggplot2 library 
plot1 <- ggplot() +
		geom_line(aes(y=DV,x=TIME,  color=as.factor(GRP), group=ID), data=myEDAData) +
		geom_point(aes(y=DV,x=TIME,  color=as.factor(GRP), group=ID), data=myEDAData) + 
		labs(y=" Tumour volume (mm3)" ,x="Time (days)")
print(plot1)


plot2 <- ggplot(aes(y=DV,x=TIME, group=GRP), data=myEDAData) +
		geom_line()+
		geom_point() + 
		facet_wrap(~ID) +
		labs(y=" Tumour volume (mm3)" ,x="Time (days)")
print(plot2)		

#' 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 MDL file to do specific tasks.
#'  
#' To define the new MOG for estimation in Monolix using the ddmore createMogObj function:
MLX.MOG <- createMogObj(dataObj = myDataObj, 
		parObj = myParObj, 
		mdlObj = myModObj, 
		taskObj = myTaskPropertiesObj_MLX)

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

#' The ddmore "estimate" function translates the contents of an .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) 

#' The ddmore "LoadSOObj" reads and parsed existing Standardise Output Objects:
# mlx <- LoadSOObject("Monolix/UseCase1.SO.xml")

#' 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.  
parameters_mlx <- getPopulationParameters(mlx, what="estimates")$MLE
parameters_mlx

#' Currently, we need a workaround to rename the "b" parameter 
#' from Monolix (b=proportional error parameter) 
#' to "CV" to match NONMEM. 
parameters_mlx <- reshape::rename(parameters_mlx, c(b="CV"))

#' We can also return the MLEs with associated precision / variances
getPopulationParameters(mlx, what="precisions") 

#' We can also return information about the model fit including LogLikelihood, AIC, BIC
getEstimationInfo(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.
mlx.xpdb<-as.xpdb(mlx,datafile)

#+ Basic_Goodness_of_Fit_following_Monolix_Estimation
basic.gof(mlx.xpdb)
#+ Individual_Plots_following_Monolix_Estimation
ind.plots(mlx.xpdb) 

#'We can export graphs to PDF file
#+ Xpose_Monolix, eval=FALSE
pdf("GOF_MLX.pdf")
print(basic.gof(mlx.xpdb))
print(ind.plots(mlx.xpdb))
dev.off()

#' 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. Note that the NONMEM estimation
#' method specified is FOCE. We will look at how to change this via R later.
NM.MOG <- createMogObj(dataObj = myDataObj, 
		parObj = myParObj, 
		mdlObj = myModObj, 
		taskObj = myTaskPropertiesObj_NM)

#' 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")

#' Load previous results
#NM <- LoadSOObject("NONMEM/Simeoni_PAGE_Estimation_POP_NM.SO.xml")
 
#' Results from NONMEM should be comparable with results from MONOLIX
parameters_nm <- getPopulationParameters(NM, what="estimates")$MLE

#' We can compare the results across tools. We sort on the names of parameters from 
#' NONMEM output
estPars <- names(parameters_nm)
cbind(  NM=parameters_nm[estPars],
		MLX=parameters_mlx[estPars]) 

#' Covariance step cannot be requested in the current version of the Framework
#getPopulationParameters(NM, what="precisions") 

#' As before we can request additional information about the NONMEM estimation process.
print(getEstimationInfo(NM))

#' Again we convert the SO object to an Xpose format database and perform model diagnostics.
nm.xpdb<-as.xpdb(NM,datafile)

#' Perform some basic goofness of fit
basic.gof(nm.xpdb)
ind.plots(nm.xpdb) 

#'Export graphs to PDF file
#+ Xpose_NONMEM, eval=FALSE
pdf("GOF_NM.pdf")
print(basic.gof(mlx.xpdb))
print(ind.plots(mlx.xpdb))
dev.off()


#' Change estimation algorithm to SAEM
#' --------------------------------
#' Let's change the estimation algorithm from FOCE to SAEM.
myTaskPropertiesObj_NM_SAEM <- myTaskPropertiesObj_NM
myTaskPropertiesObj_NM_SAEM@ESTIMATE$algo <- "saem"

#' create new MOG, write to file and estimate.
NM.MOG <- createMogObj(dataObj = myDataObj, 
		parObj = myParObj, 
		mdlObj = myModObj, 
		taskObj = myTaskPropertiesObj_NM_SAEM)

mdlfile.NM.SAEM <- paste0(model,"_NM_SAEM.mdl")
writeMogObj(NM.MOG,mdlfile.NM.SAEM)

#+ NONMEM_SAEM, cache=TRUE
NM.SAEM <- estimate(mdlfile.NM.SAEM, target="NONMEM", subfolder="Nonmem_SAEM")

parameters_nm_saem <- getPopulationParameters(NM.SAEM, what="estimates")$MLE

#' We can compare the results across tools. We sort on the names of parameters from 
#' NONMEM output
estPars <- names(parameters_nm)
cbind(  NM=parameters_nm[estPars],
		NM.SAEM=parameters_nm_saem[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.   
myParObj <- getParameterObjects(mdlfile)[[1]]
myParObjUpdated <- updateParObj(myParObj,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_ALL <- createMogObj(dataObj = myDataObj, 
		parObj = myParObjUpdated, 
		mdlObj = myModObj, 
		taskObj = myTaskPropertiesObj_NM)

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

#' 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_ALL,samples=20, seed=12345,
		vpcOptions ="-n_simulation=10  -auto_bin=3  -stratify_on=GRP" ,
		subfolder="VPC", plot=FALSE) 

#' To replay the visualisation using information from the VPC SO file
print(xpose.VPC(vpc.info= file.path("./VPC",vpcFiles@RawResults@DataFiles$PsN_VPC_results$path),
				vpctab= file.path("./VPC",vpcFiles@RawResults@DataFiles$PsN_VPC_vpctab$path),
				by="GRP",
				main="VPC SIMEONI 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.

#' Use parameter values from the FOCE estimation
p <- parameters_nm

#' adm : dose administration times
adm <- list( time  = c(8, 12, 16), 
		amount = 30, 
		target = 'Q1')

#' f : what to sample ("observe") and at what times?
f <- list( name='WTOT', time=seq(0,45,by=0.1))
y <- list( name='Y', time=c(0.001, 8, 12, 16, 20, 24, 28, 32, 35))

#' g : define group sizes and assign treatments.
g0 <- list(size=8, level='longitudinal')
g1 <- list(size=8, 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 = list(g0,g1),
		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"))

#' ESTIMATE model parameters in WINBUGS
#' -------------------------
#' BUGS dataset conversion currently has problems with categorical covariates.
#' So we read the Data Object in MDL and change the use of GRP variable to "ignore".

BUGSDataObj <- myDataObj
BUGSDataObj@DATA_INPUT_VARIABLES$GRP$use <- "ignore"

BUGSModObj <- myModObj 
slot(BUGSModObj,"COVARIATES") <- list()

#' 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 = BUGSDataObj, 
		priorObj = getPriorObjects(mdlfile)[[1]], 
		mdlObj = BUGSModObj, 
		taskObj = getTaskPropertiesObjects(mdlfile)[["simeoni2004_BUGS_task"]])

#' 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"
file.exists(mdlfile.BUGS)
#' 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.
source(file.path(ddmore:::DDMORE.checkConfiguration(),"BlackBoxWinBUGS","ddmorewinbugs.R"))
#+ BUGS, cache=TRUE
BUGS<-run.winbugs(mdlfile.BUGS,subfolder="WinBUGS")

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


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

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

#' Summary statistics of MCMC samples
summary(coda_out)
## trace and density plots
#pdf("MCMC_diagnostics_BUGS.pdf")
#plot(coda_out, ask=F)
#dev.off()
## gelman-rubin diagnostics (if at least 2 chains)
#gelman.diag(coda_out)
#gelman.plot(coda_out)

#' The SO conversion for BUGS output is not yet complete, but eventually we should be able to 
#' perform Xpose diagnostics after estimation with BUGS as we have with Monolix and NONMEM.
# bugs<-as.xpdb(BUGS,datafile)

#' Compare estimated parameters in the 3 tools
#' ----------------------------------
#' Note that in this example, we may need to run the MCMC sampler in BUGS to get better
#' estimates of model parameters (see the MCMC diagnostics above), but for expediency
#' we have limited the sampler to 5000 samples from the posterior distribution. 
estPars <- intersect(rownames(parameters_BUGS), names(parameters_nm))
cbind(BUGS=parameters_BUGS[estPars,"Median"],
		NM=parameters_nm[estPars],
		MLX=parameters_mlx[estPars]) 

