powell.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 for the Powel's conjugate gradient class 
// Powell's derivative-free conjugate gradient class
// author:  Wenceslau Gouveia
// modified: H. Lydia Deng, 02/23/94, 03/15/94
//=========================

#include <Powell.hh>
#include <DensMatrix.hh>

static const char*  myNameIs =  "Powell Method";

const char* PowellOptima::className() const { 
	return (myNameIs);
}

PowellOptima::PowellOptima(LineSearch *p,  int it, double eps, double chang, double delt)
: LineSearchOptima(p)
{
    	iterMax 	= 	it;
    	tol 		= 	eps;
	change 		= 	chang;
	delta 		= 	delt;
    	iterNum 	= 	0;
}

PowellOptima::PowellOptima(LineSearch *p,  int it, double eps, double chang, double delt, int verb)
: LineSearchOptima(p,verb)
{
    	iterMax 	= 	it;
    	tol 		= 	eps;
	change 		= 	chang;
	delta 		= 	delt;
    	iterNum 	= 	0;
}

Model<double> PowellOptima::optimizer(Model<double>& model0)
{
    int n= model0.modSize();
    DensMatrix<double> direction(n,n);
    Model<double> model(n);
    Vector<double> values(n+2);			         // OF evaluations
    Model<double> model_init(n), model_end(n);   	// input to line search
    Vector<double> s(n);				// possible new search direction
    Vector<double> p(n);				// possible new search direction
    Vector<double> distance(model0.modParam());	   // used in the stopping criterion
    Vector<double> diff_of(n);			   // difference in OF evaluations
    Vector<double> temp(n);
    double deltak;					// used to select new search
    int i, j;

	// initializing Matrix direction
	direction = 0.;
	for (i=0; i<n; i++) direction[i][i] = 1.;

	// beginning of the Powel's method
	do
	{			// line search for the initial guess
		values[0] = ls->evaluate(model0); 
		iterNum++;
		p = direction.rowVector(0);
		model_init = ls->search(model0, p, tol, delta);

		values[1] = ls->currentValue(); 

		for (i = 2; i <= n; i++)
		{
			p = direction.rowVector(i-1); 
			model =	ls->search(model_init, p, tol, delta);
			values[i] = ls->currentValue();
			model_init = model;
		}

		// model stores the last outcome of the line search
		model_end = 2. * model - model0;
		values[n+1] = ls->evaluate(model_end);
		iterNum++;

		// computing differences in the OF evaluations
		int index_max;
		for (i = 1; i <= n; i++)
			diff_of[i-1] = values[i-1] - values[i];
		index_max = diff_of.indexMax();
		deltak = diff_of.max();

		// decide on new direction of search
		if (values[n+1] >= values[0] ||
		    ((values[0]-2.*values[n]+values[n+1]) *
		     (values[0]-values[n]-deltak) * 
		     (values[0]-values[n]-deltak) >= 
		     .5*deltak*(values[0]-values[n+1]) *
		     (values[0]-values[n+1])))
		{	// keep all search directions

			// deciding the new starting point
			if (values[n] >= values[n+1])
 				model0 = model_end;
			else
				model0 = model;
		}
		else
		{	// just one search direction will be replaced
			// this is the one pointed by index_max
			s = model.modParam() - model0.modParam();

			// the next initial point model0 is defined as
			model0 = ls->search(model, s, tol, delta);

			// and replacing the row index_max in matrix directions
			for (i = index_max + 1; i < n; i++)
				for (j = 0; j < n; j++)
					direction[i-1][j] = direction[i][j];

			// and adding the new search direction
			for (j = 0; j < n; j++)	direction[n-1][j] = s[j];
		}

		// checking the convergence of the method
		temp	=	model0.modParam();
		distance -= temp;

		double err =  (distance*distance) / (temp*temp);
		if (isVerbose) cerr << "the "<<iterNum<<"th iteration residue: "<<err<<endl;
		NonQuadraticOptima::appendResidue (err);
		distance = model0.modParam();

	} while ( residue->last() > change && iterNum < iterMax);

	return(model0);
}


Model<long> PowellOptima::optimizer(Model<long>& model0)
{
    Model<double> temp(model0);
    temp = optimizer(temp);
    Model<long> optm(temp);
    return optm;
}