[Lme4-commits] r1488 - in pkg/lme4Eigen: . R src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Dec 22 22:56:24 CET 2011
Author: dmbates
Date: 2011-12-22 22:56:23 +0100 (Thu, 22 Dec 2011)
New Revision: 1488
Modified:
pkg/lme4Eigen/NAMESPACE
pkg/lme4Eigen/R/AllClass.R
pkg/lme4Eigen/src/external.cpp
pkg/lme4Eigen/src/optimizer.cpp
pkg/lme4Eigen/src/optimizer.h
Log:
Tuned the Nelder-Mead optimizer a bit and got rid of some of the noise.
Modified: pkg/lme4Eigen/NAMESPACE
===================================================================
--- pkg/lme4Eigen/NAMESPACE 2011-12-22 19:51:52 UTC (rev 1487)
+++ pkg/lme4Eigen/NAMESPACE 2011-12-22 21:56:23 UTC (rev 1488)
@@ -85,6 +85,7 @@
# and the rest (S3 generics; regular functions):
export(GHrule,
+ NelderMead,
VarCorr,
bootMer,
devcomp,
@@ -111,7 +112,7 @@
subbars
)
-exportClasses(
+exportClasses(NelderMead,
glmerMod,
glmFamily,
glmResp,
Modified: pkg/lme4Eigen/R/AllClass.R
===================================================================
--- pkg/lme4Eigen/R/AllClass.R 2011-12-22 19:51:52 UTC (rev 1487)
+++ pkg/lme4Eigen/R/AllClass.R 2011-12-22 21:56:23 UTC (rev 1488)
@@ -457,6 +457,52 @@
)
)
+NelderMead <-
+ setRefClass("NelderMead", # Reverse communication implementation of Nelder-Mead simplex optimizer
+ fields =
+ list(
+ Ptr = "externalptr",
+ lowerbd = "numeric",
+ upperbd = "numeric",
+ xstep = "numeric",
+ xtol = "numeric"
+ ),
+ methods =
+ list(
+ initialize = function(lower, upper, xst, x0, xt, ...) {
+ stopifnot((n <- length(lower <- as.numeric(lower))) > 0L,
+ length(upper <- as.numeric(upper)) == n,
+ all(lower < upper),
+ length(xst <- as.numeric(xst)) == n,
+ all(xst != 0),
+ length(x0 <- as.numeric(x0)) == n,
+ all(x0 >= lower),
+ all(x0 <= upper),
+ all(is.finite(x0)),
+ length(xt <- as.numeric(xt)) == n,
+ all(xt > 0))
+ lowerbd <<- lower
+ upperbd <<- upper
+ xstep <<- xst
+ xtol <<- xt
+ Ptr <<- .Call(NelderMead_Create, lowerbd, upperbd, xstep, x0, xtol)
+ },
+ ptr = function() {
+ if (length(lowerbd))
+ if (.Call(isNullExtPtr, Ptr))
+ Ptr <<- .Call(NelderMead_Create, lowerbd, upperbd, xstep, x0, xtol)
+ Ptr
+ },
+ newf = function(value) {
+ stopifnot(length(value <- as.numeric(value)) == 1L)
+ .Call(NelderMead_newf, ptr(), value)
+ },
+ value = function() .Call(NelderMead_value, ptr()),
+ xeval = function() .Call(NelderMead_xeval, ptr()),
+ xpos = function() .Call(NelderMead_xpos, ptr())
+ )
+ )
+
setClass("merMod",
representation(Gp = "integer",
call = "call",
Modified: pkg/lme4Eigen/src/external.cpp
===================================================================
--- pkg/lme4Eigen/src/external.cpp 2011-12-22 19:51:52 UTC (rev 1487)
+++ pkg/lme4Eigen/src/external.cpp 2011-12-22 21:56:23 UTC (rev 1488)
@@ -596,11 +596,11 @@
SEXP NelderMead_Create(SEXP lb_, SEXP ub_, SEXP xstep0_, SEXP x_, SEXP xtol_) {
BEGIN_RCPP;
MVec lb(as<MVec>(lb_)), ub(as<MVec>(ub_)), xstep0(as<MVec>(xstep0_)), x(as<MVec>(x_)), xtol(as<MVec>(xtol_));
- Rcpp::Rcout << "lb: " << lb.adjoint() << std::endl;
- Rcpp::Rcout << "ub: " << ub.adjoint() << std::endl;
- Rcpp::Rcout << "xstep0: " << xstep0.adjoint() << std::endl;
- Rcpp::Rcout << "x: " << x.adjoint() << std::endl;
- Rcpp::Rcout << "xtol: " << xtol.adjoint() << std::endl;
+ // Rcpp::Rcout << "lb: " << lb.adjoint() << std::endl;
+ // Rcpp::Rcout << "ub: " << ub.adjoint() << std::endl;
+ // Rcpp::Rcout << "xstep0: " << xstep0.adjoint() << std::endl;
+ // Rcpp::Rcout << "x: " << x.adjoint() << std::endl;
+ // Rcpp::Rcout << "xtol: " << xtol.adjoint() << std::endl;
Nelder_Mead *ans =
new Nelder_Mead(lb, ub, xstep0, x, optimizer::nl_stop(as<MVec>(xtol_)));
return wrap(XPtr<Nelder_Mead>(ans, true));
Modified: pkg/lme4Eigen/src/optimizer.cpp
===================================================================
--- pkg/lme4Eigen/src/optimizer.cpp 2011-12-22 19:51:52 UTC (rev 1487)
+++ pkg/lme4Eigen/src/optimizer.cpp 2011-12-22 21:56:23 UTC (rev 1488)
@@ -55,7 +55,8 @@
d_xeval( x),
d_minf( std::numeric_limits<Scalar>::infinity()),
d_stage( nm_restart),
- d_stop( stp)
+ d_stop( stp),
+ d_verb( 10)
{
if (!d_n || d_lb.size() != d_n ||
d_ub.size() != d_n || d_xstep.size() != d_n)
@@ -95,9 +96,9 @@
* @return status
*/
nm_status Nelder_Mead::newf(const Scalar& f) {
- Rcpp::Rcout << "f = " << f
- << " at " << d_xeval.adjoint() << std::endl;
d_stop.incrEvals();
+ if (d_verb > 0 && (d_stop.ev() % d_verb) == 0)
+ Rcpp::Rcout << "f = " << value() << " at " << d_x.adjoint() << std::endl;
if (d_stop.forced()) return nm_forced;
if (f < d_minf) {
d_minf = f;
@@ -142,44 +143,45 @@
nm_status Nelder_Mead::restart(const Scalar& f) {
d_fl = d_vals.minCoeff(&d_il);
d_fh = d_vals.maxCoeff(&d_ih);
- Rcpp::Rcout << "restart: neval = " << d_stop.ev()
- << ", d_fl = " << d_fl
- << ", d_fh = " << d_fh
- << ", d_il = " << d_il
- << ", d_ih = " << d_ih << std::endl;
+ // Rcpp::Rcout << "restart: neval = " << d_stop.ev()
+ // << ", d_fl = " << d_fl
+ // << ", d_fh = " << d_fh
+ // << ", d_il = " << d_il
+ // << ", d_ih = " << d_ih << std::endl;
d_c = (d_pts.rowwise().sum() - d_pts.col(d_ih)) / d_n; // compute centroid
- Rcpp::Rcout << "d_c: " << d_c.adjoint() << std::endl;
+ // Rcpp::Rcout << "centroid: " << d_c.adjoint() << std::endl;
// Check if the simplex has gotten to be too small. First
// calculate the coefficient-wise maximum abs deviation for
// the centroid.
d_xcur = (d_pts.colwise() - d_c).array().abs().rowwise().maxCoeff();
- if (d_stop.x(VectorXd::Zero(d_n), d_xcur)) return nm_xcvg;
+ // Rcpp::Rcout << "Radius: " << d_xcur.adjoint() << std::endl;
+ if (d_stop.x(VectorXd::Constant(d_n, 0.), d_xcur)) return nm_xcvg;
if (!reflectpt(d_xcur, d_c, alpha, d_pts.col(d_ih))) return nm_xcvg;
- Rcpp::Rcout << "d_xcur (after reflectpt): " << d_xcur.adjoint() << std::endl;
+ // Rcpp::Rcout << "d_xcur (after reflectpt): " << d_xcur.adjoint() << std::endl;
d_xeval = d_xcur;
d_stage = nm_postreflect;
return nm_active;
}
nm_status Nelder_Mead::postreflect(const Scalar& f) {
- Rcpp::Rcout << "postreflect: ";
+ // Rcpp::Rcout << "postreflect: ";
if (f < d_fl) { // new best point, try to expand
if (!reflectpt(d_xeval, d_c, gamm, d_pts.col(d_ih))) return nm_xcvg;
- Rcpp::Rcout << "New best point" << std::endl;
+ // Rcpp::Rcout << "New best point" << std::endl;
d_stage = nm_postexpand;
f_old = f;
return nm_active;
}
if (f < d_fh) { // accept new point
- Rcpp::Rcout << "Accept new point" << std::endl;
+ // Rcpp::Rcout << "Accept new point" << std::endl;
d_vals[d_ih] = f;
d_pts.col(d_ih) = d_xeval;
return restart(f);
}
// new worst point, contract
- Rcpp::Rcout << "New worst point" << std::endl;
+ // Rcpp::Rcout << "New worst point" << std::endl;
if (!reflectpt(d_xcur, d_c, d_fh <= f ? -beta : beta, d_pts.col(d_ih))) return nm_xcvg;
f_old = f;
d_xeval = d_xcur;
@@ -189,11 +191,11 @@
nm_status Nelder_Mead::postexpand(const Scalar& f) {
if (f < d_vals[d_ih]) { // expanding improved
- Rcpp::Rcout << "successful expand" << std::endl;
+ // Rcpp::Rcout << "successful expand" << std::endl;
d_pts.col(d_ih) = d_xeval;
d_vals[d_ih] = f;
} else {
- Rcpp::Rcout << "unsuccessful expand" << std::endl;
+ // Rcpp::Rcout << "unsuccessful expand" << std::endl;
d_pts.col(d_ih) = d_xcur;
d_vals[d_ih] = f_old;
}
@@ -202,12 +204,12 @@
nm_status Nelder_Mead::postcontract(const Scalar& f) {
if (f < f_old && f < d_fh) {
- Rcpp::Rcout << "successful contraction:" << std::endl;
+ // Rcpp::Rcout << "successful contraction:" << std::endl;
d_pts.col(d_ih) = d_xeval;
d_vals[d_ih] = f;
return restart(f);
}
- Rcpp::Rcout << "unsuccessful contraction, shrink simplex" << std::endl;
+ // Rcpp::Rcout << "unsuccessful contraction, shrink simplex" << std::endl;
for (Index i = 0; i <= d_n; ++i) {
if (i != d_il) {
if (!reflectpt(d_xeval, d_pts.col(d_il), -delta, d_pts.col(i))) return nm_xcvg;
Modified: pkg/lme4Eigen/src/optimizer.h
===================================================================
--- pkg/lme4Eigen/src/optimizer.h 2011-12-22 19:51:52 UTC (rev 1487)
+++ pkg/lme4Eigen/src/optimizer.h 2011-12-22 21:56:23 UTC (rev 1488)
@@ -102,6 +102,7 @@
nm_status d_stat;
nm_stage d_stage;
nl_stop d_stop;
+ Index d_verb; /**< verbosity, if > 0 results are displayed every d_verb evaluations */
public:
Nelder_Mead(const VectorXd&, const VectorXd&, const VectorXd&, const VectorXd&,
const nl_stop&);
More information about the Lme4-commits
mailing list