> # test.mem.R: test earth C code memory usage under both normal and error conditions
> 
> source("test.prolog.R")
> library(earth)
Loading required package: Formula
Loading required package: plotmo
Loading required package: plotrix
Loading required package: TeachingDemos
> 
> # the data we will build the models on
> ncases <- 10
> x <- matrix(1:ncases, ncol=1)
> colnames(x) <- "x"
> max <- max(x)
> y <- sin(3 * x / max(x))
> colnames(y) <- "y"
> 
> nmodels <- 5
> nlm            <- double(length=nmodels) # mem used for each lm model
> nstandardearth <- double(length=nmodels) # mem used for each earth model
> ngoodallowed   <- double(length=nmodels)
> nbadallowed    <- double(length=nmodels)
> nbadendspan    <- double(length=nmodels)
> 
> max.mem.change <- function(mem.start, gc.start)
+ {
+     mem <- memory.size()
+     gc <- gc(full=TRUE) # returns cells left after garbage collection
+     max(abs(mem - mem.start),
+         abs(gc[1,1] - gc.start[1,1]),  # Ncells
+         abs(gc[2,1] - gc.start[2,1]))  # Vcells
+ }
> plotmem <- function(nlm, nstandardearth, ngoodallowed, nbadallowed, nbadendspan)
+ {
+     min <- min(nlm, nstandardearth, ngoodallowed, nbadallowed, nbadendspan)
+     max <- max(nlm, nstandardearth, ngoodallowed, nbadallowed, nbadendspan)
+     min <- min - 1
+     max <- max + 3
+     yjitter <- (max - min) / 130 # minimize overplotting
+ 
+     # in the graphs, lines should be horizontal (at least after the first iter)
+     # if a line increases after the first iter, it means that memory is not being released
+     plot( 1:nmodels, nlm, type="l", main="memory used by each model",
+                           xlab="nmodels", ylab="memory change", ylim=c(min, max))
+     lines(1:nmodels, nstandardearth + 1 * yjitter, col=2)
+     lines(1:nmodels, ngoodallowed   + 2 * yjitter, col=3)
+     lines(1:nmodels, nbadallowed    + 3 * yjitter, col=1, lty=2)
+     lines(1:nmodels, nbadendspan    + 4 * yjitter, col=2, lty=2)
+ 
+     legend(x="topright", bg="white",
+            legend=c("lm", "standardearth", "goodallowed", "badallowed", "badendspan"),
+            lty=c(1,1,1,2,2),
+            col=c(1,2,3,1,2))
+ }
> good.allowedfunc  <- function(degree, pred, parents, namesx, first)
+ {
+     pred != 999
+ }
> bad.allowedfunc  <- function(degree, pred, parents, namesx, first)
+ {
+     # this stop is silent because call earth using try(..., silent=TRUE)
+     stop("early exit from bad.allowedfunc")
+ }
> cat("initial redundant run of lm\n") # else initial nlm very large (why?)
initial redundant run of lm
> print(summary(lm(y~x)))

Call:
lm(formula = y ~ x)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.45497 -0.17062  0.06399  0.23525  0.32264 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)  0.76940    0.21312    3.61  0.00688 **
x           -0.01891    0.03435   -0.55  0.59701   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.312 on 8 degrees of freedom
Multiple R-squared:  0.0365,	Adjusted R-squared:  -0.08394 
F-statistic: 0.303 on 1 and 8 DF,  p-value: 0.597

> for(i in 0:nmodels) {
+     try(lm(y~x), silent=FALSE)
+     gc <- gc(full=TRUE)
+     if(i <= 0) {
+         mem.start <- memory.size()
+         gc.start <- gc(full=TRUE)
+     } else
+         nlm[i] <- max.mem.change(mem.start, gc.start)
+ }
> cat("actual run of lm\n")
actual run of lm
> # We use 0:nmodels, because we build the first model at iter 0,
> # but don't save results from iter 0 (i.e. we the ignore first model).
> # This is because the first model sometimes leaves some memory allocated (why?).
> print(summary(lm(y~x)))

Call:
lm(formula = y ~ x)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.45497 -0.17062  0.06399  0.23525  0.32264 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)  0.76940    0.21312    3.61  0.00688 **
x           -0.01891    0.03435   -0.55  0.59701   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.312 on 8 degrees of freedom
Multiple R-squared:  0.0365,	Adjusted R-squared:  -0.08394 
F-statistic: 0.303 on 1 and 8 DF,  p-value: 0.597

> for(i in 0:nmodels) {
+     try(lm(y~x), silent=FALSE)
+     gc <- gc(full=TRUE)
+     if(i <= 0) {
+         mem.start <- memory.size()
+         gc.start <- gc(full=TRUE)
+     } else
+         nlm[i] <- max.mem.change(mem.start, gc.start)
+ }
> # standard earth model
> cat("earth(y~x)\n")
earth(y~x)
> print(summary(earth(y~x)))
Call: earth(formula=y~x)

            coefficients
(Intercept)    1.1315971
h(5-x)        -0.1980423
h(x-5)        -0.1787669

Selected 3 of 3 terms, and 1 of 1 predictors
Termination condition: RSq changed by less than 0.001 at 3 terms
Importance: x
Number of terms at each degree of interaction: 1 2 (additive model)
GCV 0.01888444    RSS 0.0472111    GRSq 0.8107224    RSq 0.941581
> for(i in 0:nmodels) {
+     try(earth(y~x), silent=FALSE)
+     gc <- gc(full=TRUE)
+     if(i <= 0) {
+         mem.start <- memory.size()
+         gc.start <- gc(full=TRUE)
+     } else
+         nstandardearth[i] <- max.mem.change(mem.start, gc.start)
+ }
> # earth model with an allowed func
> cat("earth(y~x, allowed = good.allowedfunc)\n")
earth(y~x, allowed = good.allowedfunc)
> print(summary(earth(y~x, allowed = good.allowedfunc)))
Call: earth(formula=y~x, allowed=good.allowedfunc)

            coefficients
(Intercept)    1.1315971
h(5-x)        -0.1980423
h(x-5)        -0.1787669

Selected 3 of 3 terms, and 1 of 1 predictors
Termination condition: RSq changed by less than 0.001 at 3 terms
Importance: x
Number of terms at each degree of interaction: 1 2 (additive model)
GCV 0.01888444    RSS 0.0472111    GRSq 0.8107224    RSq 0.941581
> for(i in 0:nmodels) {
+     try(earth(y~x, allowed = good.allowedfunc), silent=FALSE)
+     gc <- gc(full=TRUE)
+     if(i <= 0) {
+         mem.start <- memory.size()
+         gc.start <- gc(full=TRUE)
+     } else
+         ngoodallowed[i] <- max.mem.change(mem.start, gc.start)
+ }
> # try earth model with an allowed func which causes an error
> cat("earth(y~x, allowed = bad.allowedfunc)\n")
earth(y~x, allowed = bad.allowedfunc)
> expect.err(try(earth(y~x, allowed = bad.allowedfunc), silent=FALSE), "early exit from bad.allowedfunc")
Error in (function (degree, pred, parents, namesx, first)  : 
  early exit from bad.allowedfunc
Got error as expected from try(earth(y ~ x, allowed = bad.allowedfunc), silent = FALSE)
> for(i in 0:nmodels) {
+     try(earth(y~x, allowed = bad.allowedfunc), silent=TRUE)
+     gc <- gc(full=TRUE)
+     if(i <= 0) {
+         mem.start <- memory.size()
+         gc.start <- gc(full=TRUE)
+     } else
+         nbadallowed[i] <- max.mem.change(mem.start, gc.start)
+ }
> # try earth model with an arg that causes error in ForwardPass in earth.c
> cat("earth(y~x, Adjust.endspan = -999\n")
earth(y~x, Adjust.endspan = -999
> expect.err(try(earth(y~x, Adjust.endspan = -999), silent=FALSE), "Adjust.endspan is -999 but should be between 0 and 10")
Error in forward.pass(x, y, yw, weights, trace, degree, penalty, nk, thresh,  : 
  Adjust.endspan is -999 but should be between 0 and 10
Got error as expected from try(earth(y ~ x, Adjust.endspan = -999), silent = FALSE)
> for(i in 0:nmodels) {
+     try(earth(y~x, Adjust.endspan = -999), silent=TRUE)
+     gc <- gc(full=TRUE)
+     if(i <= 0) {
+         mem.start <- memory.size()
+         gc.start <- gc(full=TRUE)
+     } else
+         nbadendspan[i] <- max.mem.change(mem.start, gc.start)
+ }
> print(nlm)
[1] 9 9 9 9 9
> print(nstandardearth)
[1] 9 9 9 9 9
> print(ngoodallowed)
[1] 9 9 9 9 9
> print(nbadallowed)
[1] 9 9 9 9 9
> print(nbadendspan)
[1] 9 9 9 9 9
> 
> # printf("\n                             Min       1stQ     Median       Mean      3rdQ        Max\n")
> # printf("lm                    %s\n", paste0(sprintf("% 10.3f", summary(nlm)), collapse=" "))
> # printf("standardearth         %s\n", paste0(sprintf("% 10.3f", summary(nstandardearth)), collapse=" "))
> # printf("goodallowed           %s\n", paste0(sprintf("% 10.3f", summary(ngoodallowed)), collapse=" "))
> # printf("badallowed            %s\n", paste0(sprintf("% 10.3f", summary(nbadallowed)), collapse=" "))
> # printf("badendspan            %s\n", paste0(sprintf("% 10.3f", summary(nbadendspan)), collapse=" "))
> 
> # plot the data we are modeling
> plot(1:nrow(x), y, type="b", pch=20, xlab="x", main="the data we are modeling")
> 
> # plot memory used for each model
> plotmem(nlm, nstandardearth, ngoodallowed, nbadallowed, nbadendspan)
> 
> source("test.epilog.R")
