hardy_1codedlg.cpp

利用多面函数拟合法来拟合似大地水准面

// Hardy_1CodeDlg.cpp : implementation file
//

#include "stdafx.h"
#include "GeoinModel.h"
#include "Hardy_1CodeDlg.h"
#include "math.h"

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

/////////////////////////////////////////////////////////////////////////////
// CHardy_1CodeDlg dialog


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

extern double CenterL(double dl);
extern  double CenterB(double db);
extern  double ReturnCenterL(double dl);
extern  double ReturnCenterB(double db);
extern  double Arc2DMS(double dd);
extern  double int_radian(double intdeg);
extern  double deg_int(double gms);
extern double myu(double a, double b);

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

CHardy_1CodeDlg::CHardy_1CodeDlg(CWnd* pParent /*=NULL*/)
	: CDialog(CHardy_1CodeDlg::IDD, pParent)
{
	//{{AFX_DATA_INIT(CHardy_1CodeDlg)
	m_dCrosser = 0.3;
	m_strKnownFilePath = _T("");
	m_dMutation = 0.2;
	m_strParamFilePath = _T("");
	m_dq = 0.0;
	m_db = -1;
	m_iPopSize = 200;
	m_iMaxCalcNum = 5000;
	m_bOverwrite = true;
	m_strOutputFilePath = _T("");
	//}}AFX_DATA_INIT
}


void CHardy_1CodeDlg::DoDataExchange(CDataExchange* pDX)
{
	CDialog::DoDataExchange(pDX);
	//{{AFX_DATA_MAP(CHardy_1CodeDlg)
	DDX_Text(pDX, IDC_EDTCrosser, m_dCrosser);
	DDX_Text(pDX, IDC_EDTKnownFilePath, m_strKnownFilePath);
	DDX_Text(pDX, IDC_EDTMutation, m_dMutation);
	DDX_Text(pDX, IDC_EDTParamFilePath, m_strParamFilePath);
	DDX_Text(pDX, IDC_EDTPopSize, m_iPopSize);
	DDX_Text(pDX, IDC_EDTMaxCalcNum, m_iMaxCalcNum);
	DDX_Check(pDX, IDC_CHECKOverwrite, m_bOverwrite);
	DDX_Text(pDX, IDC_EDTOutputFilePath, m_strOutputFilePath);
	//}}AFX_DATA_MAP
}


BEGIN_MESSAGE_MAP(CHardy_1CodeDlg, CDialog)
	//{{AFX_MSG_MAP(CHardy_1CodeDlg)
	ON_BN_CLICKED(ID_Calc, OnCalc)
	ON_BN_CLICKED(IDC_BTNOpenKnownFile, OnBTNOpenKnownFile)
	ON_BN_CLICKED(IDC_BTNOpenParamFile, OnBTNOpenParamFile)
	ON_BN_CLICKED(IDC_BTNOpenOutputFile, OnBTNOpenOutputFile)
	ON_EN_CHANGE(IDC_EDTPopSize, OnChangeEDTPopSize)
	ON_EN_CHANGE(IDC_EDTMaxCalcNum, OnChangeEDTMaxCalcNum)
	ON_EN_CHANGE(IDC_EDTCrosser, OnChangeEDTCrosser)
	ON_EN_CHANGE(IDC_EDTMutation, OnChangeEDTMutation)
	//}}AFX_MSG_MAP
END_MESSAGE_MAP()

/////////////////////////////////////////////////////////////////////////////
// CHardy_1CodeDlg message handlers

//计算
void CHardy_1CodeDlg::OnCalc() 
{
	double dTemp;
	CString strPtname,strB,strL,strH,strh,strType,strEnd;
	int i,j;

	if(m_strKnownFilePath=="")
	{
		MessageBox("请先指定已知数据文件。");
		return;
	}
	if(m_strOutputFilePath=="")
	{
		MessageBox("请先指定结果输出文件。");
		return;
	}

	//◎◎读入已知坐标文件
	char buffer[256];
	FILE *File=fopen(m_strKnownFilePath,"r");

	//读取已知坐标文件,点名,B,L,H,h,i
	m_iPtNum=0;m_iCodeNum=1;

	//移到文件头
//	fseek(File,0,0);

	i=0;j=0;
	while(fgets(buffer,256,File)!=NULL)	
	{
		strPtname=strtok(buffer,",， ");			
		strB=strtok('\0',",， ");
		strL=strtok(  '\0',",， ");
		strH=strtok('\0',",， ");
		strh=strtok('\0',",， ");
		strType=strtok('\0',",， ");
		strEnd=strtok('\0',",， \r\n");
			
		if(atoi(strType)<2)  //0-拟和点，1-节点，2-不用
		{
			i++;	

			dTemp =atof(strB);       //B格式：DD.MMSSS
			m_dB[i]=int_radian(deg_int(dTemp));
			m_dB[i]=CenterB(m_dB[i]);    //中心化
			
			dTemp=atof(strL);          //L格式：DD.MMSSS
			m_dL[i]=int_radian(deg_int(dTemp));    //弧度
			m_dL[i]=CenterL(m_dL[i]);    //中心化

			m_dE[i]=atof(strH)-atof(strh);    //高程异常
		
			if (atoi(strType)==1)  //节点数据
			{
				j++;
				if(j==m_iCodeIndex)
				{
					m_dCodeB[1]=m_dB[i];
					m_dCodeL[1]=m_dL[i];
					m_dCodeE[1]=m_dE[i];
				}
			}
		}
	}
	fclose(File);
	m_iPtNum=i;		
	m_iParamNum=m_iCodeNum*3;

	if(m_iPtNum<=0)
	{
		MessageBox("已知数据文件格式有误，请查看。");
		return;
	}
	if(m_iCodeNum<=0)
	{
		MessageBox("已知数据文件中没有节点（点类型＝1），无法进行多面函数拟合，请查看。");
		return;
	}



	//◎◎遗传算法

	//产生初始种群
	GA_Initialization();
	//求目标函数值、亲和力
	GA_Evaluation(0);

	double a;
	q[0]=0.05; a=0.05;
	for(i=1; i<=m_iPopSize; i++) {a=a*0.95; q[i]=q[i-1]+a;}

	//◎◎◎遗传算法
	int iGen=0;
	//while((m_iCurrentGen<=m_iMaxCalcNum) && (m_dAbsV>0.01))
	while(iGen<=m_iMaxCalcNum)
	{
		  iGen++;
		  GA_Selection();
		  GA_Crossover();
		  GA_Mutation();
		  GA_Evaluation(iGen);
		
	}



	CString str;
	if(m_bOverwrite==0)
	{
		str="at";
	}
	else
	{
		str="wt";
	}

	if((File=fopen(m_strOutputFilePath,str))!=NULL)
	{
		//核函数类型，权函数类型
	//	fprintf(fp,"%s,%s\n",m_strCoreFunction,m_strP);
	//	fprintf(File,"-9999.9999,%d\n",m_iCodeNum);
		//转换参数a，平滑参数q，曲面参数b（0.5/-0.5），大地经度，大地纬度，高程异常
		for(j=1;j<=m_iCodeNum;j++) fprintf(File,"%.15f,%.3f,%.1f,%.9f,%.9f,%.4f\n",
			                   m_dAntibody[0][j],m_dAntibody[0][m_iCodeNum+j],m_dAntibody[0][2*m_iCodeNum+j],
							   m_dCodeB[j],m_dCodeL[j],m_dCodeE[j]);
		//fprintf(fp,"\n");
		fclose(File);
		str.Format("计算完成，共迭代了%d次，请查看结果文件Hardy_GA", iGen);   
		MessageBox(str);
		//打开文本文件
		CString Command;
		Command = "notepad ";
		Command += m_strOutputFilePath;
		WinExec(Command,SW_SHOWMAXIMIZED);

	}
	else
	{
		MessageBox("写入结果文件出错，请检查。。。");
	}
}

//限制条件
int CHardy_1CodeDlg::GA_Constraint_check(double x[])
{
	int i;
	//未知参数顺序：转移参数；平滑参数；正负双曲
	
	//正负双曲
	for(i=(2*m_iCodeNum+1);i++;i<=m_iCodeNum*3)
	{
		if(x[i]>1.0 || x[i]<0)	return 0;
	}

	return 1;
}


//初始化
void CHardy_1CodeDlg::GA_Initialization(void)
{
	int iNum,i,j,k;
	//打开参数文件	
	CString stra,strq,strb,strB,strL,strH,strE;

	//◎◎读入参数坐标文件

	iNum=0;  //文件中的参数组数
	
	if(m_strParamFilePath!="")
	{

		char buffer[256];
		FILE *File=fopen(m_strParamFilePath,"r");	
		while(fgets(buffer,256,File)!=NULL)
		{
			stra=strtok(buffer,",， ");
			strq=strtok('\0',",， ");
			strb=strtok(  '\0',",， ");	

			if(atof(stra)==-9999.9999)   //读入的是第一行：-9999.9999，核函数，权函数
			{
				iNum++;
				j=0;
			}
			else   //读入的是普通的参数行
			{
				j++;			
                //转换参数a，平滑参数q，曲面参数b（0.5/-0.5），大地经度，大地纬度，高程异常
				if(j==m_iCodeIndex)
				{	
					//要的参数是a,q,b,其他的不是未知参数，是作为大地高转换时用的
				
					//前m_iCodeNum个是转换参数a
					m_dAntibody[iNum][1]=atof(stra);
					//m_iCodeNum＋1——m_iCodeNum×2：平滑参数q
					m_dAntibody[iNum][2]=atof(strq);
					//m_iCodeNum×2＋1——m_iCodeNum×3：曲面参数b
					m_dAntibody[iNum][3]=atof(strb);
				}
				
			}
		}
		fclose(File);
	}

	//◎◎◎初始化

	//◎细胞克隆生产部分初始值
	for(i=1;i<=iNum;i++)   // 每个在【－1，1】间克隆10个
		for(k=1;k<=10;k++)
			for(j=1;j<=m_iParamNum;j++)
				 if((i-1)*10+iNum+k<=m_iPopSize) 
					 m_dAntibody[(i-1)*10+iNum+k][j]=m_dAntibody[i][j]+myu(-1.0,1.0);

	//◎随机生产部分初始值
	for(i=iNum*11+1;i<=m_iPopSize;i++)
		if(i<=m_iPopSize)
			for(j=1;j<=m_iCodeNum;j++) 
			{
				m_dAntibody[i][j]=myu(-20.0,20.0);           //a
				m_dAntibody[i][m_iCodeNum+j]=myu(-1.0,1.0);  //q
				m_dAntibody[i][2*m_iCodeNum+j]=myu(0.0,1.0); //b
			}

	//第0个赋初始值
	for(j=1;j<=m_iParamNum;j++)		
	     m_dAntibody[0][j]=m_dAntibody[1][j]+myu(-1.0,1.0);
	
}

//遗传算法求目标函数值
void CHardy_1CodeDlg::GA_Objective(void)
{	
	int i,j,k;
	double dE,dS;
	m_dAbsV=0;
	for(i=0;i<=m_iPopSize;i++)  //
	{
		m_dObjective[i]=0;	
		for(j=1;j<=m_iPtNum;j++)  //m_iPtNum — 参与拟和的已知点个数
		{
			dE=0;
			for(k=1;k<=m_iCodeNum;k++)
			{
				dS=sqrt(pow(m_dB[j]-m_dCodeB[k],2)+pow(m_dL[j]-m_dCodeL[k],2)+m_dAntibody[i][m_iCodeNum+k]);
				if(dS>999999.0 ||dS<-0.000000001) dS=0;

				if(m_dAntibody[i][2*m_iCodeNum+k]<0.5 && dS>0.0) dS=1.0/dS;  //倒双曲				
				dE+=m_dAntibody[i][k]*dS;
			}
			dE=dE-m_dE[j];
			//m_dObjective[i]=m_dObjective[i]+exp(fabs(dV));
			m_dObjective[i]=m_dObjective[i]+dE*dE;
			//m_dObjective[i]=m_dObjective[i]+exp(-fabs(dE));

		}				
		if(i==0)
		{
			m_dAbsV=m_dAbsV+fabs(dE);
			
		}
	}
	m_dAbsV=m_dAbsV/m_iPopSize;
}

//遗传算法求目标函数，并将抗体冒泡排序
void CHardy_1CodeDlg::GA_Evaluation(int gen)
{
  double a;
  int   i, j, label;
  GA_Objective();

  //抗体冒泡排序
  for(i=0; i<m_iPopSize; i++)
  {
	  label=0;  a=m_dObjective[i];
	  for(j=i+1; j<=m_iPopSize; j++)
		 if(a>m_dObjective[j])
		 {
			 a=m_dObjective[j];
			 label=j;
		 }
	  if(label!=0) 
	  {
		 a=m_dObjective[i];
		 m_dObjective[i]=m_dObjective[label];
		 m_dObjective[label]=a;

		 for(j=1; j<=m_iParamNum; j++)
		 {
			 a=m_dAntibody[i][j];
			 m_dAntibody[i][j]=m_dAntibody[label][j];
			 m_dAntibody[label][j]=a;
		 }
	  }
  }
}


//遗传算法的选择操作
void CHardy_1CodeDlg::GA_Selection()
{
  int   i, j, k;
  double r;

  //double **temp = new double* [m_iPopSize+1]; 
  //for(i = 0 ; i <=m_iPopSize+1; i++) 	temp[i] = new double[ m_iParamNum+1 ];

  double temp[POP_SIZE+1][PARAM_NUM+1];
  
  for(i=1; i<=m_iPopSize; i++)
  {
	  r=myu(0, q[m_iPopSize]);
	  for(j=0; j<=m_iPopSize; j++) 
	  {
		  if(r<=q[j]) 
		  {
			  for(k=1; k<=m_iParamNum; k++) temp[i][k]=m_dAntibody[j][k];
			  break;
		  }
	  }
  }
  for(i=1; i<=m_iPopSize; i++)
	 for(k=1; k<=m_iParamNum; k++)
		 m_dAntibody[i][k]=temp[i][k];


	/*for(i = 0 ; i <= m_iPopSize+1 ; i++) 
	delete[] temp[i]; 
	delete[] temp;*/
}


//遗传算法的交叉操作
void CHardy_1CodeDlg::GA_Crossover()
{
  int   i, j, jj, k, pop;
  double r;

  //double *x=new double[m_iParamNum+1];
  //double *y=new double[m_iParamNum+1];
  double x[PARAM_NUM+1],y[PARAM_NUM+1];

  pop=m_iPopSize/2;
  for(i=1; i<=pop; i++) 
  {
	 if(myu(0,1)>m_dCrosser) continue;
	 j=(int)myu(1,m_iPopSize);
	 jj=(int)myu(1,m_iPopSize);
	 r=myu(0,1);
	 for(k=1; k<=m_iParamNum; k++) 
	 {
		 x[k]=r*m_dAntibody[j][k]+(1-r)*m_dAntibody[jj][k];
		 y[k]=r*m_dAntibody[jj][k]+(1-r)*m_dAntibody[j][k];
	 }
	 if(GA_Constraint_check(x)==1)
		 for(k=1; k<=m_iParamNum; k++) m_dAntibody[j][k]=x[k];
	 if(GA_Constraint_check(y)==1)
		 for(k=1; k<=m_iParamNum; k++) m_dAntibody[jj][k]=y[k];
  }

  //delete[] x; 
  //delete[] y; 
}

//遗传算法的变异操作
void CHardy_1CodeDlg::GA_Mutation(void)
{
  int i, j, k;
  double infty;
  //double *x = new double[m_iParamNum+1];
  //double *y = new double[m_iParamNum+1];
  //double *direction = new double[m_iParamNum+1];
  double x[PARAM_NUM+1],y[PARAM_NUM+1],direction[PARAM_NUM+1];
  

  double INFTY=10, precision=0.0001;
  for(i=1; i<=m_iPopSize; i++)
  {
	  if(myu(0,1)>m_dMutation) continue;
	  for(k=1; k<=m_iParamNum; k++) x[k] = m_dAntibody[i][k];
	  for(k=1; k<=m_iParamNum; k++)
		  if(myu(0,1)<0.5) direction[k]=myu(-1,1);
		  else direction[k]=0;
	  infty=myu(0,INFTY);
	  while(infty>precision)
	  {
		  for(j=1; j<=m_iParamNum; j++) y[j]=x[j]+infty*direction[j];
		  if(GA_Constraint_check(y)==1)
		  {
			 for(k=1; k<=m_iParamNum; k++) m_dAntibody[i][k]=y[k];
			 break;
		  }
		  infty=myu(0,infty);
	  }
  }

 // delete[] x; 
 //delete[] y; 
 // delete[] direction; 


}

void CHardy_1CodeDlg::OnBTNOpenKnownFile() 
{
	CFileDialog file(TRUE,NULL,NULL,OFN_HIDEREADONLY|OFN_OVERWRITEPROMPT,
		         TEXT("Chart Files (*.txt)|*.txt"),NULL);
	//打开文件
	if(file.DoModal()==IDOK)
		m_strKnownFilePath=file.m_ofn.lpstrFile;
	UpdateData(false);		
}

void CHardy_1CodeDlg::OnBTNOpenParamFile() 
{
	CFileDialog file(TRUE,NULL,NULL,OFN_HIDEREADONLY|OFN_OVERWRITEPROMPT,
		         TEXT("Chart Files (*.txt)|*.txt"),NULL);
	//打开文件
	if(file.DoModal()==IDOK)
		m_strParamFilePath=file.m_ofn.lpstrFile;
	UpdateData(false);		
}

void CHardy_1CodeDlg::OnBTNOpenOutputFile() 
{
	CFileDialog file(TRUE,NULL,NULL,OFN_HIDEREADONLY|OFN_OVERWRITEPROMPT,
		         TEXT("Chart Files (*.txt)|*.txt"),NULL);
	//打开文件
	if(file.DoModal()==IDOK)
		m_strOutputFilePath=file.m_ofn.lpstrFile;
	UpdateData(false);	
	
}

void CHardy_1CodeDlg::OnChangeEDTPopSize() 
{
	UpdateData(true);
	
}

void CHardy_1CodeDlg::OnChangeEDTMaxCalcNum() 
{
	UpdateData(true);
	
}

void CHardy_1CodeDlg::OnChangeEDTCrosser() 
{
	UpdateData(true);
	
}

void CHardy_1CodeDlg::OnChangeEDTMutation() 
{
	UpdateData(true);
	
}

void CHardy_1CodeDlg::OnCancel() 
{
	CDialog::OnCancel();
}