📄 optimi.h
字号:
#include<math.h>
#include<STDIO.h>
#include<STDLIB.h>
#define SHFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d);
#define SIGN(a,b) ((b)>=0.0?fabs(a):-fabs(a))
#define FMAX(a,b) (a)>(b)?(a):(b)
#define SQR(a) ((a)==0.0?0.0:a*a)
#define ZEPS 1.0e-10
#define NR_END 1
#define FREE_ARG char*
#define TOL 0.001//1.0
#define ITMAX 100 //减小了
#define GLIMIT 100.0 //GLIMIT is the maximum magni cation allowed for a parabolic- t step.
#define GOLD 1.618034 //Here GOLD is the default ratio by which successive intervals are magni ed;
#define CGOLD 0.3819660
#define ITMAX2 200
#define TINY 1.0e-25
int ncom; // De ned in linmin.
double *pcom,*xicom,(*nrfunc)(double[]);
void free_vector(double *v, long nl, long nh)
{
free((FREE_ARG)(v+nl-NR_END));
}
double *vector(long nl,long nh)
{
double*v;
v=(double*)malloc((size_t)((nh-nl+1+NR_END)*sizeof(double)));
return v-nl+NR_END;
}
/*******************************************************************************************
Given a function func, and given distinct initial points ax and bx, this routine searches in
the downhill direction (de ned by the function as evaluated at the initial points) and returns
new points ax, bx, cx that bracket a minimum of the function. Also returned are the function
values at the three points, fa, fb, and fc.
*******************************************************************************************/
void mnbrak(double *ax,double*bx,double*cx,double*fa,double*fb,double*fc,double(*func)(double))
{
double ulim,u,r,q,fu,dum;
*fa=(*func)(*ax);
*fb=(*func)(*bx);
//Switch roles of a and b so that we can go downhill in the direction from a to b.
if(*fb>*fa)
{
SHFT(dum,*ax,*bx,dum)
SHFT(dum,*fb,*fa,dum)
}
//First guess for c.
*cx=(*bx)+GOLD*(*bx-*ax);
*fc=(*func)(*cx);
//Keep returning here until we bracket.
while(*fb>*fc)
{
//Compute u by parabolic extrapolation from a; b; c.
//TINY is used to prevent any possible division by zero.
r=(*bx-*ax)*(*fb-*fc);
q=(*bx-*cx)*(*fb-*fa);
u=(*bx)-((*bx-*cx)*q-(*bx-*ax)*r)/(2.0*SIGN(FMAX(fabs(q-r),TINY),q-r));
ulim=(*bx)+GLIMIT*(*cx-*bx);
//We won't go farther than this. Test various possibilities:
if((*bx-u)*(u-*cx)>0.0)
{
//Parabolic u is between b and c: try it.
fu=(*func)(u);
//Got a minimum between b and c.
if(fu<*fc)
{
*ax=(*bx);
*bx=u;
*fa=(*fb);
*fb=fu;
return;
}
else if (fu>*fb)
{
//Got a minimum between between a and u.
*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)
{
SHFT(*bx,*cx,u,*cx+GOLD*(*cx-*bx))
SHFT(*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);
}
//Eliminate oldest point and continue.
SHFT(*ax,*bx,*cx,u)
SHFT(*fa,*fb,*fc,fu)
}
}
////////////////////////////////////////////////////////////////////////////
double f1dim(double x)
{
int j;
double f,*xt;
xt=vector(1,ncom);
for(j=1;j<=ncom;j++)
xt[j]=pcom[j]+x*xicom[j];
f=(*nrfunc)(xt);
free_vector(xt,1,ncom);
return f;
}
/*******************************************************************************************
double brent(double ax,double bx,double cx,double(*f)(double),double tol,double*xmin)
Given a function f, and given a bracketing triplet of abscissas ax, bx, cx (such that bx is
between ax and cx, and f(bx) is less than both f(ax) and f(cx)), this routine isolates
the minimum to a fractional precision of about tol using Brent's method. The abscissa of
the minimum is returned as xmin, and the minimum function value is returned as brent, the
returned function value.
*******************************************************************************************/
double brent(double ax,double bx,double cx,double(*f)(double),double tol,double*xmin)
{
int iter;
double a,b,d,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,x,xm;
//This will be the distance moved on the step before last.
double e=0.0;
//a and b must be in ascending order,but input abscissas need not be Initializations...
a=(ax<cx?ax:cx);
b=(ax>cx?ax:cx);
x=w=v=bx;
fw=fv=fx=(*f)(x);
//Main program loop.
for(iter=1;iter<=ITMAX;iter++)
{
xm=0.5*(a+b);
tol2=2.0*(tol1=tol*fabs(x)+ZEPS);
//Test for done here.
if(fabs(x-xm)<=(tol2-0.5*(b-a)))
{
*xmin=x;
return fx;
}
//Construct a trial parabolic t.
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));
//The above conditions determine the acceptability of the parabolic t.
//Here we take the golden section step into the larger of the two segments.
else
{
//Take the parabolic step.
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);
//Now decide what to do with our function evaluation. Housekeeping follows:
if(fu<=fx)
{
if(u>=x)
a=x;
else
b=x;
SHFT(v,w,x,u)
SHFT(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;
}
}
} //Done with housekeeping. Back for another iteration.
//Never get here.
*xmin=x;
return fx;
}
/*******************************************************************************************
void linmin(double p[], double xi[], int n, double *fret, double (*func)(double []))
Given an n-dimensional point p[1..n] and an n-dimensional direction xi[1..n], moves and
resets p to where the function func(p) takes on a minimum along the direction xi from p,
and replaces xi by the actual vector displacement that p was moved. Also returns as fret
the value of func at the returned location p. This is actually all accomplished by calling the
routines mnbrak and brent.
*******************************************************************************************/
extern int ncom;
extern double *pcom,*xicom,(*nrfunc)(double[]);
void linmin(double p[],double xi[],int n,double*fret,double(*func)(double[]))
{
/* Use Function
double brent(double ax, double bx, double cx,
double (*f)(double), double tol, double *xmin);
double f1dim(double x);
void mnbrak(double *ax, double *bx, double *cx, double *fa, double *fb,double *fc, double (*func)(double));
*/
int j;
double xx,xmin,fx,fb,fa,bx,ax;
//De ne the global variables.
ncom=n;
pcom=vector(1,n);
xicom=vector(1,n);
nrfunc=func;
for(j=1;j<=n;j++)
{
pcom[j]=p[j];
xicom[j]=xi[j];
}
//Initial guess for brackets.
ax=0.0;
xx=1.0;
mnbrak(&ax,&xx,&bx,&fa,&fx,&fb,f1dim);
*fret=brent(ax,xx,bx,f1dim,TOL,&xmin);
//Construct the vector results to return.
for(j=1;j<=n;j++)
{
xi[j]*=xmin;
p[j]+=xi[j];
}
free_vector(xicom,1,n);
free_vector(pcom,1,n);
}
/*******************************************************************************************
void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret,double (*func)(double []))
Minimization of a function func of n variables.
Input:
consists of an initial starting point p[1..n];
an initial matrix xi[1..n][1..n], whose columns contain the initial set of directions (usually the n unit vectors);
ftol, the fractional tolerance in the function value such that failure to decrease by more than this amount on one iteration signals doneness.
Output:
p is set to the best point found,
xi is the then-current direction set,
fret is the returned function value at p,
and iter is the number of iterations taken.
The routine linmin is used.
/*******************************************************************************************/
void powell(double p[],double**xi,int n,double ftol,int*iter,double*fret,double(*func)(double[]))
{
/* Use Function
void linmin(double p[], double xi[], int n, double *fret,
double (*func)(double []));
int i,ibig,j;
double del,fp,fptt,t,*pt,*ptt,*xit;
*/
int i,ibig,j;
double del,fp,fptt,t,*pt,*ptt,*xit;
int flag=0;
pt=vector(1,n);
ptt=vector(1,n);
xit=vector(1,n);
*fret=(*func)(p);
//save inital point
for(j=1;j<=n;j++)
pt[j]=p[j];
for(*iter=1;;++(*iter))
{
fp=(*fret);
ibig=0;
//Will be the biggest function decrease.
del=0.0;
//In each iteration, loop over all directions in the set.
for(i=1;i<=n;i++)
{
//Copy the direction,
for(j=1;j<=n;j++)
xit[j]=xi[j][i];
fptt=(*fret);
///////////////////////////////////////////
printf("i=%d \n",i);
///////////////////////////////////////////
//minimize along it,
linmin(p,xit,n,fret,func);
//and record it if it is the largest decrease so far
if(fptt-(*fret)>del)
{
del=fptt-(*fret);
ibig=i;
}
}
//Termination criterion.
if(2.0*(fp-(*fret))<=ftol*(fabs(fp)+fabs(*fret))+TINY)
{
free_vector(xit,1,n);
free_vector(ptt,1,n);
free_vector(pt,1,n);
return;
}
if(*iter==ITMAX2)
printf("powell exceeding maximum iterations.");
//Construct the extrapolated point and the average direction moved.
//Save the old starting point.
for(j=1;j<=n;j++)
{
ptt[j]=2.0*p[j]-pt[j];
xit[j]=p[j]-pt[j];
pt[j]=p[j];
}
//Function value at extrapolated point.
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)
{
//Move to the minimum of the new direction,and save the new direction.
linmin(p,xit,n,fret,func);
for(j=1;j<=n;j++)
{
xi[j][ibig]=xi[j][n];
xi[j][n]=xit[j];
}
}
}
} //Back for another iteration.
}
double ** matrix(long nrl, long nrh, long ncl, long nch)
{
long i,nrow=nrh-nrl+1,ncol=nch-ncl+1;
double**m;
// allocate pointer storows
m=(double**)malloc((size_t)((nrow+NR_END)*sizeof(double*)));
if(!m) printf("allocationfailure1inmatrix()");
m+=NR_END;
m-=nrl;
//allocate rows and set pointer stothem
m[nrl]=(double*)malloc((size_t)((nrow*ncol+NR_END)*sizeof(double)));
if(!m[nrl])
printf("allocationfailure2inmatrix()");
m[nrl]+=NR_END;
m[nrl]-=ncl;
for(i=nrl+1;i<=nrh;i++)
m[i]=m[i-1]+ncol;
return m;
}
void free_matrix(double **m, long nrl, long nrh, long ncl, long nch)
{
free((FREE_ARG)(m[nrl]+ncl-NR_END));
free((FREE_ARG)(m+nrl-NR_END));
}⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -