/
foceiFit.R
2932 lines (2850 loc) · 100 KB
/
foceiFit.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
.regFloat1 <- rex::rex(
or(
group(some_of("0":"9"), ".", any_of("0":"9")),
group(any_of("0":"9"), ".", some_of("0":"9"))
),
maybe(group(one_of("E", "e"), maybe(one_of("+", "-")), some_of("0":"9")))
)
.regFloat2 <- rex::rex(some_of("0":"9"), one_of("E", "e"), maybe(one_of("-", "+")), some_of("0":"9"))
.regDecimalint <- rex::rex(or("0", group("1":"9", any_of("0":"9"))))
.regNum <- rex::rex(maybe("-"), or(.regDecimalint, .regFloat1, .regFloat2))
use.utf <- function() {
opt <- getOption("cli.unicode", NULL)
if (!is.null(opt)) {
isTRUE(opt)
} else {
l10n_info()$`UTF-8` && !is.latex()
}
}
is.latex <- function() {
if (!("knitr" %in% loadedNamespaces())) {
return(FALSE)
}
get("is_latex_output", asNamespace("knitr"))()
}
##' Control Options for FOCEi
##'
##' @param sigdig Optimization significant digits. This controls:
##'
##' \itemize{
##'
##' \item The tolerance of the inner and outer optimization is \code{10^-sigdig}
##'
##' \item The tolerance of the ODE solvers is
##' \code{0.5*10^(-sigdig-2)}; For the sensitivity equations and
##' steady-state solutions the default is \code{0.5*10^(-sigdig-1.5)}
##' (sensitivity changes only applicable for liblsoda)
##'
##' \item The tolerance of the boundary check is \code{5 * 10 ^ (-sigdig + 1)}
##'
##' \item The significant figures that some tables are rounded to.
##' }
##'
##' @param atolSens Sensitivity atol, can be different than atol with
##' liblsoda. This allows a less accurate solve for gradients (if desired)
##'
##' @param rtolSens Sensitivity rtol, can be different than rtol with
##' liblsoda. This allows a less accurate solve for gradients (if desired)
##'
##' @param ssAtol Steady state absolute tolerance (atol) for calculating if steady-state
##' has been archived.
##'
##' @param ssRtol Steady state relative tolerance (rtol) for
##' calculating if steady-state has been achieved.
##'
##' @param ssAtolSens Sensitivity absolute tolerance (atol) for
##' calculating if steady state has been achieved for sensitivity compartments.
##'
##' @param ssRtolSens Sensitivity relative tolerance (rtol) for
##' calculating if steady state has been achieved for sensitivity compartments.
##'
##' @param epsilon Precision of estimate for n1qn1 optimization.
##'
##' @param maxstepsOde Maximum number of steps for ODE solver.
##'
##' @param print Integer representing when the outer step is
##' printed. When this is 0 or do not print the iterations. 1 is
##' print every function evaluation (default), 5 is print every 5
##' evaluations.
##'
##' @param scaleTo Scale the initial parameter estimate to this value.
##' By default this is 1. When zero or below, no scaling is performed.
##'
##' @param scaleObjective Scale the initial objective function to this
##' value. By default this is 1.
##'
##' @param derivEps Forward difference tolerances, which is a
##' vector of relative difference and absolute difference. The
##' central/forward difference step size h is calculated as:
##'
##' \code{h = abs(x)*derivEps[1] + derivEps[2]}
##'
##' @param derivMethod indicates the method for calculating
##' derivatives of the outer problem. Currently supports
##' "switch", "central" and "forward" difference methods. Switch
##' starts with forward differences. This will switch to central
##' differences when abs(delta(OFV)) <= derivSwitchTol and switch
##' back to forward differences when abs(delta(OFV)) >
##' derivSwitchTol.
##'
##' @param derivSwitchTol The tolerance to switch forward to central
##' differences.
##'
##' @param covDerivMethod indicates the method for calculating the
##' derivatives while calculating the covariance components
##' (Hessian and S).
##'
##' @param covMethod Method for calculating covariance. In this
##' discussion, R is the Hessian matrix of the objective
##' function. The S matrix is the sum of individual
##' gradient cross-product (evaluated at the individual empirical
##' Bayes estimates).
##'
##' \itemize{
##'
##' \item "\code{r,s}" Uses the sandwich matrix to calculate the
##' covariance, that is: \code{solve(R) \%*\% S \%*\% solve(R)}
##'
##' \item "\code{r}" Uses the Hessian matrix to calculate the
##' covariance as \code{2 \%*\% solve(R)}
##'
##' \item "\code{s}" Uses the cross-product matrix to calculate the
##' covariance as \code{4 \%*\% solve(S)}
##'
##' \item "" Does not calculate the covariance step.
##' }
##'
##' @param covTryHarder If the R matrix is non-positive definite and
##' cannot be corrected to be non-positive definite try estimating
##' the Hessian on the unscaled parameter space.
##'
##' @param hessEps is a double value representing the epsilon for the Hessian calculation.
##'
##' @param centralDerivEps Central difference tolerances. This is a
##' numeric vector of relative difference and absolute difference.
##' The central/forward difference step size h is calculated as:
##'
##' \code{h = abs(x)*derivEps[1] + derivEps[2]}
##'
##' @param lbfgsLmm An integer giving the number of BFGS updates
##' retained in the "L-BFGS-B" method, It defaults to 7.
##'
##' @param lbfgsPgtol is a double precision variable.
##'
##' On entry pgtol >= 0 is specified by the user. The iteration
##' will stop when:
##'
##' \code{max(\| proj g_i \| i = 1, ..., n) <= lbfgsPgtol}
##'
##' where pg_i is the ith component of the projected gradient.
##'
##' On exit pgtol is unchanged. This defaults to zero, when the
##' check is suppressed.
##'
##' @param lbfgsFactr Controls the convergence of the "L-BFGS-B"
##' method. Convergence occurs when the reduction in the
##' objective is within this factor of the machine
##' tolerance. Default is 1e10, which gives a tolerance of about
##' \code{2e-6}, approximately 4 sigdigs. You can check your
##' exact tolerance by multiplying this value by
##' \code{.Machine$double.eps}
##'
##' @param diagXform This is the transformation used on the diagonal
##' of the \code{chol(solve(omega))}. This matrix and values are the
##' parameters estimated in FOCEi. The possibilities are:
##'
##' \itemize{
##' \item \code{sqrt} Estimates the sqrt of the diagonal elements of \code{chol(solve(omega))}. This is the default method.
##'
##' \item \code{log} Estimates the log of the diagonal elements of \code{chol(solve(omega))}
##'
##' \item \code{identity} Estimates the diagonal elements without any transformations
##' }
##' @param sumProd Is a boolean indicating if the model should change
##' multiplication to high precision multiplication and sums to
##' high precision sums using the PreciseSums package. By default
##' this is \code{FALSE}.
##'
##'
##' @param optExpression Optimize the RxODE expression to speed up
##' calculation. By default this is turned on.
##'
##' @param ci Confidence level for some tables. By default this is
##' 0.95 or 95\% confidence.
##'
##' @param useColor Boolean indicating if focei can use ASCII color codes
##'
##' @param boundTol Tolerance for boundary issues.
##'
##' @param calcTables This boolean is to determine if the foceiFit
##' will calculate tables. By default this is \code{TRUE}
##'
##' @param ... Ignored parameters
##'
##' @param maxInnerIterations Number of iterations for n1qn1
##' optimization.
##'
##' @param maxOuterIterations Maximum number of L-BFGS-B optimization
##' for outer problem.
##'
##' @param n1qn1nsim Number of function evaluations for n1qn1
##' optimization.
##'
##' @param eigen A boolean indicating if eigenvectors are calculated
##' to include a condition number calculation.
##'
##' @param addPosthoc Boolean indicating if posthoc parameters are
##' added to the table output.
##'
##' @param printNcol Number of columns to printout before wrapping
##' parameter estimates/gradient
##'
##' @param noAbort Boolean to indicate if you should abort the FOCEi
##' evaluation if it runs into troubles. (default TRUE)
##'
##' @param interaction Boolean indicate FOCEi should be used (TRUE)
##' instead of FOCE (FALSE)
##'
##' @param cholSEOpt Boolean indicating if the generalized Cholesky
##' should be used while optimizing.
##'
##' @param cholSECov Boolean indicating if the generalized Cholesky
##' should be used while calculating the Covariance Matrix.
##'
##' @param fo is a boolean indicating if this is a FO approximation routine.
##'
##' @param cholSEtol tolerance for Generalized Cholesky
##' Decomposition. Defaults to suggested (.Machine$double.eps)^(1/3)
##'
##' @param cholAccept Tolerance to accept a Generalized Cholesky
##' Decomposition for a R or S matrix.
##'
##' @param outerOpt optimization method for the outer problem
##'
##' @param innerOpt optimization method for the inner problem (not
##' implemented yet.)
##'
##' @param stateTrim Trim state amounts/concentrations to this value.
##'
##' @param resetEtaP represents the p-value for reseting the
##' individual ETA to 0 during optimization (instead of the saved
##' value). The two test statistics used in the z-test are either
##' chol(omega^-1) \%*\% eta or eta/sd(allEtas). A p-value of 0
##' indicates the ETAs never reset. A p-value of 1 indicates the
##' ETAs always reset.
##'
##' @param resetThetaP represents the p-value for reseting the
##' population mu-referenced THETA parameters based on ETA drift
##' during optimization, and resetting the optimization. A
##' p-value of 0 indicates the THETAs never reset. A p-value of 1
##' indicates the THETAs always reset and is not allowed. The
##' theta reset is checked at the beginning and when nearing a
##' local minima. The percent change in objective function where
##' a theta reset check is initiated is controlled in
##' \code{resetThetaCheckPer}.
##'
##' @param resetThetaCheckPer represents objective function
##' \% percentage below which resetThetaP is checked.
##'
##' @param resetThetaFinalP represents the p-value for reseting the
##' population mu-referenced THETA parameters based on ETA drift
##' during optimization, and resetting the optimization one final time.
##'
##' @param resetHessianAndEta is a boolean representing if the
##' individual Hessian is reset when ETAs are reset using the
##' option \code{resetEtaP}.
##'
##' @param diagOmegaBoundUpper This represents the upper bound of the
##' diagonal omega matrix. The upper bound is given by
##' diag(omega)*diagOmegaBoundUpper. If
##' \code{diagOmegaBoundUpper} is 1, there is no upper bound on
##' Omega.
##'
##' @param diagOmegaBoundLower This represents the lower bound of the
##' diagonal omega matrix. The lower bound is given by
##' diag(omega)/diagOmegaBoundUpper. If
##' \code{diagOmegaBoundLower} is 1, there is no lower bound on
##' Omega.
##'
##' @param rhobeg Beginning change in parameters for bobyqa algorithm
##' (trust region). By default this is 0.2 or 20% of the initial
##' parameters when the parameters are scaled to 1. rhobeg and
##' rhoend must be set to the initial and final values of a trust
##' region radius, so both must be positive with 0 < rhoend <
##' rhobeg. Typically rhobeg should be about one tenth of the
##' greatest expected change to a variable. Note also that
##' smallest difference abs(upper-lower) should be greater than or
##' equal to rhobeg*2. If this is not the case then rhobeg will be
##' adjusted.
##'
##' @param rhoend The smallest value of the trust region radius that
##' is allowed. If not defined, then 10^(-sigdig-1) will be used.
##'
##' @param npt The number of points used to approximate the objective
##' function via a quadratic approximation for bobyqa. The value
##' of npt must be in the interval [n+2,(n+1)(n+2)/2] where n is
##' the number of parameters in par. Choices that exceed 2*n+1 are
##' not recommended. If not defined, it will be set to 2*n + 1
##' @param eval.max Number of maximum evaluations of the objective function
##'
##' @param iter.max Maximum number of iterations allowed.
##'
##' @param rel.tol Relative tolerance before nlminb stops.
##'
##' @param x.tol X tolerance for nlmixr optimizers
##'
##' @param abstol Absolute tolerance for nlmixr optimizer
##'
##' @param reltol tolerance for nlmixr
##'
##' @param gillK The total number of possible steps to determine the
##' optimal forward/central difference step size per parameter (by
##' the Gill 1983 method). If 0, no optimal step size is
##' determined. Otherwise this is the optimal step size
##' determined.
##'
##' @param gillRtol The relative tolerance used for Gill 1983
##' determination of optimal step size.
##'
##' @param scaleType The scaling scheme for nlmixr. The supported types are:
##'
##' \itemize{
##' \item \code{nlmixr} In this approach the scaling is performed by the following equation:
##'
##' v_{scaled} = (v_{current} - v_{init})/scaleC[i] + scaleTo
##'
##' The \code{scaleTo} parameter is specified by the \code{normType},
##' and the scales are specified by \code{scaleC}.
##'
##' \item \code{norm} This approach uses the simple scaling provided
##' by the \code{normType} argument.
##'
##' \item \code{mult} This approach does not use the data
##' normalization provided by \code{normType}, but rather uses
##' multiplicative scaling to a constant provided by the \code{scaleTo}
##' argument.
##'
##' In this case:
##'
##' v_{scaled} = v_{current}/v_{init}*scaleTo
##'
##' \item \code{multAdd} This approach changes the scaling based on
##' the parameter being specified. If a parameter is defined in an
##' exponential block (ie exp(theta)), then it is scaled on a
##' linearly, that is:
##'
##' v_{scaled} = (v_{current}-v_{init}) + scaleTo
##'
##' Otherwise the parameter is scaled multiplicatively.
##'
##' v_{scaled} = v_{current}/v_{init}*scaleTo
##'
##' }
##'
##' @param scaleC The scaling constant used with
##' \code{scaleType=nlmixr}. When not specified, it is based on
##' the type of parameter that is estimated. The idea is to keep
##' the derivatives similar on a log scale to have similar
##' gradient sizes. Hence parameters like log(exp(theta)) would
##' have a scaling factor of 1 and log(theta) would have a scaling
##' factor of ini_value (to scale by 1/value; ie
##' d/dt(log(ini_value)) = 1/ini_value or scaleC=ini_value)
##'
##' \itemize{
##'
##' \item For parameters in an exponential (ie exp(theta)) or
##' parameters specifying powers, boxCox or yeoJohnson
##' transformations , this is 1.
##'
##' \item For additive, proportional, lognormal error structures,
##' these are given by 0.5*abs(initial_estimate)
##'
##' \item Factorials are scaled by abs(1/digamma(inital_estimate+1))
##'
##' \item parameters in a log scale (ie log(theta)) are transformed
##' by log(abs(initial_estimate))*abs(initial_estimate)
##'
##' }
##'
##' These parameter scaling coefficients are chose to try to keep
##' similar slopes among parameters. That is they all follow the
##' slopes approximately on a log-scale.
##'
##' While these are chosen in a logical manner, they may not always
##' apply. You can specify each parameters scaling factor by this
##' parameter if you wish.
##'
##' @param scaleC0 Number to adjust the scaling factor by if the initial
##' gradient is zero.
##'
##' @param scaleCmax Maximum value of the scaleC to prevent overflow.
##'
##' @param scaleCmin Minimum value of the scaleC to prevent underflow.
##'
##' @param normType This is the type of parameter
##' normalization/scaling used to get the scaled initial values
##' for nlmixr. These are used with \code{scaleType} of.
##'
##' With the exception of \code{rescale2}, these come
##' from
##' \href{https://en.wikipedia.org/wiki/Feature_scaling}{Feature
##' Scaling}. The \code{rescale2} The rescaling is the same type
##' described in the
##' \href{http://apmonitor.com/me575/uploads/Main/optimization_book.pdf}{OptdesX}
##' software manual.
##'
##' In general, all all scaling formula can be described by:
##'
##' v_{scaled} = (v_{unscaled}-C_{1})/C_{2}
##'
##' Where
##'
##'
##' The other data normalization approaches follow the following formula
##'
##' v_{scaled} = (v_{unscaled}-C_{1})/C_{2};
##'
##' \itemize{
##'
##' \item \code{rescale2} This scales all parameters from (-1 to 1).
##' The relative differences between the parameters are preserved
##' with this approach and the constants are:
##'
##' C_{1} = (max(all unscaled values)+min(all unscaled values))/2
##'
##' C_{2} = (max(all unscaled values) - min(all unscaled values))/2
##'
##'
##' \item \code{rescale} or min-max normalization. This rescales all
##' parameters from (0 to 1). As in the \code{rescale2} the
##' relative differences are preserved. In this approach:
##'
##' C_{1} = min(all unscaled values)
##'
##' C_{2} = max(all unscaled values) - min(all unscaled values)
##'
##'
##' \item \code{mean} or mean normalization. This rescales to center
##' the parameters around the mean but the parameters are from 0
##' to 1. In this approach:
##'
##' C_{1} = mean(all unscaled values)
##'
##' C_{2} = max(all unscaled values) - min(all unscaled values)
##'
##' \item \code{std} or standardization. This standardizes by the mean
##' and standard deviation. In this approach:
##'
##' C_{1} = mean(all unscaled values)
##'
##' C_{2} = sd(all unscaled values)
##'
##' \item \code{len} or unit length scaling. This scales the
##' parameters to the unit length. For this approach we use the Euclidean length, that
##' is:
##'
##' C_{1} = 0
##'
##' C_{2} = sqrt(v_1^2 + v_2^2 + ... + v_n^2)
##'
##'
##' \item \code{constant} which does not perform data normalization. That is
##'
##' C_{1} = 0
##'
##' C_{2} = 1
##'
##' }
##'
##' @param gillStep When looking for the optimal forward difference
##' step size, this is This is the step size to increase the
##' initial estimate by. So each iteration the new step size =
##' (prior step size)*gillStep
##'
##' @param gillFtol The gillFtol is the gradient error tolerance that
##' is acceptable before issuing a warning/error about the gradient estimates.
##'
##' @param gillKcov The total number of possible steps to determine
##' the optimal forward/central difference step size per parameter
##' (by the Gill 1983 method) during the covariance step. If 0,
##' no optimal step size is determined. Otherwise this is the
##' optimal step size determined.
##'
##' @param gillStepCov When looking for the optimal forward difference
##' step size, this is This is the step size to increase the
##' initial estimate by. So each iteration during the covariance
##' step is equal to the new step size = (prior step size)*gillStepCov
##'
##' @param gillFtolCov The gillFtol is the gradient error tolerance
##' that is acceptable before issuing a warning/error about the
##' gradient estimates during the covariance step.
##'
##' @param rmatNorm A parameter to normalize gradient step size by the
##' parameter value during the calculation of the R matrix
##'
##' @param smatNorm A parameter to normalize gradient step size by the
##' parameter value during the calculation of the S matrix
##'
##' @param covGillF Use the Gill calculated optimal Forward difference
##' step size for the instead of the central difference step size
##' during the central difference gradient calculation.
##'
##' @param optGillF Use the Gill calculated optimal Forward difference
##' step size for the instead of the central difference step size
##' during the central differences for optimization.
##'
##' @param covSmall The covSmall is the small number to compare
##' covariance numbers before rejecting an estimate of the
##' covariance as the final estimate (when comparing sandwich vs
##' R/S matrix estimates of the covariance). This number controls
##' how small the variance is before the covariance matrix is
##' rejected.
##'
##' @param adjLik In nlmixr, the objective function matches NONMEM's
##' objective function, which removes a 2*pi constant from the
##' likelihood calculation. If this is TRUE, the likelihood
##' function is adjusted by this 2*pi factor. When adjusted this
##' number more closely matches the likelihood approximations of
##' nlme, and SAS approximations. Regardless of if this is turned
##' on or off the objective function matches NONMEM's objective
##' function.
##'
##' @param gradTrim The parameter to adjust the gradient to if the
##' |gradient| is very large.
##'
##' @param gradCalcCentralSmall A small number that represents the value
##' where |grad| < gradCalcCentralSmall where forward differences
##' switch to central differences.
##'
##' @param gradCalcCentralLarge A large number that represents the value
##' where |grad| > gradCalcCentralLarge where forward differences
##' switch to central differences.
##'
##' @param etaNudge By default initial ETA estimates start at zero;
##' Sometimes this doesn't optimize appropriately. If this value is
##' non-zero, when the n1qn1 optimization didn't perform
##' appropriately, reset the Hessian, and nudge the ETA up by this
##' value; If the ETA still doesn't move, nudge the ETA down by this
##' value. By default this value is qnorm(1-0.05/2)*1/sqrt(3), the
##' first of the Gauss Quadrature numbers times by the 0.95\% normal
##' region. If this is not successful try the second eta nudge
##' number (below). If +-etaNudge2 is not successful, then assign
##' to zero and do not optimize any longer
##'
##' @param etaNudge2 This is the second eta nudge. By default it is
##' qnorm(1-0.05/2)*sqrt(3/5), which is the n=3 quadrature point
##' (excluding zero) times by the 0.95\% normal region
##'
##' @param maxOdeRecalc Maximum number of times to reduce the ODE
##' tolerances and try to resolve the system if there was a bad
##' ODE solve.
##'
##' @param odeRecalcFactor The factor to increase the rtol/atol with
##' bad ODE solving.
##'
##' @param repeatGillMax If the tolerances were reduced when
##' calculating the initial Gill differences, the Gill difference
##' is repeated up to a maximum number of times defined by this
##' parameter.
##'
##' @param stickyRecalcN The number of bad ODE solves before reducing
##' the atol/rtol for the rest of the problem.
##'
##' @param nRetries If FOCEi doesn't fit with the current parameter
##' estimates, randomly sample new parameter estimates and restart
##' the problem. This is similar to 'PsN' resampling.
##'
##' @param eventFD Finite difference step for forward or central
##' difference estimation of event-based gradients
##'
##' @param eventType Event gradient type for dosing events; Can be
##' "gill", "central" or "forward"
##'
##' @param gradProgressOfvTime This is the time for a single objective
##' function evaluation (in seconds) to start progress bars on gradient evaluations
##'
##' @param singleOde This option allows a single ode model to include
##' the PK parameter information instead of splitting it into a
##' function and a RxODE model
##'
##' @param badSolveObjfAdj The objective function adjustment when the
##' ODE system cannot be solved. It is based on each individual bad
##' solve.
##'
##' @inheritParams configsaem
##' @inheritParams RxODE::rxSolve
##' @inheritParams minqa::bobyqa
##' @inheritParams foceiFit
##'
##' @details
##'
##' Note this uses the R's L-BFGS-B in \code{\link{optim}} for the
##' outer problem and the BFGS \code{\link[n1qn1]{n1qn1}} with that
##' allows restoring the prior individual Hessian (for faster
##' optimization speed).
##'
##' However the inner problem is not scaled. Since most eta estimates
##' start near zero, scaling for these parameters do not make sense.
##'
##' This process of scaling can fix some ill conditioning for the
##' unscaled problem. The covariance step is performed on the
##' unscaled problem, so the condition number of that matrix may not
##' be reflective of the scaled problem's condition-number.
##'
##' @author Matthew L. Fidler
##'
##' @return The control object that changes the options for the FOCEi
##' family of estimation methods
##'
##' @seealso \code{\link{optim}}
##' @seealso \code{\link[n1qn1]{n1qn1}}
##' @seealso \code{\link[RxODE]{rxSolve}}
##' @export
foceiControl <- function(sigdig = 3, ...,
epsilon = NULL, # 1e-4,
maxInnerIterations = 1000,
maxOuterIterations = 5000,
n1qn1nsim = NULL,
method = c("liblsoda", "lsoda", "dop853"),
transitAbs = NULL, atol = NULL, rtol = NULL,
atolSens = NULL, rtolSens = NULL,
ssAtol = NULL, ssRtol = NULL, ssAtolSens = NULL, ssRtolSens = NULL,
minSS = 10L, maxSS = 1000L,
maxstepsOde = 500000L, hmin = 0L, hmax = NA_real_, hini = 0, maxordn = 12L, maxords = 5L, cores,
covsInterpolation = c("locf", "linear", "nocb", "midpoint"),
print = 1L,
printNcol = floor((getOption("width") - 23) / 12),
scaleTo = 1.0,
scaleObjective = 0,
normType = c("rescale2", "mean", "rescale", "std", "len", "constant"),
scaleType = c("nlmixr", "norm", "mult", "multAdd"),
scaleCmax = 1e5,
scaleCmin = 1e-5,
scaleC = NULL,
scaleC0 = 1e5,
derivEps = rep(20 * sqrt(.Machine$double.eps), 2),
derivMethod = c("switch", "forward", "central"),
derivSwitchTol = NULL,
covDerivMethod = c("central", "forward"),
covMethod = c("r,s", "r", "s", ""),
hessEps = (.Machine$double.eps)^(1 / 3),
eventFD = sqrt(.Machine$double.eps),
eventType = c("gill", "central", "forward"),
centralDerivEps = rep(20 * sqrt(.Machine$double.eps), 2),
lbfgsLmm = 7L,
lbfgsPgtol = 0,
lbfgsFactr = NULL,
eigen = TRUE,
addPosthoc = TRUE,
diagXform = c("sqrt", "log", "identity"),
sumProd = FALSE,
optExpression = TRUE,
ci = 0.95,
useColor = crayon::has_color(),
boundTol = NULL,
calcTables = TRUE,
noAbort = TRUE,
interaction = TRUE,
cholSEtol = (.Machine$double.eps)^(1 / 3),
cholAccept = 1e-3,
resetEtaP = 0.15,
resetThetaP = 0.05,
resetThetaFinalP = 0.15,
diagOmegaBoundUpper = 5, # diag(omega) = diag(omega)*diagOmegaBoundUpper; =1 no upper
diagOmegaBoundLower = 100, # diag(omega) = diag(omega)/diagOmegaBoundLower; = 1 no lower
cholSEOpt = FALSE,
cholSECov = FALSE,
fo = FALSE,
covTryHarder = FALSE,
## Ranking based on run 025
## L-BFGS-B: 20970.53 (2094.004 429.535)
## bobyqa: 21082.34 (338.677 420.754)
## lbfgsb3* (modified for tolerances):
## nlminb: 20973.468 (755.821 458.343)
## mma: 20974.20 (Time: Opt: 3000.501 Cov: 467.287)
## slsqp: 21023.89 (Time: Opt: 460.099; Cov: 488.921)
## lbfgsbLG: 20974.74 (Time: Opt: 946.463; Cov:397.537)
outerOpt = c("nlminb", "bobyqa", "lbfgsb3c", "L-BFGS-B", "mma", "lbfgsbLG", "slsqp", "Rvmmin"),
innerOpt = c("n1qn1", "BFGS"),
##
rhobeg = .2,
rhoend = NULL,
npt = NULL,
## nlminb
rel.tol = NULL,
x.tol = NULL,
eval.max = 4000,
iter.max = 2000,
abstol = NULL,
reltol = NULL,
resetHessianAndEta = FALSE,
stateTrim = Inf,
gillK = 10L,
gillStep = 4,
gillFtol = 0,
gillRtol = sqrt(.Machine$double.eps),
gillKcov = 10L,
gillStepCov = 2,
gillFtolCov = 0,
rmatNorm = TRUE,
smatNorm = TRUE,
covGillF = TRUE,
optGillF = TRUE,
covSmall = 1e-5,
adjLik = TRUE, ## Adjust likelihood by 2pi for FOCEi methods
gradTrim = Inf,
maxOdeRecalc = 5,
odeRecalcFactor = 10^(0.5),
gradCalcCentralSmall = 1e-4,
gradCalcCentralLarge = 1e4,
etaNudge = qnorm(1-0.05/2)/sqrt(3),
etaNudge2=qnorm(1-0.05/2) * sqrt(3/5),
stiff,
nRetries = 3,
seed = 42,
resetThetaCheckPer = 0.1,
etaMat = NULL,
repeatGillMax = 3,
stickyRecalcN = 5,
gradProgressOfvTime = 10,
addProp = c("combined2", "combined1"),
singleOde = TRUE,
badSolveObjfAdj=100) {
if (is.null(boundTol)) {
boundTol <- 5 * 10^(-sigdig + 1)
}
if (is.null(epsilon)) {
epsilon <- 10^(-sigdig - 1)
}
if (is.null(abstol)) {
abstol <- 10^(-sigdig - 1)
}
if (is.null(reltol)) {
reltol <- 10^(-sigdig - 1)
}
if (is.null(rhoend)) {
rhoend <- 10^(-sigdig - 1)
}
if (is.null(lbfgsFactr)) {
lbfgsFactr <- 10^(-sigdig - 1) / .Machine$double.eps
}
if (is.null(atol)) {
atol <- 0.5 * 10^(-sigdig - 2)
}
if (is.null(rtol)) {
rtol <- 0.5 * 10^(-sigdig - 2)
}
if (is.null(atolSens)) {
atolSens <- 0.5 * 10^(-sigdig - 1.5)
}
if (is.null(rtolSens)) {
rtolSens <- 0.5 * 10^(-sigdig - 1.5)
}
if (is.null(ssAtol)) {
ssAtol <- 0.5 * 10^(-sigdig - 2)
}
if (is.null(ssRtol)) {
ssRtol <- 0.5 * 10^(-sigdig - 2)
}
if (is.null(ssAtolSens)) {
ssAtolSens <- 0.5 * 10^(-sigdig - 1.5)
}
if (is.null(ssRtolSens)) {
ssRtolSens <- 0.5 * 10^(-sigdig - 1.5)
}
if (is.null(rel.tol)) {
rel.tol <- 10^(-sigdig - 1)
}
if (is.null(x.tol)) {
x.tol <- 10^(-sigdig - 1)
}
if (is.null(derivSwitchTol)) {
derivSwitchTol <- 2 * 10^(-sigdig - 1)
}
## if (is.null(gillRtol)){
## ## FIXME: there is a way to calculate this according to the
## ## Gill paper but it is buried in their optimization book.
## gillRtol <- 10 ^ (-sigdig - 1);
## }
.xtra <- list(...)
if (is.null(transitAbs) && !is.null(.xtra$transit_abs)) { # nolint
transitAbs <- .xtra$transit_abs # nolint
}
if (missing(covsInterpolation) && !is.null(.xtra$covs_interpolation)) { # nolint
covsInterpolation <- .xtra$covs_interpolation # nolint
}
if (missing(maxInnerIterations) && !is.null(.xtra$max_iterations)) { # nolint
maxInnerIterations <- .xtra$max_iterations # nolint
}
if (!missing(stiff) && missing(method)) {
if (RxODE::rxIs(stiff, "logical")) {
if (stiff) {
method <- "lsoda"
warning("stiff=TRUE has been replaced with method = \"lsoda\".")
} else {
method <- "dop853"
warning("stiff=FALSE has been replaced with method = \"dop853\".")
}
}
} else {
if (inherits(method, "numeric")) {
method <- as.integer(method)
}
if (!RxODE::rxIs(method, "integer")) {
if (inherits(method, "character")) {
method <- match.arg(method)
} else {
method <- "liblsoda"
warning("could not figure out method, using 'liblsoda'")
}
}
}
## .methodIdx <- c("lsoda"=1L, "dop853"=0L, "liblsoda"=2L);
## method <- as.integer(.methodIdx[method]);
if (RxODE::rxIs(scaleType, "character")) {
.scaleTypeIdx <- c("norm" = 1L, "nlmixr" = 2L, "mult" = 3L, "multAdd" = 4L)
scaleType <- as.integer(.scaleTypeIdx[match.arg(scaleType)])
}
if (RxODE::rxIs(eventType, "character")) {
.eventTypeIdx <- c("gill" = 1L, "central" = 2L, "forward" = 3L)
eventType <- as.integer(.eventTypeIdx[match.arg(eventType)])
}
if (RxODE::rxIs(normType, "character")) {
.normTypeIdx <- c("rescale2" = 1L, "rescale" = 2L, "mean" = 3L, "std" = 4L, "len" = 5L, "constant" = 6)
normType <- as.integer(.normTypeIdx[match.arg(normType)])
}
derivMethod <- match.arg(derivMethod)
.methodIdx <- c("forward" = 0L, "central" = 1L, "switch" = 3L)
derivMethod <- as.integer(.methodIdx[derivMethod])
covDerivMethod <- .methodIdx[match.arg(covDerivMethod)]
if (length(covsInterpolation) > 1) covsInterpolation <- covsInterpolation[1]
if (!RxODE::rxIs(covsInterpolation, "integer")) {
covsInterpolation <- tolower(match.arg(
covsInterpolation,
c("linear", "locf", "LOCF", "constant", "nocb", "NOCB", "midpoint")
))
}
## if (covsInterpolation == "constant") covsInterpolation <- "locf";
## covsInterpolation <- as.integer(which(covsInterpolation == c("linear", "locf", "nocb", "midpoint")) - 1);
if (missing(cores)) {
cores <- RxODE::rxCores()
}
if (missing(n1qn1nsim)) {
n1qn1nsim <- 10 * maxInnerIterations + 1
}
if (length(covMethod) == 1) {
if (covMethod == "") {
covMethod <- 0L
}
}
if (RxODE::rxIs(covMethod, "character")) {
covMethod <- match.arg(covMethod)
.covMethodIdx <- c("r,s" = 1L, "r" = 2L, "s" = 3L)
covMethod <- .covMethodIdx[match.arg(covMethod)]
}
.outerOptTxt <- "custom"
if (RxODE::rxIs(outerOpt, "character")) {
outerOpt <- match.arg(outerOpt)
.outerOptTxt <- outerOpt
if (outerOpt == "bobyqa") {
RxODE::rxReq("minqa")
outerOptFun <- .bobyqa
outerOpt <- -1L
} else if (outerOpt == "nlminb") {
outerOptFun <- .nlminb
outerOpt <- -1L
} else if (outerOpt == "mma") {
outerOptFun <- .nloptr
outerOpt <- -1L
} else if (outerOpt == "slsqp") {
outerOptFun <- .slsqp
outerOpt <- -1L
} else if (outerOpt == "lbfgsbLG") {
outerOptFun <- .lbfgsbLG
outerOpt <- -1L
} else if (outerOpt == "Rvmmin") {
outerOptFun <- .Rvmmin
outerOpt <- -1L
} else {
.outerOptIdx <- c("L-BFGS-B" = 0L, "lbfgsb3c" = 1L)
outerOpt <- .outerOptIdx[outerOpt]
if (outerOpt == 1L) {
RxODE::rxReq("lbfgsb3c")
}
outerOptFun <- NULL
}
} else if (is(outerOpt, "function")) {
outerOptFun <- outerOpt
outerOpt <- -1L
}
if (RxODE::rxIs(innerOpt, "character")) {
.innerOptFun <- c("n1qn1" = 1L, "BFGS" = 2L)
innerOpt <- setNames(.innerOptFun[match.arg(innerOpt)], NULL)
}
if (resetEtaP > 0 & resetEtaP < 1) {
.resetEtaSize <- qnorm(1 - (resetEtaP / 2))
} else if (resetEtaP <= 0) {
.resetEtaSize <- Inf
} else {
.resetEtaSize <- 0
}
if (resetThetaP > 0 & resetThetaP < 1) {
.resetThetaSize <- qnorm(1 - (resetThetaP / 2))
} else if (resetThetaP <= 0) {
.resetThetaSize <- Inf
} else {
stop("Cannot always reset THETAs")
}
if (resetThetaFinalP > 0 & resetThetaFinalP < 1) {
.resetThetaFinalSize <- qnorm(1 - (resetThetaFinalP / 2))
} else if (resetThetaP <= 0) {
.resetThetaFinalSize <- Inf
} else {
stop("Cannot always reset THETAs")
}
if (inherits(addProp, "numeric")) {
if (addProp == 1) {
addProp <- "combined1"
} else if (addProp == 2) {
addProp <- "combined2"
} else {
stop("addProp must be 1, 2, \"combined1\" or \"combined2\"", call.=FALSE)
}
} else {
addProp <- match.arg(addProp)
}
.ret <- list(
maxOuterIterations = as.integer(maxOuterIterations),
maxInnerIterations = as.integer(maxInnerIterations),
method = method,
transitAbs = transitAbs,
atol = atol,
rtol = rtol,
atolSens = atolSens,
rtolSens = rtolSens,
ssAtol = ssAtol,
ssRtol = ssRtol,
ssAtolSens = ssAtolSens,
ssRtolSens = ssRtolSens,
minSS = minSS, maxSS = maxSS,
maxstepsOde = maxstepsOde,
hmin = hmin,
hmax = hmax,
hini = hini,
maxordn = maxordn,
maxords = maxords,
cores = cores,
covsInterpolation = covsInterpolation,
n1qn1nsim = as.integer(n1qn1nsim),
print = as.integer(print),
lbfgsLmm = as.integer(lbfgsLmm),
lbfgsPgtol = as.double(lbfgsPgtol),
lbfgsFactr = as.double(lbfgsFactr),
scaleTo = scaleTo,
epsilon = epsilon,
derivEps = derivEps,
derivMethod = derivMethod,
covDerivMethod = covDerivMethod,
covMethod = covMethod,
centralDerivEps = centralDerivEps,
eigen = as.integer(eigen),
addPosthoc = as.integer(addPosthoc),
diagXform = match.arg(diagXform),
sumProd = sumProd,
optExpression = optExpression,
outerOpt = as.integer(outerOpt),
ci = as.double(ci),
sigdig = as.double(sigdig),
scaleObjective = as.double(scaleObjective),
useColor = as.integer(useColor),
boundTol = as.double(boundTol),
calcTables = calcTables,
printNcol = as.integer(printNcol),
noAbort = as.integer(noAbort),
interaction = as.integer(interaction),
cholSEtol = as.double(cholSEtol),
hessEps = as.double(hessEps),
cholAccept = as.double(cholAccept),
resetEtaSize = as.double(.resetEtaSize),
resetThetaSize = as.double(.resetThetaSize),
resetThetaFinalSize = as.double(.resetThetaFinalSize),
diagOmegaBoundUpper = diagOmegaBoundUpper,
diagOmegaBoundLower = diagOmegaBoundLower,
cholSEOpt = as.integer(cholSEOpt),
cholSECov = as.integer(cholSECov),
fo = as.integer(fo),
covTryHarder = as.integer(covTryHarder),
outerOptFun = outerOptFun,
## bobyqa
rhobeg = as.double(rhobeg),
rhoend = as.double(rhoend),
npt = npt,
## nlminb
rel.tol = as.double(rel.tol),
x.tol = as.double(x.tol),
eval.max = eval.max,
iter.max = iter.max,
innerOpt = innerOpt,
## BFGS
abstol = abstol,
reltol = reltol,
derivSwitchTol = derivSwitchTol,
resetHessianAndEta = as.integer(resetHessianAndEta),
stateTrim = as.double(stateTrim),
gillK = as.integer(gillK),