?? extremum.inl
字號:
for(i=1; i<2*n; i++)
if(xx(n,i)>fr)
{
r = i;
fr = xx(n,i);
}
g = 0;
j = 0;
fg = xx(n,0);
if(r==0)
{
g = 1;
j = 1;
fg = xx(n,1);
}
for(i=j+1; i<2*n; i++)
if(i!=r)
if(xx(n,i)>fg)
{
g = i;
fg = xx(n,i);
}
for(i=0; i<n; i++)
{
xf[i] = 0.0;
for(j=0; j<2*n; j++)
if(j!=r) xf[i]=xf[i]+xx(i,j)/(2.0*n-1.0);
xt[i]=(1.0+alpha)*xf[i]-alpha*xx(i,r);
}
jt = 1;
while(jt==1)
{
jt = 0;
z = FunctionValueTarget(xt,n);
while(z>fg)
{
for(i=0; i<n; i++) xt[i]=(xt[i]+xf[i])/2.0;
z = FunctionValueTarget(xt,n);
}
j=0;
for(i=0; i<n; i++)
{
if(a[i]>xt[i])
{
xt[i]=xt[i]+0.000001;
j=1;
}
if(b[i]<xt[i])
{
xt[i]=xt[i]-0.000001;
j=1;
}
}
if(j!=0) jt=1;
else
{
ConditionValue(n,m,xt,c,d,w);
j=0;
kt=1;
while((kt==1)&&(j<m))
{
if((c[j]<=w[j])&&(d[j]>=w[j])) j=j+1;
else kt=0;
}
if(j<m)
{
for(i=0; i<n; i++) xt[i]=(xt[i]+xf[i])/2.0;
jt=1;
}
}
}
for(i=0; i<n; i++) xx(i,r)=xt[i];
xx(n,r) = z;
fr = 0.0;
fg = 0.0;
for(j=0; j<2*n; j++)
{
fj = xx(n,j);
fr = fr + fj / (2.0*n);
fg = fg + fj * fj;
}
fr = (fg - 2.0 * n * fr * fr) / (2.0 * n - 1.0);
if(fr>=eps)
{
kk = kk + 1;
if(kk<k) it = 1;
}
}
for(i=0; i<n; i++)
{
x[i] = 0.0;
for(j=0; j<2*n; j++) x[i]=x[i]+xx(i,j)/(2.0*n);
}
z = FunctionValueTarget(x,n);
x[n] = z;
return(kk);
}
//確定1維函數(shù)極小值點所在區(qū)間
template <class _Ty>
void MinimizerInterval(_Ty& ax, _Ty& bx, _Ty& cx, _Ty& fa, _Ty& fb, _Ty& fc)
{
_Ty r,q,dum,gold = 1.0 + GoldNo;
_Ty u,ulim,fu,tiny = 1e-20;
int glimit = 100;
fa = FunctionValue(ax);
fb = FunctionValue(bx);
if (fb > fa)
{
dum = ax;
ax = bx;
bx = dum;
dum = fb;
fb = fa;
fa = dum;
}
cx = bx + gold * (bx - ax);
fc = FunctionValue(cx);
while (fb >= fc)
{
r = (bx - ax) * (fb - fc);
q = (bx - cx) * (fb - fa);
dum = q - r;
if (Abs(dum) < tiny)
{
dum = tiny;
}
u = bx - ((bx - cx) * q - (bx - ax) * r) / (2 * dum);
ulim = bx + glimit * (cx - bx);
if ((bx - u) * (u - cx) > 0)
{
fu = FunctionValue(u);
if (fu < fc)
{
ax = bx;
fa = fb;
bx = u;
fb = fu;
return;
}
else
{
if (fu > fb)
{
cx = u;
fc = fu;
return;
}
}
u = cx + gold * (cx - bx);
fu = FunctionValue(u);
}
else
{
if ((cx - u) * (u - ulim) > 0)
{
fu = FunctionValue(u);
if (fu < fc)
{
bx = cx;
cx = u;
u = cx + gold * (cx - bx);
fb = fc;
fc = fu;
fu = FunctionValue(u);
}
}
else
{
if ((u - ulim) * (ulim - cx) > 0 || FloatEqual((u - ulim) * (ulim - cx),0))
{
u = ulim;
fu = FunctionValue(u);
}
else
{
u = cx + gold * (cx - bx);
fu = FunctionValue(u);
}
}
}
ax = bx;
bx = cx;
cx = u;
fa = fb;
fb = fc;
fc = fu;
}
}
//確定一維函數(shù)極小值點所在區(qū)間(重載)
template <class _Ty>
int MinimizerInterval(_Ty& ax, _Ty& bx, _Ty& cx, int& sign,
_Ty b, _Ty& fa, _Ty& fb, _Ty& fc)
{
_Ty r,q,dum,gold = GoldNo + 1.0;
_Ty u,ulim,fu,tiny = 1e-20;
int glimit = 100;
fa = FunctionValue(ax);
fb = FunctionValue(bx);
if (fb > fa)
{
sign = -1;
fa *= sign;
fb *= sign;
}
cx = bx + gold * (bx - ax);
fc = sign * FunctionValue(cx);
while (fb >= fc)
{
r = (bx - ax) * (fb - fc);
q = (bx - cx) * (fb - fa);
dum = q - r;
if (Abs(dum) < tiny)
{
dum = tiny;
}
u = bx - ((bx - cx) * q - (bx - ax) * r) / (2 * dum);
ulim = bx + glimit * (cx - bx);
if(u > b) //超出f(x)右邊界
{
cout<<" The MinimizerInterval exceed right border. Program is over!"<<endl;
return 0;
}
if ((bx - u) * (u - cx) > 0)
{
fu = sign * FunctionValue(u);
if (fu < fc)
{
ax = bx;
fa = fb;
bx = u;
fb = fu;
return 1;
}
else
{
if (fu > fb)
{
cx = u;
fc = fu;
return 1;
}
}
u = cx + gold * (cx - bx);
fu = sign * FunctionValue(u);
}
else
{
if ((cx - u) * (u - ulim) > 0)
{
fu = sign * FunctionValue(u);
if (fu < fc)
{
bx = cx;
cx = u;
u = cx + gold * (cx - bx);
fb = fc;
fc = fu;
fu = sign * FunctionValue(u);
}
}
else
{
if ((u - ulim) * (ulim - cx) > 0 || FloatEqual((u - ulim) * (ulim - cx),0))
{
u = ulim;
fu = sign * FunctionValue(u);
}
else
{
u = cx + gold * (cx - bx);
fu = sign * FunctionValue(u);
}
}
}
ax = bx;
bx = cx;
cx = u;
fa = fb;
fb = fc;
fc = fu;
}
return 1;
}
//1維函數(shù)極小值的黃金分割法
template <class _Ty>
_Ty ExtremumGold1D(_Ty ax, _Ty bx, _Ty cx, int sign, _Ty tol, _Ty& xmin)
{
_Ty x0,x1,x2,x3,f0,f1,f2,f3,r = GoldNo;
_Ty c = 1.0 - GoldNo;
x0 = ax;
x3 = cx;
if (Abs(cx - bx) > Abs(bx - ax))
{
x1 = bx;
x2 = bx + c * (cx - bx);
}
else
{
x2 = bx;
x1 = bx - c * (bx - ax);
}
f1 = sign * FunctionValue(x1);
f2 = sign * FunctionValue(x2);
while (Abs(x3 - x0) > tol * (Abs(x1) + Abs(x2)))
{
if (f2 < f1)
{
x0 = x1;
x1 = x2;
x2 = r * x1 + c * x3;
f0 = f1;
f1 = f2;
f2 = sign * FunctionValue(x2);
}
else
{
x3 = x2;
x2 = x1;
x1 = r * x2 + c * x0;
f3 = f2;
f2 = f1;
f1 = sign * FunctionValue(x1);
}
}
if (f1 < f2)
{
xmin = x1;
if (sign == -1)
f1 = sign * f1;
return f1;
}
else
{
xmin = x2;
if (sign == -1)
f2 = sign * f2;
return f2;
}
}
//1維函數(shù)極小值點的不用導數(shù)布倫特(Brent)法
template <class _Ty>
_Ty ExtremumBrentNonDerivative1D(_Ty ax, _Ty bx, _Ty cx, int sign,
_Ty tol, _Ty& xmin)
{
int done,iter,itmax = 100;
_Ty d,fu,r,q,p,xm,tol1,tol2,a,b,cgold = 1.0 - GoldNo;
_Ty u,etemp,dum,v,w,x,e,fx,fv1,fw,zeps = DOUBLEERROR;
a = ax;
if (cx < ax)
{
a = cx;
}
b = ax;
if (cx > ax)
{
b = cx;
}
v = bx;
w = v;
x = v;
e = 0.0;
fx = sign * FunctionValue(x);
fv1 = fx;
fw = fx;
for (iter = 1; iter<=itmax; iter++)
{
xm = 0.5 * (a + b);
tol1 = tol * Abs(x) + zeps;
tol2 = 2.0 * tol1;
if (Abs(x - xm) <= tol2 - 0.5 * (b - a))
{
break;
}
done = -1;
if (Abs(e) > tol1)
{
r = (x - w) * (fx - fv1);
q = (x - v) * (fx - fw);
p = (x - v) * q - (x - w) * r;
q = 2.0 * (q - r);
if (q > 0.0)
{
p = -p;
}
q = Abs(q);
etemp = e;
e = d;
dum = Abs(0.5 * q * etemp);
if (Abs(p) < dum && p > q * (a - x) && p < q * (b - x))
{
d = p / q;
u = x + d;
if (u - a < tol2 || b - u < tol2)
{
d = Abs(tol1) * Sgn(xm - x);
}
done = 0;
}
}
if (done)
{
if (x > xm || FloatEqual(x,xm))
{
e = a - x;
}
else
{
e = b - x;
}
d = cgold * e;
}
if (Abs(d) > tol1 || FloatEqual(d,tol1))
{
u = x + d;
}
else
{
u = x + Abs(tol1) * Sgn(d);
}
fu = sign * FunctionValue(u);
if (fu < fx || FloatEqual(fu,fx))
{
if (u > x || FloatEqual(u,x))
{
a = x;
}
else
{
b = x;
}
v = w;
fv1 = fw;
w = x;
fw = fx;
x = u;
fx = fu;
}
else
{
if (u < x)
{
a = u;
}
else
{
b = u;
}
if (fu < fw || FloatEqual(fu,fw) || FloatEqual(w,x))
{
v = w;
fv1 = fw;
w = u;
fw = fu;
}
else
{
if (fu < fv1 || FloatEqual(fu,fv1) || FloatEqual(v,x) || FloatEqual(v,w))
{
v = u;
fv1 = fu;
}
}
}
}
if (iter > itmax)
{
cout<<" ExtremumBrentNonDerivative1D exceed maximum iterations."<<endl;
}
xmin = x;
if(sign == -1)
{
fx *= sign;
}
return fx;
}
#endif // _EXTREMUM_INL
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -