// Extremum.inl		計算極值函數(方法)定義頭文件
// Ver 1.0.0.0
// 版權所有(C) 何渝, 2002
// 最后修改: 2002.5.31.

#ifndef _EXTREMUM_INL		//避免多次編譯
#define _EXTREMUM_INL

//一維極值連分式法
template <class _Ty>
void ExtremumFraction1D(valarray<_Ty>& x, _Ty eps, 
										int k, valarray<int>& js)
{	
    _Ty xx,h1,dx,y[10],b[10],h2(0);
	valarray<_Ty> z(2);
    int jt = 1;
	js[0] = 0;
    while(jt == 1)
    {
		int j=0;
        while(j<=7)
        {
			if(j <=2) 
				xx=x[0]+0.01*j;
            else xx=h2;
            FunctionValueEF1D(xx,z);
            if(Abs(z[1]) < eps)
            {
				jt=0; 
				j=10;
			}
            else
            {
				h1=z[1]; 
				h2=xx;
                if(j==0)
				{
					y[0]=h1;
					b[0]=h2;
				}
                else
                {
					y[j]=h1;
					int m=0;
					int i=0;
                    while((m==0)&&(i<=j-1))
					{
						if(FloatEqual((h2-b[i]),0)) 
							m=1;
						else
							h2=(h1-y[i])/(h2-b[i]);
                        i=i+1;
					}
                    b[j]=h2;
                    if(m!=0)
						b[j]=1.0e+35;
                    h2=0.0;
                    for(i=j-1; i>=0; i--)
						h2=-y[i]/(b[i+1]+h2);
                    h2=h2+b[0];
                }
                j=j+1;
            }
        }
        x[0]=h2;
        js[0]=js[0]+1;
        if(js[0]==k)  jt=0;
    }
    xx = x[0];
    FunctionValueEF1D(xx,z);
	x[1] = z[0];
    if(Abs(x[0]) < 1.0 || FloatEqual(Abs(x[0]), 1)) dx=1.0e-05;
    else	dx = Abs(x[0] * 1.0e-05);
    xx = x[0] - dx;
    FunctionValueEF1D(xx,z);
	h1 = z[0];
    xx = x[0] + dx;
    FunctionValueEF1D(xx,z);
	h2 = z[0];
    js[1] = -1;
    if((h1 + h2 - 2.0 * x[1]) < 0.0 || FloatEqual((h1 + h2 - 2.0 * x[1]),0)) js[1]=1;
}

//n維極值連分式法
template <class _Ty>
void ExtremumFractionND(valarray<_Ty>& x, _Ty eps, 
										int k, valarray<int>& js)
{
    int i,j,m,l,jt(1),il;
    _Ty y[10],b[10],p,z,t,h1,h2(0),f,dx;

	int n = x.size()-1;			//自變量個數

    js[0] = 0;
	//jt = 1;
	//h2 = 0;
    while(jt == 1)
    {
		t=0.0; 
		js[0]=js[0]+1;
        for(i=1; i<=n; i++)
        {
			f=FunctionValueEFND(x,i);
            t=t+Abs(f);
        }
        if(t<eps) jt=0;
        else
        {
			for(i = 0; i < n; i++)
              { 
				il=5;
                while(il!=0)
                {
					j=0;
					t=x[i];
					il=il-1;
                    while(j<=7)
                    {
						if(j<=2) z=t+j*0.01;
                        else z=h2;
                        x[i]=z;
                        f=FunctionValueEFND(x,i+1);
                        if(FloatEqual(f,0))
                        {
							j=10;
							il=0;
						}
                        else
                        {
							h1=f;
							h2=z;
                            if(j==0)
                            {
								y[0]=h1;
								b[0]=h2;
							}
                            else
                            {
								y[j]=h1;
								m=0;
								l=0;
                                while((m==0)&&(l<=j-1))
                                {
									p=h2-b[l];
                                    if(FloatEqual(p,0)) m=1;
                                    else h2=(h1-y[l])/p;
                                    l=l+1;
                                }
                                b[j]=h2;
                                if(m!=0) b[j]=1.0e+35;
                                h2=0.0;
                                for(l=j-1; l>=0; l--)
                                h2=-y[l]/(b[l+1]+h2);
                                h2=h2+b[0];
                            }
                            j=j+1;
                        }
                    }
                    x[i]=h2; 
                }
                x[i]=z;
            }
            if(js[0]==k) jt=0;
        }
    }
    js[1]=1;
    f=FunctionValueEFND(x,0); 
	x[n]=f;
    dx=0.00001;
	t=x[0];
    x[0]=t+dx;
	h1=FunctionValueEFND(x,0);
    x[0]=t-dx;
	h2=FunctionValueEFND(x,0);
    x[0]=t;
    t=h1+h2-2.0*f;
    if(t>0.0) js[1]=-1;
    j=1; jt=1;
    while(jt==1)
    {
		j=j+1;
		dx=0.00001;
		jt=0;
        t=x[j-1];
        x[j-1]=t+dx;
		h2=FunctionValueEFND(x,0);
        x[j-1]=t-dx;
		h1=FunctionValueEFND(x,0);
        x[j-1]=t;
		t=h1+h2-2.0*f;
        if((t*js[1]<0.0)&&(j<n)) jt=1;
    }
    if(t*js[1]>0.0) js[1]=0;
}

//線性規劃
template <class _Ty>
int ExtremumLinePrograming(matrix<_Ty>& a, valarray<_Ty>& b, 
								valarray<_Ty>& c, valarray<_Ty>& x)
{
    _Ty y;
	
	int m = b.size();		//不等式約束條件的個數
	int n = c.size() - m;	//變量個數

	valarray<int> js(m);
	matrix<_Ty> p(m, m);
	matrix<_Ty> d(m, (m+n));
	for(int i = 0; i < m; i++) js[i] = n + i;
    int mn = m + n;
	_Ty s(0);
    while(1==1)
    {
		for(i = 0; i < m; i++)
          for(int j = 0; j < m; j++)
			  p(i,j) = a(i, js[j]);
        int l=MatrixInversionGS(p);
        if(l==0)
        {
			x[n]=s;
            return(l);
        }

		d = p * a;			//矩陣相乘
        
        for(i=0; i<mn; i++) x[i]=0.0;
        for(i=0; i<m; i++)
        {
			s=0.0;
            for(int j=0; j<m; j++)
              s=s + p(i,j) * b[j];
            x[js[i]]=s;
        }
        int k = -1;
		_Ty dd=1.0e-35;
        for(int j = 0; j < mn; j++)
        {
			_Ty z(0);
            for(i=0; i<m; i++)
              z = z + c[js[i]] * d(i, j);
            z = z - c[j];
            if(z > dd)
			{
				dd = z;
				k = j;
			}
        }
        if(k == -1)
          {
			s=0.0;
            for(j=0; j<n; j++)
				s = s + c[j] * x[j];
			x[n] = s;
            return(1);
        }
        j = -1;
        dd = 1.0e+20;
        for(i=0; i<m; i++)
			if(d(i, k) >= 1.0e-20 )
            { 
			  y=x[js[i]] / d(i, k);
              if(y < dd)
			  {
				  dd = y; 
				  j = i;
			  }
			}
        if(j==-1) 
        {
			x[n]=s;
            return(0);
        }
        js[j]=k;
    }
}

//n維極值單形調優法
template <class _Ty>
int ExtremumSimplexND(_Ty d, _Ty u, _Ty v, valarray<_Ty>& x,  
				_Ty eps, int k, matrix<_Ty>& xx, valarray<_Ty>& f)
{
	int r,g,l;
	_Ty fe,ft,ff;
	
	int n = x.size() - 1;	//變量個數

	valarray<_Ty> xt(n);
	valarray<_Ty> xf(n);
	valarray<_Ty> xe(n);

    int kk=0;
	_Ty nn=1.0*n;
    _Ty fr=sqrt(nn+1.0);
    _Ty fl=d*(fr-1.0)/(1.414*nn);
    _Ty fg=d*(fr+nn-1.0)/(1.414*nn);
    for(int i=0; i<n; i++)
    for(int j=0; j<=n; j++)
       xx(i,j) = 0;
    for(i=1; i<=n; i++)
    for(int j=0; j<n; j++)
       xx(j,i) = fl;
    for(i=1; i<=n; i++)
       xx((i-1),i) = fg;
    for(i=0; i<=n; i++)
    { 
		for(int j=0; j<n; j++)	xt[j]=xx(j,i);
        f[i]=FunctionValueESND(xt,n);
    }
    ft=1.0+eps;
    while((kk<k)&&(ft>eps))
    {
		kk=kk+1;
        fr=f[0]; fl=f[0]; r=0; l=0;
        for(i=1; i<=n; i++)
        {
			if(f[i]>fr)
			{
				r=i;
				fr=f[i];
			}
            if(f[i]<fl)
			{
				l=i;
				fl=f[i];
			}
        }
        g=0;
		fg=f[0];
        int j=0;
        if(r==0)
		{
			g=1;
			fg=f[1];
			j=1;
		}
        for(i=j+1; i<=n; i++)
			if((i!=r)&&(f[i]>fg))
            {
				g=i;
				fg=f[i];
			}
        for(j=0; j<n; j++)
        {
			xf[j]=0.0;
            for(i=0; i<=n; i++)
				if(i!=r)	xf[j]=xf[j]+xx(j,i)/nn;
            xt[j]=2.0*xf[j]-xx(j,r);
        }
        ft=FunctionValueESND(xt,n);
        if(ft<f[l])
        {
			for(j=0; j<n; j++)	xf[j]=(1.0+u)*xt[j]-u*xf[j];
            ff=FunctionValueESND(xf,n);
            if(ff<f[l])
            {
				for(j=0; j<n; j++)	xx(j,r)=xf[j];
                f[r]=ff;
            }
            else
			{
				for(j=0; j<n; j++)	xx(j,r)=xt[j];
                f[r]=ft;
            }
        }
        else
        {
			if(ft<=f[g])
            {
				for(j=0; j<=n-1; j++)	xx(j,r)=xt[j];
                f[r]=ft;
            }
            else 
            {
				if(ft<=f[r])
                {
					for(j=0; j<=n-1; j++)	xx(j,r)=xt[j];
                    f[r]=ft;
                }
                for(j=0; j<n; j++)	xf[j]=v*xx(j,r)+(1.0-v)*xf[j];
                ff=FunctionValueESND(xf,n);
                if(ff>f[r])
                  for(i=0; i<=n; i++)
                  {
					  for(j=0; j<=n-1; j++)
					  {
						  xx(j,i)=(xx(j,i)+xx(j,l))/2.0;
                          x[j]=xx(j,i);
						  xe[j]=x[j];
                      }
                      fe=FunctionValueESND(xe,n);
					  f[i]=fe;
                  }
                else
                {
					for(j=0; j<=n-1; j++)	xx(j,r)=xf[j];
                    f[r]=ff;
                }
            }
        }
        ff=0.0;
		ft=0.0;
        for(i=0; i<=n; i++)
        {
			ff=ff+f[i]/(1.0+nn);
            ft=ft+f[i]*f[i];
        }
        ft=(ft-(1.0+n)*ff*ff)/nn;
    }
    for(int j=0; j<n; j++)
    {
		x[j] = 0.0;
        for(i=0; i<=n; i++)	x[j] = x[j] + xx(j,i) / (1.0+nn);
        xe[j] = x[j];
    }
    fe = FunctionValueESND(xe,n);
	x[n] = fe;
    return(kk);
}

//n維極值復形調優法
template <class _Ty>
int ExtremumComplexND(int m, valarray<_Ty>& a, valarray<_Ty>& b,
		_Ty alpha, _Ty eps, valarray<_Ty>& x, matrix<_Ty>& xx, int k)
{
    int r, g, kt, jt, kk;
    _Ty fj, fr, fg, z;
	double rr = 0;				//隨機數種子
	int n = a.size();			//變量個數

	valarray<_Ty> c(m);
	valarray<_Ty> d(m);
	valarray<_Ty> w(m);
	valarray<_Ty> xt(m);
	valarray<_Ty> xf(m);

    for(int i=0; i<n; i++)	xx(i,0) = x[i];
    xx(n,0) = FunctionValueTarget(x,n);
    for(int j=1; j<2*n; j++)
    {
		for(i=0; i<n; i++)
		{
			xx(i,j) = a[i] + (b[i] - a[i]) * rand_01_One(rr);
            x[i] = xx(i,j);
        }
        int it = 1;
        while(it==1)
        {
			it = 0;
			r = 0;
			g = 0;
            while((r<n)&&(g==0))
            {
				if((a[r]<=x[r])&&(b[r]>=x[r])) r = r + 1;
                else g = 1;
            }
            if(g==0)
            {
				ConditionValue(n,m,x,c,d,w);
                r = 0;
                while((r<m)&&(g==0))
                {
					if((c[r]<=w[r])&&(d[r]>=w[r])) r=r+1;
                    else g=1;
                }
            }
            if(g!=0)
            {
				for(r=0; r<n; r++)
                {
					z = 0.0;
                    for(g=0; g<j; g++)	z = z + xx(r,g) / (1.0*j);
                    xx(r,j) = (xx(r,j) + z) / 2.0;
                    x[r]=xx(r,j);
                }
                it = 1;
            }
            else xx(n,j) = FunctionValueTarget(x,n);
        }
    }
    kk = 1;
	int it = 1;
    while(it==1)
    {
		it = 0;
        fr = xx(n,0);
		r = 0;