cubicls.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 cubic line search class
// Armijo and Goldstein's line search algorithm
// author:  Doug Hart, Adapted to the COOOL library by Wenceslau Gouveia
// Modified to fit into new classes.  H. Lydia Deng, 02/21/94, 03/15/94
//========================================================================


#include <limits.h>
#include <float.h>
#include <CubicLS.hh>

CubicLineSearch::CubicLineSearch(ObjectiveFunction* f, int it) 
: LineSearch(f)
{	iterMax	=	it;
}

CubicLineSearch::CubicLineSearch(ObjectiveFunction* f, int it, Vector<double>* interval)
: LineSearch(f,interval)
{	iterMax	=	it;
}

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

// Code for the Cubic Line Search
Model<double> CubicLineSearch::search(Model<double>& current_solution, Vector<double>& p, double slope, double lambda)
{
        int tst = 0;					// useful vars
	double alpha2=0, alpha_tmp=0, alpha_prev=0;
	double alpha_prev2=0, alpha=0;
	double f1=0, f2=0, fprev=0;
	double a=0, b=0;
	double c=0, cm11=0, cm12=0, cm21=0, cm22=0;
	double disc=0;
	double new_m=0, old_m=0;
	Model<double> new_solution(current_solution);

	old_m = fp->performance(current_solution);	// Evaluation at the  
							// current solution
	iterNum = 0; iterNum++;				// iteration counter
	alpha = 1.;					// updating step

	new_solution = current_solution.update(1,alpha,p);	// updating
	new_m = fp->performance(new_solution);     	// Evaluation at the 
							// new solution
	iterNum++;

	// Implementing Goldstein's test for alpha too small
	while (new_m < old_m + (1. - lambda)*alpha*slope && iterNum< iterMax)
	{
		alpha *= 3;
		new_solution = current_solution.update(1, alpha, p);
		new_m = fp->performance(new_solution);
		iterNum++;
	}
	if (iterNum == iterMax)
		cerr << "Alpha over flowed! \n";

	// Armijo's test for alpha too large
	alpha_prev = alpha; // H.L. Deng, 6/13/95
	while (new_m > old_m + lambda*alpha*slope && iterNum < iterMax)
	{
		alpha2 = alpha * alpha;
		f1 = new_m - old_m - slope * alpha;

		if (tst == 0)
		{
			alpha_tmp = -slope * alpha2 / (f1 * 2.);
						// tentative alpha

			tst = 1;
		}
		else
		{
			alpha_prev2 = alpha_prev * alpha_prev;
			f2 = fprev - old_m - alpha_prev * slope;

			c = 1. / (alpha - alpha_prev);
			cm11 = 1. / alpha2;
			cm12 = -1. / alpha_prev2;
			cm21 = -alpha_prev / alpha2;
			cm22 = alpha / alpha_prev2;

			a = c * (cm11 * f1 + cm12 * f2);
			b = c * (cm21 * f1 + cm22 * f2);
			disc = b * b - 3. * a * slope;

			if ((Abs(a) > FLT_MIN) && (disc > FLT_MIN))
				alpha_tmp = (-b + sqrt(disc)) / (3. * a);
			else
				alpha_tmp = slope * alpha2 / (2. * f1);

			if (alpha_tmp >= .5 * alpha)
				alpha_tmp = .5 * alpha;
		}
		alpha_prev = alpha;
		fprev = new_m;

		if (alpha_tmp < .1 * alpha)
			alpha *= .1;
		else
			alpha = alpha_tmp;

		new_solution = current_solution.update(1, alpha, p);
		new_m = fp->performance(new_solution);
		iterNum++;
	}
	if (iterNum == iterMax)
		cerr << "Alpha under flowed! \n";
	
        value = new_m;
	appendSearchNumber();
	return (new_solution);				// # of iterations
}

Model<long> CubicLineSearch::search(Model<long>& current_solution, Vector<double>& p, double slope, double lambda)
{
    Model<double> temp(current_solution);
    temp = search(temp, p, slope, lambda);
    Model<long> optm(temp);
    return optm;
}