[Lme4-commits] r1782 - in pkg/lme4: R src

[noreply at r-forge.r-project.org]

Author: bbolker
Date: 2012-06-27 11:48:31 +0200 (Wed, 27 Jun 2012)
New Revision: 1782

Modified:
   pkg/lme4/R/lmer.R
   pkg/lme4/R/optimizer.R
   pkg/lme4/src/optimizer.cpp
Log:

   added more verbosity to Nelder-Mead
   *changed* MinfMax (minimum allowed objective function value) from double.xmin to -double.xmax (!)
   start values no longer ignored by lmer (but FIXME: still ignored by [gn]lmer?)
   added restart from boundary for lmer (ditto)



Modified: pkg/lme4/R/lmer.R
===================================================================
--- pkg/lme4/R/lmer.R	2012-06-27 08:38:59 UTC (rev 1781)
+++ pkg/lme4/R/lmer.R	2012-06-27 09:48:31 UTC (rev 1782)
@@ -81,6 +81,12 @@
                  contrasts = NULL, devFunOnly=FALSE,
                  optimizer="Nelder_Mead", ...)
 {
+    verbose <- as.integer(verbose)
+    restart <- TRUE
+    if (!is.null(control$restart)) {
+        restart <- control$restart
+        control$restart <- NULL
+    }
     if (sparseX) warning("sparseX = TRUE has no effect at present")
     mf <- mc <- match.call()
     ## '...' handling up front, safe-guarding against typos ("familiy") :
@@ -144,14 +150,27 @@
     rho$resp <- mkRespMod(fr, if(REML) p else 0L)
 
     devfun <- mkdevfun(rho, 0L)
-    devfun(reTrms$theta) # one evaluation to ensure all values are set
 
+    ## FIXME: should we apply start elsewhere? what about reTrms$theta?
+    if (!is.null(start)) {
+        rho$pp$setTheta(unlist(start))
+    }
+
+    devfun(rho$pp$theta) # one evaluation to ensure all values are set
+
     if (devFunOnly) return(devfun)
 
     opt <- optwrap(optimizer,
                    devfun, rho$pp$theta, lower=reTrms$lower, control=control,
-                   adj=FALSE)
+                   adj=FALSE, verbose=verbose)
 
+    if (restart && any(rho$pp$theta==0)) {
+        if (verbose) message("some theta parameters on the boundary, restarting")
+        opt <- optwrap(optimizer,
+                       devfun, rho$pp$theta, lower=reTrms$lower, control=control,
+                       adj=FALSE, verbose=verbose)
+    }
+        
     mkMerMod(environment(devfun), opt, reTrms, fr, mc)
 }## { lmer }
 
@@ -281,6 +300,7 @@
                   contrasts = NULL, mustart, etastart, devFunOnly = FALSE,
                   tolPwrss = 1e-7, optimizer=c("bobyqa","Nelder_Mead"), ...)
 {
+    ## FIXME: does start= do anything? test & fix
     verbose <- as.integer(verbose)
     mf <- mc <- match.call()
                                         # extract family, call lmer for gaussian
@@ -2092,9 +2112,9 @@
                        xst <- c(rep.int(0.1, length(environment(fn)$pp$theta)),  ## theta parameters
                                 sqrt(diag(environment(fn)$pp$unsc())))
                    control$xst <- 0.2*xst
-                   control$verbose <- verbose
                    if (is.null(control$xt)) control$xt <- control$xst*5e-4
                })
+    if (optimizer=="Nelder_Mead") control$verbose <- verbose
     if (optimizer=="bobyqa" && all(par==0)) par[] <- 0.001  ## minor kluge
     arglist <- list(fn=fn, par=par, lower=lower, upper=upper, control=control)
     ## optimx: must pass method in control (?) because 'method' was previously

Modified: pkg/lme4/R/optimizer.R
===================================================================
--- pkg/lme4/R/optimizer.R	2012-06-27 08:38:59 UTC (rev 1781)
+++ pkg/lme4/R/optimizer.R	2012-06-27 09:48:31 UTC (rev 1782)
@@ -62,7 +62,7 @@
     nM <- NelderMead$new(lower=lower, upper=upper, x0=par, xst=xst, xt=xt)
     cc <- do.call(function(iprint=0L, maxfun=10000L, FtolAbs=1e-5,
                            FtolRel=1e-15, XtolRel=1e-7,
-                           MinfMax=.Machine$double.xmin, ...) {
+                           MinfMax=-.Machine$double.xmax, ...) {
         if (length(list(...))>0) warning("unused control arguments ignored")
         list(iprint=iprint, maxfun=maxfun, FtolAbs=FtolAbs, FtolRel=FtolRel,
              XtolRel=XtolRel, MinfMax=MinfMax)

Modified: pkg/lme4/src/optimizer.cpp
===================================================================
--- pkg/lme4/src/optimizer.cpp	2012-06-27 08:38:59 UTC (rev 1781)
+++ pkg/lme4/src/optimizer.cpp	2012-06-27 09:48:31 UTC (rev 1782)
@@ -101,14 +101,35 @@
 	d_stop.incrEvals();
 	if (d_verb > 0 && (d_stop.ev() % d_verb) == 0)
 	    Rcpp::Rcout << "(NM) " << d_stop.ev() << ": " << "f = " << value() << " at " << d_x.adjoint() << std::endl;
-	if (d_stop.forced()) return nm_forced;
+	if (d_stop.forced()) {
+	    if (d_verb==1) {
+		Rcpp::Rcout << "(NM) stop_forced" << std::endl;
+	    }
+	    return nm_forced;
+	}
 	if (f < d_minf) {
 	    d_minf = f;
 	    d_x = d_xeval;	// save the value generating current minimum
-	    if (d_minf < d_stop.minfMax()) return nm_minf_max;
+	    if (d_minf < d_stop.minfMax()) {
+		if (d_verb==1) {
+		    Rcpp::Rcout << "(NM) nm_minf_max: " << d_minf << ", " 
+				<< d_stop.minfMax() << ", " << d_x << std::endl;
+		}
+		return nm_minf_max;
+	    }
 	}
-	if (d_stop.evals()) return nm_evals;
-	if (init_pos <= d_n) return init(f);
+	if (d_stop.evals()) {
+	    if (d_verb==1) {
+		Rcpp::Rcout << "(NM) nm_evals" << std::endl;
+	    }
+	    return nm_evals;
+	}
+	if (init_pos <= d_n) {
+	    if (d_verb==1) {
+		Rcpp::Rcout << "(NM) init_pos <= d_n" << std::endl;
+	    }
+	    return init(f);
+	}
 	switch (d_stage) {
 	case nm_restart:      return restart(f);
 	case nm_postreflect:  return postreflect(f);