simplex.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 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();
	    }
	}
    }
}