#include <SequentialMinimizer.h>

Inheritance diagram for cmsmath::OneDimMinimizer:

Public Types
enum	ImproveRet { Unchanged = 2, Done = 1, NotDone = 0 }
Public Member Functions
const char *	cname () const
ImproveRet	improve (int steps=1, double ytol=0, double xtol=0, bool force=true)
void	init (const MinimizerContext &ctx, unsigned int idx, double xmin, double xmax, double xstep, const std::string &name)
void	initDefault (const MinimizerContext &ctx, unsigned int idx)
void	initUnbound (const MinimizerContext &ctx, unsigned int idx, double xstep, const std::string &name)
bool	isInit () const
double	max () const
double	min () const
bool	minimize (int steps=1, double ytol=0, double xtol=0)
void	moveTo (double x)
const std::string &	name () const
	OneDimMinimizer (const MinimizerContext &ctx, unsigned int idx, double xmin, double xmax, double xstep, const std::string &name)
	OneDimMinimizer ()
	OneDimMinimizer (const MinimizerContext &ctx, unsigned int idx)
double	step () const
Private Member Functions
void	assign (int to, int from)
bool	doloop (int steps, double ytol, double xtol)
double	eval (double x)
double	eval ()
void	goldenBisection ()
	do the golden bisection
double	parabolaFit ()
	do the parabola fit
bool	parabolaStep ()
	do the parabola fit
void	reseek ()
void	seek ()
	search for a triplet of points bracketing the maximum. return false in case of errors
double &	x ()
	evaluate function at x
Private Attributes
const MinimizerContext *	f_
unsigned int	idx_
std::string	name_
double	xi_ [3]
double	xmax_
double	xmin_
double	xstep_
double	yi_ [3]

Detailed Description

Definition at line 26 of file SequentialMinimizer.h.

Member Enumeration Documentation

enum cmsmath::OneDimMinimizer::ImproveRet

improve minimum, re-evaluating points (other parameters of the function might have changed) return value is: 0: minimum changed, and did not converge there yet 1: minimum changed, but did converge to it 2: minimum did not change significantly in case 2, then the final position is NOT changed at all. force = true will instead force the update even if it's trivial

Enumerator:

Unchanged
Done
NotDone

Definition at line 63 of file SequentialMinimizer.h.

{ Unchanged = 2, Done = 1, NotDone = 0 };

Constructor & Destructor Documentation

cmsmath::OneDimMinimizer::OneDimMinimizer ( ) [inline]

Definition at line 28 of file SequentialMinimizer.h.

: f_(0), idx_(0) {}

cmsmath::OneDimMinimizer::OneDimMinimizer	(	const MinimizerContext &	ctx,
		unsigned int	idx
	)		`[inline]`

Definition at line 29 of file SequentialMinimizer.h.

                                                                           :
                f_(&ctx), idx_(idx) {}

cmsmath::OneDimMinimizer::OneDimMinimizer	(	const MinimizerContext &	ctx,
		unsigned int	idx,
		double	xmin,
		double	xmax,
		double	xstep,
		const std::string &	name
	)		`[inline]`

Definition at line 31 of file SequentialMinimizer.h.

                                                                                                                                          : 
                f_(&ctx), idx_(idx), name_(name), xmin_(xmin), xmax_(xmax), xstep_(xstep) {}

Member Function Documentation

void cmsmath::OneDimMinimizer::assign	(	int	to,
		int	from
	)		`[inline, private]`

Definition at line 106 of file SequentialMinimizer.h.

References Capri::details::from(), xi_, and yi_.

{ xi_[to] = xi_[from]; yi_[to] = yi_[from]; }

const char* cmsmath::OneDimMinimizer::cname ( ) const [inline]

Definition at line 35 of file SequentialMinimizer.h.

References name_.

Referenced by cmsmath::SequentialMinimizer::improve().

{ return name_.c_str(); }

bool cmsmath::OneDimMinimizer::doloop	(	int	steps,
		double	ytol,
		double	xtol
	)		`[private]`

basic loop return false if ended steps, true if reached tolerance

Definition at line 113 of file SequentialMinimizer.cc.

References DEBUG_ODM_printf, i, max(), and relval_steps::steps.

{
    if (steps <= 0) steps = 100;
    for (int i = 0; i < steps; ++i) {
        if (xtol > 0 && (xi_[2] - xi_[0]) < xtol) {
            return true;
        }
        goldenBisection();
        if (ytol > 0 && (std::max(yi_[2],yi_[0]) - yi_[1]) < ytol) {
            DEBUG_ODM_printf("ODM: intermediate ytol for %s: ymin %.8f, ymax %.8f, diff %.8f\n", name_.c_str(), yi_[1], std::max(yi_[2],yi_[0]), std::max(yi_[2],yi_[0]) - yi_[1]);
            return true;
        }
        DEBUG_ODM_printf("ODM: step %d/%d done for %s: ymin %.8f, ymax %.8f, diff %.8f\n", i+1, steps, name_.c_str(), yi_[1], std::max(yi_[2],yi_[0]), std::max(yi_[2],yi_[0]) - yi_[1]);
        DEBUG_ODM_printf("ODM: %s x = [%.4f, %.4f, %.4f], y = [%.4f, %.4f, %.4f]\n", name_.c_str(), xi_[0], xi_[1], xi_[2], yi_[0], yi_[1], yi_[2]);
    }

    return false;
}

double cmsmath::OneDimMinimizer::eval ( ) [inline, private]

Definition at line 103 of file SequentialMinimizer.h.

References cmsmath::MinimizerContext::eval(), and f_.

{ return f_->eval(); }

double cmsmath::OneDimMinimizer::eval ( double x ) [inline, private]

Definition at line 104 of file SequentialMinimizer.h.

References cmsmath::MinimizerContext::cleanEval(), f_, and idx_.

{ return f_->cleanEval(idx_, x); }

void cmsmath::OneDimMinimizer::goldenBisection ( ) [private]

do the golden bisection

pre-condition: the endpoint equal to x1 is not the highest

Definition at line 199 of file SequentialMinimizer.cc.

References assign.

{
    //pre-condition: always a sorted interval
    assert(xi_[0] < xi_[2]);
    assert(xi_[0] <= xi_[1] && xi_[1] <= xi_[2]);
    // if midpoint is not not one of the extremes, it's not higher than that extreme
    assert(xi_[1] == xi_[0] || yi_[1] <= yi_[0]); 
    assert(xi_[1] == xi_[2] || yi_[1] <= yi_[2]);
    if (xi_[0] == xi_[1] || xi_[1] == xi_[2]) {
        int isame = (xi_[0] == xi_[1] ? 0 : 2);
        assert(yi_[isame] <= yi_[2-isame]);
        xi_[1] = 0.5*(xi_[0]+xi_[2]);
        yi_[1] = eval(xi_[1]);
        if (yi_[1] < yi_[isame]) {
            // maximum is in the interval-
            // leave as is, next bisection steps will find it
        } else {
            // maximum remains on the boundary, leave both points there
            assign(2-isame, 1);
            assign(1, isame); 
        }
    } else {
        int inear = 2, ifar = 0;
        if (xi_[2]-xi_[1] > xi_[1] - xi_[0]) {
            inear = 0, ifar = 2;
        } else {
            inear = 2, ifar = 0;
        }
        double xc = xi_[1]*GOLD_R1 + xi_[ifar]*GOLD_R2;
        double yc = eval(xc);
        //DEBUG_ODM_printf("ODM: goldenBisection:\n\t\tfar = (%.2f,%.8f)\n\t\tnear = (%.2f,%.8f)\n\t\tcenter  = (%.2f,%.8f)\n\t\tnew  = (%.2f,%.8f)\n",
        //            xi_[ifar],  yi_[ifar], xi_[inear], yi_[inear], xi_[1], yi_[1], xc, yc);
        if (yc < yi_[1]) {   // then use 1, c, far
            assign(inear, 1);
            xi_[1] = xc; yi_[1] = yc;
        } else {  // then use c, 1, near
            xi_[ifar] = xc; yi_[ifar] = yc;
        }
    }
    //post-condition: always a sorted interval
    assert(xi_[0] < xi_[2]);
    assert(xi_[0] <= xi_[1] && xi_[1] <= xi_[2]);
    // if midpoint is not not one of the extremes, it's not higher than that extreme
    assert(xi_[1] == xi_[0] || yi_[1] <= yi_[0]); 
    assert(xi_[1] == xi_[2] || yi_[1] <= yi_[2]);

}

cmsmath::OneDimMinimizer::ImproveRet cmsmath::OneDimMinimizer::improve	(	int	steps = `1`,
		double	ytol = `0`,
		double	xtol = `0`,
		bool	force = `true`
	)

Definition at line 47 of file SequentialMinimizer.cc.

References DEBUG_ODM_printf, run_regression::done, max(), and x.

Referenced by cmsmath::SequentialMinimizer::improve().

{
    double x0 = x();
    if (x0 < xi_[0] || x0 > xi_[2]) {
        // could happen if somebody outside this loop modifies some parameters
        DEBUG_ODM_printf("ODM: ALERT: variable %s outside bounds x = [%.4f, %.4f, %.4f], x0 = %.4f\n", name_.c_str(), xi_[0], xi_[1], xi_[2], x0);
        x0 = x() = xi_[1]; 
    } else {
        xi_[1] = x0;
    }
    double y0 = eval();
    yi_[1] = y0;
    yi_[0] = eval(xi_[0]);
    yi_[2] = eval(xi_[2]);
    if (xtol == 0) xtol = (fabs(xi_[1])+XTOL)*XTOL;

    DEBUG_ODM_printf("ODM: start of improve %s x = [%.4f, %.4f, %.4f], y = [%.4f, %.4f, %.4f]\n", name_.c_str(), xi_[0], xi_[1], xi_[2], yi_[0], yi_[1], yi_[2]);

    if (yi_[1] <= yi_[0] && yi_[1] <= yi_[2]) {
        if (ytol > 0 && (std::max(yi_[2],yi_[0]) - yi_[1]) < ytol) {
            DEBUG_ODM_printf("ODM: immediate ytol for %s: ymin %.8f, ymax %.8f, diff %.8f\n", name_.c_str(), yi_[1], std::max(yi_[2],yi_[0]), std::max(yi_[2],yi_[0]) - yi_[1]);
            if (!force || parabolaStep()) return Unchanged;
        }
        if (xtol > 0 && (xi_[2] - xi_[0]) < xtol) {
            DEBUG_ODM_printf("ODM: immediate xtol for %s: xmin %.8f, xmax %.8f, diff %.8f\n", name_.c_str(), xi_[0], xi_[2], xi_[2] - xi_[0]);
            if (!force || parabolaStep()) return Unchanged;
        }
    } else {
        reseek();
    }

    //post-condition: always a sorted interval
    assert(xi_[0] < xi_[2]);
    assert(xi_[0] <= xi_[1] && xi_[1] <= xi_[2]);
    // if midpoint is not not one of the extremes, it's not higher than that extreme
    assert(xi_[1] == xi_[0] || yi_[1] <= yi_[0]); 
    assert(xi_[1] == xi_[2] || yi_[1] <= yi_[2]);

    bool done = doloop(steps,ytol,xtol);
    parabolaStep(); 

    //post-condition: always a sorted interval
    assert(xi_[0] < xi_[2]);
    assert(xi_[0] <= xi_[1] && xi_[1] <= xi_[2]);
    // if midpoint is not not one of the extremes, it's not higher than that extreme
    assert(xi_[1] == xi_[0] || yi_[1] <= yi_[0]); 
    assert(xi_[1] == xi_[2] || yi_[1] <= yi_[2]);


    if (ytol > 0 && fabs(yi_[1] - y0) < ytol) {
        DEBUG_ODM_printf("ODM: final ytol for %s: ymin(old) %.8f, ymin(new) %.8f, diff %.8f\n", name_.c_str(), y0, yi_[1], y0 -yi_[1]);
        if (!force) x() = x0;
        return Unchanged;
    } 

    if (xtol > 0 && fabs(xi_[1] - x0) < xtol) {
        x() = (force ? xi_[1] : x0);
        return Unchanged;
    }
    DEBUG_ODM_printf("ODM: doloop for %s is %s\n", name_.c_str(), done ? "Done" : "NotDone");
    DEBUG_ODM_printf("ODM: end of improve %s x = [%.4f, %.4f, %.4f], y = [%.4f, %.4f, %.4f]\n", name_.c_str(), xi_[0], xi_[1], xi_[2], yi_[0], yi_[1], yi_[2]);
    x() = xi_[1];
    xstep_ = xi_[2] - xi_[0];
    return done ? Done : NotDone; 
}

void cmsmath::OneDimMinimizer::init	(	const MinimizerContext &	ctx,
		unsigned int	idx,
		double	xmin,
		double	xmax,
		double	xstep,
		const std::string &	name
	)		`[inline]`

Definition at line 41 of file SequentialMinimizer.h.

References f_, idx_, name(), name_, xmax_, xmin_, and xstep_.

Referenced by initDefault(), and initUnbound().

                                                                                                                                    {
                f_ = &ctx; idx_ = idx; 
                xmin_ = xmin; xmax_ = xmax; xstep_ = xstep;
                name_ = name;
            }

void cmsmath::OneDimMinimizer::initDefault	(	const MinimizerContext &	ctx,
		unsigned int	idx
	)

Definition at line 30 of file SequentialMinimizer.cc.

References infinity, and init().

                                                                                      {
    init(ctx, idx, -std::numeric_limits<double>::infinity(), std::numeric_limits<double>::infinity(), 1.0, Form("#%d", idx));
}

void cmsmath::OneDimMinimizer::initUnbound	(	const MinimizerContext &	ctx,
		unsigned int	idx,
		double	xstep,
		const std::string &	name
	)		`[inline]`

Definition at line 46 of file SequentialMinimizer.h.

References infinity, and init().

                                                                                                                 {
                init(ctx, idx, -std::numeric_limits<double>::infinity(), std::numeric_limits<double>::infinity(), xstep, name);
            }

bool cmsmath::OneDimMinimizer::isInit ( ) const [inline]

Definition at line 40 of file SequentialMinimizer.h.

References f_.

Referenced by cmsmath::SequentialMinimizer::minimize().

{ return f_ != 0; }

double cmsmath::OneDimMinimizer::max ( ) const [inline]

Definition at line 36 of file SequentialMinimizer.h.

References xmax_.

Referenced by cmsmath::SequentialMinimizer::doFullMinim().

{ return  xmax_; }

double cmsmath::OneDimMinimizer::min ( ) const [inline]

Definition at line 37 of file SequentialMinimizer.h.

References xmin_.

Referenced by cmsmath::SequentialMinimizer::doFullMinim().

{ return  xmin_; }

bool cmsmath::OneDimMinimizer::minimize	(	int	steps = `1`,
		double	ytol = `0`,
		double	xtol = `0`
	)

do N steps of golden bisection (or until tolerance in x, y is reached) return true if converged, false if finished steps

Definition at line 35 of file SequentialMinimizer.cc.

References run_regression::done, and x.

Referenced by cmsmath::SequentialMinimizer::minimize().

{
    // get initial bracket
    seek();
    
    bool done = doloop(steps,ytol,xtol);
    parabolaStep();
    x() = xi_[1];
    xstep_ = xi_[2] - xi_[0];
    return done;
}

void cmsmath::OneDimMinimizer::moveTo ( double x )

Definition at line 285 of file SequentialMinimizer.cc.

References max(), min, and x.

Referenced by cmsmath::SequentialMinimizer::doFullMinim().

                                            {
    if (x == xmax_) {
        xi_[0] = xmax_ - (xi_[2]-xi_[0]); yi_[0] = eval(xi_[0]);
        xi_[1] = xmax_; yi_[1] = eval(xmax_);
        xi_[2] = xmax_; yi_[2] = yi_[1];
    } else if (x == xmin_) {
        xi_[2] = xmin_ + (xi_[2]-xi_[0]); yi_[2] = eval(xi_[0]);
        xi_[1] = xmin_; yi_[1] = eval(xmin_);
        xi_[0] = xmin_; yi_[0] = yi_[1];
    } else {
        double dx = xi_[2] - xi_[0];
        xi_[1] = x; yi_[0] = eval(x);
        xi_[0] = std::max(xmin_, x-dx); yi_[0] = eval(xi_[0]);
        xi_[2] = std::min(xmax_, x+dx); yi_[2] = eval(xi_[2]);
    }
}

const std::string& cmsmath::OneDimMinimizer::name ( void ) const [inline]

Definition at line 34 of file SequentialMinimizer.h.

References name_.

Referenced by cmsmath::SequentialMinimizer::doFullMinim(), and init().

{ return name_; }

double cmsmath::OneDimMinimizer::parabolaFit ( ) [private]

do the parabola fit

Definition at line 248 of file SequentialMinimizer.cc.

References x.

{
    if (xi_[0] == xi_[1] || xi_[1] == xi_[2]) { 
        return xi_[1]; 
    }
    double dx0 = xi_[1] - xi_[0], dx2 = xi_[1] - xi_[2];
    double dy0 = yi_[1] - yi_[0], dy2 = yi_[1] - yi_[2];
    double num = dx0*dx0*dy2 - dx2*dx2*dy0;
    double den = dx0*dy2     - dx2*dy0;
    if (den != 0) {
        double x = xi_[1] - 0.5*num/den;
        if (xi_[0] < x && x < xi_[2]) {
            return x;
        }
    } 
    return xi_[1];
}

bool cmsmath::OneDimMinimizer::parabolaStep ( ) [private]

do the parabola fit

Definition at line 266 of file SequentialMinimizer.cc.

                                          {
    double xc = parabolaFit();
    if (xc != xi_[1]) {
        double yc = eval(xc);
        if (yc < yi_[1]) {
            xi_[1] = xc; 
            yi_[1] = yc;
            //post-condition: always a sorted interval
            assert(xi_[0] < xi_[2]);
            assert(xi_[0] <= xi_[1] && xi_[1] <= xi_[2]);
            // if midpoint is not not one of the extremes, it's not higher than that extreme
            assert(xi_[1] == xi_[0] || yi_[1] <= yi_[0]); 
            assert(xi_[1] == xi_[2] || yi_[1] <= yi_[2]);
            return true;
        }
    }
    return false;
}

void cmsmath::OneDimMinimizer::reseek ( ) [private]

re-search for a triplet of points bracketing the maximum. return false in case of errors assume that value and error make sense

Definition at line 141 of file SequentialMinimizer.cc.

References assign, i, max(), min, and x.

{
    double xtol2 = 2*(fabs(xi_[1])+XTOL)*XTOL;
    if (xstep_ < xtol2) xstep_ = xtol2;
    xi_[1] = x(); 
    yi_[1] = eval(xi_[1]);
    xi_[0] = std::max(xmin_, xi_[1]-xstep_);
    yi_[0] = (xi_[0] == xi_[1] ? yi_[1] : eval(xi_[0]));
    xi_[2] = std::min(xmax_, xi_[1]+xstep_);
    yi_[2] = (xi_[2] == xi_[1] ? yi_[1] : eval(xi_[2]));
    assert(xi_[0] < xi_[2]);
    assert(xi_[0] <= xi_[1] && xi_[1] <= xi_[2]);

    for (;;) {
        //DEBUG_ODM_printf("ODM: bracketing %s in x = [%.4f, %.4f, %.4f], y = [%.4f, %.4f, %.4f]\n", name_.c_str(), xi_[0], xi_[1], xi_[2], yi_[0], yi_[1], yi_[2]);
        if (yi_[0] < yi_[1]) {
            assign(2,1); // 2:=1
            assign(1,0); // 1:=0
            xi_[0] = std::max(xmin_, xi_[1]-xstep_);
            yi_[0] = (xi_[0] == xi_[1] ? yi_[1] : eval(xi_[0]));
        } else if (yi_[2]  < yi_[1]) {
            assign(0,1); // 0:=1
            assign(1,2); // 1:=2
            xi_[2] = std::min(xmax_, xi_[1]+xstep_);
            yi_[2] = (xi_[2] == xi_[1] ? yi_[1] : eval(xi_[2]));
        } else if (yi_[2] == yi_[1] && yi_[1] == yi_[0]) {
            // function is identical in three points --> constant?
            // try a scan
            const int nPoints = 20;
            double xi[nPoints], yi[nPoints];
            double dx = (xmax_-xmin_)/nPoints, x = xmin_ - 0.5*dx;
            bool isConst = true;
            int iFound = 0; 
            for (int i = 0; i < nPoints; ++i, x += dx) {
                xi[i] = x; yi[i] = eval(x);
                if (yi[i] != yi_[1]) isConst = false;
                if (yi[i] < yi[iFound]) { iFound = i; }
            }
            if (isConst) break;
            xi_[0] = (iFound == 0        ? xmin_ : xi[iFound-1]);
            xi_[2] = (iFound > nPoints-1 ? xmax_ : xi[iFound+1]);
            xi_[1] = iFound; yi_[1] = yi_[iFound];
            break;
        } else {
            xstep_ /= 2;
            break;
        }
        xstep_ *= 2;
    }
    //DEBUG_ODM_printf("ODM: bracketed minimum of %s in [%.4f, %.4f, %.4f]\n", name_.c_str(), xi_[0], xi_[1], xi_[2]);
    //post-condition: always a sorted interval
    assert(xi_[0] < xi_[2]);
    assert(xi_[0] <= xi_[1] && xi_[1] <= xi_[2]);
    // if midpoint is not not one of the extremes, it's not higher than that extreme
    assert(xi_[1] == xi_[0] || yi_[1] <= yi_[0]); 
    assert(xi_[1] == xi_[2] || yi_[1] <= yi_[2]);
}

void cmsmath::OneDimMinimizer::seek ( ) [private]

search for a triplet of points bracketing the maximum. return false in case of errors

Definition at line 132 of file SequentialMinimizer.cc.

References max().

{
    if (std::isfinite(xmax_-xmin_)) {
        xstep_ = std::max(xstep_, 0.2*(xmax_-xmin_));
    } else {
        xstep_ = 1.0;
    }
    reseek();
}