Skip to content

Commit

Permalink
Merge pull request #486 from nlmixr2/485-nls
Browse files Browse the repository at this point in the history 
Scaling now falls back to unscaled
  • Loading branch information
@mattfidler
mattfidler authored Oct 16, 2024
2 parents 562dd28 + ae64f2e commit d33f9e4
Show file tree
Hide file tree
Showing 2 changed files with 84 additions and 26 deletions.
88 changes: 62 additions & 26 deletions src/scale.h
View file
Original file line number Diff line number Diff line change
@@ -1,3 +1,4 @@
#include <math.h>
#define min2( a , b ) ( (a) < (b) ? (a) : (b) )
#define max2( a , b ) ( (a) > (b) ? (a) : (b) )
#define expit(alpha, low, high) _powerDi(alpha, 1.0, 4, low, high)
Expand Down Expand Up @@ -123,14 +124,21 @@ static inline void scaleSetup(scaling *scale,
mx = max2(scale->initPar[k],mx);
scale->scaleC[k] = NA_REAL;
}
if (mx == mn) {
warning(_("all parameters are the same value, switch to length normType"));
if (fabs(mx-mn) < DBL_EPSILON) {
for (unsigned int k = scale->npars-1; k--;){
len += scale->initPar[k]*scale->initPar[k];
}
scale->c1 = 0;
scale->c2 = sqrt(len);
scale->normType = 5;
if (len < DBL_EPSILON){
warning(_("all parameters are zero, cannot scale, run unscaled"));
scale->c1 = 0;
scale->c2 = 1;
scale->normType = normTypeConstant;
} else {
warning(_("all parameters are the same value, switch to length normType"));
scale->c1 = 0;
scale->c2 = sqrt(len);
scale->normType = 5;
}
} else {
scale->c1 = (mx+mn)/2;
scale->c2 = (mx-mn)/2;
Expand All @@ -142,14 +150,21 @@ static inline void scaleSetup(scaling *scale,
mx = max2(scale->initPar[k],mx);
scale->scaleC[k] = NA_REAL;
}
if (mx == mn) {
warning(_("all parameters are the same value, switch to length normType"));
if (fabs(mx-mn) < DBL_EPSILON) {
for (unsigned int k = scale->npars-1; k--;){
len += scale->initPar[k]*scale->initPar[k];
}
scale->c1 = 0;
scale->c2 = sqrt(len);
scale->normType = 5;
if (len < DBL_EPSILON){
warning(_("all parameters are zero, cannot scale, run unscaled"));
scale->c1 = 0;
scale->c2 = 1;
scale->normType = normTypeConstant;
} else {
warning(_("all parameters are the same value, switch to length normType"));
scale->c1 = 0;
scale->c2 = sqrt(len);
scale->normType = 5;
}
} else {
scale->c1 = mn;
scale->c2 = (mx-mn);
Expand All @@ -163,14 +178,21 @@ static inline void scaleSetup(scaling *scale,
mean += (scale->initPar[k]-mean)/oN;
scale->scaleC[k] = NA_REAL;
}
if (mx == mn) {
warning(_("all parameters are the same value, switch to length normType"));
if (fabs(mx-mn) < DBL_EPSILON) {
for (unsigned int k = scale->npars-1; k--;){
len += scale->initPar[k]*scale->initPar[k];
}
scale->c1 = 0;
scale->c2 = sqrt(len);
scale->normType = 5;
if (len < DBL_EPSILON){
warning(_("all parameters are zero, cannot scale, run unscaled"));
scale->c1 = 0;
scale->c2 = 1;
scale->normType = normTypeConstant;
} else {
warning(_("all parameters are the same value, switch to length normType"));
scale->c1 = 0;
scale->c2 = sqrt(len);
scale->normType = 5;
}
} else {
scale->c1 = mean;
scale->c2 = (mx-mn);
Expand All @@ -186,14 +208,21 @@ static inline void scaleSetup(scaling *scale,
s += (scale->initPar[k]-mean)*(scale->initPar[k]-oM);
scale->scaleC[k] = NA_REAL;
}
if (mx == mn) {
warning("all parameters are the same value, switch to length norm type");
if (fabs(mx-mn) < DBL_EPSILON) {
for (unsigned int k = scale->npars-1; k--;){
len += scale->initPar[k]*scale->initPar[k];
}
scale->c1 = 0;
scale->c2 = sqrt(len);
scale->normType = 5;
if (len < DBL_EPSILON){
warning(_("all parameters are zero, cannot scale, run unscaled"));
scale->c1 = 0;
scale->c2 = 1;
scale->normType = normTypeConstant;
} else {
warning(_("all parameters are the same value, switch to length normType"));
scale->c1 = 0;
scale->c2 = sqrt(len);
scale->normType = 5;
}
} else {
scale->c1 = mean;
scale->c2 = sqrt(s/(oN-1));
Expand All @@ -204,8 +233,15 @@ static inline void scaleSetup(scaling *scale,
len += scale->initPar[k]*scale->initPar[k];
scale->scaleC[k] = NA_REAL;
}
scale->c1 = 0;
scale->c2 = sqrt(len);
if (len < DBL_EPSILON){
warning(_("all parameters are zero, cannot scale, run unscaled"));
scale->c1 = 0;
scale->c2 = 1;
scale->normType = normTypeConstant;
} else {
scale->c1 = 0;
scale->c2 = sqrt(len);
}
break;
case normTypeConstant:
// No Normalization
Expand All @@ -222,23 +258,23 @@ static inline void scaleSetup(scaling *scale,


static inline double scaleGetScaleC(scaling *scale, int i){
if (ISNA(scale->scaleC[i])){
if (ISNA(scale->scaleC[i]) || isnan(scale->scaleC[i])) {
switch (scale->xPar[i]){
case 1: // log
scale->scaleC[i]=1.0;
break;
case 2: // diag^2
scale->scaleC[i]=1.0/fabs(scale->initPar[i]);
scale->scaleC[i]=1.0/max2(fabs(scale->initPar[i]), scale->scaleCmin);
break;
case 3: // exp(diag)
scale->scaleC[i] = 1.0/2.0;
break;
case 4: // Identity diagonal chol(Omega ^-1)
case 5: // off diagonal chol(Omega^-1)
scale->scaleC[i] = 1.0/(2.0*fabs(scale->initPar[i]));
scale->scaleC[i] = 1.0/max2(2.0*fabs(scale->initPar[i]), scale->scaleCmin);
break;
default:
scale->scaleC[i]= 1.0/fabs(scale->initPar[i]);
scale->scaleC[i]= 1.0/max2(fabs(scale->initPar[i]), scale->scaleCmin);
break;
}
}
Expand Down
22 changes: 22 additions & 0 deletions tests/testthat/test-nls.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -48,6 +48,28 @@ nmTest({

})

test_that("nls all 1 issue", {

pheno <- function() {
ini({
tcl <- log(1) # typical value of clearance
tv <- log(1) # typical value of volume
add.err <- 0.1 # residual variability
})
model({
cl <- exp(tcl ) # individual value of clearance
v <- exp(tv) # individual value of volume
ke <- cl / v # elimination rate constant
d/dt(A1) = - ke * A1 # model differential equation
cp = A1 / v # concentration in plasma
cp ~ add(add.err) # define error model
})
}

expect_error(.nlmixr(pheno, nlmixr2data::pheno_sd, est="nls", nlsControl(algorithm="LM",print=0L)), NA)

})

test_that("nls makes sense", {

d <- nlmixr2data::theo_sd
Expand Down

0 comments on commit d33f9e4

Please sign in to comment.