// Basic.cpp: implementation of the Basic class.
//  隨機射線法,用于解決大型優化問題的最優解.
//  2006年1月由本人編制,如有任何問題,請聯系xiaofc4395@sina.com

//  這個程序目前只支持變量有界約束,不支持其它類型的約束.
//  函數本身是求極大值的,如果需要求極小值,需要做簡單處理(前面*負號 ^_^ ),例程中就是求極小值的.

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

#include "stdafx.h"
#include "RanMethod.h"
#include "Basic.h"
#include "mat.h"

#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

Basic::Basic()
{

}

Basic::~Basic()
{

}

mat Basic::RandMethod(mat xmin, mat xmax,mat xstart,double eps)
{
	//xmin,xmax,xstart---必須是列數相等的列向量
	int n=xstart.row;
	ASSERT( xmin.row==xmax.row && xmin.col==xmax.col);
	ASSERT( xstart.row==xmax.row && xstart.col==xmax.col);
	ASSERT( xstart.col==1);
//-----------------------------------------
	double A,M,OPT,Y;
	mat Z,E1,E0,D,X,XX,XOPT,XTEMP;
	int N,NP,Q;
	A=0.2;
	N=1;
	NP=1;
	OPT=minf(xstart);
	E0=norm(xmin,xmax,xstart);
step6:
	Q=1;
	M=1;
	D=GenerateDX(n);
step8:
	Z=E0+D*A*M;
	E1=norm11(Z);
	X=E1;
	XTEMP=invnorm(xmin,xmax,E1);
	Y=minf(XTEMP);
	N=N+1;
	if(Y>OPT)
		NP=N;
	if(Y<=OPT && Q==1)
		goto step10;
	goto step14;
step10:
	if(N-NP >= n*10)
		goto step19;
	Q=2;
	D=D*(-1.0);
	goto step8;
step14:
	Q=2;
	if(Y>OPT)
		goto step15;
	goto step17;
step15:
	E0=E1;
	XOPT=E0;
	OPT=Y;
	M=2*M;
	goto step8;
step17:
	if(A<=eps)
		goto step20;
	if(N-NP >= n*10)
		goto step19;
	goto step6;
step19:
	A=0.1*A;
	NP=N;
	goto step6;
step20:
	mat out;
	if(XOPT.p==NULL)
		out=xstart;
	else
		out=invnorm(xmin,xmax,XOPT);
	result=OPT;
	return out;
}

double Basic::Random()
{
	double out=rand()/32768.0;
    return out;
}

mat Basic::GenerateDX(int m)
{
	//m---向量的列數
	mat out;out.zero(m,1);
	for(int i=1;i<=m;i++)
	{
		double C=Random();
		double flag=ADDSUB(C);	
		double DX=flag*C*C;
		out.set(i,1,DX);
	}
	return out;
}

mat Basic::invnorm(mat xmin, mat xmax, mat xnow)
{
	int m=xnow.row;
	int n=xnow.col;
	ASSERT( xmin.row==xmax.row && xmin.row==xnow.row);
	ASSERT( xmin.col==xmax.col && xmin.col==xnow.col);
	mat deta=xmax-xmin;
	mat out;out.zero(m,n);
	for (int i=1;i<=m;i++)
	{
		for(int j=1;j<=n;j++)
		{
			out.set(i,j,xmin.r(i,j)+xnow.r(i,j)*deta.r(i,j));
		}
	}
	return out;
}

mat Basic::norm(mat xmin, mat xmax, mat xnow)
{
	int m=xnow.row;
	int n=xnow.col;
	ASSERT( xmin.row==xmax.row && xmin.row==xnow.row);
	ASSERT( xmin.col==xmax.col && xmin.col==xnow.col);
	mat deta=xmax-xmin;
	mat out;out.zero(m,n);
	for (int i=1;i<=m;i++)
	{
		for(int j=1;j<=n;j++)
		{
			double temp=(xmax.r(i,j)-xnow.r(i,j))/(xmax.r(i,j)-xmin.r(i,j));
			out.set(i,j,temp);
		}
	}
	return out;	
}

double Basic::ADDSUB(double c)
{
	double out;
	if( c >= 0.5)
		out=1.0;
	else
		out=-1.0;
	return out;
}

mat Basic::norm11(mat x)
{
	for (int i=1;i<=x.row;i++)
	{
		for(int j=1;j<=x.col;j++)
		{
			double Z=x.r(i,j);
			if (Z < 0.0)
				Z=-Z;
			if(Z > 1.0) 
				Z=2.0-Z;
			if (Z< 0.0) 
				Z=-Z;
			x.set(i,j,Z);
		}
	}
	return x;
}