## ----setup, include=FALSE---------------------------------------------------------------------------------------------
backup_options = options()
options(width = 120)

knitr::opts_chunk$set(
echo       = TRUE,
eval       = FALSE,          # chunks do NOT execute by default (CRAN compliance)
warning    = FALSE,
message    = FALSE,
comment    = "#>",
cache      = FALSE,
fig.width  = 8,
fig.height = 5,
fig.align  = "center",
out.width  = "90%"
)

library(PFIM)

# paths$data    : pre-computed RDS files shipped in inst/extdata/
# paths$outputs : text outputs  -> tempdir() (CRAN-compliant)
# paths$reports : HTML reports  -> tempdir()
paths = list(
data    = system.file("extdata", package = "PFIM"),
outputs = file.path(tempdir(), "pfim_out"),
reports = file.path(tempdir(), "pfim_rep")
)
dir.create(paths$outputs, showWarnings = FALSE, recursive = TRUE)
dir.create(paths$reports, showWarnings = FALSE, recursive = TRUE)

# Shared plot options
plotOptions = list(unitTime = "hour", unitOutcomes = c("mcg/mL", "DI%"))

## ----model-library, echo = TRUE, eval = FALSE, comment=''-------------------------------------------------------------
# modelFromLibrary = list("PKModel" = "Linear1InfusionSingleDose_ClV")

## ----model-parameters-------------------------------------------------------------------------------------------------
# modelParameters = list(
#   ModelParameter(name = "V",  distribution = LogNormal(mu = 50, omega = sqrt(0.26))),
#   ModelParameter(name = "Cl", distribution = LogNormal(mu = 5,  omega = sqrt(0.34)))
# )

## ----error-model, echo = TRUE, eval = FALSE, comment=''---------------------------------------------------------------
# errorModelRespPK = Combined1( output = "RespPK", sigmaInter = 0.5, sigmaSlope = sqrt( 0.15 ) )
# modelError = list( errorModelRespPK )

## ----administration, echo = TRUE, eval = FALSE, comment=''------------------------------------------------------------
# administrationRespPK = Administration(
#   outcome  = "RespPK",
#   Tinf     = rep(1, 5),           # 1-hour infusion for each of the 5 doses
#   timeDose = seq(0, 96, 24),      # dosing at t = 0, 24, 48, 72, 96 h
#   dose     = c(400, rep(200, 4))  # 400 mg loading + 4 × 200 mg maintenance
# )

## ----sampling-times-eval, echo = TRUE, eval = FALSE, comment=''-------------------------------------------------------
# samplingTimesRespPK = SamplingTimes(
#   outcome   = "RespPK",
#   samplings = c(1, 12, 24, 44, 72, 120)
# )

## ----arm-design-eval, echo = TRUE, eval = FALSE, comment=''-----------------------------------------------------------
# arm1 = Arm(
#   name            = "arm1",
#   size            = 150,
#   administrations = list(administrationRespPK),
#   samplingTimes   = list(samplingTimesRespPK)
# )
# 
# design1 = Design(name = "design1", arms = list(arm1))

## ----evaluation-run, echo = TRUE, eval = FALSE, comment=''------------------------------------------------------------
# # --- Population FIM ---
# evaluationPop = Evaluation(
#   name                = "evaluationPop",
#   modelFromLibrary    = modelFromLibrary,
#   modelParameters     = modelParameters,
#   modelError          = modelError,
#   outputs             = list("RespPK"),
#   designs             = list(design1),
#   fimType             = "population",
#   odeSolverParameters = list(atol = 1e-8, rtol = 1e-8)
# )
# evaluationPopResults = run(evaluationPop)
# 
# # --- Individual FIM ---
# evaluationInd = Evaluation(
#   name                = "evaluationInd",
#   modelFromLibrary    = modelFromLibrary,
#   modelParameters     = modelParameters,
#   modelError          = modelError,
#   outputs             = list("RespPK"),
#   designs             = list(design1),
#   fimType             = "individual",
#   odeSolverParameters = list(atol = 1e-8, rtol = 1e-8)
# )
# evaluationIndResults = run(evaluationInd)
# 
# # --- Bayesian FIM ---
# evaluationBay = Evaluation(
#   name                = "evaluationBay",
#   modelFromLibrary    = modelFromLibrary,
#   modelParameters     = modelParameters,
#   modelError          = modelError,
#   outputs             = list("RespPK"),
#   designs             = list(design1),
#   fimType             = "Bayesian",
#   odeSolverParameters = list(atol = 1e-8, rtol = 1e-8)
# )
# evaluationBayResults = run(evaluationBay)

## ----results-show-pop, echo = TRUE, eval = FALSE, comment=''----------------------------------------------------------
# show(evaluationPopResults)

## ----results-pop_from_outputs, echo = FALSE,  eval=TRUE, comment=''---------------------------------------------------
lines = readLines("outputs/vignette2_evaluation_PopFIM_show.txt")
cat(paste(lines, collapse = "\n"))

## ----results-show-ind, echo = TRUE, eval = FALSE, comment=''----------------------------------------------------------
# show(evaluationIndResults)

## ----results-ind_from_outputs, echo = FALSE,  eval=TRUE, comment=''---------------------------------------------------
lines = readLines("outputs/vignette2_evaluation_IndFIM_show.txt")
cat(paste(lines, collapse = "\n"))

## ----results-show-bay, echo = TRUE, eval = FALSE, comment=''----------------------------------------------------------
# show(evaluationBayResults)

## ----results-bay_from_outputs, echo = FALSE,  eval=TRUE, comment=''---------------------------------------------------
lines = readLines("outputs/vignette2_evaluation_BayFIM_show.txt")
cat(paste(lines, collapse = "\n"))

## ----eval-accessors, echo = TRUE, eval = FALSE, comment=''------------------------------------------------------------
# cat("--- Fisher Information Matrix (population FIM) ---\n")
# print(getFisherMatrix(evaluationPopResults))
# 
# cat("--- Correlation Matrix (population FIM) ---\n")
# print(getCorrelationMatrix(evaluationPopResults))
# 
# cat("--- Standard Errors ---\n")
# print(getSE(evaluationPopResults))
# 
# cat("--- Relative Standard Errors ---\n")
# print(getRSE(evaluationPopResults))
# 
# cat("--- Shrinkage ---\n")
# print(getShrinkage(evaluationPopResults))
# 
# cat("--- D-Criterion ---\n")
# print(getDcriterion(evaluationPopResults))

## ----eval-plots, echo = TRUE, eval = FALSE, comment=''----------------------------------------------------------------
# plotOptions = list(unitTime = "hour", unitOutcomes = "mcg/mL")
# # plot() is the unified OO entry point -- dispatches on class, returns a named list:
# #   $evaluation         -> nested [["design"]][["arm"]][["outcome"]]
# #   $sensitivityIndices -> nested [["design"]][["arm"]][["outcome"]][["param"]]
# #   $SE / $RSE          -> ggplot2 bar charts
# # Predicted concentration profile with sampling times overlaid
# plotsEval = plot(evaluationPopResults,
#                  plotOptions = plotOptions,
#                  which       = c("evaluation", "sensitivityIndices", "SE", "RSE"))
# plotOutcomesRespPK = plotsEval$evaluation[["design1"]][["arm1"]][["RespPK"]]
# print(plotOutcomesRespPK)

## ----plot-eval-pk_from_figures, echo = FALSE,  eval=TRUE, out.width="50%",comment='', fig.align = "center", fig.cap="PK response profile -- arm1 (population FIM)"----
knitr::include_graphics("figures/vignette2_evaluation_popFim_design1_arm1_RespPK.png")

## ----eval-plots-si-V-Cl, echo = TRUE, eval = FALSE, comment=''--------------------------------------------------------
# # $sensitivityIndices is computed in the same plot() call above
# print(plotsEval$sensitivityIndices[["design1"]][["arm1"]][["RespPK"]][["V"]])
# print(plotsEval$sensitivityIndices[["design1"]][["arm1"]][["RespPK"]][["Cl"]])

## ----plot-si-V_from_figures, echo = FALSE,  eval=TRUE, out.width="50%",comment='', fig.align = "center", fig.cap="Sensitivity index for V -- RespPK, arm1"----
knitr::include_graphics("figures/vignette2_evaluation_popFim_design1_SI_RespPK_V.png")

## ----plot-si-cl_from_figures, echo = FALSE,  eval=TRUE, out.width="50%",comment='', fig.align = "center", fig.cap="Sensitivity index for Cl -- RespPK, arm1"----
knitr::include_graphics("figures/vignette2_evaluation_popFim_design1_SI_RespPK_Cl.png")

## ----plot-se, echo = TRUE, eval= FALSE,   comment=''------------------------------------------------------------------
# # Standard error and RSE bar charts -- computed in the same plot() call above
# print(plotsEval$SE)
# print(plotsEval$RSE)

## ----plot-se_figures, out.width="33%", echo = FALSE, eval= TRUE,  comment='', fig.align = "center",fig.cap="Standard Errors (SE)"----
knitr::include_graphics("figures/vignette2_evaluation_popFim_design1_SE.png")

## ----plot-rse_figures, out.width="33%", echo = FALSE, eval= TRUE,  comment='', fig.align = "center",fig.cap="Relative Standard Errors (RSE)"----
knitr::include_graphics("figures/vignette2_evaluation_popFim_design1_RSE.png")

## ----optim-shared-setup, echo = TRUE, eval = FALSE, comment=''--------------------------------------------------------
# samplingTimesOptim = SamplingTimes(
#   outcome   = "RespPK",
#   samplings = c(1, 48, 72, 120)   # initial guess: window boundaries
# )
# 
# samplingConstraintsRespPK = SamplingTimeConstraints(
#   outcome                = "RespPK",
#   initialSamplings       = c(1, 48, 72, 120),
#   numberOfTimesByWindows = c(2, 2),
#   samplingsWindows       = list(c(1,  48),    # Window 1: post-dose 1 to pre-dose 3
#                                 c(72, 120)),  # Window 2: near-SS to post-last-dose
#   minSampling            = 5                  # min 5 h between consecutive samples
# )
# 
# arm2 = Arm(
#   name                     = "arm2",
#   size                     = 150,
#   administrations          = list(administrationRespPK),
#   samplingTimes            = list(samplingTimesOptim),
#   samplingTimesConstraints = list(samplingConstraintsRespPK)
# )
# 
# # numberOfArms = 150: upper bound on distinct elementary protocols
# design2 = Design(name = "design2", arms = list(arm2), numberOfArms = 150)

## ----pso-run, echo = TRUE, eval = FALSE, comment=''-------------------------------------------------------------------
# optimizationPSO = Optimization(
#   name                = "PSO",
#   modelFromLibrary    = modelFromLibrary,
#   modelParameters     = modelParameters,
#   modelError          = modelError,
#   optimizer           = "PSOAlgorithm",
#   optimizerParameters = list(
#     maxIteration                = 100,
#     populationSize              = 50,
#     personalLearningCoefficient = 2.05,
#     globalLearningCoefficient   = 2.05,
#     seed                        = 42,
#     showProcess                 = FALSE
#   ),
#   designs  = list(design2),
#   fimType  = "population",
#   outputs  = list("RespPK")
# )
# optimizationPSO = run(optimizationPSO)

## ----pso-load, echo = TRUE, eval = FALSE, comment=''------------------------------------------------------------------
# show(optimizationPSO)

## ----results-pso-popFIM_from_outputs, echo = FALSE,  eval=TRUE, comment=''--------------------------------------------
lines = readLines("outputs/vignette2_optimization_PSO_populationFIM_show.txt")
cat(paste(lines, collapse = "\n"))

## ----pgbo-run, echo = TRUE, eval = FALSE, comment=''------------------------------------------------------------------
# optimizationPGBO = Optimization(
#   name                = "PGBO",
#   modelFromLibrary    = modelFromLibrary,
#   modelParameters     = modelParameters,
#   modelError          = modelError,
#   optimizer           = "PGBOAlgorithm",
#   optimizerParameters = list(
#     N              = 30,
#     muteEffect     = 0.65,
#     maxIteration   = 1000,
#     purgeIteration = 200,
#     seed           = 42,
#     showProcess    = FALSE
#   ),
#   designs  = list(design2),
#   fimType  = "population",
#   outputs  = list("RespPK")
# )
# optimizationPGBO = run(optimizationPGBO)

## ----pgbo-load, echo = TRUE, eval = FALSE, comment=''-----------------------------------------------------------------
# show(optimizationPGBO)

## ----results-pgbo-popFIM_from_outputs, echo = FALSE,  eval=TRUE, comment=''-------------------------------------------
lines = readLines("outputs/vignette2_optimization_PGBO_populationFIM_show.txt")
cat(paste(lines, collapse = "\n"))

## ----simplex-run, echo = TRUE, eval = FALSE, comment=''---------------------------------------------------------------
# optimizationSimplex = Optimization(
#   name                = "Simplex",
#   modelFromLibrary    = modelFromLibrary,
#   modelParameters     = modelParameters,
#   modelError          = modelError,
#   optimizer           = "SimplexAlgorithm",
#   optimizerParameters = list(
#     pctInitialSimplexBuilding = 10,
#     maxIteration              = 1000,
#     tolerance                 = 1e-10,
#     showProcess               = FALSE
#   ),
#   designs  = list(design2),
#   fimType  = "population",
#   outputs  = list("RespPK")
# )
# optimizationSimplex = run(optimizationSimplex)

## ----simplex-load, echo = TRUE, eval = FALSE, comment=''--------------------------------------------------------------
# show(optimizationSimplex)

## ----results-simplex-popFIM_from_outputs, echo = FALSE,  eval=TRUE, comment=''----------------------------------------
lines = readLines("outputs/vignette2_optimization_Simplex_populationFIM_show.txt")
cat(paste(lines, collapse = "\n"))

## ----teardown, echo=FALSE, include=FALSE------------------------------------------------------------------------------
# options(backup_options)