//============================================================// 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 Simplex class// The flexible polyhedron algorithm// author:  Wenceslau Gouveia// modified:  H. Lydia Deng, 02/23/94,  /03/14/94//=============================#include <Simplex.hh>void Simplex::reset(Vector<double>& lambda){    int i, imin = 0;    for(i = 1; i < nd + 1; i++)	if((*fv)[i] < (*fv)[imin])       	    imin = i;    for(i = 0; i < nd + 1; i++) {	if(i == imin)	    continue;	int j = (i > imin) ? (i - 1) : i;	models[i][j] = models[imin][j] + lambda[i];	(*fv)[i] = fp->performance(models[i]);    }    formPsum();}Simplex::Simplex(ObjectiveFunction* f, Model<double>* m, int iter, double a, double b, double y): NonQuadraticOptima(){    iterMax	=	iter;    alpha 	= 	a;    beta 	=	b;    gamma 	=	y;    fp		= 	f;    models	=	m;    nd		=	models[0].modSize();    psum 	= 	new 	Vector<double>(nd);    fv 		= 	new 	Vector<double>(nd+1);    for (int i = 0; i < nd+1; i++) (*fv)[i] = fp->performance(models[i]);}Simplex::Simplex(ObjectiveFunction* f, Model<double>* m, int iter, double a, double b, 	double y, int verb): NonQuadraticOptima(verb){    iterMax	=	iter;    alpha 	= 	a;    beta 	=	b;    gamma 	=	y;    fp		= 	f;    models	=	m;    nd		=	models[0].modSize();    psum 	= 	new 	Vector<double>(nd);    fv 		= 	new 	Vector<double>(nd+1);    for (int i = 0; i < nd+1; i++) (*fv)[i] = fp->performance(models[i]);}Simplex::~Simplex(){	delete  fv;  	delete  psum;}void Simplex::formPsum(){    for(int i = 0; i < nd; i++) 	for (int j = 0; j<nd+1; j++)	    (*psum)[i] += models[j][i];}double Simplex::tryNewPoint(int ihigh, const double lever){    Model<double> trypoint(nd);    double f1 	= 	(1.0 - lever)/nd;    double f2 	= 	f1 - lever;    trypoint	=	models[ihigh].update(-f2, f1, psum[0]);    double yt 	= 	fp->performance(trypoint);    if(yt < (*fv)[ihigh]) {	(*fv)[ihigh] 	= 	yt;	psum[0]		+=	(trypoint - models[ihigh]).modParam();	models[ihigh] 	= 	trypoint;    }    return yt;}Model<double> Simplex::optimizer(const double atol){    Vector<double> atl(nd);    for(int i = 0; i < nd; i++)	atl[i] = atol;    return optimizer(atl);}Model<double> Simplex::optimizer(Vector<double>& atollist){    Model<double>	bestModel(nd);    double 		fmax;    int 		imax = 0;    int 		imin = 0;    int			i, inext;    formPsum();    for(;;) {	bestModel	=	models[0];	value		=	(*fv)[0];	imin 		= 	0;	imax 		= 	0;	inext 		= 	1;	if((*fv)[0] < (*fv)[1]) {	    imax = 1;	    inext = 0;	}	for(i = 1; i < nd + 1; i++){ 		if((*fv)[i] < value) { 			value 		= 	(*fv)[i]; 			bestModel 	= 	models[i]; 			imin 		= 	i; 		}	    	if((*fv)[i] > (*fv)[imax]) {		    inext 		= 	imax;		    imax 		= 	i;	    	}	    }	fmax = (*fv)[imax];	double xerror = 0.0;	double cmax, cmin;	for(i = 0; i < nd; i++) {	    cmax	=	models[0][i];	    cmin	=	models[0][i];	    for (int j = 1; j < nd+1; j++) {		 cmax 	= 	Max(cmax, models[j][i]);		 cmin 	= 	Min(cmin, models[j][i]);	     }	    xerror 	= 	Max(xerror, Abs(cmax - cmin)/atollist[i]);	}	if (isVerbose) cerr << "error: "<<xerror<<endl;	if(xerror < 1.0 || fp->iterations() > iterMax)		    return bestModel;	double ytry 	= 	tryNewPoint(imax, -alpha);	if(ytry < value)	    ytry = tryNewPoint(imax, gamma);	else if(ytry > (*fv)[inext]) {	    double save = (*fv)[imax];	    ytry = tryNewPoint(imax, beta);	    if(ytry > save) {		for(i = 0; i < nd + 1; i++)		    if(i != imin)			models[i] = 0.5*(models[i]+models[imin]).modParam();		formPsum();	    }	}    }}