#include #include #include #include #include /* DBL_EPSILON */ #include "c_optim.h" // #include "uncmin.cpp" namespace rstpm2 { using namespace Rcpp; double min(double a, double b) { return a < b ? a : b; } double max(double a, double b) { return a < b ? b : a; } double bound(double x, double lower, double upper) { return x < lower ? lower : (x > upper ? upper : x); } // void nmmin(int n, double *Bvec, double *X, double *Fmin, optimfn fminfn, // int *fail, double abstol, double intol, void *ex, // double alpha, double bet, double gamm, int trace, // int *fncount, int maxit) NelderMead::NelderMead(int trace, int maxit, double abstol, double reltol, double alpha, double beta, double gamma, double epshess, bool hessianp) : trace(trace), maxit(maxit), abstol(abstol), reltol(reltol), alpha(alpha), beta(beta), gamma(gamma), epshess(epshess), hessianp(hessianp) { } void NelderMead::optim(optimfn fn, NumericVector init, void * ex) { n = init.size(); coef = clone(init); nmmin(n, &init[0], &coef[0], &Fmin, fn, &fail, abstol, reltol, ex, alpha, beta, gamma, trace, &fncount, maxit); if (hessianp) hessian = calc_hessian(fn, ex); } NumericMatrix NelderMead::calc_hessian(optimfn fn, void * ex) { int n = coef.size(); NumericMatrix hess(n,n); double tmpi,tmpj,f1,f0,fm1,hi,hj,fij,fimj,fmij,fmimj; f0 = fn(n,&coef[0],ex); for(int i=0; i mask(n,1); vmmin(n, &init[0], &Fmin, fn, gr, maxit, trace, &mask[0], abstol, reltol, report, ex, &fncount, &grcount, &fail); coef = clone(init); if (hessianp) hessian = calc_hessian(gr, ex); } double BFGS::calc_objective(optimfn fn, NumericVector coef, void * ex) { return fn(coef.size(), &coef[0], ex); } double BFGS::calc_objective(optimfn fn, void * ex) { return fn(coef.size(), &coef[0], ex); } NumericMatrix BFGS::calc_hessian(optimgr gr, void * ex) { int n = coef.size(); NumericVector df1(clone(coef)); NumericVector df2(clone(coef)); NumericMatrix hess(n,n); double tmp; for(int i=0; i typsize(n,1.0), gpls(n,0.0), a(n*n,0.0), wrk(n*8,0.0); double norm, fpls; NumericVector xpls(n); // stepmax calculations if (stepmx == -1.0) { norm = 0.0; for (int i=0; i= 2) msg = 17; } double Brent_fmin(double ax, double bx, double (*f)(double, void *), void *info, double tol) { /* c is the squared inverse of the golden ratio */ const double c = (3. - sqrt(5.)) * .5; /* Local variables */ double a, b, d, e, p, q, r, u, v, w, x; double t2, fu, fv, fw, fx, xm, eps, tol1, tol3; /* eps is approximately the square root of the relative machine precision. */ eps = DBL_EPSILON; tol1 = eps + 1.;/* the smallest 1.000... > 1 */ eps = sqrt(eps); a = ax; b = bx; v = a + c * (b - a); w = v; x = v; d = 0.;/* -Wall */ e = 0.; fx = (*f)(x, info); fv = fx; fw = fx; tol3 = tol / 3.; /* main loop starts here ----------------------------------- */ for(;;) { xm = (a + b) * .5; tol1 = eps * fabs(x) + tol3; t2 = tol1 * 2.; /* check stopping criterion */ if (fabs(x - xm) <= t2 - (b - a) * .5) break; p = 0.; q = 0.; r = 0.; if (fabs(e) > tol1) { /* fit parabola */ r = (x - w) * (fx - fv); q = (x - v) * (fx - fw); p = (x - v) * q - (x - w) * r; q = (q - r) * 2.; if (q > 0.) p = -p; else q = -q; r = e; e = d; } if (fabs(p) >= fabs(q * .5 * r) || p <= q * (a - x) || p >= q * (b - x)) { /* a golden-section step */ if (x < xm) e = b - x; else e = a - x; d = c * e; } else { /* a parabolic-interpolation step */ d = p / q; u = x + d; /* f must not be evaluated too close to ax or bx */ if (u - a < t2 || b - u < t2) { d = tol1; if (x >= xm) d = -d; } } /* f must not be evaluated too close to x */ if (fabs(d) >= tol1) u = x + d; else if (d > 0.) u = x + tol1; else u = x - tol1; fu = (*f)(u, info); /* update a, b, v, w, and x */ if (fu <= fx) { if (u < x) b = x; else a = x; v = w; w = x; x = u; fv = fw; fw = fx; fx = fu; } else { if (u < x) a = u; else b = u; if (fu <= fw || w == x) { v = w; fv = fw; w = u; fw = fu; } else if (fu <= fv || v == x || v == w) { v = u; fv = fu; } } } /* end of main loop */ return x; } // R CMD INSTALL ~/src/R/microsimulation // R -q -e "require(microsimulation); .Call('test_nmmin',1:2,PACKAGE='microsimulation')" // .Call("optim_stpm2",init,X,XD,rep(bhazard,nrow(X)),wt,ifelse(event,1,0),package="rstpm2") } // namespace rstpm2