powell .cpp

高斯赛德尔迭代法

//Numerical Recipes in C 
//一书及所附程序，有完整程序。不过我封装了它的C++版本，可以对但参数或多参数求极值，完整的头文件为： 

#include <iostream.h>
#include <math.h>
#ifndef __OPTIMIZATION_H__ 
#define __OPTIMIZATION_H__ 

////////////////////////////////////////////////////////////////////////// 
// class TOptimization 
// 
// $ 求函数一个或多个参数的最小值 
// 
// 该类默认对一个参数优化，一般只要输入优化参数数目 
// 优化目标函数就可以使用。 
// 
// ...主要代码来自： 
// Numerical Recipes in C++ 
// The Art of Scientific Computing 
// Second Edition 
// William H. Press Saul A. Teukolsky 
// William T. Vetterling Brian P. Flannery 
// 
// 中译本： 
// C++ 数值算法(第二版) 胡健伟 赵志勇 薛运华 等译 
// 电子工业出版社 北京 (2005) 
// 
// Author: Jian Feng 
// Email: [e-mail] 
// Dec. 9, 2006 
// 
////////////////////////////////////////////////////////////////////////// 
// 
// 输入函数： 
// 
// @MaxIterationStep: 最大迭代次数， 默认 1000 
// @ParameterNumbers: 优化参数数目， 默认 1 
// @InitMatrix: 初始化矩阵参数(N*N)， 默认 
// @Tolerance: 容许公差， 默认 1E-7 
// 
// 执行函数： 
// 
// @ExecutePowell: 利用powell方法进行多参数优化 
// @ExecuteBrent: 利用brent方法进行单参数优化 
// 
// 输出函数： 
// 
// @OptimizatedParameters: 优化结果数据 
// @ObjectiveFunctionValue: 目标函数在优化值处的值 
// 
// 使用示例： 
// 
// 1. 单参数 
// double objfun(double a){ 
// double sum = 0; 
// for(int i = 0; i < DataPoints; ++i) 
// sum += SQR(Exps[i] - Theo(a)); 
// } 
// double value 
// TOptimization opt; 
// if(opt.ExecuteBrent(objfun, -10, -1)) opt.OptimizatedParameters(&value); 
// 
// 2. 多参数 
// double objfun(double *a){ 
// double sum = 0; 
// for(int i = 0; i < DataPoints; ++i) 
// sum += SQR(Exps[i] - Theo(a)); 
// } 
// double value[3] 
// TOptimization opt(3); 
// double ival[3] = {-1, 0, 1}; 
// if(opt.ExecutePowell(objfun, ival)) opt.OptimizatedParameters(value); 
// 
namespace{ 
	static int ncom; //公用变量 
	static double *pcom_p; //公用变量 
	static double *xicom_p; //公用变量 
	static double (*nrfunc)(double*); //公用函数指针 
} 

class TOptimization 
{ 
	private: 
	typedef double (*Reff)(double *); 
	typedef double (*Ptrf)(double ); 

	public: 
	TOptimization(int n = 1); 
	~TOptimization(){ FreeMemory(); } 

	//主要方法 
	void ParameterNumbers(int n){ FreeMemory(); num = n; AllocateMemory(); } 

	//利用powell方法对一个或多个参数优化 
	bool ExecutePowell(Reff obj, double *a = 0); 

	//利用brent方法对一个参数优化，需给出参数所在的区间 
	bool ExecuteBrent(Ptrf obj, double vFrom = 0, double vTo = 1); 

	void OptimizatedParameters(double *a){ for(int i=0; i<num; ++i) a[i]=coef[i];} 

	void OptimizatedParameters(double &a){ a = vmin; } 

	//void OptimizatedParameters(double *a){ 
	// if(method) for(int i=0; i<num; ++i) a[i]=coef[i]; 
	// else *a = vmin; 
	//} 

	//其它方法 
	void InitMatrix(double **m) 
	{ 
		for(int i=0; i<num; ++i) 
		for(int j = 0; j<num; ++j) 
		matx[i][j]=m[i][j]; 
		setm = true; 
	} 

	void MaxIterationStep(int s){ ITMAX = s; } 

	void Tolerance(double eps){ ftol = eps; } 

	double ObjectiveFunctionValue()const{ return fret; } 


	private: 

	double brent(double ax, double bx, double cx, Ptrf f, double tol, double &xmin, int &flag); 
	void mnbrak(double &ax, double &bx, double &cx, double &fa, double &fb, double &fc, Ptrf func); 
	void linmin(double *p, double *xi, double &fret, Reff func); 
	bool powell(double *p, double **xi, double ftol, int &iter, double &fret, Reff func); 

	void shft2(double &a, double &b, const double c){ a=b; b=c; } 

	void shft3(double &a, double &b, double &c, const double d){ a=b; b=c; c=d; } 

	double SQR(double x){ return x * x; } 

	void SWAP(double &a, double &b){ double dum=a; a=b; b=dum; } 

	double SIGN(const double &a, const double &b){return b >= 0?(a>=0?a:-a):(a>=0?-a:a);} 

	double MAX(const double &a, const double &b){return b > a ? (b) : (a);} 


	void AllocateMemory(); 

	void FreeMemory(); 

	static double f1dim(double x) 
	{ 
		int j; 
		double *xt = new double [ncom]; 
		//Vec_Dp &pcom=*pcom_p,&xicom=*xicom_p; 
		double *pcom = pcom_p, *xicom = xicom_p; 
		for (j=0;j<ncom;j++) 
		xt[j]=pcom[j]+x*xicom[j]; 
		//delete []xt; 
		double val = nrfunc(xt); 
		delete []xt; 
		return val; 
	} 

	bool setm; //是否设置优化方向初始矩阵 

	int num; //优化参数 
	int ITMAX; //最大迭代数 
	int iter; //实际迭代步数 
	int method; //优化方法 0: 1-D brent, 2: M-D Powell 

	double vmin; //一维优化参数 
	double ftol; //容许差 
	double fret; //目标函数值 


	double *coef; //多维优化参数值 
	double **matx; //多维优化参数方向的初始值 

}; 

////////////////////////////////////////////////////////////////////////// 

inline TOptimization::TOptimization(int n ) 
{ 
	num = n; 
	ftol = 1e-7; 
	ITMAX = 1000; 
	iter = 0; 
	fret = 0.; 
	vmin = 0.; 
	method = 0; 
	setm = false; 

	AllocateMemory(); 
} 


inline void TOptimization::AllocateMemory() 
{ 
	pcom_p = new double [num]; 
	xicom_p = new double [num]; 
	coef = new double [num]; 
	matx = new double *[num]; 
	for(int i = 0; i < num; ++i) 
	{ 
		coef[i] = 0.; 
		matx[i] = new double [num]; 
		for(int j = 0; j < num; ++j) 
		matx[i][j]=(i == j ? 1.0 : 0.0); 
	} 
} 


inline void TOptimization::FreeMemory() 
{ 
	for(int i = 0; i < num; ++i) 
	{ 
		delete []matx[i]; 
	} 
	delete []matx; 
	delete []pcom_p; 
	delete []xicom_p; 
	delete []coef; 
} 

inline bool TOptimization::ExecutePowell(Reff obj, double *a) 
{ 
	method = 1; 
	if(a) 
	for(int i = 0; i < num; ++i) coef[i] = a[i]; 
	return powell(coef, matx, ftol, iter, fret, obj); 
} 


inline bool TOptimization::ExecuteBrent(Ptrf obj, double vFrom, double vTo) 
{ 
	method = 0; 
	int flag; 
	double cx, fa, fb, fc; 
	mnbrak(vFrom,vTo,cx,fa,fb,fc,obj); 
	fret = brent(vFrom,vTo,cx,obj, ftol,vmin, flag); 
	return flag ? true : false; 
} 


inline void TOptimization::mnbrak(double &ax, double &bx, double &cx, double &fa, double &fb, double &fc, Ptrf func) 
{
	const double GOLD=1.618034,GLIMIT=100.0,TINY=1.0e-20; 
	double ulim,u,r,q,fu; 

	fa=func(ax); 
	fb=func(bx); 
	if (fb > fa) 
	{ 
		SWAP(ax,bx); 
		SWAP(fb,fa); 
	} 
	cx=bx+GOLD*(bx-ax); 
	fc=func(cx); 
	while (fb > fc) 
	{ 
		r=(bx-ax)*(fb-fc); 
		q=(bx-cx)*(fb-fa); 
		u=bx-((bx-cx)*q-(bx-ax)*r)/ 
		(2.0*SIGN(MAX(fabs(q-r),TINY),q-r)); 
		ulim=bx+GLIMIT*(cx-bx); 
		if ((bx-u)*(u-cx) > 0.0) 
		{ 
			fu=func(u); 
			if (fu < fc) 
			{ 
				ax=bx; 
				bx=u; 
				fa=fb; 
				fb=fu; 
				return; 
			}
			else if (fu > fb) 
			{ 
				cx=u; 
				fc=fu; 
				return; 
			} 
			u=cx+GOLD*(cx-bx); 
			fu=func(u); 
		}
		else if ((cx-u)*(u-ulim) > 0.0) 
		{ 
			fu=func(u); 
			if (fu < fc) 
			{ 
				shft3(bx,cx,u,cx+GOLD*(cx-bx)); 
				shft3(fb,fc,fu,func(u)); 
			} 
		} 
		else if ((u-ulim)*(ulim-cx) >= 0.0) 
		{ 
			u=ulim; 
			fu=func(u); 
		} 
		else 
		{ 
			u=cx+GOLD*(cx-bx); 
			fu=func(u); 
		} 
		shft3(ax,bx,cx,u); 
		shft3(fa,fb,fc,fu); 
	} 
} 


inline double TOptimization::brent(double ax, double bx, double cx, Ptrf f, double tol, double &xmin, int &flag) 
{ 
	flag = 1; 
	const double CGOLD=0.3819660; 
	const double ZEPS=1.0e-20; 
	int iter; 
	double a,b,d=0.0,etemp,fu,fv,fw,fx; 
	double p,q,r,tol1,tol2,u,v,w,x,xm; 
	double e=0.0; 

	a=(ax < cx ? ax : cx); 
	b=(ax > cx ? ax : cx); 
	x=w=v=bx; 
	fw=fv=fx=f(x); 
	for (iter=0;iter<ITMAX;iter++) 
	{ 
		xm=0.5*(a+b); 
		tol2=2.0*(tol1=tol*fabs(x)+ZEPS); 
		if (fabs(x-xm) <= (tol2-0.5*(b-a))) 
		{ 
			xmin=x; 
			return fx; 
		} 
		if (fabs(e) > tol1) 
		{ 
			r=(x-w)*(fx-fv); 
			q=(x-v)*(fx-fw); 
			p=(x-v)*q-(x-w)*r; 
			q=2.0*(q-r); 
			if(q > 0.0)
			p = -p; 
			q=fabs(q); 
			etemp=e; 
			e=d; 
			if (fabs(p) >= fabs(0.5*q*etemp) || p <= q*(a-x) || p >= q*(b-x)) 
				d=CGOLD*(e=(x >= xm ? a-x : b-x)); 
			else
			{ 
				d=p/q; 
				u=x+d; 
				if (u-a < tol2 || b-u < tol2) 
				d=SIGN(tol1,xm-x); 
			} 
		} 
		else 
		{ 
			d=CGOLD*(e=(x >= xm ? a-x : b-x)); 
		} 
		u=(fabs(d) >= tol1 ? x+d : x+SIGN(tol1,d)); 
		fu=f(u); 
		if (fu <= fx) 
		{ 
			if (u >= x) a=x; else b=x; 
			shft3(v,w,x,u); 
			shft3(fv,fw,fx,fu); 
		} 
		else 
		{ 
			if (u < x) a=u; else b=u; 
			if (fu <= fw || w == x) 
			{ 
				v=w; 
				w=u; 
				fv=fw; 
				fw=fu; 
			} 
			else if (fu <= fv || v == x || v == w) 
			{ 
				v=u; 
				fv=fu; 
			} 
		} 
	} 
	flag = 0; 
	xmin=x; 
	return fx; 
} 


inline void TOptimization::linmin(double *p, double *xi, double &fret, Reff func) 
{ 
	int j, flag; 
	const double TOL=1.0e-8; 
	double xx,xmin,fx,fb,fa,bx,ax; 

	int n=num; 
	ncom=n; 
	//pcom_p=new Vec_Dp(n); 
	//xicom_p=new Vec_Dp(n); 
	nrfunc=func; 
	//Vec_Dp &pcom=*pcom_p,&xicom=*xicom_p; 
	double *pcom = pcom_p, *xicom = xicom_p; 
	for (j=0;j<n;j++) 
	{ 
		pcom[j]=p[j]; 
		xicom[j]=xi[j]; 
	} 
	ax=0.0; 
	xx=1.0; 
	mnbrak(ax,xx,bx,fa,fx,fb,f1dim); 
	fret=brent(ax,xx,bx,f1dim,TOL,xmin, flag); 
	for (j=0;j<n;j++)
	{ 
		xi[j] *= xmin; 
		p[j] += xi[j]; 
	} 
	//delete xicom_p; 
	//delete pcom_p; 
} 


inline bool TOptimization::powell(double *p, double **xi, double ftol, int &iter, double &fret, Reff func) 
{ 
	const int ITMAX=500; 
	const double TINY=1.0e-20; 
	int i,j,ibig; 
	double del,fp,fptt,t; 

	int n=num; 
	//Vec_Dp pt(n),ptt(n),xit(n); 
	double *pt, *ptt, *xit; 
	for(i = 0; i < n; ++i) 
	{ 
		pt = new double [n]; 
		ptt = new double [n]; 
		xit = new double [n]; 
	} 
	fret=func(p); 
	for (j=0;j<n;j++) pt[j]=p[j]; 
	for (iter=0;;++iter)
	{ 
		fp=fret; 
		ibig=0; 
		del=0.0; 
		for (i=0;i<n;i++)
		{ 
			for (j=0;j<n;j++) xit[j]=xi[j][i]; 
			fptt=fret; 
			linmin(p,xit,fret,func); 
			if (fptt-fret > del)
			{ 
				del=fptt-fret; 
				ibig=i+1; 
			} 
		} 
		if (2.0*(fp-fret) <= ftol*(fabs(fp)+fabs(fret))+TINY) 
		{ 
			delete []pt; 
			delete []ptt; 
			delete []xit; 
			return true; 
		} 
		if (iter == ITMAX) 
		{ 
			delete []pt; 
			delete []ptt; 
			delete []xit; 
			return false; 
			//cerr<<"powell exceeding maximum iterations."; 
		} 
		for (j=0;j<n;j++) 
		{ 
			ptt[j]=2.0*p[j]-pt[j]; 
			xit[j]=p[j]-pt[j]; 
			pt[j]=p[j]; 
			
		} 
		fptt=func(ptt); 
		if (fptt < fp) 
		{ 
			t=2.0*(fp-2.0*fret+fptt)*SQR(fp-fret-del)-del*SQR(fp-fptt); 
			if (t < 0.0) 
			{ 
				linmin(p,xit,fret,func); 
				for (j=0;j<n;j++) 
				{ 
					xi[j][ibig-1]=xi[j][n-1]; 
					xi[j][n-1]=xit[j]; 
				} 
			} 
		} 
	} 
} 

#endif