###########################################
# Example: Theophylline project
# Script to initialize population parameters in Monolix and run scenarios
# 
# Monolix files for the fit without IIV and for the fit with IIV are required
# with all relevant fitting settings set:
#   - normal distribution for parameters (instead of log-normal)
#   -'add decreasing variability' for fit without IIV
# Here, THEOPHYLLINE_project_pop_train.mlxtran is the file without IIV
# and THEOPHYLLINE_project_ind_train.mlxtran is the file with IIV
#
# The structural model reads
# dA/dt = - NNa(A) - D*NNat(t)
# dC/dt = NNa(A) + NNat(t) + NNc(C) + D*NNct(t)
###########################################

#remove variables from the global environment
rm(list=ls())

#load packages
library(tidyverse) #for easier data manipulation
library(lixoftConnectors) #for connection to Monolix

#initialize lixoftconnecotr with Monolix
initializeLixoftConnectors(software = "monolix", path = "C:/ProgramData/Lixoft/MonolixSuite2023R1")

#set your working directory
setwd("~/your_working_directory")

#define ranges for activation points
acmin <- 0 #minimal expected concentration in the absorption compartment
acmax <- 4 #maximal expected concentration in the absorption compartment (~maximal dose)
atmin <- 0.5 #minimal expected time for absorption process
atmax <- 3.5 #maximal expected time for absorption process

ccmin <- 0.5 #minimal expected concentration in the central compartment
ccmax <- 10 #maximal expected concentration in the central compartment (~maximal observed concentration)
ctmin <- 0.5 #minimal expected time for distribution process
ctmax <- 20 #maximal expected time for distribution process

#seed for random number generator
set.seed(1234)

#sampling weights from a discrete pseudo-normal distribution (higher probability for values closer to 0)
#weights for the time-dependent absorption NN
wat1_s  <- sample(c(-0.3,-0.2,-0.2,-0.1,-0.1,-0.1,0.1,0.1,0.1,0.2,0.2,0.3),5,replace = F)
wat2_s <- sample(c(-0.3,-0.2,-0.2,-0.1,-0.1,-0.1,0.1,0.1,0.1,0.2,0.2,0.3),5,replace = F)

#weights and bias for the amount-dependent absorption NN
wa1_s  <- sample(c(-0.3,-0.2,-0.2,-0.1,-0.1,-0.1,0.1,0.1,0.1,0.2,0.2,0.3),5,replace = F)
wa2_s <- sample(c(-0.3,-0.2,-0.2,-0.1,-0.1,-0.1,0.1,0.1,0.1,0.2,0.2,0.3),5,replace = F)
ba2_s <- sample(c(-0.3,-0.2,-0.2,-0.1,-0.1,-0.1,0.1,0.1,0.1,0.2,0.2,0.3),1,replace = F)

#weights for the time-dependent central NN
wct1_s  <- sample(c(-0.3,-0.2,-0.2,-0.1,-0.1,-0.1,0.1,0.1,0.1,0.2,0.2,0.3),5,replace = F)
wct2_s <- sample(c(-0.3,-0.2,-0.2,-0.1,-0.1,-0.1,0.1,0.1,0.1,0.2,0.2,0.3),5,replace = F)

#weights and bias for the concentration-dependent central NN
wc1_s  <- sample(c(-0.3,-0.2,-0.2,-0.1,-0.1,-0.1,0.1,0.1,0.1,0.2,0.2,0.3),5,replace = F)
wc2_s <- sample(c(-0.3,-0.2,-0.2,-0.1,-0.1,-0.1,0.1,0.1,0.1,0.2,0.2,0.3),5,replace = F)
bc2_s <- sample(c(-0.3,-0.2,-0.2,-0.1,-0.1,-0.1,0.1,0.1,0.1,0.2,0.2,0.3),1,replace = F)

#sampling random activation points for the absorption NNs
atacs <- runif(5,atmin,atmax) #time activation points
acacs <- runif(5,acmin,acmax) #amount activation points

#sampling random activation points for the central NNs
ctacs <- runif(5,ccmin,ccmax) #time activation points
ccacs <- runif(5,ctmin,ctmax) #concentration activation points

#calculating biases from weights (w) and activation points (ap) according to b = -w * ap
bat1_s <- -(wat1_s^2) * atacs #time dependent absorption biases (w' = w^2)
ba1_s <- -wa1_s * acacs #amount dependent absorption biases

bct1_s <- -(wct1_s^2) * ctacs #time dependent central biases (w' = w^2)
bc1_s <- -wc1_s * ccacs #concentration dependent central biases

#vector for initial parameters (in same order as they are in the mlxtran file)
ini_parms <- c(wat1_s,bat1_s,wat2_s,wa1_s,ba1_s,wa2_s,ba2_s,wc1_s,bc1_s,wc2_s,bc2_s,wct1_s,bct1_s,wct2_s)

#load the mlx file without IIV
loadProject("THEOPHYLLINE_project_pop_train.mlxtran")

#get the original initial values data frame
og_ini <- getPopulationParameterInformation()

#create a new intial values data frame with previously sampled/calculated values
new_ini <- og_ini %>%
  mutate(initialValue = ifelse(grepl("pop",name),ini_parms,initialValue))

#set the new initial values
setPopulationParameterInformation(new_ini)

#save the project and run the population estimation without IIV
saveProject("THEOPHYLLINE_project_pop_train_newini1.mlxtran")
runScenario()

#get the population estimates
est_pop_parms <- getEstimatedPopulationParameters()

#load the mlx file with IIV
loadProject("THEOPHYLLINE_project_ind_train.mlxtran")

#get the original initial values data frame
og_ind_ini <- getPopulationParameterInformation()

#create a new initial values data frame with previously estimated population values and omegas of 0.1
new_ind_ini <- og_ind_ini %>%
  mutate(initialValue = ifelse(grepl("pop",name),est_pop_parms,initialValue)) %>%
  mutate(initialValue = ifelse(grepl("omega",name),0.1,initialValue))

#set the new intial values
setPopulationParameterInformation(new_ind_ini)

#save the project an run the population and individual estimation with IIV
saveProject("THEOPHYLLINE_project_ind_train_newini1.mlxtran")
runScenario()