📄 hardy_onecodedlg.cpp

📁 利用多面函数拟合法来拟合似大地水准面
💻 CPP
字号:
// Hardy_OneCodeDlg.cpp : implementation file
//

#include "stdafx.h"
#include "GeoinModel.h"
#include "Hardy_OneCodeDlg.h"
#include <math.h>

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

/////////////////////////////////////////////////////////////////////////////
// CHardy_OneCodeDlg 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_OneCodeDlg::CHardy_OneCodeDlg(CWnd* pParent /*=NULL*/)
	: CDialog(CHardy_OneCodeDlg::IDD, pParent)
{
	//{{AFX_DATA_INIT(CHardy_OneCodeDlg)
	m_dCrosser = 0.3;
	m_strKnownFilePath = _T("");
	m_dMutation = 0.2;
	m_strParamFilePath = _T("");
	m_dq = 0.0;
	m_db = -1;
	m_iPopSize = 120;
	m_iMaxCalcNum = 5000;
	m_bOverwrite = true;
	m_strOutputFilePath = _T("");
	m_iCodeIndex = 1;
	m_iOpenOutputFile = true;
	//}}AFX_DATA_INIT
}


void CHardy_OneCodeDlg::DoDataExchange(CDataExchange* pDX)
{
	CDialog::DoDataExchange(pDX);
	//{{AFX_DATA_MAP(CHardy_OneCodeDlg)
	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);
	DDX_Text(pDX, IDC_EDTCodeIndex, m_iCodeIndex);
	DDX_Check(pDX, IDC_CHKOpenOutputFile, m_iOpenOutputFile);
	//}}AFX_DATA_MAP
}


BEGIN_MESSAGE_MAP(CHardy_OneCodeDlg, CDialog)
	//{{AFX_MSG_MAP(CHardy_OneCodeDlg)
	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)
	ON_EN_CHANGE(IDC_EDTB, OnChangeEdtb)
	//}}AFX_MSG_MAP
END_MESSAGE_MAP()

/////////////////////////////////////////////////////////////////////////////
// CHardy_OneCodeDlg message handlers

//计算
void CHardy_OneCodeDlg::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=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++;
				m_dCodeB[j]=m_dB[i];
				m_dCodeL[j]=m_dL[i];
				m_dCodeE[j]=m_dE[i];
			}
		}
	}
	fclose(File);
	m_iPtNum=i;	
	m_iCodeNum=j;	
	m_iParamNum=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;
	m_dAbsV=9999.9999;
	m_iRepeatNum=0;
	while((iGen<=m_iMaxCalcNum) && (m_dAbsV>0.005) && m_iRepeatNum<100)
	//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)
	{
		if(m_iCodeIndex==1) fprintf(File,"-9999.9999,%d\n",m_iCodeNum);	
	
		//转换参数a，平滑参数q，曲面参数b（0.5/-0.5），大地经度，大地纬度，高程异常
		fprintf(File,"%.15f,%.3f,%.1f,%.9f,%.9f,%.4f\n",
			          m_dAntibody[0][1],m_dAntibody[0][2], 
					  m_dAntibody[0][3],
					  m_dCodeB[m_iCodeIndex],m_dCodeL[m_iCodeIndex],m_dCodeE[m_iCodeIndex]);
		
		fclose(File);
		str.Format("计算完成，共迭代了%d次，请查看结果文件Hardy_GA", iGen);   
		MessageBox(str);

		if(m_iOpenOutputFile==true)
		{
			//打开文本文件
			CString Command;
			Command = "notepad ";
			Command += m_strOutputFilePath;
			WinExec(Command,SW_SHOWMAXIMIZED);
		}

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

//限制条件
int CHardy_OneCodeDlg::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_OneCodeDlg::GA_Initialization(void)
{
	int i,j,k;
	//打开参数文件	
	CString stra,strq,strb,strB,strL,strH,strE;

	//◎◎读入参数坐标文件
	if(m_strParamFilePath!="")
	{
		char buffer[256];
		FILE *File=fopen(m_strParamFilePath,"r");
		
		if(fgets(buffer,256,File)!=NULL)
		{
			for(i=1;i<=m_iCodeNum;i++)
			{	
				if(fgets(buffer,256,File)!=NULL)
				{
					stra=strtok(buffer,",， ");
					strq=strtok('\0',",， ");
					strb=strtok(  '\0',",， ");	
					
					m_dParama[i]=atof(stra);
					m_dParamq[i]=atof(strq);
					m_dParamb[i]=atof(strb);
				}
			}
		}
		fclose(File);
	}

	//◎◎◎初始化
	m_dAntibody[0][1]=m_dParama[m_iCodeIndex];
	m_dAntibody[0][2]=m_dParamq[m_iCodeIndex];
	m_dAntibody[0][3]=m_dParamb[m_iCodeIndex];

	for(i=1;i<=m_iPopSize;i++)
		for(j=1;j<=m_iCodeNum;j++)
		{
			m_dAntibody[i][1]=m_dAntibody[0][1]+myu(-0.1,0.1);
			m_dAntibody[i][2]=m_dAntibody[0][2]+myu(-0.1,0.1);
			if (myu(0.0,1)>0.5) 
				m_dAntibody[i][3]=-0.5;
			else
				m_dAntibody[i][3]=0.5;

		}

	//◎◎◎求出全局变量m_dXigmaE，除了指定节点外，其他节点的E-Q(x,y)
		double de,ds;
		for(i=1;i<=m_iPtNum;i++)
		{
			de=0;
			for(j=1;j<=m_iCodeNum;j++)
			{
				if(j!=m_iCodeIndex)
				{
					ds=0.0;
					ds=sqrt(pow(m_dB[i]-m_dCodeB[j],2)+pow(m_dL[i]-m_dCodeL[j],2)+m_dParamq[j]);

					if(m_dParamb[j]>0.5 && ds>0.0) ds=1.0/ds;  //倒双曲		
					de+=m_dParama[j]*ds;
				}
			}
			m_dXigmaE[i]=m_dE[i]-de;
		}
}

//遗传算法求目标函数值
void CHardy_OneCodeDlg::GA_Objective(void)
{	
	int i,j,k;
	double de,ds,dv;
	dv=m_dAbsV;
	m_dAbsV=0;
	for(i=0;i<=m_iPopSize;i++)  //
	{
		m_dObjective[i]=0;
		de=0;
		for(j=1;j<=m_iPtNum;j++)
		{
			ds=0.0;
			ds=sqrt(pow(m_dB[j]-m_dCodeB[m_iCodeIndex],2)+pow(m_dL[j]-m_dCodeL[m_iCodeIndex],2)+m_dAntibody[i][2]);

			if(m_dAntibody[i][3]>0.5 && ds>0.0) ds=1.0/ds;  //倒双曲	
			de=m_dAntibody[i][1]*ds;

			m_dObjective[i]+=pow(m_dXigmaE[j]-de,2);
		//	m_dObjective[i]+=exp(fabs(m_dXigmaE[i]-de));

			if(i==0)
			{				
				m_dAbsV=m_dAbsV+fabs(m_dXigmaE[j]-de);
				
			}

		}




	}	
	m_dAbsV=m_dAbsV/m_iPopSize;

	if(dv<=m_dAbsV) 
	{
		m_iRepeatNum++;
	}
	else 
		m_iRepeatNum=0;

}

//遗传算法求目标函数，并将抗体冒泡排序
void CHardy_OneCodeDlg::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_OneCodeDlg::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[POPSIZE+1][PARAMNUM+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_OneCodeDlg::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[PARAMNUM+1],y[PARAMNUM+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_OneCodeDlg::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[PARAMNUM+1],y[PARAMNUM+1],direction[PARAMNUM+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_OneCodeDlg::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_OneCodeDlg::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_OneCodeDlg::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_OneCodeDlg::OnChangeEDTPopSize() 
{
	UpdateData(true);
	
}

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

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

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

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

void CHardy_OneCodeDlg::OnChangeEdtb() 
{
	UpdateData(true);
	
}
💿 文件大小 143 K
👤 上传用户 zxj334
📂 所属分类交通/航空行业
🏷️ 相关标签

#函数
⌨️ 快捷键说明

复制代码 Ctrl + C
搜索代码 Ctrl + F
全屏模式 F11
切换主题 Ctrl + Shift + D
显示快捷键 ?
增大字号 Ctrl + =
减小字号 Ctrl + -