diff --git a/.gitignore b/.gitignore index 6cd3c721e..d6cce0004 100755 --- a/.gitignore +++ b/.gitignore @@ -64,7 +64,6 @@ tests/testthat/Runs/ .vscode/ .Rd2pdf65996/ other/ -Examples Experimental .httr-oauth inst/rust/template/* @@ -73,8 +72,6 @@ inst/options/PMoptions.json inst/doc target/ -*_test.R -*.csv *.txt 1 2 diff --git a/Cargo.lock b/Cargo.lock index 0a6d1402e..a4064bec1 100755 --- a/Cargo.lock +++ b/Cargo.lock @@ -1286,9 +1286,9 @@ dependencies = [ [[package]] name = "pharmsol" -version = "0.20.0" +version = "0.21.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "98b8e2ab3a0e91cd4b20c28544cb3676e8df31aa490cf5680ec0531259b5fa4e" +checksum = "2fc25564d039d0cd5701013aa3785a339b14cf0b51409d7b817320bc360dc944" dependencies = [ "argmin", "argmin-math", @@ -1320,6 +1320,7 @@ version = "0.1.0" dependencies = [ "anyhow", "extendr-api", + "libloading", "pmcore", "rayon", "tracing", @@ -1328,13 +1329,12 @@ dependencies = [ [[package]] name = "pmcore" -version = "0.21.0" +version = "0.22.1" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "74646008c207f79cdfd54152cefae76ebcba895fd6005b23de42dfb838a2ca95" +checksum = "3866100507aa3bcba475381af3102d84b5e503bce82cb82f56cf0fa46cc1e408" dependencies = [ "anyhow", "argmin", - "argmin-math", "csv", "faer", "faer-ext", diff --git a/Examples/Rscript/examples.R b/Examples/Rscript/examples.R new file mode 100644 index 000000000..15ccf6921 --- /dev/null +++ b/Examples/Rscript/examples.R @@ -0,0 +1,778 @@ +# INTRODUCTION ------------------------------------------------------------ + +# Lines that start with "#" are comments and ignored by R. Follow the +# directions in them. Execute each non-comment line in this script by +# putting your cursor on it and sending it to the R console. +# You can do this in several ways: +# Windows +# R-studio +# 1) The Run button at the top +# 2) Ctrl-Enter +# R GUI - when the script window is active +# 1) The Run line or selection button at the top +# 2) Ctrl-R +# Mac +# R-studio +# 1) The Run button at the top +# 2) Command-Enter +# R GUI +# 1) Command-Enter + +# This script also serves to introduce several R programming functions +# and techniques. For any function, you can get help by typing ?function_name +# in the R console (the lower left window pane in RStudio). + +# Load Pmetrics into memory. You must include this line at the beginning +# of every script. + +library(Pmetrics) + +# EXERCISE 1 - NPAG RUN ------------------------------------------------ + +# EXAMPLE NPAG RUN - tlag, ka, kel, vol + +# It is useful to annotate your runs as above, so that you can remember +# what you did later! + + +# Tell R where your working directory is going to be. +# Windows users: Make sure that you separate directories with a +# forward slash "/" or double backslashes "\\". Unfortunately, Windows is the only OS that uses +# backslashes "\", so R conforms to Unix/Linux style. + +wd <- "C:/Users/siel3/code/LAPKB/Pmetrics_rust/Examples" + +wd <- glue::glue("{getwd()}/inst/Examples/Runs") + +# change to the working directory to the Examples folder +setwd(wd) + +# DATA OBJECT + +# Pmetrics always needs data and a model to run +# create our first data object + +# create a new data object by reading a file +# set the limit of quantification (loq) to 1: see ?PM_data for help +exData <- PM_data$new(data = "../src/ex.csv", loq = 1) + +# you can look at this file directly by opening it in +# a spreadsheet program like Excel, or a text editor + +# exData is an R6 object, which means that contains both data and methods to +# process that data, for example: +exData$data # contains your original datafile +exData$standard_data # contains the standardized and validated data, +exData$summary() # prints the summary of the data to the terminal, or + +# another way to do that is using the more common S3 framework in R: +summary(exData) + +# To look at the contents of an object: +names(exData) + +# other examples of things that can be done with this object are +exData # view the original data in the viewer +exData$print(standard = TRUE) # view the standardized data in the viewer +exData$print(viewer = FALSE) # view original data in console +exData$plot() #plot the raw data; more on that later + +# MODEL OBJECT +# You can specify a model by reading a file or directly as an object. We'll do both. +# The following code creates the same model as in /src/model.txt file. +# See PMmanual() for details on creating models in R compared to text files. +# The advantage of creating them in R is that one does not need to copy model +# files into folders to provide necessary inputs. + +mod1 <- PM_model$new( + pri = list( + Ka = ab(0.1, 0.9), + Ke = ab(0.001, 0.1), + V = ab(30, 120), + lag1 = ab(0, 4) + ), + cov = list( + wt = interp(), + africa = interp("none"), + age = interp(), + gender = interp("none"), + height = interp() + ), + eqn = function(){ + two_comp_bolus + }, + lag = function(){ + lag[1] = lag1 + }, + out = function(){ + Y[1] = X[2]/V + }, + err = list( + proportional(5, c(0.02, 0.05, -0.0002, 0)) + ) +) + + +# look at it +mod1 + +#plot it +mod1$plot() + + +# in the working directory we have another file "model.txt" that contains the old +# representation of the same model we previously presented, let's take a look at it. +system("cat ../src/model.txt") + +# PM_model$new() also accepts the path to a model file +# create the same model using this file +mod1b <- PM_model$new("../src/model.txt") +mod1b + +# PM_model provides a method to update the different elements of a model, for example: +mod1b$update( + pri = list( + ka = ab(0.001, 5) +)) + +# It is case sensitive, so ka is different from Ka. To remove a parameter, set it to NULL. + + +mod1b$arg_list$pri + +# to copy a model use the $clone() method. +mod1b <- mod1$clone() + +# simply using mod1b <- mod1 will cause mod1b to be changed if mod1 is changed, +# as R6 objects use reference semantics. For more details you can refer to +# https://adv-r.hadley.nz/r6.html, Section 14.4. + +#lastly, use the app! PMmanual() and the article on models for details on this. +build_model() #start from scratch +build_model(exData) #start with data to match covariates +build_model(mod1) #start with a model and update it + + + +# To keep everything tidy, we are working in a folder specific to store the runs + + +run1 <- mod1$fit(data = exData, run = 1, overwrite = TRUE) # execute the fit and return the results to run1 + + +# +# After the run is complete the results are returned to the assigned object, here 'run1' + +# you need get the extracted information back into R. +# They will be sequentially numbered as /1, /2, /3,... in your working directory. + +# One benefit of having this fit object is that it is possible to run multiple +# fittings without needing to move data files around + +getwd() +list.files() + +# Result Object - the result is already returned to run1 above, but if you need to load it later, +# you can use PM_load() +run1 <- PM_load(1) + +# Create a PM_result object by reading a run folder. The "1" in the parentheses tells Pmetrics to +# look in the /1 folder. + +# Plot the raw data using R6 with various options. Type ?plot.PM_data in the R console for help. +run1$data$plot() +run1$data$plot(overlay = FALSE, xlim = c(119, 145)) +run1$data$plot(xlim = c(119, 146), group = "gender", group_names = c("Male", "Female"), + marker = list(color = c("blue","red"), symbol = c("circle","triangle-up"))) + +run1$data$plot(xlim = c(119, 146), group = "gender", group_names = c("Male", "Female"), marker = list(color = "Set2")) + + +# The following are the older S3 method with plot(...) for the first two examples +# You can use R6 or S3 for any Pmetrics object +# We will focus on R6 as the more modern way. +plot(run1$data) + +# here's a summary of the original data file; ?summary.PMmatrix for help +run1$data$summary() + +# Plot some observed vs. predicted data. Type ?plot.PMop in the R console for help. +# Click on points to highlight all points from the same subject. +run1$op$plot() +run1$op$plot(pred.type = "pop") +run1$op$plot(line = list(lm = list(ci = 0, color = "red"), loess = FALSE)) + +# the original op object data can be accessed via +run1$op$data + +# get a summary with bias and imprecision of the population predictions; +# ?summary.PMop for help +run1$op$summary(pred.type = "pop") + +# the S3 way +summary(run1$op, pred.type = "pop") + +# look at the summary for the posterior predictions (default pred.type) based +# on means of parameter values +run1$op$summary(icen = "mean") + +# The OP plot can be disaggregated into a Tidy compatible format from $data (see https://www.tidyverse.org/) +# This allow pre processing in ways more flexible than the default plot method. +library(tidyverse) +run1$op$data %>% plot() +run1$op$data %>% + filter(pred > 5) %>% + filter(pred < 10) %>% + summary() + +# see a header with the first 10 rows of the op object +head(run1$op$data, 10) + + +# Plot final population joint density information. Type ?plot.PMfinal in the R console for help. +run1$final$plot() + +# add a kernel density curve +run1$final$plot(line = list(density = list(color = "red"))) +run1$final$data %>% plot() + + +# A bivariate plot. Plotting formulae in R are of the form 'y~x' +run1$final$plot(ke ~ v, + marker = list(color = "red", symbol = "diamond")) + + + +# or the S3 way +plot(run1$final) + +# The original final object can be accessed using +run1$final$data +names(run1$final$data) + +# see the population points +run1$final$popPoints + +# or +run1$final$data$popPoints + +# see the population mean parameter values +run1$final$popMean + +# see a summary with confidence intervals around the medians +# and the Median Absolute Weighted Difference (MAWD); +# ?summary.PMfinal for help +run1$final$summary() + +# summarize the cycle information; ?summary.PM_cycle for help +run1$cycle$summary() +run1$cycle$data %>% summary() + + +# Plot cycle information +# Type ?plot.PM_cycle in the R console for help. +run1$cycle$plot() + +# names of the cycle object; ?makeCycle for help +names(run1$cycle$data) + +# gamma/lamda value on last 6 cycles +tail(run1$cycle$data$gamlam) + +# Plot covariate information. Type ?plot.PMcov in the R console for help. +# Recall that plotting formulae in R are of the form 'y~x' +run1$cov$plot(v ~ wt) + +run1$cov$data %>% + filter(age > 25) %>% + plot(v ~ wt) + +# Plot +run1$cov$plot(ke ~ age, line = list(loess = FALSE, lm = TRUE), + marker = list(symbol = 3)) + +# Another plot with mean Bayesian posterior parameter and covariate values... +# Remember the 'icen' argument? +run1$cov$plot(v ~ wt, icen = "mean") + +# When time is the x variable, the y variable is aggregated by subject. +# In R plot formulae, calculations on the fly can be included using the I() function +run1$cov$plot(I(v * wt) ~ time) + +# The previous cov object can be seen via: +run1$cov + +# but to access individual elements, use: +run1$cov$data[, 1:3] # for example +names(run1$cov) + +# summarize with mean covariates; ?summary.PMcov for help +run1$cov$summary(icen = "mean") + + +# Look at all possible covariate-parameter relationships by multiple linear regression with forward +# and backward elimination - type ?PM_step in the R console for help. +run1$step() +# or on the cov object directly +run1$cov$step() +# icen works here too.... +run1$step(icen = "mean") +# forward elimination only +run1$step(direction = "forward") + + +# EXERCISE 2 - NPAG WITH COVARIATES --------------------------------------- + +# Again, without copying files, let's create another run object, this time using +# a model that include covariates + +# First clone mod1 +mod2 <- mod1$clone() + +# Then update it +mod2 <- mod2$update( + pri = list( + V0 = ab(30, 120), + V = NULL + ), + sec = function(x) { + V = V0*(WT/55) + }, + err = list( + proportional(2.39, c(0.02, 0.05, -0.0002, 0), fixed = TRUE) + ) +) +# we can also make a model object by loading a file +mod2b <- PM_model$new("../src/model2.txt") + + +run2 <- mod2$fit(data = exData, run = 2, overwrite = TRUE) + +# for future loading +run2 <- PM_load(2) + + + + +# EXERCISE 3 - COMPARING MODELS ------------------------------------------- + + +# Let's compare model 1 and model 2. You can compare any number of models. +# Type ?PM_compare for help. +PM_compare(run1, run2) + + + +# EXERCISE 4 - MODEL VALIDATION ------------------------------------------- + +# MODEL VALIDATION EXAMPLES +# Example of Pmetrics visual predictive check and prediction-corrected visual predictive check +# for model validation - be sure to have executed the NPAG run above +# Type ?makeValid in the R console for help. +# Choose wt as the covariate to bin. Accept all default bin sizes. +run2$validate(limits = c(0, 3)) + +# To see what it contains, use: +run2$valid + +# Default visual predictive check; ?plot.PM_valid for help +run2$valid$plot() + +# or old S3 +plot(run2$valid) + + +# Generate a prediction-corrected visual predictive check; type ?plot.PMvalid in the R console for help. +run2$valid$plot(type = "pcvpc") + +# Create an npde plot +run2$valid$plot(type = "npde") + +# Here is another way to generate a visual predicive check... +npc_2 <- run2$valid$simdata$plot(obs = run2$op, log = FALSE, binSize = 0.5) + +# The jagged appearance of the plot when binSize=0 is because different subjects have +# different doses, covariates, and observation times, which are all combined in one simulation. +# Collapsing simulation times within 1 hour bins (binSize=1) smooths +# the plot, but can change the P-values in the numerical predictive check below. + +npc_2 +# ...and here is a numerical predictive check +# P-values are binomial test of proportion of observations less than +# the respective quantile + + +# EXERCISE 5 - SIMULATOR RUN ---------------------------------------------- + +setwd(wd) +dir.create("../Sim") +setwd("../Sim") + +# The following will simulate 100 sets of parameters/concentrations using the +# first subject in the data file as a template. +# Limits are put on the simulated parameter ranges to be the same as in the model. +# The population parameter values from the NPAG run in exercise 2 are used for the Monte Carlo Simulation. +simdata <- run2$sim(include = 1, limits = NA, nsim = 100) + +# Below is the alternate way to simulate, which is particularly useful if you define +# your own population parameters. See ?SIMrun for details on this as well as +# the article on simulation linked by PMmanual(). +poppar <- list( + wt = 1, + mean = c(0.6, 0.05, 77.8, 1.2), + cov = diag(c(0.07, 0.0004, 830, 0.45)) +) + +simOther <- PM_sim$new(poppar = poppar, data = exData, model = mod1, + include = 1, limits = NA, nsim = 100) + + +# simulate from a model with new data +sim_new <- run2$sim( + data = "../src/ptaex1.csv", + include = 2, limits = NA, + predInt = c(120, 144, 0.5) +) + + + +# Plot them; ?plot.PM_sim for help +simdata$plot() +simOther$plot() +sim_new$plot(log = FALSE) + +# Simulate using multiple subjects as templates +simdata <- run2$sim(include = 1:4, limits = NA, nsim = 100) + +# Plot the third simulation; use include/exclude to specify the ID numbers, +# which are the same as the ID numbers in the template data file +simdata$plot(include = 2) + +# or in S3 +plot(simdata$data, include = 3) + +# Parse and combine multiple files and plot them. Note that combining simulations from templates +# with different simulated observation times can lead to unpredictable plots +simdata2 <- run2$sim(include = 1:4, limits = NA, nsim = 100) +simdata2$plot() + +# simulate with covariates +# in this case we use the covariate-parameter correlations from run 2, which +# are found in the cov.2 object; we re-define the mean weight to be 50 with +# SD of 20, and limits of 10 to 70 kg. We fix africa, gender and height covariates, +# but allow age (the last covariate) to be simulated, using the mean, sd, and +# limits in the original population, since we didn't specify them. +# See ?SIMrun for more help on this and the Pmetrics manual. + +covariate <- list( + cov = run2$cov, + mean = list(wt = 50), + sd = list(wt = 20), + limits = list(wt = c(10, 70)), + fix = c("africa", "gender", "height") +) + +# now simulate with this covariate list object +simdata3 <- run2$sim(include = 1:4, limits = NA, nsim = 100, covariate = covariate) + +# compare difference in simulations without covariates simulated... +# PM_sim's plot function defaults to the first simulation if there +# are multiple simulations and "at" is not specified. +simdata$plot() + +# ...and with covariates simulated +simdata3$plot() + +# Here are the simulated parameters and covariates for the first subject's +# template; note that both wt and age are simulated, using proper covariances +# with simulated PK parameters +simdata3$data$parValues + +# We can summarize simulations too. See ?summary.PM_sim for help. +simdata3$summary(field = "obs", include = 3) + +# look in the working directory and find the "c_simdata.csv" and "c_simmodel.txt" files +# which were made when you simulated with covariates. Compare to original +# "simdata.csv" and "simmoddel.txt" files to note that simulated covariates become +# Primary block variables, and are removed from the template data file. + +# EXERCISE 6 - SAVING PMETRICS OBJECTS ------------------------------------ + +setwd(wd) + +# The following objects have methods to save them to or load them from files: +# PM_fit +# PM_result +# PM_sim +# PM_pta + +# Example - save the PM_result (run2) to the "2" folder +run2$save(file = "2/outputs/run2.rds") # rds is the recommended file extension +list.files("2/outputs") +copy_run2 <- PM_load(file = "2/outputs/run2.rds") +copy_run2 + +# If you want to save multiple objects into one single file, R provides the +# following functionality + +save(exData, mod1, run1, simdata, file = "2/test_drug.Rdata") +list.files("2") +load("2/test_drug.Rdata") + +# or +save.image("2/workspace.Rdata") # This will save all variables in your environment +list.files("2") +load("2/workspace.Rdata") + + +# EXERCISE 7 - CONTINUING RUNS OR EXTERNAL VALIDATIONS -------------------- + +# Example of a run with a non-uniform density +# This is a good way to continue a previous run, +# in this case it continues where run 1 left off + +# note that we can supply a run number to model, data, and prior arguments. The numbers do not +# have to be the same. This will copy the appropriate files from the specified run to be used +# in the current run. By specifying a prior, we are starting with the non-uniform density from the +# end of the specified fun. +run3 <- mod2$fit(data = exData, prior = 2) +run3 <- PM_load(3) + +# We could also generate Bayesian posterior parameter estimates for a new population this way: +# run3 <- mod2$fit(data=PM_data("newPop.csv"), prior = 2, cycles = 0) +# This won't run because we don't have a newPop.csv file, +# but shows you how it could be done. + + + +# EXERCISE 8 - PROBABILITY OF TARGET ATTAINMENT --------------------------- + +# Note: these can be computationally intense and take some time. + +# Examples of probability of target attainment analysis +# Be sure to have executed the NPAG run above and used PM_load(2) in EXERCISE 2 +# Type ?PM_sim, ?PM_pta, or ?plot.PM_pta into the R console for help. + + +# simulate with the template data file that contains different doses +# Look at PM_sim for help on arguments to this function, including predInt, +# seed, limits, nsim. + +simlist1 <- run1$sim( + limits = c(0, 3), data = "../src/ptaex1.csv", + predInt = c(120, 144, 0.5), seed = rep(-17, 4) +) + +# now simulate with covariates; make sure that you defined the covariate +# object first in Exercise 5 above and have loaded the results of Exercise 2 +# with PM_load(2) +simlist2 <- run2$sim( + limits = 5, data = "../src/ptaex1.csv", + predInt = c(120, 144, 0.5), seed = rep(-17, 4), + covariate = covariate +) + +# make the first PMpta object to calculate the time above each target for at +# least 60% of the dosing +# interval from 120 to 144 hours. Include labels for the simulations. +# ?makePTA for help +# define simulation labels first + +simlabels <- c("600 mg daily", "1200 mg daily", "300 mg bid", "600 mg bid") + +pta1_2 <- PM_pta$new( + simdata = simlist1, + target = c(0.25, 0.5, 1, 2, 4, 8, 16, 32), target_type = "time", + success = 0.6, start = 120, end = 144 +) + +pta1b_2 <- PM_pta$new( + simdata = simlist2, + simlabels = simlabels, + target = c(0.25, 0.5, 1, 2, 4, 8, 16, 32), target_type = "time", + success = 0.6, start = 120, end = 144 +) + +# summarize the results +pta1_2$summary() + +# in the summary()$pta, reg_num is the simulation template ID number; +# target in this case is the MIC; prop_success is the proportion of the simulated +# profiles for each dose/MIC that are above the success threshold (0.6); pdi.mean and pdi.sd +# are the mean and standard deviation of the pharmacodynamic index (PDI), in this case proportion of the interval > MIC. +# In the $pdi, target and simnum are the same, but now the median and confidence +# intervals (default 95%) PDI are shown. +# ?summary.PMpta for help + +# Plot the first without covariates. We didn't include simulation +# labels in the makePTA command, so generics are used here, but we move it to +# the bottom left; ?legend for help on arguments to supply to the +# legend list argument to plot.PMpta. +pta1_2$plot(ylab = "Proportion with %T>MIC of at least 60%", grid = TRUE, legend = list(x = "bottomleft")) + +pta1b_2$summary() + +# Plot the second with covariates simulated. Note the regimen labels are included, but we move +# the legend to the bottom left. +pta1b_2$plot( + ylab = "Proportion with %T>MIC of at least 60%", grid = TRUE, + legend = list(x = "bottomleft") +) + +# Now we'll define success as free auc:mic > 100 with a free drug fraction of 50% +pta2_2 <- PM_pta$new( + simdata = simlist2, + simlabels = simlabels, target = c(0.25, 0.5, 1, 2, 4, 8, 16, 32), + free_fraction = 0.7, + target_type = "auc", success = 100, start = 120, end = 144 +) +summary(pta2_2) +pta2_2$plot( + ylab = "Proportion with AUC/MIC of at least 100", grid = TRUE, + legend = list(x = "bottomleft") +) + +# success is Cmax/MIC >=10 +pta3_2 <- PM_pta$new( + simdata = simlist2, + simlabels = simlabels, + target = c(0.25, 0.5, 1, 2, 4, 8, 16, 32), + target_type = "peak", success = 10, start = 120, end = 144 +) +pta3_2$summary() +pta3_2$plot(ylab = "Proportion with peak/MIC of at least 10", grid = TRUE) + +# success = Cmin:MIC > 1 +pta4_2 <- PM_pta$new( + simdata = simlist2, + simlabels = simlabels, + target = c(0.25, 0.5, 1, 2, 4, 8, 16, 32), + target_type = "min", success = 1, start = 120, end = 144 +) +pta4_2$summary() +pta4_2$plot(ylab = "Proportion with Cmin/MIC of at least 1", grid = TRUE, legend = list(x = "bottomleft")) + +# now plot the PDI (pharmacodynamic index) of each regimen, rather than the proportion +# of successful profiles. A PDI plot is always available for PMpta objects. +pta4_2$plot(at = 1, type = "pdi", ylab = "Cmin:MIC", grid = TRUE) + +# Each regimen has the 90% confidence interval PDI around the median curve, +# in the corresponding, semi-transparent color. Make the CI much narrower... +pta4_2$plot(at = 1, type = "pdi", ci = 0.1) + +# ...or gone altogether, put back the grid, redefine the colors, and make lines narrower +pta4_2$plot( + + at = 1, type = "pdi", ci = 0, grid = TRUE, + line = list( + color = c("blue", "purple", "black", "brown"), + width = 1 + ) +) + +# now let's repeat the analysis but simulate the distribution of MICs +# using susceptibility of Staphylococcus aureus to vancomycin contained +# in the mic1 dataset within Pmetrics (?mic1) + +# see the source with ?mic1 +pta4b_2 <- PM_pta$new( + simdata = simlist2, + simlabels = c("600 mg daily", "1200 mg daily", "300 mg bid", "600 mg bid"), + target = makePTAtarget(mic1), target_type = "min", success = 1, start = 120, end = 144 +) + +pta4b_2$summary() +# plot it +pta4b_2$plot( + grid = TRUE, ylab = "Proportion with Cmin/MIC of at least 1", + marker = list(color = "red"), line = list(color = "black") +) +pta4b_2$plot(type = "pdi", grid = TRUE, ylab = "Proportion with Cmin/MIC of at least 1") + +# note that the plot changes since target MICs are no longer discrete +# since most of the MICs are very low, the regimens all look very similar + +# success = concentration at time 3 hours:MIC > 2 +pta5_2 <- PM_pta$new( + simdata = simlist2, + simlabels = simlabels, + target = c(0.25, 0.5, 1, 2, 4, 8, 16, 32), target_type = 123, success = 2, start = 120, end = 144 +) +pta5_2$summary() +pta5_2$plot(ylab = "Proportion with C3/MIC of at least 1", grid = TRUE, legend = list(x = .3, y = 0.1)) + + +# success is trough >10 +pta6_2 <- PM_pta$new( + simdata = simlist2, + simlabels = simlabels, + target = 10, target_type = 144, success = 1, start = 120, end = 144 +) +plot(pta6_2) +pta6_2$summary() + +# EXERCISE 10 - OPTIMAL SAMPLE TIMES -------------------------------------- + +setwd(wd) +dir.create("../MMopt") +setwd("../MMopt") + +# calculate MM-optimal sample times for Run 2, and the 1200 mg once daily dose in the PTA +# By specifying the predInt to start and stop at 120 and 144 hours, with an interval of 1 hour, +# we are sampling at steady state. Including "subject 2", means only the 1200 mg once daily dose +# will serve as a simulation template. + +run2$opt( + data = "../src/ptaex1.csv", + nsamp = 2, predInt = c(120, 140, 1), + limits = NA, + include = 2 +) +# see the optimal sample times and the Bayes Risk of misclassification, +# which is only useful to compare optimal sampling regimens, i.e. the +# absolute value is less helpful, but is the statistic minimized by the +# selected optimal sample times for a given model + +mmopt_2 <- PM_opt$new( + poppar = run2$final, + model = run2$model, + data = "../src/ptaex1.csv", + nsamp = 2, predInt = c(120, 140, 1), + include = 2 +) + +# plot it, with the red lines indicating the optimal sample times. +# see ?plot.MMopt for help +mmopt_2$plot() +plot(mmopt_2) +plot(mmopt_2, line = list(color = "slategrey"), times = list(color = "orange")) + +# EXERCISE 11 - ASSAY ERROR ----------------------------------------------- +# see ?makeErrorPoly for more help +# This will let you choose the best set of C0, C1, C2, C3 for your modeling, +# based on assay validation data which includes the "obs", which are the +# nominal concentrations of the standards, and "sd", which is the standard +# deviation of replicate measurements of each of the standards, i.e. the +# inter-day and/or intra-day standard deviation + +obs <- c(0, 25, 50, 100, 250, 500, 1000, 2000, 5000) +sd <- c(0.5, 6.4, 8.6, 12, 8.6, 37.2, 60.1, 165.7, 483) + +# See plots.pdf, page 50 +makeErrorPoly(obs = obs, sd = sd) + +# choose the one with the best R-squared that will never result in a +# negative value for the SD + + + +# Ancillary functions ----------------------------------------------------- + +# Be sure to check out the help files for the following functions: +# +# makeAUC() - calculate AUC from a variety of inputs +# makeNCA() - non-compartmental analysis +# NM2PM() - convert NONMEM data files to Pmetrics data files +# qgrowth() - CDC growth charts +# zBMI() - CDC Pediatric BMI z-scores and percentiles +# ss.PK() - sample size for Phase 1 PK studies diff --git a/Examples/Runs/.gitkeep b/Examples/Runs/.gitkeep new file mode 100644 index 000000000..e69de29bb diff --git a/Examples/src/bad.csv b/Examples/src/bad.csv new file mode 100644 index 000000000..0d644a159 --- /dev/null +++ b/Examples/src/bad.csv @@ -0,0 +1 @@ +POPDATA DEC_11,,,,,,,,,,,,,,,,,, #ID,EVID,TIME,DUR,DOSE,ADDL,II,INPUT,OUT,OUTEQ,C0,C1,C2,C3,WT,AFRICA,AGE,GENDER,HEIGHT 1,1,0,1,600,.,.,1,.,..,.,.,.,.,46.7,1,21,1,160 1,.,24,0,600,.,.,1,.,.,.,.,.,.,46.7,1,21,1,160 1,1,48,0,600,.,.,1,.,.,.,.,.,.,46.7,1,21,1,160 1,1,72,0,600,.,.,1,.,.,.,.,.,.,46.7,1,21,1,160 1,1,96,0,600,.,.,1,.,.,.,.,.,.,46.7,1,21,1,160 1,0,120,.,.,.,.,.,10.44,1,0.02,0.0506,-0.0002,0,46.7,1,21,1,160 1,1,120,0,600,.,.,1,.,.,.,.,.,.,46.7,1,21,1,160 1,0,121,.,.,.,.,.,12.89,1,0.02,0.0506,-0.0002,0,46.7,1,21,1,160 1,0,122,.,.,.,.,.,14.98,1,0.02,0.0506,-0.0002,0,46.7,1,21,1,160 1,0,125.99,.,.,.,.,.,16.69,1,0.02,0.0506,-0.0002,0,46.7,1,21,1,160 1,0,129,.,.,.,.,.,20.15,1,0.02,0.0506,-0.0002,0,46.7,1,21,1,160 1,0,132,.,.,.,.,.,14.97,1,0.02,0.0506,-0.0002,0,46.7,1,21,1,160 1,0,143.98,.,.,.,.,.,12.57,1,0.02,0.0506,-0.0002,0,46.7,1,21,1,160 2,1,0,0,600,.,.,1,.,.,.,.,.,.,66.5,1,30,1,174 2,1,24,0,600,.,.,1,.,.,.,.,.,.,66.5,1,30,1,174 2,1,48,0,600,.,.,1,.,.,.,.,.,.,66.5,1,30,1,174 2,1,72,0,600,.,.,1,.,.,.,.,.,.,66.5,1,30,1,174 2,1,96,0,600,.,.,1,.,.,.,.,.,.,66.5,1,30,1,174 2,0,120,.,.,.,.,.,3.56,1,0.02,0.0506,-0.0002,0,66.5,1,30,1,174 2,1,120,0,600,.,.,1,.,.,.,.,.,.,66.5,1,30,1,174 2,0,120.98,.,.,.,.,.,5.84,1,0.02,0.0506,-0.0002,0,66.5,1,30,1,174 2,0,121.98,.,.,.,.,.,6.54,1,0.02,0.0506,-0.0002,0,66.5,1,30,1,174 2,0,126,.,.,.,.,.,6.14,1,0.02,0.0506,-0.0002,0,66.5,1,30,1,174 2,0,129.02,.,.,.,.,.,6.56,1,0.02,0.0506,-0.0002,0,66.5,1,30,1,174 2,0,132.02,.,.,.,.,.,4.44,1,0.02,0.0506,-0.0002,0,66.5,1,30,1,174 2,0,144,.,.,.,.,.,3.76,1,0.02,0.0506,-0.0002,0,66.5,1,30,1,174 3,1,0,0,600,.,.,1,.,.,.,.,.,.,46.7,1,24,0,164 3,1,24,0,600,.,.,1,.,.,.,.,.,.,46.7,1,24,0,164 3,1,48,0,600,.,.,1,.,.,.,.,.,.,46.7,1,24,0,164 3,1,72,0,600,.,.,1,.,.,.,.,.,.,46.7,1,24,0,164 3,1,96,0,600,.,.,1,.,.,.,.,.,.,46.7,1,24,0,164 3,1,120,0,600,.,.,1,.,.,.,.,.,.,46.7,1,24,0,164 3,0,120.08,.,.,.,.,.,4.06,1,0.02,0.0506,-0.0002,0,46.7,1,24,0,164 3,0,121.07,.,.,.,.,.,3.24,1,0.02,0.0506,-0.0002,0,46.7,1,24,0,164 3,0,122.08,.,.,.,.,.,3.09,1,0.02,0.0506,-0.0002,0,46.7,1,24,0,164 3,0,126.08,.,.,.,.,.,7.98,1,0.02,0.0506,-0.0002,0,46.7,1,24,0,164 3,0,129.05,.,.,.,.,.,7.23,1,0.02,0.0506,-0.0002,0,46.7,1,24,0,164 3,0,132.1,.,.,.,.,.,4.71,1,0.02,0.0506,-0.0002,0,46.7,1,24,0,164 3,0,144.08,.,.,.,.,.,3.82,1,0.02,0.0506,-0.0002,0,46.7,1,24,0,164 4,1,0,0,600,.,.,1,.,.,.,.,.,.,50.8,1,25,1,165 4,1,24,0,600,.,.,1,.,.,.,.,.,.,50.8,1,25,1,165 4,1,48,0,600,.,.,1,.,.,.,.,.,.,50.8,1,25,1,165 4,1,72,0,600,.,.,1,.,.,.,.,.,.,50.8,1,25,1,165 4,1,96,0,600,.,.,1,.,.,.,.,.,.,50.8,1,25,1,165 4,0,120,.,.,.,.,.,2.1,1,0.02,0.0506,-0.0002,0,50.8,1,25,1,165 4,1,120,0,600,.,.,1,.,.,.,.,.,.,50.8,1,25,1,165 4,0,121,.,.,.,.,.,3.05,1,0.02,0.0506,-0.0002,0,50.8,1,25,1,165 4,0,122.02,.,.,.,.,.,5.21,1,0.02,0.0506,-0.0002,0,50.8,1,25,1,165 4,0,126,.,.,.,.,.,5.09,1,0.02,0.0506,-0.0002,0,50.8,1,25,1,165 4,0,129.03,.,.,.,.,.,4.24,1,0.02,0.0506,-0.0002,0,50.8,1,25,1,165 4,0,132,.,.,.,.,.,3.69,1,0.02,0.0506,-0.0002,0,50.8,1,25,1,165 4,0,144.02,.,.,.,.,.,1.96,1,0.02,0.0506,-0.0002,0,50.8,1,25,1,165 5,1,0,0,600,.,.,1,.,.,.,.,.,.,65.8,1,22,1,181 5,1,24,0,600,.,.,1,.,.,.,.,.,.,65.8,1,22,1,181 5,1,48,0,600,.,.,1,.,.,.,.,.,.,65.8,1,22,1,181 5,1,72,0,600,.,.,1,.,.,.,.,.,.,65.8,1,22,1,181 5,1,96,0,600,.,.,1,.,.,.,.,.,.,65.8,1,22,1,181 5,0,120,.,.,.,.,.,2.93,1,0.02,0.0506,-0.0002,0,65.8,1,22,1,181 5,1,120,0,600,.,.,1,.,.,.,.,.,.,65.8,1,22,1,181 5,0,121,.,.,.,.,.,2.64,1,0.02,0.0506,-0.0002,0,65.8,1,22,1,181 5,0,122,.,.,.,.,.,4.8,1,0.02,0.0506,-0.0002,0,65.8,1,22,1,181 5,0,126,.,.,.,.,.,3.7,1,0.02,0.0506,-0.0002,0,65.8,1,22,1,181 5,0,129.02,.,.,.,.,.,4.13,1,0.02,0.0506,-0.0002,0,65.8,1,22,1,181 5,0,132,.,.,.,.,.,2.81,1,0.02,0.0506,-0.0002,0,65.8,1,22,1,181 5,0,144,.,.,.,.,.,2.21,1,0.02,0.0506,-0.0002,0,65.8,1,22,1,181 6,1,0,0,600,.,.,1,.,.,.,.,.,.,65,1,23,1,177 6,1,24,0,600,.,.,1,.,.,.,.,.,.,65,1,23,1,177 6,1,48,0,600,.,.,1,.,.,.,.,.,.,65,1,23,1,177 6,1,72,0,600,.,.,1,.,.,.,.,.,.,65,1,23,1,177 6,1,96,0,600,.,.,1,.,.,.,.,.,.,65,1,23,1,177 6,0,120,.,.,.,.,.,6.92,1,0.02,0.0506,-0.0002,0,65,1,23,1,177 6,1,120,0,600,.,.,1,.,.,.,.,.,.,65,1,23,1,177 6,0,121,.,.,.,.,.,6.89,1,0.02,0.0506,-0.0002,0,65,1,23,1,177 6,0,121.98,.,.,.,.,.,6.64,1,0.02,0.0506,-0.0002,0,65,1,23,1,177 6,0,126,.,.,.,.,.,13.72,1,0.02,0.0506,-0.0002,0,65,1,23,1,177 6,0,129,.,.,.,.,.,12.69,1,0.02,0.0506,-0.0002,0,65,1,23,1,177 6,0,131.98,.,.,.,.,.,10.58,1,0.02,0.0506,-0.0002,0,65,1,23,1,177 6,0,144.98,.,.,.,.,.,6.62,1,0.02,0.0506,-0.0002,0,65,1,23,1,177 7,1,0,0,600,.,.,1,.,.,.,.,.,.,51.7,1,27,0,161 7,1,24,0,600,.,.,1,.,.,.,.,.,.,51.7,1,27,0,161 7,1,48,0,600,.,.,1,.,.,.,.,.,.,51.7,1,27,0,161 7,1,72,0,600,.,.,1,.,.,.,.,.,.,51.7,1,27,0,161 7,1,96,0,600,.,.,1,.,.,.,.,.,.,51.7,1,27,0,161 7,0,120,.,.,.,.,.,5.41,1,0.02,0.0506,-0.0002,0,51.7,1,27,0,161 7,1,120,0,600,.,.,1,.,.,.,.,.,.,51.7,1,27,0,161 7,0,121.03,.,.,.,.,.,4.46,1,0.02,0.0506,-0.0002,0,51.7,1,27,0,161 7,0,122.03,.,.,.,.,.,4.54,1,0.02,0.0506,-0.0002,0,51.7,1,27,0,161 7,0,126.02,.,.,.,.,.,12.19,1,0.02,0.0506,-0.0002,0,51.7,1,27,0,161 7,0,129.08,.,.,.,.,.,12.1,1,0.02,0.0506,-0.0002,0,51.7,1,27,0,161 7,0,132.03,.,.,.,.,.,8.61,1,0.02,0.0506,-0.0002,0,51.7,1,27,0,161 7,0,144.03,.,.,.,.,.,6.37,1,0.02,0.0506,-0.0002,0,51.7,1,27,0,161 8,1,0,0,600,.,.,1,.,.,.,.,.,.,51.2,1,22,1,163 8,1,24,0,600,.,.,1,.,.,.,.,.,.,51.2,1,22,1,163 8,1,48,0,600,.,.,1,.,.,.,.,.,.,51.2,1,22,1,163 8,1,72,0,600,.,.,1,.,.,.,.,.,.,51.2,1,22,1,163 8,1,96,0,600,.,.,1,.,.,.,.,.,.,51.2,1,22,1,163 8,0,120,.,.,.,.,.,6.19,1,0.02,0.0506,-0.0002,0,51.2,1,22,1,163 8,1,120,0,600,.,.,1,.,.,.,.,.,.,51.2,1,22,1,163 8,0,121.03,.,.,.,.,.,6.33,1,0.02,0.0506,-0.0002,0,51.2,1,22,1,163 8,0,122,.,.,.,.,.,6.24,1,0.02,0.0506,-0.0002,0,51.2,1,22,1,163 8,0,125.98,.,.,.,.,.,13.03,1,0.02,0.0506,-0.0002,0,51.2,1,22,1,163 8,0,128.98,.,.,.,.,.,11.86,1,0.02,0.0506,-0.0002,0,51.2,1,22,1,163 8,0,132,.,.,.,.,.,11.45,1,0.02,0.0506,-0.0002,0,51.2,1,22,1,163 8,0,143.98,.,.,.,.,.,7.83,1,0.02,0.0506,-0.0002,0,51.2,1,22,1,163 9,1,0,0,600,.,.,1,.,.,.,.,.,.,55,1,23,1,174 9,1,24,0,600,.,.,1,.,.,.,.,.,.,55,1,23,1,174 9,1,48,0,600,.,.,1,.,.,.,.,.,.,55,1,23,1,174 9,1,72,0,600,.,.,1,.,.,.,.,.,.,55,1,23,1,174 9,1,96,0,600,.,.,1,.,.,.,.,.,.,55,1,23,1,174 9,0,120,.,.,.,.,.,2.85,1,0.02,0.0506,-0.0002,0,55,1,23,1,174 9,1,120,0,600,.,.,1,.,.,.,.,.,.,55,1,23,1,174 9,0,120.97,.,.,.,.,.,3.7,1,0.02,0.0506,-0.0002,0,55,1,23,1,174 9,0,122,.,.,.,.,.,6.65,1,0.02,0.0506,-0.0002,0,55,1,23,1,174 9,0,125.98,.,.,.,.,.,6.81,1,0.02,0.0506,-0.0002,0,55,1,23,1,174 9,0,128.98,.,.,.,.,.,6.51,1,0.02,0.0506,-0.0002,0,55,1,23,1,174 9,0,132,.,.,.,.,.,7.48,1,0.02,0.0506,-0.0002,0,55,1,23,1,174 9,0,143.98,.,.,.,.,.,4.51,1,0.02,0.0506,-0.0002,0,55,1,23,1,174 10,1,0,0,600,.,.,1,.,.,.,.,.,.,52.1,1,32,1,163 10,1,24,0,600,.,.,1,.,.,.,.,.,.,52.1,1,32,1,163 10,1,48,0,600,.,.,1,.,.,.,.,.,.,52.1,1,32,1,163 10,1,72,0,600,.,.,1,.,.,.,.,.,.,52.1,1,32,1,163 10,1,96,0,600,.,.,1,.,.,.,.,.,.,52.1,1,32,1,163 10,0,120,.,.,.,.,.,2.93,1,0.02,0.0506,-0.0002,0,52.1,1,32,1,163 10,1,120,0,600,.,.,1,.,.,.,.,.,.,52.1,1,32,1,163 10,0,121,.,.,.,.,.,4.36,1,0.02,0.0506,-0.0002,0,52.1,1,32,1,163 10,0,122.02,.,.,.,.,.,7.79,1,0.02,0.0506,-0.0002,0,52.1,1,32,1,163 10,0,126,.,.,.,.,.,11.02,1,0.02,0.0506,-0.0002,0,52.1,1,32,1,163 10,0,129,.,.,.,.,.,8.86,1,0.02,0.0506,-0.0002,0,52.1,1,32,1,163 10,0,131.97,.,.,.,.,.,6.09,1,0.02,0.0506,-0.0002,0,52.1,1,32,1,163 10,0,144,.,.,.,.,.,4.15,1,0.02,0.0506,-0.0002,0,52.1,1,32,1,163 11,1,0,0,600,.,.,1,.,.,.,.,.,.,56.5,1,34,1,165 11,1,24,0,600,.,.,1,.,.,.,.,.,.,56.5,1,34,1,165 11,1,48,0,600,.,.,1,.,.,.,.,.,.,56.5,1,34,1,165 11,1,72,0,600,.,.,1,.,.,.,.,.,.,56.5,1,34,1,165 11,1,96,0,600,.,.,1,.,.,.,.,.,.,56.5,1,34,1,165 11,0,120,.,.,.,.,.,2.09,1,0.02,0.0506,-0.0002,0,56.5,1,34,1,165 11,1,120,0,600,.,.,1,.,.,.,.,.,.,56.5,1,34,1,165 11,0,121.03,.,.,.,.,.,2.68,1,0.02,0.0506,-0.0002,0,56.5,1,34,1,165 11,0,122,.,.,.,.,.,4.71,1,0.02,0.0506,-0.0002,0,56.5,1,34,1,165 11,0,125.98,.,.,.,.,.,7.71,1,0.02,0.0506,-0.0002,0,56.5,1,34,1,165 11,0,129,.,.,.,.,.,6.31,1,0.02,0.0506,-0.0002,0,56.5,1,34,1,165 11,0,132,.,.,.,.,.,5.82,1,0.02,0.0506,-0.0002,0,56.5,1,34,1,165 11,0,144.13,.,.,.,.,.,2.63,1,0.02,0.0506,-0.0002,0,56.5,1,34,1,165 12,1,0,0,600,.,.,1,.,.,.,.,.,.,47.9,1,54,0,160 12,1,24,0,600,.,.,1,.,.,.,.,.,.,47.9,1,54,0,160 12,1,48,0,600,.,.,1,.,.,.,.,.,.,47.9,1,54,0,160 12,1,72,0,600,.,.,1,.,.,.,.,.,.,47.9,1,54,0,160 12,1,96,0,600,.,.,1,.,.,.,.,.,.,47.9,1,54,0,160 12,0,120,.,.,.,.,.,7.09,1,0.02,0.0506,-0.0002,0,47.9,1,54,0,160 12,1,120,0,600,.,.,1,.,.,.,.,.,.,47.9,1,54,0,160 12,0,121.03,.,.,.,.,.,6.18,1,0.02,0.0506,-0.0002,0,47.9,1,54,0,160 12,0,122.13,.,.,.,.,.,8.66,1,0.02,0.0506,-0.0002,0,47.9,1,54,0,160 12,0,126,.,.,.,.,.,11.16,1,0.02,0.0506,-0.0002,0,47.9,1,54,0,160 12,0,129,.,.,.,.,.,9.51,1,0.02,0.0506,-0.0002,0,47.9,1,54,0,160 12,0,132,.,.,.,.,.,8.14,1,0.02,0.0506,-0.0002,0,47.9,1,54,0,160 12,0,144,.,.,.,.,.,7.89,1,0.02,0.0506,-0.0002,0,47.9,1,54,0,160 13,1,0,0,600,.,.,1,.,.,.,.,.,.,60.5,1,24,1,180 13,1,24,0,600,.,.,1,.,.,.,.,.,.,60.5,1,24,1,180 13,1,48,0,600,.,.,1,.,.,.,.,.,.,60.5,1,24,1,180 13,1,72,0,600,.,.,1,.,.,.,.,.,.,60.5,1,24,1,180 13,1,96,0,600,.,.,1,.,.,.,.,.,.,60.5,1,24,1,180 13,0,120,.,.,.,.,.,6.62,1,0.02,0.0506,-0.0002,0,60.5,1,24,1,180 13,1,120,0,600,.,.,1,.,.,.,.,.,.,60.5,1,24,1,180 13,0,121,.,.,.,.,.,3.18,1,0.02,0.0506,-0.0002,0,60.5,1,24,1,180 13,0,122,.,.,.,.,.,5.41,1,0.02,0.0506,-0.0002,0,60.5,1,24,1,180 13,0,126,.,.,.,.,.,10.18,1,0.02,0.0506,-0.0002,0,60.5,1,24,1,180 13,0,129.02,.,.,.,.,.,12.84,1,0.02,0.0506,-0.0002,0,60.5,1,24,1,180 13,0,132,.,.,.,.,.,12.35,1,0.02,0.0506,-0.0002,0,60.5,1,24,1,180 13,0,144,.,.,.,.,.,8.06,1,0.02,0.0506,-0.0002,0,60.5,1,24,1,180 14,1,0,0,600,.,.,1,.,.,.,.,.,.,59.2,1,26,1,174 14,1,24,0,600,.,.,1,.,.,.,.,.,.,59.2,1,26,1,174 14,1,48,0,600,.,.,1,.,.,.,.,.,.,59.2,1,26,1,174 14,1,72,0,600,.,.,1,.,.,.,.,.,.,59.2,1,26,1,174 14,1,96,0,600,.,.,1,.,.,.,.,.,.,59.2,1,26,1,174 14,0,120,.,.,.,.,.,3.63,1,0.02,0.0506,-0.0002,0,59.2,1,26,1,174 14,1,120,0,600,.,.,1,.,.,.,.,.,.,59.2,1,26,1,174 14,0,121,.,.,.,.,.,4.49,1,0.02,0.0506,-0.0002,0,59.2,1,26,1,174 14,0,122,.,.,.,.,.,5.5,1,0.02,0.0506,-0.0002,0,59.2,1,26,1,174 14,0,126,.,.,.,.,.,7.28,1,0.02,0.0506,-0.0002,0,59.2,1,26,1,174 14,0,129,.,.,.,.,.,5.27,1,0.02,0.0506,-0.0002,0,59.2,1,26,1,174 14,0,132,.,.,.,.,.,4.89,1,0.02,0.0506,-0.0002,0,59.2,1,26,1,174 14,0,144,.,.,.,.,.,2.68,1,0.02,0.0506,-0.0002,0,59.2,1,26,1,174 15,1,0,0,450,.,.,1,.,.,.,.,.,.,43,1,19,0,150 15,1,24,0,450,.,.,1,.,.,.,.,.,.,43,1,19,0,150 15,1,48,0,450,.,.,1,.,.,.,.,.,.,43,1,19,0,150 15,1,72,0,450,.,.,1,.,.,.,.,.,.,43,1,19,0,150 15,1,96,0,450,.,.,1,.,.,.,.,.,.,43,1,19,0,150 15,0,120,.,.,.,.,.,5.53,1,0.02,0.0506,-0.0002,0,43,1,19,0,150 15,1,120,0,450,.,.,1,.,.,.,.,.,.,43,1,19,0,150 15,0,121,.,.,.,.,.,4.81,1,0.02,0.0506,-0.0002,0,43,1,19,0,150 15,0,122,.,.,.,.,.,8.14,1,0.02,0.0506,-0.0002,0,43,1,19,0,150 15,0,126,.,.,.,.,.,9.96,1,0.02,0.0506,-0.0002,0,43,1,19,0,150 15,0,129,.,.,.,.,.,8.55,1,0.02,0.0506,-0.0002,0,43,1,19,0,150 15,0,132.05,.,.,.,.,.,7.54,1,0.02,0.0506,-0.0002,0,43,1,19,0,150 15,0,144.05,.,.,.,.,.,5.74,1,0.02,0.0506,-0.0002,0,43,1,19,0,150 16,1,0,0,600,.,.,1,.,.,.,.,.,.,64.4,1,25,1,173 16,1,24,0,600,.,.,1,.,.,.,.,.,.,64.4,1,25,1,173 16,1,48,0,600,.,.,1,.,.,.,.,.,.,64.4,1,25,1,173 16,1,72,0,600,.,.,1,.,.,.,.,.,.,64.4,1,25,1,173 16,1,96,0,600,.,.,1,.,.,.,.,.,.,64.4,1,25,1,173 16,0,120,.,.,.,.,.,5.48,1,0.02,0.0506,-0.0002,0,64.4,1,25,1,173 16,1,120,0,600,.,.,1,.,.,.,.,.,.,64.4,1,25,1,173 16,0,121,.,.,.,.,.,6.59,1,0.02,0.0506,-0.0002,0,64.4,1,25,1,173 16,0,122,.,.,.,.,.,8.91,1,0.02,0.0506,-0.0002,0,64.4,1,25,1,173 16,0,126,.,.,.,.,.,10.57,1,0.02,0.0506,-0.0002,0,64.4,1,25,1,173 16,0,129,.,.,.,.,.,9.52,1,0.02,0.0506,-0.0002,0,64.4,1,25,1,173 16,0,132,.,.,.,.,.,7.83,1,0.02,0.0506,-0.0002,0,64.4,1,25,1,173 16,0,143.97,.,.,.,.,.,4.96,1,0.02,0.0506,-0.0002,0,64.4,1,25,1,173 17,1,0,0,600,.,.,1,.,.,.,.,.,.,54.8,1,23,1,170 17,1,24,0,600,.,.,1,.,.,.,.,.,.,54.8,1,23,1,170 17,1,48,0,600,.,.,1,.,.,.,.,.,.,54.8,1,23,1,170 17,1,72,0,600,.,.,1,.,.,.,.,.,.,54.8,1,23,1,170 17,1,96,0,600,.,.,1,.,.,.,.,.,.,54.8,1,23,1,170 17,0,120,.,.,.,.,.,2.11,1,0.02,0.0506,-0.0002,0,54.8,1,23,1,170 17,1,120,0,600,.,.,1,.,.,.,.,.,.,54.8,1,23,1,170 17,0,121.02,.,.,.,.,.,1.86,1,0.02,0.0506,-0.0002,0,54.8,1,23,1,170 17,0,122.02,.,.,.,.,.,6.92,1,0.02,0.0506,-0.0002,0,54.8,1,23,1,170 17,0,126,.,.,.,.,.,9.11,1,0.02,0.0506,-0.0002,0,54.8,1,23,1,170 17,0,129,.,.,.,.,.,6.96,1,0.02,0.0506,-0.0002,0,54.8,1,23,1,170 17,0,132,.,.,.,.,.,5.64,1,0.02,0.0506,-0.0002,0,54.8,1,23,1,170 17,0,144.08,.,.,.,.,.,3.59,1,0.02,0.0506,-0.0002,0,54.8,1,23,1,170 18,1,0,0,450,.,.,1,.,.,.,.,.,.,44.3,1,20,0,164 18,1,24,0,450,.,.,1,.,.,.,.,.,.,44.3,1,20,0,164 18,1,48,0,450,.,.,1,.,.,.,.,.,.,44.3,1,20,0,164 18,1,72,0,450,.,.,1,.,.,.,.,.,.,44.3,1,20,0,164 18,1,96,0,450,.,.,1,.,.,.,.,.,.,44.3,1,20,0,164 18,0,120,.,.,.,.,.,7.95,1,0.02,0.0506,-0.0002,0,44.3,1,20,0,164 18,1,120,0,450,.,.,1,.,.,.,.,.,.,44.3,1,20,0,164 18,0,120.98,.,.,.,.,.,7.47,1,0.02,0.0506,-0.0002,0,44.3,1,20,0,164 18,0,121.98,.,.,.,.,.,8.67,1,0.02,0.0506,-0.0002,0,44.3,1,20,0,164 18,0,126,.,.,.,.,.,13.83,1,0.02,0.0506,-0.0002,0,44.3,1,20,0,164 18,0,129.17,.,.,.,.,.,14.01,1,0.02,0.0506,-0.0002,0,44.3,1,20,0,164 18,0,132.17,.,.,.,.,.,8.97,1,0.02,0.0506,-0.0002,0,44.3,1,20,0,164 18,0,143.97,.,.,.,.,.,8.4,1,0.02,0.0506,-0.0002,0,44.3,1,20,0,164 19,1,0,0,600,.,.,1,.,.,.,.,.,.,50,1,36,1,168 19,1,24,0,600,.,.,1,.,.,.,.,.,.,50,1,36,1,168 19,1,48,0,600,.,.,1,.,.,.,.,.,.,50,1,36,1,168 19,1,72,0,600,.,.,1,.,.,.,.,.,.,50,1,36,1,168 19,1,96,0,600,.,.,1,.,.,.,.,.,.,50,1,36,1,168 19,0,120,.,.,.,.,.,5.42,1,0.02,0.0506,-0.0002,0,50,1,36,1,168 19,1,120,0,600,.,.,1,.,.,.,.,.,.,50,1,36,1,168 19,0,121,.,.,.,.,.,7.08,1,0.02,0.0506,-0.0002,0,50,1,36,1,168 19,0,122,.,.,.,.,.,7.27,1,0.02,0.0506,-0.0002,0,50,1,36,1,168 19,0,125.98,.,.,.,.,.,20.07,1,0.02,0.0506,-0.0002,0,50,1,36,1,168 19,0,128.98,.,.,.,.,.,18.24,1,0.02,0.0506,-0.0002,0,50,1,36,1,168 19,0,132,.,.,.,.,.,15.36,1,0.02,0.0506,-0.0002,0,50,1,36,1,168 19,0,144,.,.,.,.,.,10.92,1,0.02,0.0506,-0.0002,0,50,1,36,1,168 20,1,0,0,600,.,.,1,.,.,.,.,.,.,59,1,31,1,170 20,1,24,0,600,.,.,1,.,.,.,.,.,.,59,1,31,1,170 20,1,48,0,600,.,.,1,.,.,.,.,.,.,59,1,31,1,170 20,1,72,0,600,.,.,1,.,.,.,.,.,.,59,1,31,1,170 20,1,96,0,600,.,.,1,.,.,.,.,.,.,59,1,31,1,170 20,0,120,.,.,.,.,.,4.71,1,0.02,0.0506,-0.0002,0,59,1,31,1,170 20,1,120,0,600,.,.,1,.,.,.,.,.,.,59,1,31,1,170 20,0,120.77,.,.,.,.,.,4.5,1,0.02,0.0506,-0.0002,0,59,1,31,1,170 20,0,121.75,.,.,.,.,.,3.35,1,0.02,0.0506,-0.0002,0,59,1,31,1,170 20,0,125.67,.,.,.,.,.,12.35,1,0.02,0.0506,-0.0002,0,59,1,31,1,170 20,0,128.67,.,.,.,.,.,11.56,1,0.02,0.0506,-0.0002,0,59,1,31,1,170 20,0,143.67,.,.,.,.,.,6.45,1,0.02,0.0506,-0.0002,0,59,1,31,1,170 \ No newline at end of file diff --git a/Examples/src/ex.csv b/Examples/src/ex.csv new file mode 100644 index 000000000..392a45499 --- /dev/null +++ b/Examples/src/ex.csv @@ -0,0 +1,260 @@ +#ID,TIME,DOSE,OUT,WT,AFRICA,AGE,GENDER,HEIGHT +1,0,600,.,46.7,1,21,1,160 +1,24,600,.,46.7,1,21,1,160 +1,48,600,.,46.7,1,21,1,160 +1,72,600,.,46.7,1,21,1,160 +1,96,600,.,46.7,1,21,1,160 +1,120,.,10.44,46.7,1,21,1,160 +1,120,600,.,46.7,1,21,1,160 +1,121,.,12.89,46.7,1,21,1,160 +1,122,.,14.98,46.7,1,21,1,160 +1,125.99,.,16.69,46.7,1,21,1,160 +1,129,.,20.15,46.7,1,21,1,160 +1,132,.,14.97,46.7,1,21,1,160 +1,143.98,.,12.57,46.7,1,21,1,160 +2,0,600,.,66.5,1,30,1,174 +2,24,600,.,66.5,1,30,1,174 +2,48,600,.,66.5,1,30,1,174 +2,72,600,.,66.5,1,30,1,174 +2,96,600,.,66.5,1,30,1,174 +2,120,.,3.56,66.5,1,30,1,174 +2,120,600,.,66.5,1,30,1,174 +2,120.98,.,5.84,66.5,1,30,1,174 +2,121.98,.,6.54,66.5,1,30,1,174 +2,126,.,6.14,66.5,1,30,1,174 +2,129.02,.,6.56,66.5,1,30,1,174 +2,132.02,.,4.44,66.5,1,30,1,174 +2,144,.,3.76,66.5,1,30,1,174 +3,0,600,.,46.7,1,24,0,164 +3,24,600,.,46.7,1,24,0,164 +3,48,600,.,46.7,1,24,0,164 +3,72,600,.,46.7,1,24,0,164 +3,96,600,.,46.7,1,24,0,164 +3,120,600,.,46.7,1,24,0,164 +3,120.08,.,4.06,46.7,1,24,0,164 +3,121.07,.,3.24,46.7,1,24,0,164 +3,122.08,.,3.09,46.7,1,24,0,164 +3,126.08,.,7.98,46.7,1,24,0,164 +3,129.05,.,7.23,46.7,1,24,0,164 +3,132.1,.,4.71,46.7,1,24,0,164 +3,144.08,.,3.82,46.7,1,24,0,164 +4,0,600,.,50.8,1,25,1,165 +4,24,600,.,50.8,1,25,1,165 +4,48,600,.,50.8,1,25,1,165 +4,72,600,.,50.8,1,25,1,165 +4,96,600,.,50.8,1,25,1,165 +4,120,.,2.1,50.8,1,25,1,165 +4,120,600,.,50.8,1,25,1,165 +4,121,.,3.05,50.8,1,25,1,165 +4,122.02,.,5.21,50.8,1,25,1,165 +4,126,.,5.09,50.8,1,25,1,165 +4,129.03,.,4.24,50.8,1,25,1,165 +4,132,.,3.69,50.8,1,25,1,165 +4,144.02,.,1.96,50.8,1,25,1,165 +5,0,600,.,65.8,1,22,1,181 +5,24,600,.,65.8,1,22,1,181 +5,48,600,.,65.8,1,22,1,181 +5,72,600,.,65.8,1,22,1,181 +5,96,600,.,65.8,1,22,1,181 +5,120,.,2.93,65.8,1,22,1,181 +5,120,600,.,65.8,1,22,1,181 +5,121,.,2.64,65.8,1,22,1,181 +5,122,.,4.8,65.8,1,22,1,181 +5,126,.,3.7,65.8,1,22,1,181 +5,129.02,.,4.13,65.8,1,22,1,181 +5,132,.,2.81,65.8,1,22,1,181 +5,144,.,2.21,65.8,1,22,1,181 +6,0,600,.,65,1,23,1,177 +6,24,600,.,65,1,23,1,177 +6,48,600,.,65,1,23,1,177 +6,72,600,.,65,1,23,1,177 +6,96,600,.,65,1,23,1,177 +6,120,.,6.92,65,1,23,1,177 +6,120,600,.,65,1,23,1,177 +6,121,.,6.89,65,1,23,1,177 +6,121.98,.,6.64,65,1,23,1,177 +6,126,.,13.72,65,1,23,1,177 +6,129,.,12.69,65,1,23,1,177 +6,131.98,.,10.58,65,1,23,1,177 +6,144.98,.,6.62,65,1,23,1,177 +7,0,600,.,51.7,1,27,0,161 +7,24,600,.,51.7,1,27,0,161 +7,48,600,.,51.7,1,27,0,161 +7,72,600,.,51.7,1,27,0,161 +7,96,600,.,51.7,1,27,0,161 +7,120,.,5.41,51.7,1,27,0,161 +7,120,600,.,51.7,1,27,0,161 +7,121.03,.,4.46,51.7,1,27,0,161 +7,122.03,.,4.54,51.7,1,27,0,161 +7,126.02,.,12.19,51.7,1,27,0,161 +7,129.08,.,12.1,51.7,1,27,0,161 +7,132.03,.,8.61,51.7,1,27,0,161 +7,144.03,.,6.37,51.7,1,27,0,161 +8,0,600,.,51.2,1,22,1,163 +8,24,600,.,51.2,1,22,1,163 +8,48,600,.,51.2,1,22,1,163 +8,72,600,.,51.2,1,22,1,163 +8,96,600,.,51.2,1,22,1,163 +8,120,.,6.19,51.2,1,22,1,163 +8,120,600,.,51.2,1,22,1,163 +8,121.03,.,6.33,51.2,1,22,1,163 +8,122,.,6.24,51.2,1,22,1,163 +8,125.98,.,13.03,51.2,1,22,1,163 +8,128.98,.,11.86,51.2,1,22,1,163 +8,132,.,11.45,51.2,1,22,1,163 +8,143.98,.,7.83,51.2,1,22,1,163 +9,0,600,.,55,1,23,1,174 +9,24,600,.,55,1,23,1,174 +9,48,600,.,55,1,23,1,174 +9,72,600,.,55,1,23,1,174 +9,96,600,.,55,1,23,1,174 +9,120,.,2.85,55,1,23,1,174 +9,120,600,.,55,1,23,1,174 +9,120.97,.,3.7,55,1,23,1,174 +9,122,.,6.65,55,1,23,1,174 +9,125.98,.,6.81,55,1,23,1,174 +9,128.98,.,6.51,55,1,23,1,174 +9,132,.,7.48,55,1,23,1,174 +9,143.98,.,4.51,55,1,23,1,174 +10,0,600,.,52.1,1,32,1,163 +10,24,600,.,52.1,1,32,1,163 +10,48,600,.,52.1,1,32,1,163 +10,72,600,.,52.1,1,32,1,163 +10,96,600,.,52.1,1,32,1,163 +10,120,.,2.93,52.1,1,32,1,163 +10,120,600,.,52.1,1,32,1,163 +10,121,.,4.36,52.1,1,32,1,163 +10,122.02,.,7.79,52.1,1,32,1,163 +10,126,.,11.02,52.1,1,32,1,163 +10,129,.,8.86,52.1,1,32,1,163 +10,131.97,.,6.09,52.1,1,32,1,163 +10,144,.,4.15,52.1,1,32,1,163 +11,0,600,.,56.5,1,34,1,165 +11,24,600,.,56.5,1,34,1,165 +11,48,600,.,56.5,1,34,1,165 +11,72,600,.,56.5,1,34,1,165 +11,96,600,.,56.5,1,34,1,165 +11,120,.,2.09,56.5,1,34,1,165 +11,120,600,.,56.5,1,34,1,165 +11,121.03,.,2.68,56.5,1,34,1,165 +11,122,.,4.71,56.5,1,34,1,165 +11,125.98,.,7.71,56.5,1,34,1,165 +11,129,.,6.31,56.5,1,34,1,165 +11,132,.,5.82,56.5,1,34,1,165 +11,144.13,.,2.63,56.5,1,34,1,165 +12,0,600,.,47.9,1,54,0,160 +12,24,600,.,47.9,1,54,0,160 +12,48,600,.,47.9,1,54,0,160 +12,72,600,.,47.9,1,54,0,160 +12,96,600,.,47.9,1,54,0,160 +12,120,.,7.09,47.9,1,54,0,160 +12,120,600,.,47.9,1,54,0,160 +12,121.03,.,6.18,47.9,1,54,0,160 +12,122.13,.,8.66,47.9,1,54,0,160 +12,126,.,11.16,47.9,1,54,0,160 +12,129,.,9.51,47.9,1,54,0,160 +12,132,.,8.14,47.9,1,54,0,160 +12,144,.,7.89,47.9,1,54,0,160 +13,0,600,.,60.5,1,24,1,180 +13,24,600,.,60.5,1,24,1,180 +13,48,600,.,60.5,1,24,1,180 +13,72,600,.,60.5,1,24,1,180 +13,96,600,.,60.5,1,24,1,180 +13,120,.,6.62,60.5,1,24,1,180 +13,120,600,.,60.5,1,24,1,180 +13,121,.,3.18,60.5,1,24,1,180 +13,122,.,5.41,60.5,1,24,1,180 +13,126,.,10.18,60.5,1,24,1,180 +13,129.02,.,12.84,60.5,1,24,1,180 +13,132,.,12.35,60.5,1,24,1,180 +13,144,.,8.06,60.5,1,24,1,180 +14,0,600,.,59.2,1,26,1,174 +14,24,600,.,59.2,1,26,1,174 +14,48,600,.,59.2,1,26,1,174 +14,72,600,.,59.2,1,26,1,174 +14,96,600,.,59.2,1,26,1,174 +14,120,.,3.63,59.2,1,26,1,174 +14,120,600,.,59.2,1,26,1,174 +14,121,.,4.49,59.2,1,26,1,174 +14,122,.,5.5,59.2,1,26,1,174 +14,126,.,7.28,59.2,1,26,1,174 +14,129,.,5.27,59.2,1,26,1,174 +14,132,.,4.89,59.2,1,26,1,174 +14,144,.,2.68,59.2,1,26,1,174 +15,0,450,.,43,1,19,0,150 +15,24,450,.,43,1,19,0,150 +15,48,450,.,43,1,19,0,150 +15,72,450,.,43,1,19,0,150 +15,96,450,.,43,1,19,0,150 +15,120,.,5.53,43,1,19,0,150 +15,120,450,.,43,1,19,0,150 +15,121,.,4.81,43,1,19,0,150 +15,122,.,8.14,43,1,19,0,150 +15,126,.,9.96,43,1,19,0,150 +15,129,.,8.55,43,1,19,0,150 +15,132.05,.,7.54,43,1,19,0,150 +15,144.05,.,5.74,43,1,19,0,150 +16,0,600,.,64.4,1,25,1,173 +16,24,600,.,64.4,1,25,1,173 +16,48,600,.,64.4,1,25,1,173 +16,72,600,.,64.4,1,25,1,173 +16,96,600,.,64.4,1,25,1,173 +16,120,.,5.48,64.4,1,25,1,173 +16,120,600,.,64.4,1,25,1,173 +16,121,.,6.59,64.4,1,25,1,173 +16,122,.,8.91,64.4,1,25,1,173 +16,126,.,10.57,64.4,1,25,1,173 +16,129,.,9.52,64.4,1,25,1,173 +16,132,.,7.83,64.4,1,25,1,173 +16,143.97,.,4.96,64.4,1,25,1,173 +17,0,600,.,54.8,1,23,1,170 +17,24,600,.,54.8,1,23,1,170 +17,48,600,.,54.8,1,23,1,170 +17,72,600,.,54.8,1,23,1,170 +17,96,600,.,54.8,1,23,1,170 +17,120,.,2.11,54.8,1,23,1,170 +17,120,600,.,54.8,1,23,1,170 +17,121.02,.,1.86,54.8,1,23,1,170 +17,122.02,.,6.92,54.8,1,23,1,170 +17,126,.,9.11,54.8,1,23,1,170 +17,129,.,6.96,54.8,1,23,1,170 +17,132,.,5.64,54.8,1,23,1,170 +17,144.08,.,3.59,54.8,1,23,1,170 +18,0,450,.,44.3,1,20,0,164 +18,24,450,.,44.3,1,20,0,164 +18,48,450,.,44.3,1,20,0,164 +18,72,450,.,44.3,1,20,0,164 +18,96,450,.,44.3,1,20,0,164 +18,120,.,7.95,44.3,1,20,0,164 +18,120,450,.,44.3,1,20,0,164 +18,120.98,.,7.47,44.3,1,20,0,164 +18,121.98,.,8.67,44.3,1,20,0,164 +18,126,.,13.83,44.3,1,20,0,164 +18,129.17,.,14.01,44.3,1,20,0,164 +18,132.17,.,8.97,44.3,1,20,0,164 +18,143.97,.,8.4,44.3,1,20,0,164 +19,0,600,.,50,1,36,1,168 +19,24,600,.,50,1,36,1,168 +19,48,600,.,50,1,36,1,168 +19,72,600,.,50,1,36,1,168 +19,96,600,.,50,1,36,1,168 +19,120,.,5.42,50,1,36,1,168 +19,120,600,.,50,1,36,1,168 +19,121,.,7.08,50,1,36,1,168 +19,122,.,7.27,50,1,36,1,168 +19,125.98,.,20.07,50,1,36,1,168 +19,128.98,.,18.24,50,1,36,1,168 +19,132,.,15.36,50,1,36,1,168 +19,144,.,10.92,50,1,36,1,168 +20,0,600,.,59,1,31,1,170 +20,24,600,.,59,1,31,1,170 +20,48,600,.,59,1,31,1,170 +20,72,600,.,59,1,31,1,170 +20,96,600,.,59,1,31,1,170 +20,120,.,4.71,59,1,31,1,170 +20,120,600,.,59,1,31,1,170 +20,120.77,.,4.5,59,1,31,1,170 +20,121.75,.,3.35,59,1,31,1,170 +20,125.67,.,12.35,59,1,31,1,170 +20,128.67,.,11.56,59,1,31,1,170 +20,143.67,.,6.45,59,1,31,1,170 \ No newline at end of file diff --git a/Examples/src/ex_full.csv b/Examples/src/ex_full.csv new file mode 100644 index 000000000..0f3b8679c --- /dev/null +++ b/Examples/src/ex_full.csv @@ -0,0 +1,261 @@ +POPDATA DEC_11,,,,,,,,,,,,,,,,,, +#ID,EVID,TIME,DUR,DOSE,ADDL,II,INPUT,OUT,OUTEQ,C0,C1,C2,C3,WT,AFRICA,AGE,GENDER,HEIGHT +1,1,0,0,600,.,.,1,.,.,.,.,.,.,46.7,1,21,1,160 +1,1,24,0,600,.,.,1,.,.,.,.,.,.,46.7,1,21,1,160 +1,1,48,0,600,.,.,1,.,.,.,.,.,.,46.7,1,21,1,160 +1,1,72,0,600,.,.,1,.,.,.,.,.,.,46.7,1,21,1,160 +1,1,96,0,600,.,.,1,.,.,.,.,.,.,46.7,1,21,1,160 +1,0,120,.,.,.,.,.,10.44,1,0.02,0.0506,-0.0002,0,46.7,1,21,1,160 +1,1,120,0,600,.,.,1,.,.,.,.,.,.,46.7,1,21,1,160 +1,0,121,.,.,.,.,.,12.89,1,0.02,0.0506,-0.0002,0,46.7,1,21,1,160 +1,0,122,.,.,.,.,.,14.98,1,0.02,0.0506,-0.0002,0,46.7,1,21,1,160 +1,0,125.99,.,.,.,.,.,16.69,1,0.02,0.0506,-0.0002,0,46.7,1,21,1,160 +1,0,129,.,.,.,.,.,20.15,1,0.02,0.0506,-0.0002,0,46.7,1,21,1,160 +1,0,132,.,.,.,.,.,14.97,1,0.02,0.0506,-0.0002,0,46.7,1,21,1,160 +1,0,143.98,.,.,.,.,.,12.57,1,0.02,0.0506,-0.0002,0,46.7,1,21,1,160 +2,1,0,0,600,.,.,1,.,.,.,.,.,.,66.5,1,30,1,174 +2,1,24,0,600,.,.,1,.,.,.,.,.,.,66.5,1,30,1,174 +2,1,48,0,600,.,.,1,.,.,.,.,.,.,66.5,1,30,1,174 +2,1,72,0,600,.,.,1,.,.,.,.,.,.,66.5,1,30,1,174 +2,1,96,0,600,.,.,1,.,.,.,.,.,.,66.5,1,30,1,174 +2,0,120,.,.,.,.,.,3.56,1,0.02,0.0506,-0.0002,0,66.5,1,30,1,174 +2,1,120,0,600,.,.,1,.,.,.,.,.,.,66.5,1,30,1,174 +2,0,120.98,.,.,.,.,.,5.84,1,0.02,0.0506,-0.0002,0,66.5,1,30,1,174 +2,0,121.98,.,.,.,.,.,6.54,1,0.02,0.0506,-0.0002,0,66.5,1,30,1,174 +2,0,126,.,.,.,.,.,6.14,1,0.02,0.0506,-0.0002,0,66.5,1,30,1,174 +2,0,129.02,.,.,.,.,.,6.56,1,0.02,0.0506,-0.0002,0,66.5,1,30,1,174 +2,0,132.02,.,.,.,.,.,4.44,1,0.02,0.0506,-0.0002,0,66.5,1,30,1,174 +2,0,144,.,.,.,.,.,3.76,1,0.02,0.0506,-0.0002,0,66.5,1,30,1,174 +3,1,0,0,600,.,.,1,.,.,.,.,.,.,46.7,1,24,0,164 +3,1,24,0,600,.,.,1,.,.,.,.,.,.,46.7,1,24,0,164 +3,1,48,0,600,.,.,1,.,.,.,.,.,.,46.7,1,24,0,164 +3,1,72,0,600,.,.,1,.,.,.,.,.,.,46.7,1,24,0,164 +3,1,96,0,600,.,.,1,.,.,.,.,.,.,46.7,1,24,0,164 +3,1,120,0,600,.,.,1,.,.,.,.,.,.,46.7,1,24,0,164 +3,0,120.08,.,.,.,.,.,4.06,1,0.02,0.0506,-0.0002,0,46.7,1,24,0,164 +3,0,121.07,.,.,.,.,.,3.24,1,0.02,0.0506,-0.0002,0,46.7,1,24,0,164 +3,0,122.08,.,.,.,.,.,3.09,1,0.02,0.0506,-0.0002,0,46.7,1,24,0,164 +3,0,126.08,.,.,.,.,.,7.98,1,0.02,0.0506,-0.0002,0,46.7,1,24,0,164 +3,0,129.05,.,.,.,.,.,7.23,1,0.02,0.0506,-0.0002,0,46.7,1,24,0,164 +3,0,132.1,.,.,.,.,.,4.71,1,0.02,0.0506,-0.0002,0,46.7,1,24,0,164 +3,0,144.08,.,.,.,.,.,3.82,1,0.02,0.0506,-0.0002,0,46.7,1,24,0,164 +4,1,0,0,600,.,.,1,.,.,.,.,.,.,50.8,1,25,1,165 +4,1,24,0,600,.,.,1,.,.,.,.,.,.,50.8,1,25,1,165 +4,1,48,0,600,.,.,1,.,.,.,.,.,.,50.8,1,25,1,165 +4,1,72,0,600,.,.,1,.,.,.,.,.,.,50.8,1,25,1,165 +4,1,96,0,600,.,.,1,.,.,.,.,.,.,50.8,1,25,1,165 +4,0,120,.,.,.,.,.,2.1,1,0.02,0.0506,-0.0002,0,50.8,1,25,1,165 +4,1,120,0,600,.,.,1,.,.,.,.,.,.,50.8,1,25,1,165 +4,0,121,.,.,.,.,.,3.05,1,0.02,0.0506,-0.0002,0,50.8,1,25,1,165 +4,0,122.02,.,.,.,.,.,5.21,1,0.02,0.0506,-0.0002,0,50.8,1,25,1,165 +4,0,126,.,.,.,.,.,5.09,1,0.02,0.0506,-0.0002,0,50.8,1,25,1,165 +4,0,129.03,.,.,.,.,.,4.24,1,0.02,0.0506,-0.0002,0,50.8,1,25,1,165 +4,0,132,.,.,.,.,.,3.69,1,0.02,0.0506,-0.0002,0,50.8,1,25,1,165 +4,0,144.02,.,.,.,.,.,1.96,1,0.02,0.0506,-0.0002,0,50.8,1,25,1,165 +5,1,0,0,600,.,.,1,.,.,.,.,.,.,65.8,1,22,1,181 +5,1,24,0,600,.,.,1,.,.,.,.,.,.,65.8,1,22,1,181 +5,1,48,0,600,.,.,1,.,.,.,.,.,.,65.8,1,22,1,181 +5,1,72,0,600,.,.,1,.,.,.,.,.,.,65.8,1,22,1,181 +5,1,96,0,600,.,.,1,.,.,.,.,.,.,65.8,1,22,1,181 +5,0,120,.,.,.,.,.,2.93,1,0.02,0.0506,-0.0002,0,65.8,1,22,1,181 +5,1,120,0,600,.,.,1,.,.,.,.,.,.,65.8,1,22,1,181 +5,0,121,.,.,.,.,.,2.64,1,0.02,0.0506,-0.0002,0,65.8,1,22,1,181 +5,0,122,.,.,.,.,.,4.8,1,0.02,0.0506,-0.0002,0,65.8,1,22,1,181 +5,0,126,.,.,.,.,.,3.7,1,0.02,0.0506,-0.0002,0,65.8,1,22,1,181 +5,0,129.02,.,.,.,.,.,4.13,1,0.02,0.0506,-0.0002,0,65.8,1,22,1,181 +5,0,132,.,.,.,.,.,2.81,1,0.02,0.0506,-0.0002,0,65.8,1,22,1,181 +5,0,144,.,.,.,.,.,2.21,1,0.02,0.0506,-0.0002,0,65.8,1,22,1,181 +6,1,0,0,600,.,.,1,.,.,.,.,.,.,65,1,23,1,177 +6,1,24,0,600,.,.,1,.,.,.,.,.,.,65,1,23,1,177 +6,1,48,0,600,.,.,1,.,.,.,.,.,.,65,1,23,1,177 +6,1,72,0,600,.,.,1,.,.,.,.,.,.,65,1,23,1,177 +6,1,96,0,600,.,.,1,.,.,.,.,.,.,65,1,23,1,177 +6,0,120,.,.,.,.,.,6.92,1,0.02,0.0506,-0.0002,0,65,1,23,1,177 +6,1,120,0,600,.,.,1,.,.,.,.,.,.,65,1,23,1,177 +6,0,121,.,.,.,.,.,6.89,1,0.02,0.0506,-0.0002,0,65,1,23,1,177 +6,0,121.98,.,.,.,.,.,6.64,1,0.02,0.0506,-0.0002,0,65,1,23,1,177 +6,0,126,.,.,.,.,.,13.72,1,0.02,0.0506,-0.0002,0,65,1,23,1,177 +6,0,129,.,.,.,.,.,12.69,1,0.02,0.0506,-0.0002,0,65,1,23,1,177 +6,0,131.98,.,.,.,.,.,10.58,1,0.02,0.0506,-0.0002,0,65,1,23,1,177 +6,0,144.98,.,.,.,.,.,6.62,1,0.02,0.0506,-0.0002,0,65,1,23,1,177 +7,1,0,0,600,.,.,1,.,.,.,.,.,.,51.7,1,27,0,161 +7,1,24,0,600,.,.,1,.,.,.,.,.,.,51.7,1,27,0,161 +7,1,48,0,600,.,.,1,.,.,.,.,.,.,51.7,1,27,0,161 +7,1,72,0,600,.,.,1,.,.,.,.,.,.,51.7,1,27,0,161 +7,1,96,0,600,.,.,1,.,.,.,.,.,.,51.7,1,27,0,161 +7,0,120,.,.,.,.,.,5.41,1,0.02,0.0506,-0.0002,0,51.7,1,27,0,161 +7,1,120,0,600,.,.,1,.,.,.,.,.,.,51.7,1,27,0,161 +7,0,121.03,.,.,.,.,.,4.46,1,0.02,0.0506,-0.0002,0,51.7,1,27,0,161 +7,0,122.03,.,.,.,.,.,4.54,1,0.02,0.0506,-0.0002,0,51.7,1,27,0,161 +7,0,126.02,.,.,.,.,.,12.19,1,0.02,0.0506,-0.0002,0,51.7,1,27,0,161 +7,0,129.08,.,.,.,.,.,12.1,1,0.02,0.0506,-0.0002,0,51.7,1,27,0,161 +7,0,132.03,.,.,.,.,.,8.61,1,0.02,0.0506,-0.0002,0,51.7,1,27,0,161 +7,0,144.03,.,.,.,.,.,6.37,1,0.02,0.0506,-0.0002,0,51.7,1,27,0,161 +8,1,0,0,600,.,.,1,.,.,.,.,.,.,51.2,1,22,1,163 +8,1,24,0,600,.,.,1,.,.,.,.,.,.,51.2,1,22,1,163 +8,1,48,0,600,.,.,1,.,.,.,.,.,.,51.2,1,22,1,163 +8,1,72,0,600,.,.,1,.,.,.,.,.,.,51.2,1,22,1,163 +8,1,96,0,600,.,.,1,.,.,.,.,.,.,51.2,1,22,1,163 +8,0,120,.,.,.,.,.,6.19,1,0.02,0.0506,-0.0002,0,51.2,1,22,1,163 +8,1,120,0,600,.,.,1,.,.,.,.,.,.,51.2,1,22,1,163 +8,0,121.03,.,.,.,.,.,6.33,1,0.02,0.0506,-0.0002,0,51.2,1,22,1,163 +8,0,122,.,.,.,.,.,6.24,1,0.02,0.0506,-0.0002,0,51.2,1,22,1,163 +8,0,125.98,.,.,.,.,.,13.03,1,0.02,0.0506,-0.0002,0,51.2,1,22,1,163 +8,0,128.98,.,.,.,.,.,11.86,1,0.02,0.0506,-0.0002,0,51.2,1,22,1,163 +8,0,132,.,.,.,.,.,11.45,1,0.02,0.0506,-0.0002,0,51.2,1,22,1,163 +8,0,143.98,.,.,.,.,.,7.83,1,0.02,0.0506,-0.0002,0,51.2,1,22,1,163 +9,1,0,0,600,.,.,1,.,.,.,.,.,.,55,1,23,1,174 +9,1,24,0,600,.,.,1,.,.,.,.,.,.,55,1,23,1,174 +9,1,48,0,600,.,.,1,.,.,.,.,.,.,55,1,23,1,174 +9,1,72,0,600,.,.,1,.,.,.,.,.,.,55,1,23,1,174 +9,1,96,0,600,.,.,1,.,.,.,.,.,.,55,1,23,1,174 +9,0,120,.,.,.,.,.,2.85,1,0.02,0.0506,-0.0002,0,55,1,23,1,174 +9,1,120,0,600,.,.,1,.,.,.,.,.,.,55,1,23,1,174 +9,0,120.97,.,.,.,.,.,3.7,1,0.02,0.0506,-0.0002,0,55,1,23,1,174 +9,0,122,.,.,.,.,.,6.65,1,0.02,0.0506,-0.0002,0,55,1,23,1,174 +9,0,125.98,.,.,.,.,.,6.81,1,0.02,0.0506,-0.0002,0,55,1,23,1,174 +9,0,128.98,.,.,.,.,.,6.51,1,0.02,0.0506,-0.0002,0,55,1,23,1,174 +9,0,132,.,.,.,.,.,7.48,1,0.02,0.0506,-0.0002,0,55,1,23,1,174 +9,0,143.98,.,.,.,.,.,4.51,1,0.02,0.0506,-0.0002,0,55,1,23,1,174 +10,1,0,0,600,.,.,1,.,.,.,.,.,.,52.1,1,32,1,163 +10,1,24,0,600,.,.,1,.,.,.,.,.,.,52.1,1,32,1,163 +10,1,48,0,600,.,.,1,.,.,.,.,.,.,52.1,1,32,1,163 +10,1,72,0,600,.,.,1,.,.,.,.,.,.,52.1,1,32,1,163 +10,1,96,0,600,.,.,1,.,.,.,.,.,.,52.1,1,32,1,163 +10,0,120,.,.,.,.,.,2.93,1,0.02,0.0506,-0.0002,0,52.1,1,32,1,163 +10,1,120,0,600,.,.,1,.,.,.,.,.,.,52.1,1,32,1,163 +10,0,121,.,.,.,.,.,4.36,1,0.02,0.0506,-0.0002,0,52.1,1,32,1,163 +10,0,122.02,.,.,.,.,.,7.79,1,0.02,0.0506,-0.0002,0,52.1,1,32,1,163 +10,0,126,.,.,.,.,.,11.02,1,0.02,0.0506,-0.0002,0,52.1,1,32,1,163 +10,0,129,.,.,.,.,.,8.86,1,0.02,0.0506,-0.0002,0,52.1,1,32,1,163 +10,0,131.97,.,.,.,.,.,6.09,1,0.02,0.0506,-0.0002,0,52.1,1,32,1,163 +10,0,144,.,.,.,.,.,4.15,1,0.02,0.0506,-0.0002,0,52.1,1,32,1,163 +11,1,0,0,600,.,.,1,.,.,.,.,.,.,56.5,1,34,1,165 +11,1,24,0,600,.,.,1,.,.,.,.,.,.,56.5,1,34,1,165 +11,1,48,0,600,.,.,1,.,.,.,.,.,.,56.5,1,34,1,165 +11,1,72,0,600,.,.,1,.,.,.,.,.,.,56.5,1,34,1,165 +11,1,96,0,600,.,.,1,.,.,.,.,.,.,56.5,1,34,1,165 +11,0,120,.,.,.,.,.,2.09,1,0.02,0.0506,-0.0002,0,56.5,1,34,1,165 +11,1,120,0,600,.,.,1,.,.,.,.,.,.,56.5,1,34,1,165 +11,0,121.03,.,.,.,.,.,2.68,1,0.02,0.0506,-0.0002,0,56.5,1,34,1,165 +11,0,122,.,.,.,.,.,4.71,1,0.02,0.0506,-0.0002,0,56.5,1,34,1,165 +11,0,125.98,.,.,.,.,.,7.71,1,0.02,0.0506,-0.0002,0,56.5,1,34,1,165 +11,0,129,.,.,.,.,.,6.31,1,0.02,0.0506,-0.0002,0,56.5,1,34,1,165 +11,0,132,.,.,.,.,.,5.82,1,0.02,0.0506,-0.0002,0,56.5,1,34,1,165 +11,0,144.13,.,.,.,.,.,2.63,1,0.02,0.0506,-0.0002,0,56.5,1,34,1,165 +12,1,0,0,600,.,.,1,.,.,.,.,.,.,47.9,1,54,0,160 +12,1,24,0,600,.,.,1,.,.,.,.,.,.,47.9,1,54,0,160 +12,1,48,0,600,.,.,1,.,.,.,.,.,.,47.9,1,54,0,160 +12,1,72,0,600,.,.,1,.,.,.,.,.,.,47.9,1,54,0,160 +12,1,96,0,600,.,.,1,.,.,.,.,.,.,47.9,1,54,0,160 +12,0,120,.,.,.,.,.,7.09,1,0.02,0.0506,-0.0002,0,47.9,1,54,0,160 +12,1,120,0,600,.,.,1,.,.,.,.,.,.,47.9,1,54,0,160 +12,0,121.03,.,.,.,.,.,6.18,1,0.02,0.0506,-0.0002,0,47.9,1,54,0,160 +12,0,122.13,.,.,.,.,.,8.66,1,0.02,0.0506,-0.0002,0,47.9,1,54,0,160 +12,0,126,.,.,.,.,.,11.16,1,0.02,0.0506,-0.0002,0,47.9,1,54,0,160 +12,0,129,.,.,.,.,.,9.51,1,0.02,0.0506,-0.0002,0,47.9,1,54,0,160 +12,0,132,.,.,.,.,.,8.14,1,0.02,0.0506,-0.0002,0,47.9,1,54,0,160 +12,0,144,.,.,.,.,.,7.89,1,0.02,0.0506,-0.0002,0,47.9,1,54,0,160 +13,1,0,0,600,.,.,1,.,.,.,.,.,.,60.5,1,24,1,180 +13,1,24,0,600,.,.,1,.,.,.,.,.,.,60.5,1,24,1,180 +13,1,48,0,600,.,.,1,.,.,.,.,.,.,60.5,1,24,1,180 +13,1,72,0,600,.,.,1,.,.,.,.,.,.,60.5,1,24,1,180 +13,1,96,0,600,.,.,1,.,.,.,.,.,.,60.5,1,24,1,180 +13,0,120,.,.,.,.,.,6.62,1,0.02,0.0506,-0.0002,0,60.5,1,24,1,180 +13,1,120,0,600,.,.,1,.,.,.,.,.,.,60.5,1,24,1,180 +13,0,121,.,.,.,.,.,3.18,1,0.02,0.0506,-0.0002,0,60.5,1,24,1,180 +13,0,122,.,.,.,.,.,5.41,1,0.02,0.0506,-0.0002,0,60.5,1,24,1,180 +13,0,126,.,.,.,.,.,10.18,1,0.02,0.0506,-0.0002,0,60.5,1,24,1,180 +13,0,129.02,.,.,.,.,.,12.84,1,0.02,0.0506,-0.0002,0,60.5,1,24,1,180 +13,0,132,.,.,.,.,.,12.35,1,0.02,0.0506,-0.0002,0,60.5,1,24,1,180 +13,0,144,.,.,.,.,.,8.06,1,0.02,0.0506,-0.0002,0,60.5,1,24,1,180 +14,1,0,0,600,.,.,1,.,.,.,.,.,.,59.2,1,26,1,174 +14,1,24,0,600,.,.,1,.,.,.,.,.,.,59.2,1,26,1,174 +14,1,48,0,600,.,.,1,.,.,.,.,.,.,59.2,1,26,1,174 +14,1,72,0,600,.,.,1,.,.,.,.,.,.,59.2,1,26,1,174 +14,1,96,0,600,.,.,1,.,.,.,.,.,.,59.2,1,26,1,174 +14,0,120,.,.,.,.,.,3.63,1,0.02,0.0506,-0.0002,0,59.2,1,26,1,174 +14,1,120,0,600,.,.,1,.,.,.,.,.,.,59.2,1,26,1,174 +14,0,121,.,.,.,.,.,4.49,1,0.02,0.0506,-0.0002,0,59.2,1,26,1,174 +14,0,122,.,.,.,.,.,5.5,1,0.02,0.0506,-0.0002,0,59.2,1,26,1,174 +14,0,126,.,.,.,.,.,7.28,1,0.02,0.0506,-0.0002,0,59.2,1,26,1,174 +14,0,129,.,.,.,.,.,5.27,1,0.02,0.0506,-0.0002,0,59.2,1,26,1,174 +14,0,132,.,.,.,.,.,4.89,1,0.02,0.0506,-0.0002,0,59.2,1,26,1,174 +14,0,144,.,.,.,.,.,2.68,1,0.02,0.0506,-0.0002,0,59.2,1,26,1,174 +15,1,0,0,450,.,.,1,.,.,.,.,.,.,43,1,19,0,150 +15,1,24,0,450,.,.,1,.,.,.,.,.,.,43,1,19,0,150 +15,1,48,0,450,.,.,1,.,.,.,.,.,.,43,1,19,0,150 +15,1,72,0,450,.,.,1,.,.,.,.,.,.,43,1,19,0,150 +15,1,96,0,450,.,.,1,.,.,.,.,.,.,43,1,19,0,150 +15,0,120,.,.,.,.,.,5.53,1,0.02,0.0506,-0.0002,0,43,1,19,0,150 +15,1,120,0,450,.,.,1,.,.,.,.,.,.,43,1,19,0,150 +15,0,121,.,.,.,.,.,4.81,1,0.02,0.0506,-0.0002,0,43,1,19,0,150 +15,0,122,.,.,.,.,.,8.14,1,0.02,0.0506,-0.0002,0,43,1,19,0,150 +15,0,126,.,.,.,.,.,9.96,1,0.02,0.0506,-0.0002,0,43,1,19,0,150 +15,0,129,.,.,.,.,.,8.55,1,0.02,0.0506,-0.0002,0,43,1,19,0,150 +15,0,132.05,.,.,.,.,.,7.54,1,0.02,0.0506,-0.0002,0,43,1,19,0,150 +15,0,144.05,.,.,.,.,.,5.74,1,0.02,0.0506,-0.0002,0,43,1,19,0,150 +16,1,0,0,600,.,.,1,.,.,.,.,.,.,64.4,1,25,1,173 +16,1,24,0,600,.,.,1,.,.,.,.,.,.,64.4,1,25,1,173 +16,1,48,0,600,.,.,1,.,.,.,.,.,.,64.4,1,25,1,173 +16,1,72,0,600,.,.,1,.,.,.,.,.,.,64.4,1,25,1,173 +16,1,96,0,600,.,.,1,.,.,.,.,.,.,64.4,1,25,1,173 +16,0,120,.,.,.,.,.,5.48,1,0.02,0.0506,-0.0002,0,64.4,1,25,1,173 +16,1,120,0,600,.,.,1,.,.,.,.,.,.,64.4,1,25,1,173 +16,0,121,.,.,.,.,.,6.59,1,0.02,0.0506,-0.0002,0,64.4,1,25,1,173 +16,0,122,.,.,.,.,.,8.91,1,0.02,0.0506,-0.0002,0,64.4,1,25,1,173 +16,0,126,.,.,.,.,.,10.57,1,0.02,0.0506,-0.0002,0,64.4,1,25,1,173 +16,0,129,.,.,.,.,.,9.52,1,0.02,0.0506,-0.0002,0,64.4,1,25,1,173 +16,0,132,.,.,.,.,.,7.83,1,0.02,0.0506,-0.0002,0,64.4,1,25,1,173 +16,0,143.97,.,.,.,.,.,4.96,1,0.02,0.0506,-0.0002,0,64.4,1,25,1,173 +17,1,0,0,600,.,.,1,.,.,.,.,.,.,54.8,1,23,1,170 +17,1,24,0,600,.,.,1,.,.,.,.,.,.,54.8,1,23,1,170 +17,1,48,0,600,.,.,1,.,.,.,.,.,.,54.8,1,23,1,170 +17,1,72,0,600,.,.,1,.,.,.,.,.,.,54.8,1,23,1,170 +17,1,96,0,600,.,.,1,.,.,.,.,.,.,54.8,1,23,1,170 +17,0,120,.,.,.,.,.,2.11,1,0.02,0.0506,-0.0002,0,54.8,1,23,1,170 +17,1,120,0,600,.,.,1,.,.,.,.,.,.,54.8,1,23,1,170 +17,0,121.02,.,.,.,.,.,1.86,1,0.02,0.0506,-0.0002,0,54.8,1,23,1,170 +17,0,122.02,.,.,.,.,.,6.92,1,0.02,0.0506,-0.0002,0,54.8,1,23,1,170 +17,0,126,.,.,.,.,.,9.11,1,0.02,0.0506,-0.0002,0,54.8,1,23,1,170 +17,0,129,.,.,.,.,.,6.96,1,0.02,0.0506,-0.0002,0,54.8,1,23,1,170 +17,0,132,.,.,.,.,.,5.64,1,0.02,0.0506,-0.0002,0,54.8,1,23,1,170 +17,0,144.08,.,.,.,.,.,3.59,1,0.02,0.0506,-0.0002,0,54.8,1,23,1,170 +18,1,0,0,450,.,.,1,.,.,.,.,.,.,44.3,1,20,0,164 +18,1,24,0,450,.,.,1,.,.,.,.,.,.,44.3,1,20,0,164 +18,1,48,0,450,.,.,1,.,.,.,.,.,.,44.3,1,20,0,164 +18,1,72,0,450,.,.,1,.,.,.,.,.,.,44.3,1,20,0,164 +18,1,96,0,450,.,.,1,.,.,.,.,.,.,44.3,1,20,0,164 +18,0,120,.,.,.,.,.,7.95,1,0.02,0.0506,-0.0002,0,44.3,1,20,0,164 +18,1,120,0,450,.,.,1,.,.,.,.,.,.,44.3,1,20,0,164 +18,0,120.98,.,.,.,.,.,7.47,1,0.02,0.0506,-0.0002,0,44.3,1,20,0,164 +18,0,121.98,.,.,.,.,.,8.67,1,0.02,0.0506,-0.0002,0,44.3,1,20,0,164 +18,0,126,.,.,.,.,.,13.83,1,0.02,0.0506,-0.0002,0,44.3,1,20,0,164 +18,0,129.17,.,.,.,.,.,14.01,1,0.02,0.0506,-0.0002,0,44.3,1,20,0,164 +18,0,132.17,.,.,.,.,.,8.97,1,0.02,0.0506,-0.0002,0,44.3,1,20,0,164 +18,0,143.97,.,.,.,.,.,8.4,1,0.02,0.0506,-0.0002,0,44.3,1,20,0,164 +19,1,0,0,600,.,.,1,.,.,.,.,.,.,50,1,36,1,168 +19,1,24,0,600,.,.,1,.,.,.,.,.,.,50,1,36,1,168 +19,1,48,0,600,.,.,1,.,.,.,.,.,.,50,1,36,1,168 +19,1,72,0,600,.,.,1,.,.,.,.,.,.,50,1,36,1,168 +19,1,96,0,600,.,.,1,.,.,.,.,.,.,50,1,36,1,168 +19,0,120,.,.,.,.,.,5.42,1,0.02,0.0506,-0.0002,0,50,1,36,1,168 +19,1,120,0,600,.,.,1,.,.,.,.,.,.,50,1,36,1,168 +19,0,121,.,.,.,.,.,7.08,1,0.02,0.0506,-0.0002,0,50,1,36,1,168 +19,0,122,.,.,.,.,.,7.27,1,0.02,0.0506,-0.0002,0,50,1,36,1,168 +19,0,125.98,.,.,.,.,.,20.07,1,0.02,0.0506,-0.0002,0,50,1,36,1,168 +19,0,128.98,.,.,.,.,.,18.24,1,0.02,0.0506,-0.0002,0,50,1,36,1,168 +19,0,132,.,.,.,.,.,15.36,1,0.02,0.0506,-0.0002,0,50,1,36,1,168 +19,0,144,.,.,.,.,.,10.92,1,0.02,0.0506,-0.0002,0,50,1,36,1,168 +20,1,0,0,600,.,.,1,.,.,.,.,.,.,59,1,31,1,170 +20,1,24,0,600,.,.,1,.,.,.,.,.,.,59,1,31,1,170 +20,1,48,0,600,.,.,1,.,.,.,.,.,.,59,1,31,1,170 +20,1,72,0,600,.,.,1,.,.,.,.,.,.,59,1,31,1,170 +20,1,96,0,600,.,.,1,.,.,.,.,.,.,59,1,31,1,170 +20,0,120,.,.,.,.,.,4.71,1,0.02,0.0506,-0.0002,0,59,1,31,1,170 +20,1,120,0,600,.,.,1,.,.,.,.,.,.,59,1,31,1,170 +20,0,120.77,.,.,.,.,.,4.5,1,0.02,0.0506,-0.0002,0,59,1,31,1,170 +20,0,121.75,.,.,.,.,.,3.35,1,0.02,0.0506,-0.0002,0,59,1,31,1,170 +20,0,125.67,.,.,.,.,.,12.35,1,0.02,0.0506,-0.0002,0,59,1,31,1,170 +20,0,128.67,.,.,.,.,.,11.56,1,0.02,0.0506,-0.0002,0,59,1,31,1,170 +20,0,143.67,.,.,.,.,.,6.45,1,0.02,0.0506,-0.0002,0,59,1,31,1,170 \ No newline at end of file diff --git a/Examples/src/ptaex1.csv b/Examples/src/ptaex1.csv new file mode 100644 index 000000000..93a207c9a --- /dev/null +++ b/Examples/src/ptaex1.csv @@ -0,0 +1,10 @@ +POPDATA DEC_11,,,,,,,,,,,,,,,,,, +#ID,EVID,TIME,DUR,DOSE,ADDL,II,INPUT,OUT,OUTEQ,C0,C1,C2,C3,WT,AFRICA,AGE,GENDER,HEIGHT +1,1,0,0,600,5,24,1,.,.,.,.,.,.,46.7,1,21,1,160 +1,0,144,.,.,.,.,.,-1,1,0.02,0.0506,-2.00E-04,0,46.7,1,21,1,160 +2,1,0,0,1200,5,24,1,.,.,.,.,.,.,46.7,1,21,1,160 +2,0,144,.,.,.,.,.,-1,1,0.02,0.0506,-2.00E-04,0,46.7,1,21,1,160 +3,1,0,0,300,11,12,1,.,.,.,.,.,.,46.7,1,21,1,160 +3,0,144,.,.,.,.,.,-1,1,0.02,0.0506,-2.00E-04,0,46.7,1,21,1,160 +4,1,0,0,600,11,12,1,.,.,.,.,.,.,46.7,1,21,1,160 +4,0,144,.,.,.,.,.,-1,1,0.02,0.0506,-2.00E-04,0,46.7,1,21,1,160 \ No newline at end of file diff --git a/Examples/src/simTemp.csv b/Examples/src/simTemp.csv new file mode 100644 index 000000000..a4d379b8a --- /dev/null +++ b/Examples/src/simTemp.csv @@ -0,0 +1,10 @@ +POPDATA DEC_11,,,,,,,,,,,,,,,,,, +#ID,EVID,TIME,DUR,DOSE,ADDL,II,INPUT,OUT,OUTEQ,C0,C1,C2,C3,WT,AFRICA,AGE,GENDER,HEIGHT +1,1,0,0,500,5,24,1,.,.,.,.,.,.,46.7,1,21,1,160 +1,0,144,.,.,.,.,.,-1,1,.,.,.,.,.,.,.,.,. +2,1,0,0,1000,5,24,1,.,.,.,.,.,.,46.7,1,21,1,160 +2,0,144,.,.,.,.,.,-1,1,.,.,.,.,.,.,.,.,. +3,1,0,0,250,11,12,1,.,.,.,.,.,.,46.7,1,21,1,160 +3,0,144,.,.,.,.,.,-1,1,.,.,.,.,.,.,.,.,. +4,1,0,0,500,11,12,1,.,.,.,.,.,.,46.7,1,21,1,160 +4,0,144,.,.,.,.,.,-1,1,.,.,.,.,.,.,.,.,. \ No newline at end of file diff --git a/NAMESPACE b/NAMESPACE index 07e9b8d43..63feb0401 100755 --- a/NAMESPACE +++ b/NAMESPACE @@ -47,6 +47,8 @@ export(NPreport) export(NPrun) export(PMFortranConfig) export(PM_batch) +export(PM_bestdose) +export(PM_bestdose_problem) export(PM_build) export(PM_compare) export(PM_cov) @@ -94,6 +96,7 @@ export(add_renal) export(add_shapes) export(add_smooth) export(additive) +export(bestdose) export(build_model) export(build_plot) export(cli_ask) diff --git a/R/PM_bestdose.R b/R/PM_bestdose.R new file mode 100644 index 000000000..684c04e02 --- /dev/null +++ b/R/PM_bestdose.R @@ -0,0 +1,381 @@ +bestdose_parse_prior <- function(prior) { + if (inherits(prior, "PM_result")) { + theta_path <- file.path(prior$rundir, "outputs", "theta.csv") + if (!file.exists(theta_path)) { + cli::cli_abort("theta.csv not found in PM_result outputs") + } + theta_path + } else if (inherits(prior, "PM_final")) { + temp_path <- tempfile(fileext = ".csv") + bestdose_write_prior_csv(prior, temp_path) + temp_path + } else if (is.character(prior)) { + if (!file.exists(prior)) { + cli::cli_abort("Prior file not found: {prior}") + } + prior + } else { + cli::cli_abort("prior must be PM_result, PM_final, or path to theta.csv") + } +} + +bestdose_write_prior_csv <- function(prior, path) { + df <- as.data.frame(prior$popPoints) + df$prob <- prior$popProb + write.csv(df, path, row.names = FALSE, quote = FALSE) +} + +bestdose_parse_model <- function(model) { + if (inherits(model, "PM_model")) { + compiled_path <- model$binary_path + if (is.null(compiled_path) || !file.exists(compiled_path)) { + cli::cli_abort("Model must be compiled first. Use model$compile()") + } + + kind <- if (!is.null(model$model_list$analytical) && model$model_list$analytical) { + "analytical" + } else { + "ode" + } + + list(path = compiled_path, kind = kind, model = model) + } else if (is.character(model)) { + if (!file.exists(model)) { + cli::cli_abort("Model file not found: {model}") + } + kind <- if (grepl("analytical", model, ignore.case = TRUE)) { + "analytical" + } else { + "ode" + } + list(path = model, kind = kind, model = NULL) + } else { + cli::cli_abort("model must be PM_model or path to compiled model") + } +} + +bestdose_parse_data <- function(data) { + if (inherits(data, "PM_data")) { + temp_path <- tempfile(fileext = ".csv") + write.csv(data$standard_data, temp_path, row.names = FALSE, quote = FALSE) + temp_path + } else if (is.character(data)) { + if (!file.exists(data)) { + cli::cli_abort("Data file not found: {data}") + } + data + } else { + cli::cli_abort("data must be PM_data or path to CSV file") + } +} + +bestdose_default_settings <- function(prior, model) { + if (inherits(prior, "PM_result")) { + return(prior$settings) + } + + param_ranges <- lapply(model$model_list$pri, function(x) { + c(x$min, x$max) + }) + names(param_ranges) <- tolower(names(param_ranges)) + + list( + algorithm = "NPAG", + ranges = param_ranges, + error_models = list( + list( + initial = 0.0, + type = "additive", + coeff = c(0.0, 0.2, 0.0, 0.0) + ) + ), + max_cycles = 500, + points = 2028, + seed = 22, + prior = "prior.csv", + idelta = 0.25, + tad = 0.0 + ) +} + +#' @title +#' Object to contain BestDose optimization results +#' +#' @description +#' `r lifecycle::badge("experimental")` +#' +#' This object is created after a successful BestDose optimization run. +#' BestDose finds optimal dosing regimens to achieve target drug concentrations +#' or AUC values using Bayesian optimization. +#' +#' @export +PM_bestdose <- R6::R6Class( + "PM_bestdose", + public = list( + result = NULL, + problem = NULL, + bias_weight = NULL, + initialize = function(prior = NULL, + model = NULL, + past_data = NULL, + target = NULL, + dose_range = list(min = 0, max = 1000), + bias_weight = 0.5, + target_type = "concentration", + time_offset = NULL, + settings = NULL, + result = NULL, + problem_obj = NULL, + bias_override = NULL) { + if (!is.null(result)) { + private$.set_result(result, problem_obj, bias_override) + return(invisible(self)) + } + + if (is.null(target)) { + cli::cli_abort("target must be supplied when computing a new BestDose result") + } + + problem <- PM_bestdose_problem$new( + prior = prior, + model = model, + past_data = past_data, + target = target, + dose_range = dose_range, + bias_weight = bias_weight, + target_type = target_type, + time_offset = time_offset, + settings = settings + ) + + raw <- problem$optimize_raw(bias_weight = bias_weight) + private$.set_result(raw$result, problem, raw$bias_weight) + invisible(self) + }, + + #' @description + #' Print summary of BestDose results + print = function() { + cat("BestDose Optimization Results\n") + cat("==============================\n\n") + cat(sprintf("Optimal doses: [%.2f, %.2f] mg\n", self$get_doses()[1], self$get_doses()[2])) + cat(sprintf("Objective function: %.10f\n", self$get_objf())) + cat(sprintf("ln(Objective): %.4f\n", log(self$get_objf()))) + cat(sprintf("Method: %s\n", self$get_method())) + cat(sprintf("Status: %s\n", self$get_status())) + if (!is.null(self$bias_weight)) { + cat(sprintf("Bias weight (lambda): %.2f\n", self$bias_weight)) + } + cat(sprintf("\nNumber of predictions: %d\n", nrow(self$result$predictions))) + if (!is.null(self$result$auc_predictions)) { + cat(sprintf("Number of AUC predictions: %d\n", nrow(self$result$auc_predictions))) + } + invisible(self) + }, + + #' @description + #' Get optimal dose values + #' @return Numeric vector of optimal doses + get_doses = function() { + self$result$doses + }, + + #' @description + #' Get concentration-time predictions + #' @return Data frame with predictions + get_predictions = function() { + self$result$predictions + }, + + #' @description + #' Get AUC predictions (if available) + #' @return Data frame with AUC predictions or NULL + get_auc_predictions = function() { + self$result$auc_predictions + }, + + #' @description + #' Get objective function value + #' @return Numeric objective function value + get_objf = function() { + self$result$objf + }, + + #' @description + #' Get optimization status + #' @return Character string with status + get_status = function() { + self$result$status + }, + + #' @description + #' Get optimization method used + #' @return Character string: "posterior" or "uniform" + get_method = function() { + self$result$method + }, + + #' @description + #' Save results to RDS file + #' @param filename Path to save file. Default: "bestdose_result.rds" + save = function(filename = "bestdose_result.rds") { + saveRDS(self, filename) + cli::cli_alert_success("Results saved to {filename}") + invisible(self) + } + ), + private = list( + .set_result = function(result, problem, bias_weight) { + self$result <- result + self$problem <- problem + self$bias_weight <- bias_weight + } + ) +) + +#' @title +#' Prepare a reusable BestDose optimization problem +#' +#' @description +#' `r lifecycle::badge("experimental")` +#' +#' Use `PM_bestdose_problem` to mirror the Rust workflow: compute posterior +#' support points once, inspect them in R, and solve for multiple bias weights +#' without repeating the expensive initialization step. +#' +#' @export +PM_bestdose_problem <- R6::R6Class( + "PM_bestdose_problem", + public = list( + handle = NULL, + theta = NULL, + theta_dim = NULL, + param_names = NULL, + posterior_weights = NULL, + population_weights = NULL, + bias_weight = NULL, + target_type = NULL, + dose_range = NULL, + model_info = NULL, + settings = NULL, + initialize = function(prior, + model, + past_data = NULL, + target, + dose_range = list(min = 0, max = 1000), + bias_weight = 0.5, + target_type = "concentration", + time_offset = NULL, + settings = NULL) { + if (!target_type %in% c("concentration", "auc_from_zero", "auc_from_last_dose")) { + cli::cli_abort("target_type must be one of: concentration, auc_from_zero, auc_from_last_dose") + } + + if (bias_weight < 0 || bias_weight > 1) { + cli::cli_abort("bias_weight must be between 0 and 1") + } + + if (is.null(dose_range$min) || is.null(dose_range$max)) { + cli::cli_abort("dose_range must have both 'min' and 'max' elements") + } + + if (dose_range$min >= dose_range$max) { + cli::cli_abort("dose_range$min must be less than dose_range$max") + } + + prior_path <- bestdose_parse_prior(prior) + model_info <- bestdose_parse_model(model) + past_data_path <- if (!is.null(past_data)) bestdose_parse_data(past_data) else NULL + target_data_path <- bestdose_parse_data(target) + + if (is.null(settings)) { + model_for_settings <- if (!is.null(model_info$model)) model_info$model else model + settings <- bestdose_default_settings(prior, model_for_settings) + } + + prep <- bestdose_prepare( + model_path = model_info$path, + prior_path = prior_path, + past_data_path = past_data_path, + target_data_path = target_data_path, + time_offset = time_offset, + dose_min = dose_range$min, + dose_max = dose_range$max, + bias_weight = bias_weight, + target_type = target_type, + params = settings, + kind = model_info$kind + ) + + if (is.character(prep)) { + cli::cli_abort(prep) + } + + dim <- as.integer(prep$theta_dim) + theta_matrix <- matrix(prep$theta_values, nrow = dim[1], ncol = dim[2]) + colnames(theta_matrix) <- prep$param_names + + self$handle <- prep$handle + self$theta <- theta_matrix + self$theta_dim <- dim + self$param_names <- prep$param_names + self$posterior_weights <- prep$posterior_weights + self$population_weights <- prep$population_weights + self$bias_weight <- prep$bias_weight + self$target_type <- prep$target_type + self$dose_range <- dose_range + self$model_info <- model_info + self$settings <- settings + + cli::cli_alert_success("BestDose problem prepared with %d support points", dim[1]) + }, + finalize = function() { + self$handle <- NULL + }, + + #' @description + #' Run optimization and return raw list (doses, objf, predictions) + optimize_raw = function(bias_weight = NULL) { + private$.run_optimize(bias_weight) + }, + + #' @description + #' Run optimization and return a `PM_bestdose` result object + optimize = function(bias_weight = NULL) { + raw <- self$optimize_raw(bias_weight) + PM_bestdose$new( + result = raw$result, + problem_obj = self, + bias_override = raw$bias_weight + ) + } + ), + private = list( + .run_optimize = function(bias_weight) { + if (is.null(self$handle)) { + cli::cli_abort("BestDose problem handle has been released") + } + + bw <- if (is.null(bias_weight)) self$bias_weight else bias_weight + + if (bw < 0 || bw > 1) { + cli::cli_abort("bias_weight must be between 0 and 1") + } + + res <- bestdose_optimize(self$handle, bw) + if (is.character(res)) { + cli::cli_abort(res) + } + + list(result = res, bias_weight = bw) + } + ) +) + +#' @export +PM_bestdose$load <- function(filename = "bestdose_result.rds") { + if (!file.exists(filename)) { + cli::cli_abort("File not found: {filename}") + } + readRDS(filename) +} diff --git a/R/PM_result.R b/R/PM_result.R index e0e324523..a2c4aaa64 100755 --- a/R/PM_result.R +++ b/R/PM_result.R @@ -54,7 +54,7 @@ PM_result <- R6::R6Class( #' Use the `$save` method on the augmented `PM_result` object to save it with the #' new optimal sampling results. opt_samp = NULL, - + #' @description #' Create new object populated with data from previous run #' @details @@ -73,9 +73,9 @@ PM_result <- R6::R6Class( if (!inherits(out[[x]], "R6")) { # older save cli::cli_abort(c("x" = "The object was saved in an older format. Please re-run the analysis.")) } else { - if(x == "model"){ - args <- list(x = out[[x]], compile = FALSE) - } else { + if (x == "model") { + args <- list(x = out[[x]], compile = FALSE) + } else { args <- list(out[[x]], path = path, quiet = TRUE) } self[[x]] <- do.call(get(paste0("PM_", x))$new, args = args) # was saved in R6 format, but remake to update if needed @@ -83,24 +83,24 @@ PM_result <- R6::R6Class( } } ) - + # these are diagnostics, not R6 self$errfile <- out$errfile self$success <- out$success - + # add the pop/post data to data if (is.null(self$data$pop) | is.null(self$data$post)) { self$data <- PM_data$new(self$data$data, quiet = TRUE) self$data$pop <- self$pop$data self$data$post <- self$post$data } - + return(self) }, #' @description #' Fit the model to the data #' #' @details - #' This method is used to fit the model in the [PM_result] object to data. + #' This method is used to fit the model in the [PM_result] object to data. #' It calls the `$fit` method of the model stored in the `model` field. #' @param data Optional data to fit. If not provided, the data stored in the #' `data` field of the [PM_result] object will be used. This can be useful to @@ -109,15 +109,15 @@ PM_result <- R6::R6Class( #' @param ... Additional arguments passed to the model's `$fit` method. #' @return Returns an invisible [PM_result]. #' @export - #' - fit = function(data, ...){ + #' + fit = function(data, ...) { if (missing(data)) { data <- self$data } res <- self$model$fit(data = data, ...) return(invisible(res)) }, - + #' @description #' Plot generic function based on type #' @param type Type of plot based on class of object @@ -129,7 +129,7 @@ PM_result <- R6::R6Class( self[[type]]$plot(...) } }, - + #' @description #' Summary generic function based on type #' @param type Type of summary based on class of object @@ -141,7 +141,7 @@ PM_result <- R6::R6Class( self[[type]]$summary(...) } }, - + #' @description #' AUC generic function based on type #' @param type Type of AUC based on class of object @@ -152,7 +152,7 @@ PM_result <- R6::R6Class( } self[[type]]$auc(...) }, - + #' @description #' Perform non-compartmental analysis #' @details @@ -180,22 +180,88 @@ PM_result <- R6::R6Class( if (!"poppar" %in% names(dots)) { dots$poppar <- self$final } - + if (!"data" %in% names(dots)) { dots$data <- self$data } - + if (!"model" %in% names(dots)) { dots$model <- self$model } - + # store copy of the final object bk_final <- self$final$clone() sim <- do.call(PM_sim$new, dots) self$final <- bk_final return(sim) }, - + + #' @description + #' Run BestDose optimization using this result as the prior + #' @details + #' BestDose finds optimal dosing regimens to achieve target drug concentrations + #' or AUC values. The algorithm uses Bayesian posterior estimation combined with + #' dual optimization to balance patient-specific adaptation and population-level + #' robustness. By default, uses the `$final`, `$model`, and `$data` objects from + #' this result. Most commonly, you will supply a different `target` data object. + #' @param target PM_data object or path to CSV with target doses/observations. + #' Required. This defines the dosing template and target values. Set dose amounts + #' to 0 for doses to be optimized. + #' @param past_data Optional: PM_data object or path to CSV with patient history. + #' If NULL (default), uses the full dataset from this result. + #' @param dose_range Named list with min and max allowable doses. + #' Default: list(min = 0, max = 1000) + #' @param bias_weight Numeric [0,1] controlling personalization level. + #' 0 = fully personalized, 1 = population-based. Default: 0.5 (balanced) + #' @param target_type One of "concentration" (default), "auc_from_zero", or + #' "auc_from_last_dose" + #' @param time_offset Optional: time offset for past/future concatenation + #' @param ... Additional arguments passed to [PM_bestdose] + #' @return A [PM_bestdose] object containing optimal doses and predictions + #' @examples + #' \dontrun{ + #' # Load NPAG result + #' result <- PM_load(1) + #' + #' # Create target data + #' target <- PM_data$new("target.csv") + #' + #' # Run BestDose optimization + #' bd <- result$bestdose( + #' target = target, + #' dose_range = list(min = 50, max = 500), + #' bias_weight = 0.5 + #' ) + #' + #' # View results + #' print(bd) + #' bd$get_doses() + #' } + bestdose = function(target, + past_data = NULL, + dose_range = list(min = 0, max = 1000), + bias_weight = 0.5, + target_type = "concentration", + time_offset = NULL, + ...) { + # Use this result's data as past_data if not specified + if (is.null(past_data)) { + past_data <- self$data + } + + PM_bestdose$new( + prior = self, + model = self$model, + past_data = past_data, + target = target, + dose_range = dose_range, + bias_weight = bias_weight, + target_type = target_type, + time_offset = time_offset, + ... + ) + }, + #' @description #' Save the current PM_result object to an .Rdata file. #' @details @@ -216,7 +282,7 @@ PM_result <- R6::R6Class( #' your current working directory, specify `run = 1` to save the result to the "outputs" #' subfolder of the "1" folder. #' @param file Custom file name. Default is "PMout.Rdata". If `run` is not specified, `file` - #' should be the full path and filename. + #' should be the full path and filename. save = function(run, file = "PMout.Rdata") { if (missing(run)) { cli::cli_inform(c( @@ -243,7 +309,7 @@ PM_result <- R6::R6Class( ) save(PMout, file = paste0(outputfolder, "/", file)) }, - + #' @description #' Validate the result by internal simulation methods. #' @param ... Arguments passed to [PM_valid]. @@ -255,9 +321,8 @@ PM_result <- R6::R6Class( " " = "For example, if your results are in {.code my_run}, use {.code my_run$save(1)} to save back to the outputs folder of run 1." )) return(invisible(self)) - }, - + #' @description #' Conduct stepwise linear regression of Bayesian posterior parameter values #' and covariates. @@ -265,7 +330,7 @@ PM_result <- R6::R6Class( step = function(...) { PM_step(self$cov$data, ...) }, - + #' @description #' Calculate optimal sampling times. #' @@ -280,7 +345,7 @@ PM_result <- R6::R6Class( }) return(invisible(self)) }, - + #' @description #' `r lifecycle::badge("deprecated")` #' @@ -344,7 +409,7 @@ PM_result$load <- function(...) { #' run4 <- PM_load(file = "Pmetrics/MyRuns/4/outputs/PMout.Rdata") # loads from Pmetrics/MyRuns/4/outputs/PMout.Rdata #' run5 <- PM_load() # loads from ./PMout.Rdata #' } -#' +#' #' @author Michael Neely and Julian Otalvaro #' @seealso [PM_final], #' [PM_cycle], [PM_op], [PM_cov], @@ -361,20 +426,22 @@ PM_load <- function(run, path = ".", file = "PMout.Rdata") { names(aux_list) <- names(Out)[i] result <- append(result, aux_list) } - + return(result) } - + found <- "" # initialize - + if (!missing(run)) { filepath <- file.path(path, run, "outputs", file) } else { filepath <- file.path(path, file) - } - - if (file.exists(filepath)) { found <- filepath } - + } + + if (file.exists(filepath)) { + found <- filepath + } + if (found != "") { result <- output2List(Out = get(load(found))) rebuild <- PM_result$new(result, path = dirname(found), quiet = TRUE) @@ -398,31 +465,32 @@ update <- function(res, found) { # start conversion n_cyc <- nrow(dat$mean) n_out <- max(res$op$outeq) - dat$gamlam <- dat$gamlam %>% select(starts_with("add")|starts_with("prop")) + dat$gamlam <- dat$gamlam %>% select(starts_with("add") | starts_with("prop")) if (ncol(gamlam) == 1 & n_out > 1) { gamlam <- cbind(gamlam, replicate((n_out - 1), gamlam[, 1])) } gamlam <- gamlam %>% - pivot_longer( - cols = everything(), - values_to = "value", names_to = c("type", "outeq"), - names_sep = "\\." - ) %>% - mutate(cycle = rep(1:n_cyc, each = n_out)) %>% - select(cycle, value, outeq, type) %>% arrange(cycle, outeq) + pivot_longer( + cols = everything(), + values_to = "value", names_to = c("type", "outeq"), + names_sep = "\\." + ) %>% + mutate(cycle = rep(1:n_cyc, each = n_out)) %>% + select(cycle, value, outeq, type) %>% + arrange(cycle, outeq) if (is.matrix(dat$mean)) { # old fortran format, but not rust format dat$mean <- tibble::tibble(cycle = 1:n_cyc) %>% - dplyr::bind_cols(tidyr::as_tibble(dat$mean)) + dplyr::bind_cols(tidyr::as_tibble(dat$mean)) dat$median <- tibble::tibble(cycle = 1:n_cyc) %>% - dplyr::bind_cols(tidyr::as_tibble(dat$median)) + dplyr::bind_cols(tidyr::as_tibble(dat$median)) dat$sd <- tibble::tibble(cycle = 1:n_cyc) %>% - dplyr::bind_cols(tidyr::as_tibble(dat$sd)) + dplyr::bind_cols(tidyr::as_tibble(dat$sd)) } msg <- c(msg, "cycle") res$cycle <- dat } } - + ####### DONE PROCESSING, INFORM ######### if (!is.null(msg)) { cat( @@ -450,6 +518,6 @@ update <- function(res, found) { cat("Results saved\n") } } - + return(res) } diff --git a/R/extendr-wrappers.R b/R/extendr-wrappers.R index 3694efead..e199aec1d 100755 --- a/R/extendr-wrappers.R +++ b/R/extendr-wrappers.R @@ -40,5 +40,13 @@ model_parameters <- function(model_path, kind) .Call(wrap__model_parameters, mod #'@export setup_logs <- function() .Call(wrap__setup_logs) +#' Run BestDose optimization to find optimal doses +#'@export +bestdose <- function(model_path, prior_path, past_data_path, target_data_path, time_offset, dose_min, dose_max, bias_weight, target_type, params, kind) .Call(wrap__bestdose, model_path, prior_path, past_data_path, target_data_path, time_offset, dose_min, dose_max, bias_weight, target_type, params, kind) + +bestdose_prepare <- function(model_path, prior_path, past_data_path, target_data_path, time_offset, dose_min, dose_max, bias_weight, target_type, params, kind) .Call(wrap__bestdose_prepare, model_path, prior_path, past_data_path, target_data_path, time_offset, dose_min, dose_max, bias_weight, target_type, params, kind) + +bestdose_optimize <- function(handle, bias_weight) .Call(wrap__bestdose_optimize, handle, bias_weight) + # nolint end diff --git a/inst/.gitignore b/inst/.gitignore index 28f483682..4472f674a 100755 --- a/inst/.gitignore +++ b/inst/.gitignore @@ -1,3 +1,5 @@ .DS_Store /.quarto/ + +template/ diff --git a/inst/Examples/Rscript/bestdose_simple_test.R b/inst/Examples/Rscript/bestdose_simple_test.R new file mode 100644 index 000000000..7cab124fa --- /dev/null +++ b/inst/Examples/Rscript/bestdose_simple_test.R @@ -0,0 +1,50 @@ +library(Pmetrics) + +setwd("inst/Examples/Runs") + +mod_onecomp <- PM_model$new( + pri = list( + ke = ab(0.001, 3.0), + v = ab(25.0, 250.0) + ), + eqn = function() { + dx[1] <- -ke * X[1] + B[1] + }, + out = function() { + Y[1] <- X[1] / v + }, + err = list( + additive(1, c(0, 0.20, 0, 0)) + ) +) + + +past_file <- "../src/bestdose_past.csv" +target_file <- "../src/bestdose_target.csv" +prior_file <- "../src/bestdose_prior.csv" + + +# Prepare the problem once (posterior + handle to optimized model) +problem <- PM_bestdose_problem$new( + prior = prior_file, + model = mod_onecomp, + past_data = past_file, + target = target_file, + dose_range = list(min = 0, max = 300), + bias_weight = 0.0, + target_type = "concentration" # "concentration", "auc_from_zero", "auc_from_last_dose" +) + +cat("\nPosterior support points:\n") +print(head(problem$theta)) + +# Reuse the same problem for different bias weights +bias_weights <- seq(0, 1, by = 0.25) +results <- lapply(bias_weights, function(lambda) { + problem$optimize(bias_weight = lambda) +}) + +for (i in seq_along(results)) { + cat("\n=== Bias weight:", bias_weights[i], "===\n") + results[[i]]$print() +} diff --git a/inst/Examples/Runs/bestdose_result.rds b/inst/Examples/Runs/bestdose_result.rds new file mode 100644 index 000000000..ead529e5b Binary files /dev/null and b/inst/Examples/Runs/bestdose_result.rds differ diff --git a/inst/Examples/src/bestdose_past.csv b/inst/Examples/src/bestdose_past.csv new file mode 100644 index 000000000..f31bbc5e9 --- /dev/null +++ b/inst/Examples/src/bestdose_past.csv @@ -0,0 +1,9 @@ +"id","evid","time","dur","dose","addl","ii","input","out","outeq","c0","c1","c2","c3" +1,1,0,0,150,0,0,1,.,.,.,.,.,. +1,0,2,.,.,.,.,.,0.759050697604428,1,.,.,.,. +1,0,4,.,.,.,.,.,0.384085169721793,1,.,.,.,. +1,0,6,.,.,.,.,.,0.194349887386702,1,.,.,.,. +1,1,12,0,75,0,0,1,.,.,.,.,.,. +1,0,14,.,.,.,.,.,0.392266577540038,1,.,.,.,. +1,0,16,.,.,.,.,.,0.198489739204705,1,.,.,.,. +1,0,18,.,.,.,.,.,0.100437250648841,1,.,.,.,. diff --git a/inst/Examples/src/bestdose_prior.csv b/inst/Examples/src/bestdose_prior.csv new file mode 100644 index 000000000..0559dea12 --- /dev/null +++ b/inst/Examples/src/bestdose_prior.csv @@ -0,0 +1,47 @@ +ke,v,prob +0.08736658442020416,104.1576635837555,0.06411826509758521 +0.1305751524925232,97.43967413902283,0.03921568628597478 +0.3540655417442322,86.90321147441864,0.03930537754808413 +0.2908282660484314,113.02810430526733,0.05369391213529181 +0.10566180405616761,140.3564077615738,0.01960784315457728 +0.9795492498397828,221.82963848114014,0.019607843139348134 +0.3214192095041275,68.54407012462616,0.019607843027666164 +0.04363290762901306,85.5460512638092,0.01960721617788871 +0.3493796042442322,75.69715678691864,0.019607073213094923 +0.09591740503311158,70.4919171333313,0.019606975069368423 +0.3318302191734314,92.54958868026733,0.019607488410309255 +0.3091094713449478,92.16600060462952,0.01969419594801783 +0.3319534166574478,118.26951622962952,0.019631823530513178 +0.06250557525157929,107.17031538486481,0.019670607599452994 +0.06885235497951508,76.57066643238068,0.01970716827239974 +0.104416601395607,91.60142242908478,0.019644091781941257 +0.11730292952060699,78.54966461658478,0.019597452684180432 +0.06250557525157929,95.70058882236481,0.019558918237699178 +0.013994081234931946,178.82485330104828,0.019739518261060587 +0.3406021231412888,99.99475717544556,0.019495129816772868 +0.07190197534561157,67.0641827583313,0.019509009848450342 +0.02050767450332642,140.61931788921356,0.019494247900270656 +0.283145404958725,125.85132718086243,0.021689589949674865 +0.28171865940093993,140.98521947860718,0.020065175979272123 +0.2902425238609314,94.48318243026733,0.02049255770567427 +0.09111055810451509,93.09410393238068,0.019720100726323925 +0.2925854926109314,128.45290899276733,0.019127402830562693 +0.2925854926109314,105.60134649276733,0.02235654024528764 +0.3019573676109314,127.61794805526733,0.01954540376390492 +0.044804392004013066,99.6524965763092,0.021082697423336363 +0.28490263152122497,92.62867093086243,0.018708872021706954 +0.08112040762901307,98.0704653263092,0.021075404370655606 +0.036757444095611574,109.3835186958313,0.018114822039674315 +0.09826037378311157,111.2731671333313,0.021727337750962878 +0.09298869409561158,113.6022686958313,0.01713443002609149 +0.3132096666574478,115.89646935462952,0.00474084568678399 +0.07658791284561157,113.3385968208313,0.002945117861059519 +0.08678084223270416,87.5463354587555,0.02933991222454541 +0.08736658442020416,128.0199682712555,0.019952065962974042 +0.05615650746822357,130.38859486579895,0.000015108929159013186 +0.056742249655723574,131.13566517829895,0.012605671648896763 +0.08736658442020416,87.5463354587555,0.007152585492468555 +0.07658791284561157,113.3825421333313,0.03210569535127533 +0.3132096666574478,115.94041466712952,0.03478763091025865 +0.05615650746822357,130.43254017829895,0.04616814785864839 +0.08736658442020416,127.9760229587555,0.00001919610085464287 diff --git a/inst/Examples/src/bestdose_target.csv b/inst/Examples/src/bestdose_target.csv new file mode 100644 index 000000000..29d8788a6 --- /dev/null +++ b/inst/Examples/src/bestdose_target.csv @@ -0,0 +1,9 @@ +"id","evid","time","dur","dose","addl","ii","input","out","outeq","c0","c1","c2","c3" +1,1,0,0,0,0,0,1,.,.,.,.,.,. +1,0,2,.,.,.,.,.,0.759050697604428,1,.,.,.,. +1,0,4,.,.,.,.,.,0.384085169721793,1,.,.,.,. +1,0,6,.,.,.,.,.,0.194349887386702,1,.,.,.,. +1,1,12,0,0,0,0,1,.,.,.,.,.,. +1,0,14,.,.,.,.,.,0.392266577540038,1,.,.,.,. +1,0,16,.,.,.,.,.,0.198489739204705,1,.,.,.,. +1,0,18,.,.,.,.,.,0.100437250648841,1,.,.,.,. diff --git a/man/PM_model.Rd b/man/PM_model.Rd index 4eb85ea03..685ee7ba1 100755 --- a/man/PM_model.Rd +++ b/man/PM_model.Rd @@ -32,41 +32,41 @@ a \code{donttest} block to avoid automatic compilation. mod_list <- list( pri = c( - CL = ab(10, 200), - V0 = ab(0, 100), - ka = ab(0, 3), - k23 = ab(0, 5), - k32 = ab(0, 5), - lag1 = ab(0, 2) - ), - cov = c( - wt = interp() - ), - sec = function() { - V = V0 * (wt/70) - ke = CL/V # define here to make eqn simpler - }, - eqn = function() { - dx[1] = -ka * x[1] - dx[2] = rateiv[1] + ka * x[1] - (ke + k23) * x[2] + k32 * x[3] - dx[3] = k23 * x[2] - k32 * x[3] - dx[4] = x[1] / V - }, - lag = function() { - tlag[1] = lag1 - }, - out = function() { - y[1] = x[1]/V - y[2] = x[4] # AUC, not fitted to any data, not required - }, - err = c( - proportional(2, c(0.1, 0.15, 0, 0)) # only applies to y[1] - ) - ) + CL = ab(10, 200), + V0 = ab(0, 100), + ka = ab(0, 3), + k23 = ab(0, 5), + k32 = ab(0, 5), + lag1 = ab(0, 2) + ), + cov = c( + wt = interp() + ), + sec = function() { + V <- V0 * (wt / 70) + ke <- CL / V # define here to make eqn simpler + }, + eqn = function() { + dx[1] <- -ka * x[1] + dx[2] <- rateiv[1] + ka * x[1] - (ke + k23) * x[2] + k32 * x[3] + dx[3] <- k23 * x[2] - k32 * x[3] + dx[4] <- x[1] / V + }, + lag = function() { + tlag[1] <- lag1 + }, + out = function() { + y[1] <- x[1] / V + y[2] <- x[4] # AUC, not fitted to any data, not required + }, + err = c( + proportional(2, c(0.1, 0.15, 0, 0)) # only applies to y[1] + ) +) \donttest{ - mod <- PM_model$new(mod_list) - } +mod <- PM_model$new(mod_list) +} } \author{ @@ -115,10 +115,10 @@ and \emph{error models}. These portions of the model have specific and defined creator functions and no additional R code is permissible. They take this form: -\if{html}{\out{
}}\preformatted{block_name = c( - var1 = creator(), +\if{html}{\out{
}}\preformatted{block_name = c( + var1 = creator(), var2 = creator() -) +) }\if{html}{\out{
}} Note the comma separating the creator functions, "\verb{c(}" to open the vector and "\verb{)}" to close the vector. @@ -128,13 +128,13 @@ equations}, \emph{model equations} (e.g. ODEs), \emph{lag time}, \emph{bioavaila and \emph{outputs}. These parts of the model are defined as R functions without arguments, but whose body contains any permissible R code. -\if{html}{\out{
}}\preformatted{block_name = function() \{ +\if{html}{\out{
}}\preformatted{block_name = function() \{ - # any valid R code + # any valid R code # can use primary or secondary parameters and covariates # lines are not separated by commas -\} +\} }\if{html}{\out{
}} Note the absence of arguments between the "\verb{()}", the opening curly brace "\verb{\{}" to start @@ -225,9 +225,9 @@ Note that \code{wt = interp()} is equivalent to \code{wt = interp("lm")}, since are not estimated for these equations but they are available to every other block in the model. For example: -\if{html}{\out{
}}\preformatted{sec = function() \{ - V = V0 * (wt/70) -\} +\if{html}{\out{
}}\preformatted{sec = function() \{ + V = V0 * (wt/70) +\} }\if{html}{\out{
}} Note that the function @@ -361,7 +361,7 @@ but will be specific to the \code{out} block. For example, -\if{html}{\out{
}}\preformatted{out = function() \{ +\if{html}{\out{
}}\preformatted{out = function() \{ V = V0 * wt # only needed if not included in sec block y[1] = x[1]/V #Vp and Vm must be defined in pri or sec blocks diff --git a/man/PM_result.Rd b/man/PM_result.Rd index c6a5e7e56..2648cc753 100755 --- a/man/PM_result.Rd +++ b/man/PM_result.Rd @@ -14,6 +14,31 @@ After a run completes, results are stored on your hard drive. They are loaded back into R with \link{PM_load} to create the \link{PM_result} object, which contains both the results and functions to analyze or plot the result. } +\examples{ + +## ------------------------------------------------ +## Method `PM_result$bestdose` +## ------------------------------------------------ + +\dontrun{ +# Load NPAG result +result <- PM_load(1) + +# Create target data +target <- PM_data$new("target.csv") + +# Run BestDose optimization +bd <- result$bestdose( + target = target, + dose_range = list(min = 50, max = 500), + bias_weight = 0.5 +) + +# View results +print(bd) +bd$get_doses() +} +} \author{ Michael Neely, Julian Otalvaro } @@ -67,6 +92,7 @@ new optimal sampling results.} \item \href{#method-PM_result-nca}{\code{PM_result$nca()}} \item \href{#method-PM_result-report}{\code{PM_result$report()}} \item \href{#method-PM_result-sim}{\code{PM_result$sim()}} +\item \href{#method-PM_result-bestdose}{\code{PM_result$bestdose()}} \item \href{#method-PM_result-save}{\code{PM_result$save()}} \item \href{#method-PM_result-validate}{\code{PM_result$validate()}} \item \href{#method-PM_result-step}{\code{PM_result$step()}} @@ -247,6 +273,85 @@ arguments, e.g. \verb{$sim(include = 1:2, predInt = 1, limits = NA)}.} } \if{html}{\out{
}} } +} +\if{html}{\out{
}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-PM_result-bestdose}{}}} +\subsection{Method \code{bestdose()}}{ +Run BestDose optimization using this result as the prior +\subsection{Usage}{ +\if{html}{\out{
}}\preformatted{PM_result$bestdose( + target, + past_data = NULL, + dose_range = list(min = 0, max = 1000), + bias_weight = 0.5, + target_type = "concentration", + time_offset = NULL, + ... +)}\if{html}{\out{
}} +} + +\subsection{Arguments}{ +\if{html}{\out{
}} +\describe{ +\item{\code{target}}{PM_data object or path to CSV with target doses/observations. +Required. This defines the dosing template and target values. Set dose amounts +to 0 for doses to be optimized.} + +\item{\code{past_data}}{Optional: PM_data object or path to CSV with patient history. +If NULL (default), uses the full dataset from this result.} + +\item{\code{dose_range}}{Named list with min and max allowable doses. +Default: list(min = 0, max = 1000)} + +\item{\code{bias_weight}}{Numeric \link{0,1} controlling personalization level. +0 = fully personalized, 1 = population-based. Default: 0.5 (balanced)} + +\item{\code{target_type}}{One of "concentration" (default), "auc_from_zero", or +"auc_from_last_dose"} + +\item{\code{time_offset}}{Optional: time offset for past/future concatenation} + +\item{\code{...}}{Additional arguments passed to \link{PM_bestdose}} +} +\if{html}{\out{
}} +} +\subsection{Details}{ +BestDose finds optimal dosing regimens to achieve target drug concentrations +or AUC values. The algorithm uses Bayesian posterior estimation combined with +dual optimization to balance patient-specific adaptation and population-level +robustness. By default, uses the \verb{$final}, \verb{$model}, and \verb{$data} objects from +this result. Most commonly, you will supply a different \code{target} data object. +} + +\subsection{Returns}{ +A \link{PM_bestdose} object containing optimal doses and predictions +} +\subsection{Examples}{ +\if{html}{\out{
}} +\preformatted{\dontrun{ +# Load NPAG result +result <- PM_load(1) + +# Create target data +target <- PM_data$new("target.csv") + +# Run BestDose optimization +bd <- result$bestdose( + target = target, + dose_range = list(min = 50, max = 500), + bias_weight = 0.5 +) + +# View results +print(bd) +bd$get_doses() +} +} +\if{html}{\out{
}} + +} + } \if{html}{\out{
}} \if{html}{\out{}} diff --git a/man/interp.Rd b/man/interp.Rd index 0d7eb549e..6a75496a5 100755 --- a/man/interp.Rd +++ b/man/interp.Rd @@ -24,7 +24,7 @@ interpolation between values or not. } \examples{ \dontrun{ -cov = c( +cov <- c( wt = interp(), # same as interp("lm") or interp("linear") visit = interp("none") ) diff --git a/man/makeNCA.Rd b/man/makeNCA.Rd index 4dc0ffcba..737122886 100755 --- a/man/makeNCA.Rd +++ b/man/makeNCA.Rd @@ -25,7 +25,7 @@ makeNCA( \itemize{ \item It can be the run number, e.g. 3, that has been previously loaded with \link{PM_load}. Either the mdata file from the run (NPAG or IT2B) can be used (default) or -the post object can be used (NPAG only) by specifying \code{postPred = T} below. If \link{x} is a run number that corresponds +the post object can be used (NPAG only) by specifying \code{postPred = T} below. If \link[ggplot2:aes_position]{ggplot2::x} is a run number that corresponds to both an NPAG and IT2B run which have been previously loaded into memory with \link{PM_load}, the NPAG run will be used. \item It can be the run number of a run that has \emph{not} been previously loaded with \link{PM_load}. In this case, the current working directory should be the Runs folder as \link{makeNCA} will call \link{PM_load}. diff --git a/man/plot.PM_cov.Rd b/man/plot.PM_cov.Rd index ea54d695a..fedc8dc5b 100755 --- a/man/plot.PM_cov.Rd +++ b/man/plot.PM_cov.Rd @@ -171,7 +171,7 @@ NPex$cov$plot(V ~ wt, } } \seealso{ -\link{PM_cov}, \link{PM_result}, \link{schema} +\link{PM_cov}, \link{PM_result}, \link[plotly:schema]{plotly::schema} Other PMplots: \code{\link{plot.PM_cycle}()}, diff --git a/man/plot.PM_cycle.Rd b/man/plot.PM_cycle.Rd index 67acb7e6b..f91068602 100755 --- a/man/plot.PM_cycle.Rd +++ b/man/plot.PM_cycle.Rd @@ -114,7 +114,7 @@ NPex$cycle$plot( } } \seealso{ -\link{PM_result}, \link{schema} +\link{PM_result}, \link[plotly:schema]{plotly::schema} Other PMplots: \code{\link{plot.PM_cov}()}, diff --git a/man/plot.PM_final.Rd b/man/plot.PM_final.Rd index 9e2374ae1..5b9129fab 100755 --- a/man/plot.PM_final.Rd +++ b/man/plot.PM_final.Rd @@ -161,7 +161,7 @@ ITex$final$plot(Ke ~ V) } } \seealso{ -\link{PM_final}, \link{schema} +\link{PM_final}, \link[plotly:schema]{plotly::schema} Other PMplots: \code{\link{plot.PM_cov}()}, diff --git a/man/plot.PM_op.Rd b/man/plot.PM_op.Rd index d3696144b..ed683317a 100755 --- a/man/plot.PM_op.Rd +++ b/man/plot.PM_op.Rd @@ -187,7 +187,7 @@ NPex$op$plot(stats = list(x = 0.5, y = 0.2, font = list(size = 7, color = "blue" } } \seealso{ -\link{PM_result}, \link{PM_op}, \link{schema} +\link{PM_result}, \link{PM_op}, \link[plotly:schema]{plotly::schema} Other PMplots: \code{\link{plot.PM_cov}()}, diff --git a/man/plot.PM_pta.Rd b/man/plot.PM_pta.Rd index 833e52de3..06168a61f 100755 --- a/man/plot.PM_pta.Rd +++ b/man/plot.PM_pta.Rd @@ -61,7 +61,7 @@ the \code{line} argument should be a list of the options for group based plottin where each group corresponds to a simulated regimen. The possible elements of the \code{line} list should be exactly named: \itemize{ -\item color Maps to the \link{plot_ly} \code{colors} argument to override default colors +\item color Maps to the \link[plotly:plot_ly]{plotly::plot_ly} \code{colors} argument to override default colors applied to the lines for each regimen. This can be a named palette, which can be obtained with \code{RColorBrewer::display.brewer.all()} or a vector of hexadecimal color names. One way to ensure reliable color palettes is to use the @@ -72,12 +72,12 @@ of JavaScript on the ColorBrewer website. The default is "Set1". Palettes with fewer colors than regimens will be recycled. A color can also be a character vector of color names, recycled as needed. For example, a print-friendly choice is \code{line = list(color = "black")}. -\item width Maps to the \link{plot_ly} \code{width} argument to override default widths +\item width Maps to the \link[plotly:plot_ly]{plotly::plot_ly} \code{width} argument to override default widths applied to the lines for each regimen. All lines will have the same width. The default value is 2. -\item dash Maps to the \link{plot_ly} \code{linetypes} argument to override default styles +\item dash Maps to the \link[plotly:plot_ly]{plotly::plot_ly} \code{linetypes} argument to override default styles applied to the lines for each regimen. If numeric, will map to \code{lty} \link{par} values. -It can also be a character vector of dash names as listed in \link{plot_ly}. +It can also be a character vector of dash names as listed in \link[plotly:plot_ly]{plotly::plot_ly}. Example: \code{line = list(color = "Blues", width = 1, dash = 2)}, whicb will result in dotted lines (dash = 2) all with width 1 but in different shades of blue. }} @@ -97,12 +97,12 @@ where each group corresponds to a simulated regimen. The possible elements of th marker color does not need to also be specified. Even if line plotting is suppressed with \code{line = F}, the default color value of "Set1" will be applied to markers, unless specified, e.g. \code{marker = list(color = "Blues")}. -\item symbol Maps to the \link{plot_ly} \code{symbols} argument to override default symbols +\item symbol Maps to the \link[plotly:plot_ly]{plotly::plot_ly} \code{symbols} argument to override default symbols applied to the markers for each regimen. If only one value is supplied for this, it will be recycled for each regimen, i.e. all will have the same symbol. See \code{plotly::schema()}, traces > scatter > attributes > marker > symbol > values for options. -\item size Maps to the \link{plot_ly} \code{size} argument to override default size +\item size Maps to the \link[plotly:plot_ly]{plotly::plot_ly} \code{size} argument to override default size applied to the markers for each regimen. All markers will have the same size. The default value is 12. }} diff --git a/man/plot.PM_sim.Rd b/man/plot.PM_sim.Rd index 2cf39f45c..7278bdb05 100755 --- a/man/plot.PM_sim.Rd +++ b/man/plot.PM_sim.Rd @@ -191,7 +191,7 @@ simEx$plot(log = FALSE, line = list(color = "orange")) } } \seealso{ -\link{PM_sim}, \link{plot_ly}, \link{schema} +\link{PM_sim}, \link[plotly:plot_ly]{plotly::plot_ly}, \link[plotly:schema]{plotly::schema} Other PMplots: \code{\link{plot.PM_cov}()}, diff --git a/src/rust/Cargo.toml b/src/rust/Cargo.toml index bca68233b..d8fe6025a 100755 --- a/src/rust/Cargo.toml +++ b/src/rust/Cargo.toml @@ -9,8 +9,8 @@ name = 'pm_rs' [dependencies] extendr-api = '*' -pmcore = {version ="=0.21.0", features = ["exa"]} -# pmcore = { path = "/Users/jotalvaro/code/LAPKB/PMcore", features = ["exa"] } +pmcore = { version = "=0.22.1", features = ["exa"] } +libloading = "0.8" rayon = "1.10.0" anyhow = "1.0.97" diff --git a/src/rust/src/bestdose_executor.rs b/src/rust/src/bestdose_executor.rs new file mode 100644 index 000000000..2d83731f6 --- /dev/null +++ b/src/rust/src/bestdose_executor.rs @@ -0,0 +1,426 @@ +use crate::{logs::RFormatLayer, settings::settings}; +use extendr_api::prelude::*; +use pmcore::bestdose::{BestDoseProblem, BestDoseResult, DoseRange, Target}; +use pmcore::prelude::{data, ODE}; +use pmcore::routines::initialization::parse_prior; +use std::path::PathBuf; + +/// Helper to parse target type from string +pub(crate) fn parse_target_type(target_str: &str) -> std::result::Result { + match target_str.to_lowercase().as_str() { + "concentration" => Ok(Target::Concentration), + "auc_from_zero" | "auc" => Ok(Target::AUCFromZero), + "auc_from_last_dose" | "auc_interval" => Ok(Target::AUCFromLastDose), + _ => Err(format!( + "Invalid target type: {}. Must be 'concentration', 'auc_from_zero', or 'auc_from_last_dose'", + target_str + )), + } +} + +/// R-compatible prediction row for BestDose output +#[derive(Debug, IntoDataFrameRow)] +pub struct BestDosePredictionRow { + id: String, + time: f64, + observed: f64, + pop_mean: f64, + pop_median: f64, + post_mean: f64, + post_median: f64, + outeq: usize, +} + +impl BestDosePredictionRow { + pub fn from_np_prediction( + pred: &pmcore::routines::output::predictions::NPPredictionRow, + id: &str, + ) -> Self { + Self { + id: id.to_string(), + time: pred.time(), + observed: pred.obs().unwrap_or(0.0), + pop_mean: pred.pop_mean(), + pop_median: pred.pop_median(), + post_mean: pred.post_mean(), + post_median: pred.post_median(), + outeq: pred.outeq(), + } + } +} + +/// R-compatible AUC prediction row +#[derive(Debug, IntoDataFrameRow)] +pub struct BestDoseAucRow { + time: f64, + auc: f64, +} + +/// Convert BestDoseResult to R-compatible list structure +pub(crate) fn convert_bestdose_result_to_r( + result: BestDoseResult, +) -> std::result::Result { + // Extract doses + let doses: Vec = result.doses(); + + // Objective function + let objf = result.objf(); + + // Status + let status_str = format!("{:?}", result.status()); + + // Predictions as data frame + let pred_rows: Vec = result + .predictions() + .predictions() + .iter() + .map(|p| BestDosePredictionRow::from_np_prediction(p, "subject_1")) + .collect(); + let pred_df = pred_rows + .into_dataframe() + .map_err(|e| format!("Failed to create predictions dataframe: {:?}", e))?; + + // AUC predictions (if available) + let auc_val = if let Some(auc_preds) = result.auc_predictions() { + let auc_rows: Vec = auc_preds + .iter() + .map(|(time, auc)| BestDoseAucRow { + time: *time, + auc: *auc, + }) + .collect(); + let auc_df = auc_rows + .into_dataframe() + .map_err(|e| format!("Failed to create AUC dataframe: {:?}", e))?; + Robj::from(auc_df) + } else { + Robj::from(()) // NULL for no AUC + }; + + // Optimization method + let method_str = format!("{}", result.optimization_method()); + + // Build the list using list! macro + let output = list!( + doses = doses, + objf = objf, + status = status_str, + predictions = pred_df, + auc_predictions = auc_val, + method = method_str + ); + + Ok(output.into()) +} +/// Opaque handle that keeps the dynamic model library alive while reusing the +/// prepared `BestDoseProblem` for multiple optimization runs. +pub struct BestDoseProblemHandle { + problem: BestDoseProblem, + #[allow(dead_code)] + library: libloading::Library, +} + +impl BestDoseProblemHandle { + #[allow(clippy::too_many_arguments)] + pub fn new( + model_path: PathBuf, + prior_path: PathBuf, + past_data_path: Option, + target_data_path: PathBuf, + time_offset: Option, + dose_min: f64, + dose_max: f64, + bias_weight: f64, + target_type: &str, + params: List, + ) -> std::result::Result { + let (library, (eq, meta)) = + unsafe { pmcore::prelude::pharmsol::exa::load::load::(model_path) }; + + let settings = settings(params, meta.get_params(), "/tmp/bestdose") + .map_err(|e| format!("Failed to parse settings: {}", e))?; + + let (population_theta, prior_weights) = + parse_prior(&prior_path.to_str().unwrap().to_string(), &settings) + .map_err(|e| format!("Failed to parse prior: {}", e))?; + + let population_weights = prior_weights + .ok_or_else(|| "Prior file must contain a 'prob' column with weights".to_string())?; + + let past_data = if let Some(path) = past_data_path { + let data = data::read_pmetrics(path.to_str().unwrap()) + .map_err(|e| format!("Failed to read past data: {}", e))?; + let subjects = data.subjects(); + if subjects.is_empty() { + return Err("Past data file contains no subjects".to_string()); + } + Some(subjects[0].clone()) + } else { + None + }; + + let target_data = { + let data = data::read_pmetrics(target_data_path.to_str().unwrap()) + .map_err(|e| format!("Failed to read target data: {}", e))?; + let subjects = data.subjects(); + if subjects.is_empty() { + return Err("Target data file contains no subjects".to_string()); + } + subjects[0].clone() + }; + + let target_enum = parse_target_type(target_type)?; + let doserange = DoseRange::new(dose_min, dose_max); + + let problem = BestDoseProblem::new( + &population_theta, + &population_weights, + past_data, + target_data, + time_offset, + eq, + doserange, + bias_weight, + settings, + target_enum, + ) + .map_err(|e| format!("Failed to create BestDose problem: {}", e))?; + + Ok(Self { problem, library }) + } + + pub fn optimize( + &self, + bias_weight: Option, + ) -> std::result::Result { + let configured_problem = match bias_weight { + Some(weight) => self.problem.clone().with_bias_weight(weight), + None => self.problem.clone(), + }; + + configured_problem + .optimize() + .map_err(|e| format!("Optimization failed: {}", e)) + } + + pub fn problem(&self) -> &BestDoseProblem { + &self.problem + } +} + +pub(crate) fn bestdose_ode( + model_path: PathBuf, + prior_path: PathBuf, + past_data_path: Option, + target_data_path: PathBuf, + time_offset: Option, + dose_min: f64, + dose_max: f64, + bias_weight: f64, + target_type: &str, + params: List, +) -> std::result::Result { + let handle = BestDoseProblemHandle::new( + model_path, + prior_path, + past_data_path, + target_data_path, + time_offset, + dose_min, + dose_max, + bias_weight, + target_type, + params, + )?; + + handle.optimize(None) +} + +/// Execute bestdose optimization for analytical models (placeholder - not yet supported) +pub(crate) fn bestdose_analytical( + _model_path: PathBuf, + _prior_path: PathBuf, + _past_data_path: Option, + _target_data_path: PathBuf, + _time_offset: Option, + _dose_min: f64, + _dose_max: f64, + _bias_weight: f64, + _target_type: &str, + _params: List, +) -> std::result::Result { + Err("BestDose for analytical models is not yet supported".to_string()) +} + +pub(crate) struct PosteriorSummary { + theta_values: Vec, + theta_dim: (i32, i32), + param_names: Vec, + posterior_weights: Vec, + population_weights: Vec, + bias_weight: f64, + target_type: Target, +} + +fn summarize_problem(problem: &BestDoseProblem) -> PosteriorSummary { + let theta = problem.posterior_theta(); + let matrix = theta.matrix(); + let nrows = matrix.nrows() as i32; + let ncols = matrix.ncols() as i32; + let mut theta_values = vec![0.0; (nrows * ncols) as usize]; + + for col in 0..ncols as usize { + for row in 0..nrows as usize { + theta_values[row + col * nrows as usize] = *matrix.get(row, col); + } + } + + PosteriorSummary { + theta_values, + theta_dim: (nrows, ncols), + param_names: theta.param_names(), + posterior_weights: problem.posterior_weights().to_vec(), + population_weights: problem.population_weights().to_vec(), + bias_weight: problem.bias_weight(), + target_type: problem.target_type(), + } +} + +fn vec_to_doubles(values: Vec, label: &str) -> std::result::Result { + Doubles::try_from(values) + .map_err(|e| format!("Failed to convert {} to doubles: {:?}", label, e)) +} + +fn dims_to_integers(dim: (i32, i32)) -> std::result::Result { + Integers::try_from(vec![dim.0, dim.1]) + .map_err(|e| format!("Failed to convert dims to integers: {:?}", e)) +} + +fn names_to_strings(names: &[String]) -> Strings { + Strings::from_values(names.iter().map(|s| s.as_str())) +} + +pub(crate) fn prepare_bestdose_problem( + model_path: PathBuf, + prior_path: PathBuf, + past_data_path: Option, + target_data_path: PathBuf, + time_offset: Option, + dose_min: f64, + dose_max: f64, + bias_weight: f64, + target_type: &str, + params: List, +) -> std::result::Result<(BestDoseProblemHandle, PosteriorSummary), String> { + let handle = BestDoseProblemHandle::new( + model_path, + prior_path, + past_data_path, + target_data_path, + time_offset, + dose_min, + dose_max, + bias_weight, + target_type, + params, + )?; + + let summary = summarize_problem(handle.problem()); + Ok((handle, summary)) +} + +pub(crate) fn bestdose_prepare_internal( + model_path: &str, + prior_path: &str, + past_data_path: Nullable, + target_data_path: &str, + time_offset: Nullable, + dose_min: f64, + dose_max: f64, + bias_weight: f64, + target_type: &str, + params: List, + kind: &str, +) -> Robj { + RFormatLayer::reset_global_timer(); + let _ = crate::setup_logs(); + + let past_path = past_data_path.into_option().map(PathBuf::from); + let time_offset = time_offset.into_option(); + + let preparation = match kind { + "ode" => prepare_bestdose_problem( + PathBuf::from(model_path), + PathBuf::from(prior_path), + past_path, + PathBuf::from(target_data_path), + time_offset, + dose_min, + dose_max, + bias_weight, + target_type, + params.clone(), + ), + "analytical" => Err("BestDose for analytical models is not yet supported".to_string()), + other => Err(format!("{} is not a supported model type", other)), + }; + + match preparation { + Ok((handle, summary)) => { + let theta_values = match vec_to_doubles(summary.theta_values, "theta_values") { + Ok(values) => values, + Err(e) => return Robj::from(e), + }; + let theta_dim = match dims_to_integers(summary.theta_dim) { + Ok(dim) => dim, + Err(e) => return Robj::from(e), + }; + let posterior_weights = + match vec_to_doubles(summary.posterior_weights, "posterior_weights") { + Ok(values) => values, + Err(e) => return Robj::from(e), + }; + let population_weights = + match vec_to_doubles(summary.population_weights, "population_weights") { + Ok(values) => values, + Err(e) => return Robj::from(e), + }; + let param_names = names_to_strings(&summary.param_names); + let handle_ptr = ExternalPtr::new(handle); + + let output = list!( + handle = handle_ptr, + theta_values = theta_values, + theta_dim = theta_dim, + param_names = param_names, + posterior_weights = posterior_weights, + population_weights = population_weights, + bias_weight = summary.bias_weight, + target_type = format!("{:?}", summary.target_type), + nspp = summary.theta_dim.0, + n_parameters = summary.theta_dim.1 + ); + + output.into() + } + Err(e) => Robj::from(format!("BestDose prepare failed: {}", e)), + } +} + +pub(crate) fn bestdose_optimize_internal( + handle: ExternalPtr, + bias_weight: Nullable, +) -> Robj { + let weight = bias_weight.into_option(); + + match handle.try_addr() { + Ok(inner) => match inner.optimize(weight) { + Ok(result) => match convert_bestdose_result_to_r(result) { + Ok(robj) => robj, + Err(e) => Robj::from(format!("Failed to convert result: {}", e)), + }, + Err(e) => Robj::from(format!("BestDose optimization failed: {}", e)), + }, + Err(e) => Robj::from(format!("Invalid BestDose handle: {}", e)), + } +} diff --git a/src/rust/src/executor.rs b/src/rust/src/executor.rs index 2be93ee94..8cf5461f0 100755 --- a/src/rust/src/executor.rs +++ b/src/rust/src/executor.rs @@ -42,7 +42,7 @@ pub(crate) fn fit( let data = data::read_pmetrics(data.to_str().unwrap()).expect("Failed to read data"); //dbg!(&data); let mut algorithm = dispatch_algorithm(settings, eq, data)?; - let result = algorithm.fit()?; + let mut result = algorithm.fit()?; result.write_outputs()?; Ok(()) } diff --git a/src/rust/src/lib.rs b/src/rust/src/lib.rs index 3ee6cd5a4..ffa43b6fc 100755 --- a/src/rust/src/lib.rs +++ b/src/rust/src/lib.rs @@ -1,4 +1,5 @@ // mod build; +mod bestdose_executor; mod executor; mod logs; mod settings; @@ -8,9 +9,11 @@ use anyhow::Result; use extendr_api::prelude::*; use pmcore::prelude::{data::read_pmetrics, pharmsol::exa::build, Analytical, ODE}; use simulation::SimulationRow; +use std::path::PathBuf; use std::process::Command; use tracing_subscriber::layer::SubscriberExt; +use crate::bestdose_executor::BestDoseProblemHandle; use crate::logs::RFormatLayer; fn validate_paths(data_path: &str, model_path: &str) { @@ -247,6 +250,110 @@ fn setup_logs() -> anyhow::Result<()> { Ok(()) } +/// Run BestDose optimization to find optimal doses +///@export +#[extendr] +fn bestdose( + model_path: &str, + prior_path: &str, + past_data_path: Nullable, + target_data_path: &str, + time_offset: Nullable, + dose_min: f64, + dose_max: f64, + bias_weight: f64, + target_type: &str, + params: List, + kind: &str, +) -> Robj { + RFormatLayer::reset_global_timer(); + let _ = setup_logs(); + + println!("Starting BestDose optimization..."); + + let past_path = match past_data_path.into_option() { + Some(p) => Some(PathBuf::from(p)), + None => None, + }; + + let time_offset_opt = time_offset.into_option(); + + let result = match kind { + "ode" => bestdose_executor::bestdose_ode( + model_path.into(), + prior_path.into(), + past_path, + target_data_path.into(), + time_offset_opt, + dose_min, + dose_max, + bias_weight, + target_type, + params.clone(), + ), + "analytical" => bestdose_executor::bestdose_analytical( + model_path.into(), + prior_path.into(), + past_path, + target_data_path.into(), + time_offset_opt, + dose_min, + dose_max, + bias_weight, + target_type, + params.clone(), + ), + _ => { + return Robj::from(format!("{} is not a supported model type", kind)); + } + }; + + match result { + Ok(bd_result) => match bestdose_executor::convert_bestdose_result_to_r(bd_result) { + Ok(r) => r, + Err(e) => Robj::from(format!("Failed to convert result: {}", e)), + }, + Err(e) => Robj::from(format!("BestDose failed: {}", e)), + } +} + +#[extendr] +fn bestdose_prepare( + model_path: &str, + prior_path: &str, + past_data_path: Nullable, + target_data_path: &str, + time_offset: Nullable, + dose_min: f64, + dose_max: f64, + bias_weight: f64, + target_type: &str, + params: List, + kind: &str, +) -> Robj { + bestdose_executor::bestdose_prepare_internal( + model_path, + prior_path, + past_data_path, + target_data_path, + time_offset, + dose_min, + dose_max, + bias_weight, + target_type, + params, + kind, + ) +} + +#[extendr] +fn bestdose_optimize( + handle: ExternalPtr, + bias_weight: Nullable, +) -> Robj { + bestdose_executor::bestdose_optimize_internal(handle, bias_weight) +} + // Macro to generate exports. // This ensures exported functions are registered with R. // See corresponding C code in `entrypoint.c`. @@ -260,6 +367,9 @@ extendr_module! { fn fit; fn model_parameters; fn setup_logs; + fn bestdose; + fn bestdose_prepare; + fn bestdose_optimize; } // To generate the exported function in R, run the following command: