@@ -30,7 +30,7 @@ PM_fit <- R6::R6Class(
30
30
model = NULL ,
31
31
# ' @field backend Backend used for calculations; default is value in PMoptions.
32
32
backend = NULL ,
33
-
33
+
34
34
# ' @description
35
35
# ' Create a new object
36
36
# ' @param data Either the name of a [PM_data]
@@ -60,45 +60,45 @@ PM_fit <- R6::R6Class(
60
60
if (! inherits(model , " PM_model" )) {
61
61
cli :: cli_abort(c(" x" = " {.code model} must be a {.cls PM_model} object" ))
62
62
}
63
-
63
+
64
64
# ### checks
65
-
65
+
66
66
# covariates
67
67
dataCov <- tolower(getCov(data )$ covnames )
68
68
modelCov <- tolower(sapply(model $ model_list $ cov , function (x ) x $ covariate ))
69
- if (length(modelCov )== 0 ) {
69
+ if (length(modelCov ) == 0 ) {
70
70
modelCov <- NA
71
71
}
72
- if (! all(is.na(dataCov )) && ! all(is.na(modelCov ))){ # if there are covariates
73
- if (! identical(dataCov , modelCov )){ # if not identical, abort
72
+ if (! all(is.na(dataCov )) && ! all(is.na(modelCov ))) { # if there are covariates
73
+ if (! identical(dataCov , modelCov )) { # if not identical, abort
74
74
msg <- glue :: glue(" Model covariates: {paste(modelCov, collapse = ', ')}; Data covariates: {paste(dataCov, collapse = ', ')}" )
75
75
cli :: cli_abort(c(
76
76
" x" = " Error: Covariates in data and model do not match." ,
77
77
" i" = msg
78
78
))
79
79
}
80
- }
81
-
82
- # output equations
83
-
80
+ }
81
+
82
+ # output equations
83
+
84
84
if (! is.null(data $ standard_data $ outeq )) {
85
85
dataOut <- max(data $ standard_data $ outeq , na.rm = TRUE )
86
86
} else {
87
87
dataOut <- 1
88
88
}
89
-
89
+
90
90
modelOut <- length(model $ model_list $ out )
91
- if (dataOut != modelOut ){
91
+ if (dataOut != modelOut ) {
92
92
cli :: cli_abort(c(
93
93
" x" = " Error: Number of output equations in data and model do not match." ,
94
94
" i" = " Check the number of output equations in the data and model."
95
95
))
96
96
}
97
-
97
+
98
98
self $ data <- data
99
99
self $ model <- model
100
100
self $ backend <- backend
101
-
101
+
102
102
if (backend == " rust" ) {
103
103
private $ setup_rust_execution()
104
104
}
@@ -111,10 +111,10 @@ PM_fit <- R6::R6Class(
111
111
# ' the simplest execution of this method is
112
112
# ' `$run()`.
113
113
# ' @param run Specify the run number of the output folder. Default if missing is the next available number.
114
- # ' @param include Vector of subject id values in the data file to include in the analysis.
114
+ # ' @param include Vector of subject id values in the data file to include in the analysis.
115
115
# ' The default (missing) is all.
116
116
# ' @param exclude A vector of subject IDs to exclude in the analysis, e.g. `c(4,6:14,16:20)`
117
- # #' @param ode Ordinary Differential Equation solver log tolerance or stiffness.
117
+ # #' @param ode Ordinary Differential Equation solver log tolerance or stiffness.
118
118
# Default is -4, i.e. 0.0001. Higher values will result in faster
119
119
# #' runs, but parameter estimates may not be as accurate.
120
120
# #' @param tol Tolerance for convergence of NPAG. Smaller numbers make it harder to converge.
@@ -126,16 +126,16 @@ PM_fit <- R6::R6Class(
126
126
# ' * The default is "sobol", which is a semi-random distribution. This is the distribution
127
127
# ' typically used when fitting a new model to the data. An example of this is
128
128
# ' on our [website](https://www.lapk.org/images/sobol_3d_plot.html).
129
- # '
129
+ # '
130
130
# ' The following all specify non-random, informative prior distributions. They
131
131
# ' are useful for either continuing a previous
132
132
# ' run which did not converge or for fitting a model to new data, whether to simply
133
133
# ' calculate Bayesian posteriors with `cycles = 0` or to revise the model to a new
134
134
# ' covergence with the new data.
135
135
# ' * The name of a suitable [PM_result] object from a prior run loaded with [PM_load].
136
- # ' This starts from the non-uniform, informative distribution obtained at the end of a prior NPAG run.
136
+ # ' This starts from the non-uniform, informative distribution obtained at the end of a prior NPAG run.
137
137
# ' Example: `run1 <- PM_load(1); fit1$run(prior = run1)`.
138
- # '
138
+ # '
139
139
# ' * A character string with the filename of a csv file containing a prior distribution with
140
140
# ' format as for 'theta.csv' in the output folder of a prior run: column headers are parameter
141
141
# ' names, and rows are the support point values. A final column with probabilities
@@ -147,21 +147,21 @@ PM_fit <- R6::R6Class(
147
147
# ' * A data frame obtained from reading an approriate file, such that the data frame
148
148
# ' is in the required format described in the filename option above. Example:
149
149
# ' `mytheta <- read_csv("mytheta.csv"); fit1$run(prior = mytheta)`.
150
- # '
151
- # ' @param density0 The proportion of the volume of the model parameter
150
+ # '
151
+ # ' @param density0 The proportion of the volume of the model parameter
152
152
# ' hyperspace used to calculate the initial number of support points if one of
153
153
# ' the semi-random, uniform distributions are selected in the `prior` argument
154
- # ' above. The initial points are
154
+ # ' above. The initial points are
155
155
# ' spread through that hyperspace and begin the search for the optimal
156
156
# ' parameter value distribution (support points) in the population.
157
157
# ' The volume of the parameter space is the product of the ranges for all parameters.
158
158
# ' For example if using two parameters `Ke` and `V`, with ranges of \[0, 5\] and \[10, 100\],
159
159
# ' the volume is (5 - 0) x (100 - 10) = 450 The default value of `density0` is 0.01, so the initial
160
160
# ' number of support points will be 0.01 x 450 = 4.5, increased to the nearest integer,
161
- # ' which is 5. The greater the initial number of points, the less chance of
162
- # ' missing the globally maximally likely parameter value distribution,
161
+ # ' which is 5. The greater the initial number of points, the less chance of
162
+ # ' missing the globally maximally likely parameter value distribution,
163
163
# ' but the slower the run.
164
- # '
164
+ # '
165
165
# #' @param indpts Index of starting grid point number. Default is missing, which allows NPAG to choose depending on the number of random parameters:
166
166
# #' 1 or 2 = index of 1; 3 = 3; 4 = 4, 5 = 6,
167
167
# #' 6 or more is 10+number of multiples for each parameter greater than 5, e.g. 6 = 101; 7 = 102, up to 108 for 13 or more parameters.
@@ -191,7 +191,7 @@ PM_fit <- R6::R6Class(
191
191
# '
192
192
# ' @author Michael Neely
193
193
# ' @export
194
-
194
+
195
195
run = function (run = NULL ,
196
196
include = NULL , exclude = NULL ,
197
197
# ode, tol, salt,
@@ -211,9 +211,9 @@ PM_fit <- R6::R6Class(
211
211
report = getPMoptions(" report_template" ),
212
212
artifacts = TRUE ) {
213
213
cwd <- getwd()
214
- intern <- TRUE # always true until (if) rust can run separately from R
215
-
216
-
214
+ intern <- TRUE # always true until (if) rust can run separately from R
215
+
216
+
217
217
# make new output directory
218
218
if (is.null(run )) {
219
219
olddir <- list.dirs(recursive = FALSE )
@@ -231,7 +231,7 @@ PM_fit <- R6::R6Class(
231
231
newdir <- as.character(run )
232
232
}
233
233
}
234
-
234
+
235
235
if (file.exists(newdir )) {
236
236
if (overwrite ) {
237
237
unlink(newdir , recursive = TRUE )
@@ -244,77 +244,84 @@ PM_fit <- R6::R6Class(
244
244
}
245
245
dir.create(newdir )
246
246
setwd(newdir )
247
-
247
+
248
248
algorithm <- tolower(algorithm )
249
-
249
+
250
250
if (getPMoptions()$ backend != " rust" ) {
251
251
setwd(cwd )
252
252
cli :: cli_abort(c(
253
253
" x" = " Error: unsupported backend." ,
254
254
" i" = " See help for {.fn setPMoptions}"
255
255
))
256
256
}
257
-
257
+
258
258
if (artifacts ) {
259
259
self $ model $ write(" model.txt" )
260
260
}
261
-
261
+
262
262
# ### Include or exclude subjects ####
263
263
if (is.null(include )) include <- unique(self $ data $ standard_data $ id )
264
264
if (is.null(exclude )) exclude <- NA
265
265
data_filtered <- self $ data $ standard_data %> % includeExclude(include , exclude )
266
-
266
+
267
267
if (nrow(data_filtered ) == 0 ) {
268
268
cli :: cli_abort(" x" = " No subjects remain after filtering." )
269
269
setwd(cwd )
270
270
return (invisible (NULL ))
271
271
}
272
-
273
-
274
-
275
-
276
-
272
+
273
+
274
+
275
+
276
+
277
277
# ### Save objects ####
278
278
self $ data <- PM_data $ new(data_filtered , quiet = TRUE )
279
279
self $ data $ write(" gendata.csv" , header = FALSE )
280
280
save(self , file = " fit.Rdata" )
281
-
281
+
282
282
# Get ranges and calculate points
283
-
283
+
284
284
ranges <- lapply(self $ model $ model_list $ pri , function (x ) {
285
285
c(x $ min , x $ max )
286
286
})
287
287
names(ranges ) <- tolower(names(ranges ))
288
-
288
+
289
289
# Set initial grid points (only applies for sobol)
290
-
290
+
291
291
vol <- prod(sapply(ranges , function (x ) x [2 ] - x [1 ]))
292
292
points <- ceiling(density0 * vol )
293
-
294
-
295
-
293
+
294
+
295
+
296
296
# set prior
297
- if (prior != " sobol" ){
298
- if (is.numeric(prior )){ # prior specified as a run number
299
- if ( ! file.exists(glue :: glue({prior }," /outputs/theta.csv" ))){
297
+ if (prior != " sobol" ) {
298
+ if (is.numeric(prior )) { # prior specified as a run number
299
+ if (! file.exists(glue :: glue(
300
+ {
301
+ prior
302
+ },
303
+ " /outputs/theta.csv"
304
+ ))) {
300
305
cli :: cli_abort(c(
301
306
" x" = " Error: {.arg prior} file does not exist." ,
302
307
" i" = " Check the file path."
303
308
))
304
- }
305
- file.copy(glue :: glue({prior }," /outputs/theta.csv" ), " theta.csv" )
309
+ }
310
+ file.copy(glue :: glue(
311
+ {
312
+ prior
313
+ },
314
+ " /outputs/theta.csv"
315
+ ), " theta.csv" )
306
316
prior <- " theta.csv"
307
-
308
-
309
- } else if (is.character(prior )) { # prior specified as a filename
317
+ } else if (is.character(prior )) { # prior specified as a filename
310
318
if (! file.exists(prior )) {
311
319
cli :: cli_abort(c(
312
320
" x" = " Error: {.arg prior} file does not exist." ,
313
321
" i" = " Check the file path."
314
322
))
315
323
}
316
- file.copy(prior , overwrite = TRUE ) # ensure in current working directory
317
-
324
+ file.copy(prior , overwrite = TRUE ) # ensure in current working directory
318
325
} else {
319
326
cli :: cli_abort(c(
320
327
" x" = " Error: {.arg prior} must be a numeric run number or character filename." ,
@@ -324,11 +331,11 @@ PM_fit <- R6::R6Class(
324
331
} else {
325
332
prior <- " sobol"
326
333
}
327
-
334
+
328
335
if (intern ) {
329
336
# ## CALL RUST
330
337
out_path <- file.path(getwd(), " outputs" )
331
-
338
+
332
339
rlang :: try_fetch(
333
340
fit(
334
341
self $ model $ binary_path ,
@@ -337,14 +344,14 @@ PM_fit <- R6::R6Class(
337
344
ranges = ranges ,
338
345
algorithm = algorithm ,
339
346
gamlam = c(self $ model $ model_list $ out $ Y1 $ err $ model $ additive , self $ model $ model_list $ out $ Y1 $ err $ model $ proportional ),
340
- error_type = c(" additive" ," proportional" )[1 + is.null(self $ model $ model_list $ out $ Y1 $ err $ model $ additive )],
347
+ error_type = c(" additive" , " proportional" )[1 + is.null(self $ model $ model_list $ out $ Y1 $ err $ model $ additive )],
341
348
error_coefficients = t(sapply(self $ model $ model_list $ out , function (x ) {
342
349
y <- x $ err $ assay $ coefficients
343
- if (length(y ) < 6 ){
344
- y <- c(y ,0 , 0 )
350
+ if (length(y ) < 6 ) {
351
+ y <- c(y , 0 , 0 )
345
352
}
346
- y }
347
- )), # matrix numeqt x 6
353
+ y
354
+ } )), # matrix numeqt x 6
348
355
max_cycles = cycles ,
349
356
prior = prior ,
350
357
ind_points = points ,
@@ -356,7 +363,7 @@ PM_fit <- R6::R6Class(
356
363
return (NULL )
357
364
}
358
365
)
359
-
366
+
360
367
PM_parse(" outputs" )
361
368
res <- PM_load(file = " PMout.Rdata" )
362
369
PM_report(res , outfile = " report.html" , template = report )
@@ -375,7 +382,7 @@ PM_fit <- R6::R6Class(
375
382
save = function (file_name = " PMfit.rds" ) {
376
383
saveRDS(self , file_name )
377
384
},
378
-
385
+
379
386
# ' @description
380
387
# ' `PM_fit` objects contain a `save` method which invokes [saveRDS] to write
381
388
# ' the object to the hard drive as an .rds file. This is the corresponding load
@@ -403,7 +410,6 @@ PM_fit <- R6::R6Class(
403
410
# check if compiled and if not, do so
404
411
self $ model $ compile()
405
412
}
406
-
407
413
) # end private
408
414
) # end PM_fit
409
415
@@ -418,6 +424,6 @@ PM_fit$load <- function(file_name = "PMfit.rds") {
418
424
if (! is.logical(bool )) {
419
425
stop(" This functions expects a logical value" )
420
426
}
421
-
427
+
422
428
rust_logical <- ifelse(bool , " true" , " false" )
423
429
}
0 commit comments