//單純形法類cpp文件
//Edit by Wang Shimin
//2004.8.25～2004.9.10
//本文件為最新文件，修改于2006.3.12

// Simplex.cpp: implementation of the Simplex class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "Simplex.h"
#include "math.h"
#include "fstream"

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

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

Simplex::Simplex()
{
	/*	初始化參數(shù)值 */
	m_precisionnum = 0; /* 連續(xù)多少步不進(jìn)化則終止 */
	m_maxnumlap = 10;   /* 最大循環(huán)次數(shù) */
	m_numlap = 1;       /* 循環(huán)次數(shù) */
	m_precision = 0.000001;  /* 收斂精度 */

//	m_sidelength = 1.0;  /* 單純形邊長 */
	m_oldfitness = 0.0; /* 上一步最優(yōu)目標(biāo)值 */
//	m_parametereq = (sqrt((double)m_numvar+1.0)+(double)m_numvar-1.0)/((double)m_numvar*sqrt(2.0))*m_sidelength;   /* 參數(shù)Q */
//	m_parameterep = m_parametereq-m_sidelength/sqrt(2.0);   /* 參數(shù)P */
	m_parameterep = 0.1; /* 參數(shù)P */
	m_parametereq = 0.2;/* 參數(shù)Q */
	m_alfa = 1.0;   /* 反射系數(shù)，大于0 */
	m_beta = 0.5;   /* 收縮系數(shù)，大于0，小于1 */
	m_gamma = 1.0;  /* 延伸系數(shù),大于0 */

	/*將指針置零*/
	p_Bndload = NULL;
	p_BAKindflag = NULL;
	p_iparameter = NULL;
	p_jparameter = NULL;
	m_constriction = NULL;
	p_prolongation = NULL;
	p_center = NULL;
	p_echo = NULL;
	p_variable = NULL;
	p_bestobjvalue = NULL;
	p_fitness = NULL;
	p_observation = NULL;
	p_calculation = NULL;

	m_blErrorExit = false;
}

Simplex::~Simplex()
{
	/*釋放內(nèi)存空間*/
	for(int i=0; i<=m_numvar; i++)
	{
		delete[]p_variable[i];
	}
	delete[]p_variable;

	delete[]m_constriction;
	delete[]p_prolongation;
	delete[]p_center;
	delete[]p_echo;
	delete[]p_bestobjvalue;
	delete[]p_fitness;
	delete[]p_observation;
	delete[]p_calculation;
	delete[]p_BAKindflag;
	delete[]p_iparameter;
	delete[]p_jparameter;
 	delete[]p_Bndload;
}

/* 計(jì)算 */
void Simplex::Compute()
{
	/* 初始化函數(shù) */
	Initialization();
	if(m_blErrorExit)return;
	
	while(Criterion())
	{
		m_oldfitness = m_lowvalue;
		/* 反射 */
		Echo();
		if(m_blErrorExit)return;

		if(fabs((m_oldfitness-m_lowvalue)/m_oldfitness)<m_precision)
			m_precisionnum++;
	
		/* 輸出反分析優(yōu)化結(jié)果 */
		WriteResult();
		if(m_blErrorExit)return;

		m_numlap++;
	}
}

/* 初始化函數(shù) */
void Simplex::Initialization()
{
	double ** pobjvalue;     /* 存儲目標(biāo)值，計(jì)算值與觀測值之間的差值  */
	pobjvalue = NULL;

	CString csBaaFileName = mGcsFileName+_T(".BAA");
	ifstream cBaaFile;
	cBaaFile.open((LPCSTR)csBaaFileName); 
	if (cBaaFile.fail())
	{
		AfxMessageBox("Cannot open input BAA file!!!",MB_OK);
		m_blErrorExit = true;
		return;
	}
	
	int itemp;
	cBaaFile >> itemp;
	cBaaFile >> m_numvar;   /* 待反演參數(shù)的個(gè)數(shù) */
	cBaaFile.close();


	CString csMEAFileName = mGcsFileName+_T(".MEA");
	ifstream cMEAFile;
	cMEAFile.open((LPCSTR)csMEAFileName); 
	if (cMEAFile.fail())
	{
		AfxMessageBox("Cannot open input MEA file!!!",MB_OK);
		m_blErrorExit = true;
		return;
	}
	
	cMEAFile >> m_numpoint;  /* 從MEA文件中讀取觀測點(diǎn)的個(gè)數(shù) */
	cMEAFile.close();

	int i,j;

	/* 為存儲識別待反演參數(shù)性質(zhì)的數(shù)組分配內(nèi)存空間 */

	if(p_BAKindflag)
		delete[]p_BAKindflag;

	if(p_iparameter)
		delete[]p_iparameter;

	if(p_jparameter)
		delete[]p_jparameter;

	p_BAKindflag = new int[m_numvar];
	p_iparameter = new int[m_numvar];
	p_jparameter = new int[m_numvar];

	/* 存儲量測值與計(jì)算值之差 */
	if(pobjvalue)
		delete[]pobjvalue;
	pobjvalue = new double *[m_numvar+1]; 
	for(i=0; i<=m_numvar; i++)
		pobjvalue[i]=new double[m_numpoint];
	/* 存儲量測值與計(jì)算值之差 */

	/* 為單純形分配動(dòng)態(tài)內(nèi)存空間 */
	if(p_variable)
		delete[]p_variable;
	p_variable = new double *[m_numvar+1];

	for(i=0; i<=m_numvar; i++)
		p_variable[i]=new double[m_numvar];

	/* *********** 建立初始單純形 ************* */
	
	/* 從BAA文件中讀入初始點(diǎn) */
	ReadBAA();
	if(m_blErrorExit)return;

	/* 形成單純形 */
	for(i=1; i<=m_numvar; i++)
	{
		for(j=0; j<m_numvar; j++)
		{
			if(j==(i-1))
				p_variable[i][j] = p_variable[0][j]*(1.0+m_parameterep);
			else
				p_variable[i][j] = p_variable[0][j]*(1.0+m_parametereq);
		}
	}

	if(m_constriction)
		delete[]m_constriction;
	m_constriction = new double[m_numvar];  /* 為收縮點(diǎn)變量分配動(dòng)態(tài)內(nèi)存空間 */

	if(p_prolongation)
		delete[]p_prolongation;
	p_prolongation = new double[m_numvar];  /* 為延伸點(diǎn)變量分配動(dòng)態(tài)內(nèi)存空間 */

	if(p_center)
		delete[]p_center;
	p_center = new double[m_numvar];  /* 為中心點(diǎn)變量分配動(dòng)態(tài)內(nèi)存空間 */

	if(p_echo)
		delete[]p_echo;
	p_echo = new double[m_numvar];   /* 為反射點(diǎn)變量分配動(dòng)態(tài)內(nèi)存空間 */ 

	if(p_bestobjvalue)
		delete[]p_bestobjvalue;
	p_bestobjvalue = new double[m_numpoint]; /* 存儲量測值與計(jì)算值之差 */
	
	if(p_fitness)
		delete[]p_fitness;
	p_fitness = new double[m_numvar+1]; /* 存儲單純形中的目標(biāo)值 */

	if(p_observation)
		delete[]p_observation;
	p_observation = new double[m_numpoint]; /* 存儲觀測值 */
	
	if(p_calculation)
		delete[]p_calculation;
	p_calculation = new double[m_numpoint]; /* 存儲計(jì)算值 */

	/* 讀入測點(diǎn)的觀測值 */
	ReadMEA();
	if(m_blErrorExit)return;

	for(i=0; i<=m_numvar; i++)
	{
		p_fitness[i] = 0.0;

		/* 將參數(shù)存儲到文件 */
		WriteToFem(p_variable[i]);
		if(m_blErrorExit)return;
		/* 有限元計(jì)算 */
		ObjFunction();
		if(m_blErrorExit) return;
		/* 讀入測點(diǎn)的計(jì)算值 */
		ReadFromFem(p_calculation);
		if(m_blErrorExit)return;
		
		/* 計(jì)算測點(diǎn)的觀測值與計(jì)算值之差的平方和，即一組參數(shù)的目標(biāo)值 */
		for(j=0; j<m_numpoint; j++)
		{
			pobjvalue[i][j] = p_observation[j]-p_calculation[j];
			p_fitness[i] += pow(pobjvalue[i][j],2.0);
		}
	}
	
	/* 求最大值和最小值 */
	m_highvalue = p_fitness[0];
	m_lowvalue = p_fitness[0];
	m_highnumber = 0;
	m_lownumber = 0;
	
	for(i=1; i<=m_numvar; i++)
	{
		if(p_fitness[i]>m_highvalue)
		{
			m_highvalue = p_fitness[i];
			m_highnumber = i;
		}
		
		if(p_fitness[i]<m_lowvalue)
		{
			m_lowvalue = p_fitness[i];
			m_lownumber = i;
		}
	}
	
	m_oldfitness = m_lowvalue;
	
	for(i=0;i<m_numpoint;i++)
	{
		p_bestobjvalue[i] = pobjvalue[m_lownumber][i];
	}
	
	delete [] pobjvalue;
}

/* 判別函數(shù) */
bool Simplex::Criterion()
{
	bool flag;
	flag = TRUE;

	if((m_numlap>m_maxnumlap)||(m_precisionnum>0))
		flag = FALSE;

	return flag;
}

/* 反射 */
void Simplex::Echo()
{
	int i,j;
	double sum;
	double * peobjvalue;
	peobjvalue = new double[m_numpoint];
	bool flag;
	flag=TRUE;

	/* 求最大值和最小值 */
	m_highvalue = p_fitness[0];
	m_lowvalue = p_fitness[0];
	m_highnumber = 0;
	m_lownumber = 0;

	for(i=1; i<=m_numvar; i++)
	{
		if(p_fitness[i]>m_highvalue)
		{
			m_highvalue = p_fitness[i];
			m_highnumber = i;
		}
		
		if(p_fitness[i]<m_lowvalue)
		{
			m_lowvalue = p_fitness[i];
			m_lownumber = i;
		}
	}
	
	for(i=0; i<m_numvar; i++)
	{
		sum = 0.0;
		
		for(j=0; j<=m_numvar; j++)
			sum += p_variable[j][i]-p_variable[m_highnumber][i];
	
		p_center[i] = fabs(sum/double(m_numvar));

		if(p_variable[m_highnumber][i]>p_variable[m_lownumber][i])
			p_echo[i] = p_variable[m_lownumber][i]-m_alfa*p_center[i];
		else
			p_echo[i] = p_variable[m_lownumber][i]+m_alfa*p_center[i];
	}

	/****************** 求中心點(diǎn)的目標(biāo)值 *******************/

	/* 將參數(shù)存儲到文件 */
	WriteToFem(p_center);
	if(m_blErrorExit)return;
	/* 有限元計(jì)算 */
	ObjFunction();
	if(m_blErrorExit) return;
	/* 讀入測點(diǎn)的計(jì)算值 */
	ReadFromFem(p_calculation);
	if(m_blErrorExit)return;
	
	m_centervalue = 0.0;

	/* 計(jì)算測點(diǎn)的觀測值與計(jì)算值之差 */
	for(i=0; i<m_numpoint; i++)
		m_centervalue += pow((p_observation[i]-p_calculation[i]),2.0);

	/****************** 求中心點(diǎn)的目標(biāo)值 *******************/

	/****************** 求反射點(diǎn)的目標(biāo)值 *******************/

	/* 將參數(shù)存儲到文件 */
	WriteToFem(p_echo);
	if(m_blErrorExit)return;
	/* 有限元計(jì)算 */
	ObjFunction();
	if(m_blErrorExit) return;
	/* 讀入測點(diǎn)的計(jì)算值 */
	ReadFromFem(p_calculation);
	if(m_blErrorExit)return;
	
	m_echovalue = 0.0;

	/* 計(jì)算測點(diǎn)的觀測值與計(jì)算值之差 */
	for(i=0; i<m_numpoint; i++)
	{
		peobjvalue[i] = p_observation[i]-p_calculation[i];
		m_echovalue += pow(peobjvalue[i],2.0);
	}

	/****************** 求反射點(diǎn)的目標(biāo)值 *******************/
	
	if(m_echovalue<=m_lowvalue)
	{
		/* 去掉原來的最大值點(diǎn)，取而代之以原來的最小點(diǎn)，
		并用反射點(diǎn)替換單純形中目標(biāo)值最小的點(diǎn) */
		for(i=0; i<m_numvar; i++)
		{
			p_variable[m_highnumber][i] = p_variable[m_lownumber][i];
			p_variable[m_lownumber][i] = p_echo[i];
		}

		p_fitness[m_highnumber] = p_fitness[m_lownumber];
		m_highvalue = m_lowvalue;

		p_fitness[m_lownumber] = m_echovalue;
		m_lowvalue = m_echovalue;
		
		for(i=0;i<m_numpoint;i++)
		{
			p_bestobjvalue[i] = peobjvalue[i];
		}
		
		/* 延伸 */
		Prolongation();
		if(m_blErrorExit)return;
	}
	else
	{
		for(i=0; i<m_numvar; i++)
		{
			if((i!=m_highnumber)&&(p_fitness[i]>m_echovalue))
				flag=FALSE;
		}
		
		if(flag=TRUE)
		{
			if(m_echovalue>=m_highvalue)
			{
				/* 收縮 */
				Constriction();
				if(m_blErrorExit)return;
			}
			else
			{
				/* 用反射點(diǎn)替換單純形中目標(biāo)值最大的點(diǎn) */
				for(i=0; i<m_numvar; i++)
					p_variable[m_highnumber][i] = p_echo[i];
				
				p_fitness[m_highnumber] = m_echovalue;
				m_highvalue = m_echovalue;
				
				/* 收縮 */
				Constriction();
				if(m_blErrorExit)return;
			}
		}
		else
		{
			/* 用反射點(diǎn)替換單純形中目標(biāo)值最大的點(diǎn) */
			for(i=0; i<m_numvar; i++)
				p_variable[m_highnumber][i] = p_echo[i];
			
			p_fitness[m_highnumber] = m_echovalue;
			m_highvalue = m_echovalue;
		}
	}

	delete [] peobjvalue;
}

/* 延伸 */
void Simplex::Prolongation()
{
	int i,flag;

	flag = 0;

	double * ppobjvalue;
	ppobjvalue = new double[m_numpoint];
	
    for(i=0; i<m_numvar; i++)
	{
		if(p_variable[m_highnumber][i]>p_variable[m_lownumber][i])
			p_prolongation[i] = p_variable[m_lownumber][i]-m_gamma*p_center[i];
		else
			p_prolongation[i] = p_variable[m_lownumber][i]+m_gamma*p_center[i];
	}
	
	/* 將參數(shù)存儲到文件 */
	WriteToFem(p_prolongation);
	if(m_blErrorExit)return;
	/* 有限元計(jì)算 */
	ObjFunction();
	if(m_blErrorExit) return;
	/* 讀入測點(diǎn)的計(jì)算值 */
	ReadFromFem(p_calculation);
	if(m_blErrorExit)return;
	
	m_prolongationvalue = 0.0;
	
	/* 計(jì)算測點(diǎn)的觀測值與計(jì)算值之差 */
	for(i=0; i<m_numpoint; i++)
	{
		ppobjvalue[i] = p_observation[i]-p_calculation[i];