brentls.cc

COOOL：CWP面向对象最优化库（CWP Object Oriented Optimization Library） COOOL是C++类的一个集合

//============================================================
// COOOL           version 1.1           ---     Nov,  1995
//   Center for Wave Phenomena, Colorado School of Mines
//============================================================
//
//   This code is part of a preliminary release of COOOL (CWP
// Object-Oriented Optimization Library) and associated class 
// libraries. 
//
// The COOOL library is a free software. You can do anything you want
// with it including make a fortune.  However, neither the authors,
// the Center for Wave Phenomena, nor anyone else you can think of
// makes any guarantees about anything in this package or any aspect
// of its functionality.
//
// Since you've got the source code you can also modify the
// library to suit your own purposes. We would appreciate it 
// if the headers that identify the authors are kept in the 
// source code.
//
//======================================================================
// Definition of the BrentLineSearch class
// Brent's line search algorithm
// author: Wenceslau Gouveia, Adapted from Numerical Recipes in C
// Modified to fit into new classes.  H. Lydia Deng, 02/21/94, 03/15/94
//========================================================================


#define CGOLD 		.3819660
#define ZEPS 		1.e-10
#define GOLD 		1.618034
#define STEP_MAX	5
#define GLIMIT 		100.
#define TINY 		1.0e-20

#include <BrentLS.hh>

BrentLineSearch::BrentLineSearch(ObjectiveFunction* f, int it) 
: LineSearch(f)
{
    iterMax 	= 	it;
    iterNum 	= 	0;
}

BrentLineSearch::~BrentLineSearch()
{
    if (fp != NULL) 	fp = NULL;
}

// Code for the BrentBrent line search
Model<double> BrentLineSearch::search(Model<double>& model0, Vector<double>& direction, 
		double tol, double delta)
{
    iterNum	=	0;
	Vector<double> steps(3);			// brackets
	Vector<double> of_values(3);			// OF evaluated inside bracket
	double dum;					// auxiliary quantity
	double r, q, ulim;				// auxiliary quantity
	double u, fu;					// define the new bracket limit

	// Variables related to Brent's algorithm
	double d, etemp, fv, fw, fx, p;			// auxiliary variables
	double tol1, tol2, v, w, xm;			// auxiliary variables
	double x, e = 0.;				// auxiliary variables
	int i;						// counter

	// Beggining of the bracketing stage

	steps[0] = 0.;
	steps[1] = STEP_MAX * delta;

	of_values[0]= fp->performance(model0);
	of_values[1] = fp->performance(model0.update(1,steps[1],direction));
	iterNum += 2;

	if (of_values[1] > of_values[0])
	{
	   dum = steps[0]; steps[0] = steps[1]; steps[1] = dum;
	   dum = of_values[0]; of_values[0] = of_values[1]; of_values[1] = dum;
	}	

	steps[2] = steps[1] + GOLD * (steps[1] - steps[0]);
	of_values[2] = fp->performance(model0.update(1,steps[2],direction));
	iterNum++;
	while (of_values[1] > of_values[2])
	{
		r = (steps[1] - steps[0]) * (of_values[1] - of_values[2]);
		q = (steps[1] - steps[2]) * (of_values[1] - of_values[0]);
		u = steps[1] - ((steps[1] - steps[2]) * q - 
				(steps[1] - steps[0]) * r) / 
		    (2. * Abs(Max(Abs(q - r), TINY)) / Sgn(q - r));
		ulim = steps[1] + GLIMIT * (steps[2] - steps[1]);

		if ((steps[1] - u) * (u - steps[2]) > 0.)
		{
			fu = fp->performance(model0.update(1,u,direction));
			iterNum++;
			if (fu < of_values[2])
			{
				steps[0] = steps[1];
				steps[1] = u;
				of_values[0] = of_values[1];
				of_values[1] = fu;
				break;
			}
			else if (fu > of_values[1])
			{
				steps[2] = u;
				of_values[2] = fu;
				break;
			}

			u = steps[2] + GOLD * (steps[2] - steps[1]);
			fu = fp->performance(model0.update(1,u,direction));
			iterNum++;

		}
		else if ((steps[2] - u) * (u - ulim) > 0.)
		{
			fu = fp->performance(model0.update(1,u,direction));
			iterNum++;
			if (fu < of_values[2])
			{
				steps[1] = steps[2];
				steps[2] = u;
				u = steps[2] + GOLD * (steps[2] - steps[1]);

				of_values[1] = of_values[2];
				of_values[2] = fu;
				fu = fp->performance(model0.update(1,u,direction));
				iterNum ++;
			}
		}
		else if ((u - ulim) * (ulim * steps[2]) >= 0.)
		{
			u = ulim;
                        fu = fp->performance(model0.update(1,u,direction));
			iterNum ++;
		}
		else
		{
			u = steps[2] + GOLD * (steps[2] - steps[1]);
			fu = fp->performance(model0.update(1,u,direction));
			iterNum ++;
		}

		steps[0] = steps[1]; steps[1] = steps[2]; steps[2] = u;
		of_values[0] = of_values[1]; of_values[1] = of_values[2]; 
		of_values[2] = fu;

	}
	
	// Sorting STEPS in ascending order
	for (i = 0; i < 2; i++)
	{
		if (steps[0] > steps[i+1])
		{
			dum = steps[0];
			steps[0] = steps[i+1];
			steps[i+1] = dum;

			dum = of_values[0];
			of_values[0] = of_values[i+1];
			of_values[i+1] = dum;
		}
	}

	if (steps[1] > steps[2])
	{
		dum = steps[1];
		steps[1] = steps[2];
		steps[2] = dum;

		dum = of_values[1];
		of_values[1] = of_values[2];
		of_values[2] = dum;
	}

	// The line minimization will be performed now using as
	// bracket the 3 first steps given in vector steps 
	// The algorithm is due to Brent

	// initializations
	x = w = v = steps[1];
	fw = fv = fx = of_values[1];

	for (; iterNum <= iterMax; iterNum++)
	{
		xm = .5 * (steps[0] + steps[2]);
		tol1 = tol * Abs(x) + ZEPS;
		tol2 = 2.0 * tol1;
		if (Abs(x - xm) <= (tol2 - .5 * (steps[2] - steps[0])))
		{
			value = fx;
			Model<double> new_model(model0.update(1,x,direction));
			
			return new_model;
		}
						// minimum along a line
		if (Abs(e) > tol1)
		{
			r = (x - w) * (fx - fv);
			q = (x - v) * (fx - fw);
			p = (x - v) * q - (x - w) * r;
			q = 2.0 * (q - r);
			if (q > 0.) p = -q;
			q = Abs(q);
			etemp = e;
			e = d;
			if (Abs(p) >= Abs(.5 * q * etemp) ||
			    p <= q * (steps[0] - x) ||
			    p >= q * (steps[2] - x))
			{	// parabolic fit
				if (x >= xm)
					e = steps[0] - x;
				else
					e = steps[2] - x;
				d = CGOLD * e;
			}
			else
			{
				d = p / q;
				u = x + d;
				if (u - steps[0] < tol2 ||
				    steps[2] - u < tol2)
					d = Abs(tol1) / Sgn(xm-x);
			}
		}
		else
		{
			if (x >= xm)
				e = steps[0] - x;
			else
				e = steps[2] - x;

			d = CGOLD * e;
		}

		if (Abs(d) >= tol1)
			u = x + d;
		else
			u = x + Abs(tol1) / Sgn(d);

		fu  = fp->performance(model0.update(1,u,direction));

		if (fu <= fx)
		{
			if (u >= x)
				steps[0] = x;
			else
				steps[2] = x;
			v = w; w = x; x = u;
			fv = fw; fw = fx; fx = fu;
		}
		else
		{
			if (u < x)
				steps[0] = u;
			else
				steps[2] = u;

			if (fu <= fw || w == x)
			{
				v = w; w = u; fv = fw; fw = fu;
			}

			else if (fu <= fw || v == x || v == w)

			{
				v = u; fv = fu;
			}
		}
	    }

    LineSearch::appendSearchNumber();
//    cout << " Maximum number of iterations reached " << endl;
    return(model0.update(1,x,direction));
}

Model<long> BrentLineSearch::search(Model<long>& model0, Vector<double>& direction, 
		double tol, double delta)
{
    Model<double> temp(model0);
    return search(temp, direction, tol, delta);
}