//============================================================// 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;}