wstspena.cpp

开发环境：Visual C++ .net2003 功能介绍：神经网络算法实验；分Console版本和MFC版本；主要用来求解TSP问题。

#include "StdAfx.h"
#include "WSTSPENA.h"
#include <math.h>
#include <fstream>
using namespace std;

#define PI 3.1415926
#define ZERO 1.0e-6
#define MINVALUE 1.0e-3
#define MBFIXPOINT 1.0e-1
#define ALPHA 0.999
#define NOTFIX false

WSTSPENA::WSTSPENA(void)
{
}

WSTSPENA::~WSTSPENA(void)
{
	delete [] m_pCity;
	delete [] m_pTempPoint;
}

//*********************************************************************************
//初始化部分
/*
功能：从指定文件中读取城市信息（城市个数、坐标）
输入：文件名p_FileName
输出：读取结果状态。true--->读取成功 false--->读取失败
*/
bool WSTSPENA::ReadCityInfoFromFile(char *p_FileName)
{	
	ifstream infile;
	infile.open(p_FileName, ios::in);
	if(!infile.is_open()) return false;
	m_FileName = p_FileName;
	
	bool needdec = false;
	CHAR l_ReadBuf[128] = {0};
	infile>>l_ReadBuf;
	if(strcmp(::strlwr(l_ReadBuf), "name") == 0)
	{
			/*
			NAME : att48
			COMMENT : 48 capitals of the US (Padberg/Rinaldi)
			TYPE : TSP
			DIMENSION : 48 <---
			EDGE_WEIGHT_TYPE : ATT
			NODE_COORD_SECTION
			*/
			infile.getline(l_ReadBuf, 255);//获得次要数据
			infile.getline(l_ReadBuf, 255);//获得次要数据
			infile.getline(l_ReadBuf, 255);//获得次要数据
			infile.getline(l_ReadBuf, 255);//获得城市个数
			int i = 0;
			while(l_ReadBuf[i++] != ':');//退出时已过了":"
			m_CityNum = atoi(l_ReadBuf+i);
			infile.getline(l_ReadBuf, 255);//获得次要数据
			infile.getline(l_ReadBuf, 255);//获得次要数据

			needdec = true;
	}
	else
	{
		m_CityNum = atoi(l_ReadBuf);
	}
	printf("城市个数：[%d]\n", m_CityNum);


	//infile>>m_CityNum;
	if( m_CityNum<= 0 ) return false;
	
	m_pCity = (WSCity*)malloc(sizeof(WSCity)*m_CityNum);
	if(!m_pCity) return false;//从全局来分配空间，以免该函数结束时将所用内存回收
	
	int i;
	int tmp;//读取城市的名字。以后可以稍加处理
	for(i = 0; i < m_CityNum; ++i)
	{
		if(!infile) return false;
		infile>>tmp>>m_pCity[i].x>>m_pCity[i].y;
		m_pCity[i].m_bFound = false;
	}

	infile.close();
	//...城市信息初始化完毕

	float maxPix = 0.0;
	for(i = 0; i < m_CityNum; ++i)
	{
		if(maxPix < m_pCity[i].x) maxPix = m_pCity[i].x;
		if(maxPix < m_pCity[i].y) maxPix = m_pCity[i].y;
	}
	m_FixTempPointDegree = maxPix/200;
	if(m_FixTempPointDegree < 0.1)m_FixTempPointDegree = 0.1;

	return true;
}
/*
功能：初始化临时点。
输入：无
输出：初始化结果状态。true--->成功 false--->失败
算法：取所有城市的重心，然后以此为圆心作圆，临时点就是均匀分布在这个圆上的。
*/
bool WSTSPENA::InitialTempPoint()
{
	if(m_TempPointFactor < 1) return false;
	
	m_TempPointNum = (int)(m_CityNum * m_TempPointFactor);
	m_pTempPoint = (WSTempPoint*)malloc( sizeof(WSTempPoint)*m_TempPointNum );

	typedef WSCity WSPointType;
	WSPointType l_Core;

	int i;
	for( i = 0; i < m_CityNum; ++i)
	{
		l_Core += m_pCity[i];
	}
	l_Core /= m_CityNum;

	//先找距重心最近的距离
	float l_MinDist = l_Core.CityDistanceSquare(m_pCity[0]);
	for(i = 1; i < m_CityNum; ++i)
	{
		if(l_MinDist > l_Core.CityDistanceSquare(m_pCity[0]))
			l_MinDist = l_Core.CityDistanceSquare(m_pCity[0]);
	}
	l_MinDist = sqrt(l_MinDist);
	if(l_MinDist < 1) l_MinDist = 1;//不要让距离太小,当然也不要太大	
	m_Temperature = l_MinDist*5;
	if(l_MinDist > 10) l_MinDist = 10;	

	printf("==========>最小路径长[%f]温度[%f]\n", l_MinDist, m_Temperature);

	//计算圆上各点
	for(i = 0; i < m_TempPointNum; ++i)
	{
		m_pTempPoint[i].m_PointPosition.x = l_Core.x + l_MinDist*sin(i*2*PI/m_TempPointNum);
		m_pTempPoint[i].m_PointPosition.y = l_Core.y + l_MinDist*cos(i*2*PI/m_TempPointNum);
		m_pTempPoint[i].m_PointPosition.m_bFound = false;
		m_pTempPoint[i].m_FoundCityNum = 0;
	}

	return true;
}

/*
功能：初始化所有参数。//参数初始化总控程序
输入：初始化需要的参数。
输出：初始化结果状态。true--->成功 false--->失败
*/
bool WSTSPENA::Initial(char *p_FileName, float p_a, float p_b, float p_Temperature,
		float p_Alpha, float p_TempPointFactor)
{
	m_a = p_a;
	m_b = p_b;
	m_Temperature = p_Temperature;
	m_Alpha = p_Alpha;
	m_TempPointFactor = p_TempPointFactor;

	//初始化各个城市的点
	ReadCityInfoFromFile(p_FileName);

	//计算临时点
	InitialTempPoint();

	return true;
}
//*********************************************************************************
/*
输入：临时点和城市距离--->p_Dist 当前温度--->p_K
输出：高斯函数计算结果
*/
float WSTSPENA::Gause(float p_Dist, float p_K)
{
	if(p_K < MINVALUE &&  p_K > -MINVALUE)
	{
		if(p_Dist < MINVALUE) return 1;
		else return 0;
	}
	return (1.0/(p_Dist*p_Dist*p_Dist*p_Dist))*exp( -(p_Dist*p_Dist) / (2*p_K*p_K) );
}
/*
输入：临时点和城市距离的平方--->p_DistSquare 当前温度的平方--->p_KSquare
输出：高斯函数计算结果
*/
float WSTSPENA::GauseSquare(float p_DistSquare, float p_KSquare)
{
	if(p_KSquare < MINVALUE)
	{
		if(p_DistSquare < MINVALUE) return 1;
		else return 0;
	}//(1.0/(p_DistSquare*p_DistSquare))*

	//return exp(-p_DistSquare)*exp( -p_DistSquare / (2*p_KSquare) );
	return ((p_KSquare*p_KSquare)/(p_DistSquare*p_DistSquare))*exp( -p_DistSquare / (2*p_KSquare) );
}
/*
输入：城市的编号--->p_CityNO 临时点的编号--->p_TempPointNO
输出：Weight计算结果的状态。true--->成功 false--->失败
问题：需要限制已找到的匹配：城市和临时点么？这里做了限制，下面的许多函数也是如是。
*/
bool WSTSPENA::Weight(unsigned int p_CityNO, unsigned int p_TempPointNO)
{
	if(p_TempPointNO >= m_TempPointNum) return false;
	if(p_CityNO >= m_CityNum) return false;

	float l_TotalValue = 0.0;
	int i;
	for(i = 0; i < m_TempPointNum; ++i) 
	{
		if(!NOTFIX) if(m_pTempPoint[i].m_PointPosition.m_bFound) continue;
		l_TotalValue += GauseSquare(
								m_pCity[p_CityNO].CityDistanceSquare(m_pTempPoint[i].m_PointPosition),
								m_Temperature * m_Temperature
								);
	}

//==================================================================
	/*
	策略：分两种情况，
	（1）所有点都不能向城市逼近时，升高温度。
	（2）当计算到一定程度后，如果某些城市不能被逼近，升高温度（缓慢升）。
	*/
	static int l_cnt = 0;//升温策略1
	static int l_cnt2 = 0;//升温策略2
	static bool l_change = false;
	if(l_TotalValue <= ZERO) 
	{
		m_Weight = 0;	
		l_cnt ++;
		l_cnt2 ++;
		if(l_change)l_cnt = 0;
		l_change = false;

		if(l_cnt >= 2) 
		{
			m_Temperature *= 20;
			l_cnt = 0;
		}
		if(l_cnt2 >= 1 && m_b < MBFIXPOINT*0.9)
		{
			m_Temperature *= 1.01;
			l_cnt2 = 0;
		}//这两个限制是为了适应一些奇异点

		//m_Weight = 0;
		//m_Temperature *= 20;
		return true;
	}
	l_change = true;
	//==================================================================


	//printf("总权[%f]\n", l_TotalValue);
	m_Weight = GauseSquare(	
						m_pCity[p_CityNO].CityDistanceSquare(m_pTempPoint[p_TempPointNO].m_PointPosition), 
						m_Temperature * m_Temperature
						) /	l_TotalValue;
	return true;
}
/*
功能：在城市点的引力和相邻临时点的张力的作用下，调整某临时点位置
输入：需要调整的临时点编号（从0起算）
输出：调整的结果状态 true--->成功 false--->失败
*/
bool WSTSPENA::AdjustTempPoint(unsigned int p_TempPointNO)
{
	if(p_TempPointNO >= m_TempPointNum) return false;

	typedef WSCity WSDeltaPoint;
	WSDeltaPoint l_DeltaYj, l_Gravitation, l_Tension ;//力的变化、引力、张力

	int i;
	//计算引力
	for(i = 0; i < m_CityNum; ++i)
	{
		if(! m_pCity[i].m_bFound)//如果城市被临时点匹配了，那么就失去了引力
		{
			if(!Weight(i, p_TempPointNO)) return false;//权重计算失败.
			l_Gravitation += (m_pCity[i] - m_pTempPoint[p_TempPointNO].m_PointPosition)*
								m_Weight;
		}
	}
	l_Gravitation *= m_a;
	
	//计算张力
	l_Tension += m_pTempPoint[(p_TempPointNO - 1 + m_TempPointNum)%m_TempPointNum ].m_PointPosition +
					m_pTempPoint[(p_TempPointNO + 1)%m_TempPointNum ].m_PointPosition -
					m_pTempPoint[ p_TempPointNO ].m_PointPosition*2;
	l_Tension *= m_b;

	l_DeltaYj = l_Gravitation + l_Tension;

	//printf("\n当前计算临时点[%d],增加值[%f,%f],权重[%f],引力[%f,%f],张力[%f,%f]\n", p_TempPointNO, l_DeltaYj.x, l_DeltaYj.y, 
	//	m_Weight, l_Gravitation.x, l_Gravitation.y, l_Tension.x, l_Tension.y);
	//printf("权重[%f]\n", m_Weight);
	//临时点调整位置
	m_pTempPoint[p_TempPointNO].m_PointPosition += l_DeltaYj;

	return true;
}

bool WSTSPENA::AdjustTempPoint(unsigned int p_TempPointNO, WSCity &p_deltaYj)
{
	if(p_TempPointNO >= m_TempPointNum) return false;

	typedef WSCity WSDeltaPoint;
	WSDeltaPoint l_Gravitation, l_Tension ;//力的变化、引力、张力

	int i;
	//计算引力
	for(i = 0; i < m_CityNum; ++i)
	{
		if(! m_pCity[i].m_bFound)//如果城市被临时点匹配了，那么就失去了引力
		{
			if(!Weight(i, p_TempPointNO)) return false;//权重计算失败.
			l_Gravitation += (m_pCity[i] - m_pTempPoint[p_TempPointNO].m_PointPosition)*
								m_Weight;
		}
	}
	l_Gravitation *= m_a;
	
	//计算张力
	l_Tension += m_pTempPoint[(p_TempPointNO - 1 + m_TempPointNum)%m_TempPointNum ].m_PointPosition +
					m_pTempPoint[(p_TempPointNO + 1)%m_TempPointNum ].m_PointPosition -
					m_pTempPoint[ p_TempPointNO ].m_PointPosition*2;
	l_Tension *= m_b;

	p_deltaYj = l_Gravitation + l_Tension;
	return true;
}

bool WSTSPENA::AdjustTempPoint()
{
	int j;
	for( j = 0; j < m_TempPointNum; ++j )
		//如果某个临时点已经被匹配了，那么就不调整它了。只调整动点
		if(! m_pTempPoint[j].m_PointPosition.m_bFound) AdjustTempPoint(j);
	return true;
}

/*
功能：同步调整临时点，效果不怎么好，所以没有使用它。
*/
bool WSTSPENA::AdjustTempPointSyn()
{
	typedef WSCity WSDeltaPoint;
	WSDeltaPoint *l_pDeltaPoint = new WSDeltaPoint[m_TempPointNum];
	int j;
	for( j = 0; j < m_TempPointNum; ++j )
		//如果某个临时点已经被匹配了，那么就不调整它了。只调整动点
	{
		if(! m_pTempPoint[j].m_PointPosition.m_bFound) AdjustTempPoint(j, l_pDeltaPoint[j]);
		else {l_pDeltaPoint[j].x = 0; l_pDeltaPoint[j].y = 0;}
	}

	for( j = 0; j < m_TempPointNum; ++j )
		//如果某个临时点已经被匹配了，那么就不调整它了。只调整动点
	{
		m_pTempPoint[j].m_PointPosition += l_pDeltaPoint[j];
	}

	delete []l_pDeltaPoint;
		
	return true;
}
/*
功能：对临时点和城市进行拟合处理，但临近点和城市足够接近时便匹配。
输入：无
输出：是否所有城市都已被匹配。true--->都匹配 false--->还有未被匹配的
*/
bool WSTSPENA::FitCity()
{
	int i, j;
	int l_FindNum = 0;
	bool l_bAllFit = true;

	for(i = 0; i < m_CityNum; ++i)
	{
		if(m_pCity[i].m_bFound)
		{
			l_FindNum++;
			continue;//如果城市被匹配了，就处理下一个
		}
		l_bAllFit = false;
		for(j = 0; j < m_TempPointNum; ++j)
		{			
			//先设置为距离平方为1的为匹配
			if(m_pCity[i].CityDistanceSquare(m_pTempPoint[j].m_PointPosition) < m_FixTempPointDegree)
			{
				m_pCity[i].m_bFound = true;
				m_pTempPoint[j].m_PointPosition = m_pCity[i];
				m_pTempPoint[j].m_FoundCityNum = i;//重新设置城市属性和临时点的属性
				break;
			}
		}
	}

	printf("已经找到的对应点城市数目：%d 温度：%f\n", l_FindNum, m_Temperature);
	return l_bAllFit;
}

bool WSTSPENA::RunENA()
{

	do{
		AdjustTempPoint();//第一次利用初始参数进行调整各临时点的位置
		AdjustParament();//参数调整		
	}while(!FitCity());

	printf("找到的城市路径总长度为：%f\n路径为:", ComputeCityRouteDist());

	int i;
	int beginCityNum, curCityNum, nextCityNum;
	
	for(i = 0; i < m_TempPointNum; ++i)
	{
		if(m_pTempPoint[i].m_PointPosition.m_bFound)
		{
			printf("%d ", m_pTempPoint[i].m_FoundCityNum);
		}
	}
	printf("\n");

	return true;
}

bool WSTSPENA::AdjustParament()
{
	m_Temperature *= m_Alpha;
	m_b *= ALPHA;//可以固定ALPHA值
	
	static float l_sDegree = m_FixTempPointDegree;

	if(m_b < MBFIXPOINT){ m_FixTempPointDegree = l_sDegree*l_sDegree;}
	else m_FixTempPointDegree = 0.01;

	if(m_Temperature < ZERO) m_Temperature = 0;
	return true;
}

float WSTSPENA::ComputeCityRouteDist()
{
	int i;
	int beginCityNum, curCityNum, nextCityNum;
	float l_CityRouteDist = 0.0;
	
	for(i = 0; i < m_TempPointNum; ++i)
	{
		if(m_pTempPoint[i].m_PointPosition.m_bFound)
		{
			beginCityNum = m_pTempPoint[i].m_FoundCityNum;
			curCityNum = beginCityNum;
			break;
		}
	}
	++i;
	for(i; i< m_TempPointNum; ++i)
	{
		if(m_pTempPoint[i].m_PointPosition.m_bFound)
		{
			nextCityNum = m_pTempPoint[i].m_FoundCityNum;
			l_CityRouteDist += m_pCity[curCityNum].CityDistance(m_pCity[nextCityNum]);
			curCityNum = nextCityNum;
		}
	}

	l_CityRouteDist += m_pCity[nextCityNum].CityDistance(m_pCity[beginCityNum]);
	return l_CityRouteDist;
}

void WSTSPENA::ShowResultRoute()
{
	ofstream outfile;
	CHAR l_OutFileName[255] = {0};
	sprintf(l_OutFileName, "%s_Route.txt", m_FileName);
	outfile.open(l_OutFileName, ios::out);
	for(int i = 0; i < m_TempPointNum; ++i)
	{
		if(m_pTempPoint[i].m_PointPosition.m_bFound)
		{
			outfile<<m_pTempPoint[i].m_PointPosition.x<<" "
				<<m_pTempPoint[i].m_PointPosition.y<<endl;
		}
	}
	outfile.close();
}