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