pso_nn.cpp

标准粒子群优化算法源码

/*********************************************************
PSO优化算法＋BP神经网络逼近二元分类问题程序
所有参数设置在"set.txt"配置文件中
*********************************************************/

#include "stdafx.h"
#include "pso_nn.h"
#include <stdlib.h>
#include <stdio.h>
#include <conio.h>
#include <math.h>
#include <Mmsystem.h>
#include "PS.h"

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

#define MAX_LINE_CHAR	1024
#define MAX_FN_LEN		50

int g_OriginalDimension=41;										//原始特征个数
char g_TrainFileName[MAX_FN_LEN],g_TestFileName[MAX_FN_LEN];	//训练数据文件名和测试数据文件名
int g_TrainDataNumber=5092,g_TestDataNumber=6890;				//训练数据和测试数据个数
int g_NormalNumber=0,g_AnormalNumber=0;							//正常、异常训练数据数目
double *g_TrainData;											//训练数据缓冲区指针

int g_InputNodeNumber=g_OriginalDimension;						//NN输入个数
int g_HideNodeNumber=8;											//NN隐含层节点个数
int g_OutputNumber=1;											//NN输出个数
int g_ParamNumber=g_InputNodeNumber*g_HideNodeNumber			//优化参数个数
	+g_HideNodeNumber*g_OutputNumber
	+g_HideNodeNumber+g_OutputNumber;
double g_WL=-.5;												//参数左边界
double g_WH=.5;													//参数右边界

int g_PopScale=20;												//粒子群规模
int g_MaxIter=200;												//最大迭代次数

double Evaluate(double *pw);									//粒子评价函数

/////////////////////////////////////////////////////////////////////////////
// The one and only application object

CWinApp theApp;

using namespace std;

/*********************************************************
从文件中读取一行数据
*********************************************************/
BOOL ReadLine(FILE *fp,char *buf)
{
	if(fgets(buf,MAX_LINE_CHAR,fp)==NULL)return FALSE;
	for(int i=0;i<MAX_LINE_CHAR;i++)
	{
		if(buf[i]==0x0a || buf[i]==0x0d)
		{
			buf[i]='\0';
			break;
		}
	}
	return TRUE;
}

/*********************************************************
按空格或制表符间隔解析一个数据元素
*********************************************************/
int ParseData(int pos,const char *buf,char *szRes)
{
	int i=pos,j=0;

	while((buf[i]==32 || buf[i]==9) && buf[i])i++;
	while(!(buf[i]==32 || buf[i]==9) && buf[i])
	{
		szRes[j]=buf[i];
		i++;
		j++;
	}
	szRes[j]=0;
	return i;
}

/*********************************************************
解析参数配置文件中的一个设定值
*********************************************************/
void GetValue(char *from,char *to)
{
	sprintf(to,strstr(from,"=")+1);
	return;
}

/*********************************************************
从配置文件中读取各参数配置
*********************************************************/
BOOL GetSet()
{
	FILE *f=fopen("set.txt","r");
	if(!f)return FALSE;
	char str[200],value[MAX_FN_LEN];
	
	ReadLine(f,str);
	GetValue(str,value);
	g_OriginalDimension=atoi(value);

	ReadLine(f,str);
	GetValue(str,g_TrainFileName);

	ReadLine(f,str);
	GetValue(str,g_TestFileName);

	ReadLine(f,str);
	GetValue(str,value);
	g_TrainDataNumber=atoi(value);

	ReadLine(f,str);
	GetValue(str,value);
	g_TestDataNumber=atoi(value);

	ReadLine(f,str);
	GetValue(str,value);
	g_HideNodeNumber=atoi(value);

	ReadLine(f,str);
	GetValue(str,value);
	g_OutputNumber=atoi(value);

	ReadLine(f,str);
	GetValue(str,value);
	g_PopScale=atoi(value);

	ReadLine(f,str);
	GetValue(str,value);
	g_MaxIter=atoi(value);

	g_InputNodeNumber=g_OriginalDimension;
	g_ParamNumber=g_InputNodeNumber*g_HideNodeNumber
		+g_HideNodeNumber*g_OutputNumber
		+g_HideNodeNumber+g_OutputNumber;

	fclose(f);
	return TRUE;
}

/*********************************************************
读取训练数据
*********************************************************/
BOOL ReadTrainData()
{
	FILE *fp=fopen(g_TrainFileName,"r");
	if(!fp)return FALSE;
	
	char buf[MAX_LINE_CHAR];
	char pData[MAX_LINE_CHAR];
	int nPos;
	int i,j;

	g_TrainData=new double[g_TrainDataNumber*(g_OriginalDimension+1)];
	for(i=0;i<g_TrainDataNumber;i++)
	{
		if(feof(fp))break;
		ReadLine(fp,buf);
		nPos=0;
		for(j=0;j<g_OriginalDimension+1;j++)
		{
			nPos=ParseData(nPos,buf,pData);
			g_TrainData[i*(g_OriginalDimension+1)+j]=atof(pData);
		}
		if(g_TrainData[i*(g_OriginalDimension+1)+g_OriginalDimension]==1.)
			g_NormalNumber++;
		else
			g_AnormalNumber++;
	}
	fclose(fp);
	return TRUE;
}

/*********************************************************
计算NN节点输出函数
*********************************************************/
double  fire(double x, double a)
{
	double e;
	e = exp(-2.0*a*x);
	return  (1.0/(1.0 + e));
}

/*********************************************************
计算NN输出函数
*********************************************************/
double *CalculateNNOutput(double *pi,double *pw)
{
	int i,j;
	double *p1=pi,*p2=pw;
	double *HideOutput=new double[g_HideNodeNumber];
	double *Output=new double[g_OutputNumber];
	int WeightNumber=g_InputNodeNumber*g_HideNodeNumber
		+g_HideNodeNumber*g_OutputNumber;

	for(i=0;i<g_HideNodeNumber;i++)
	{
		double sum=0.;
		p1=pi;
		for(j=0;j<g_InputNodeNumber;j++,p1++,p2++)
			sum+=(*p1)*(*p2);
		HideOutput[i]=fire(sum,pw[WeightNumber+i]);
	}
	for(i=0;i<g_OutputNumber;i++)
	{
		double sum=0.;
		p1=HideOutput;
		for(j=0;j<g_HideNodeNumber;j++,p1++,p2++)
			sum+=(*p1)*(*p2);
		Output[i]=fire(sum,pw[WeightNumber+g_HideNodeNumber+i]);
	}

	delete []HideOutput;
	return Output;
}

/*********************************************************
粒子外部评价函数
*********************************************************/
double Evaluate(double *pw,int &err)
{
	int i;
	double *p,*po;
	int tp=0,tn=0,fp=0,fn=0;
	double f;

	for(i=0;i<g_TrainDataNumber;i++)
	{
		p=g_TrainData+i*(g_OriginalDimension+1);
		po=CalculateNNOutput(p,pw);

		if(g_TrainData[i*(g_OriginalDimension+1)+g_OriginalDimension]==1. && po[0]>=.5)
			tp++;
		if(g_TrainData[i*(g_OriginalDimension+1)+g_OriginalDimension]==1. && po[0]<.5)
			tn++;
		if(g_TrainData[i*(g_OriginalDimension+1)+g_OriginalDimension]==0. && po[0]<.5)
			fn++;
		if(g_TrainData[i*(g_OriginalDimension+1)+g_OriginalDimension]==0. && po[0]>=.5)
			fp++;

		delete []po;
	}
//	f=(double)(tp+fn)/g_TrainDataNumber;
	f=(double)tp/g_NormalNumber+(double)fn/g_AnormalNumber;
	f/=2.;
	err=tn+fp;

	return f;
}

/*********************************************************
保存优化结果函数
*********************************************************/
void SaveWParam(double *pw)
{
	FILE *fp=fopen("weight.txt","w");
	if(!fp)return;
	for(int i=0;i<g_ParamNumber;i++)
		fprintf(fp,"%lf\n",pw[i]);
	fclose(fp);
}

void TestByChromosome(double *pw)
{
	FILE *fp=fopen(g_TestFileName,"r");
	double *x=new double[g_OriginalDimension+1];
	double *po;
	char buf[MAX_LINE_CHAR],pData[20];
	int nPos,j;
	int TP=0,TN=0,FP=0,FN=0;
	int NormalNumber=0,AnormalNumber=0;
	for(int i=0;i<g_TestDataNumber;i++)
	{
		if(feof(fp))break;
		ReadLine(fp,buf);
		nPos=0;
		for(j=0;j<g_OriginalDimension+1;j++)
		{
			nPos=ParseData(nPos,buf,pData);
			x[j]=atof(pData);
		}
		po=CalculateNNOutput(x,pw);

		if(x[g_OriginalDimension]==1.)
		{
			NormalNumber++;
			if(po[0]>=.5)
				TP++;
			else
				TN++;
		}
		else
		{
			AnormalNumber++;
			if(po[0]<.5)
				FN++;
			else
				FP++;
		}
		delete []po;
	}
	fclose(fp);
	cout<<"Error number:"<<TN+FP<<endl;
	cout<<"TP="<<TP<<"  FP="<<FP<<endl;
	cout<<"Accuracy="<<(double)(TP+FN)/g_TestDataNumber*100.<<"%  tp="
		<<(double)TP/NormalNumber*100.<<"%  fp="<<(double)FP/AnormalNumber*100.<<"%"<<endl;
	delete []x;
}

/*********************************************************
综合测试函数
*********************************************************/
void Test()
{
	FILE *fp=fopen("weight.txt","r");
	if(!fp)return;
	double *pw=new double[g_ParamNumber];
	for(int i=0;i<g_ParamNumber;i++)
		fscanf(fp,"%lf\n",pw+i);
	fclose(fp);
	fp=NULL;

	fp=fopen(g_TestFileName,"r");
	double *x=new double[g_OriginalDimension+1];
	char buf[MAX_LINE_CHAR],pData[20];
	int nPos,j;
	int TP=0,TN=0,FP=0,FN=0;
	int NormalNumber=0,AnormalNumber=0;
	for(i=0;i<g_TestDataNumber;i++)
	{
		if(feof(fp))break;
		ReadLine(fp,buf);
		nPos=0;
		for(j=0;j<g_OriginalDimension+1;j++)
		{
			nPos=ParseData(nPos,buf,pData);
			x[j]=atof(pData);
		}
		double *po=CalculateNNOutput(x,pw);
		if(x[g_OriginalDimension]==1.)
		{
			NormalNumber++;
			if(po[0]>=.5)
				TP++;
			else
				TN++;
		}
		else
		{
			AnormalNumber++;
			if(po[0]<.5)
				FN++;
			else
				FP++;
		}
		delete []po;
	}
	fclose(fp);
	cout<<"Error number:"<<TN+FP<<endl;
	cout<<"TP="<<TP<<"  FP="<<FP<<endl;
	cout<<"Accuracy="<<(double)(TP+FN)/g_TestDataNumber*100.<<"%  tp="
		<<(double)TP/NormalNumber*100.<<"%  fp="<<(double)FP/AnormalNumber*100.<<"%"<<endl;
	delete []pw;
	delete []x;
}

int _tmain(int argc, TCHAR* argv[], TCHAR* envp[])
{
	int nRetCode = 0;

	if (!AfxWinInit(::GetModuleHandle(NULL), NULL, ::GetCommandLine(), 0))
	{
		cerr << _T("Fatal Error: MFC initialization failed") << endl;
		nRetCode = 1;
	}
	else
	{
		cout<<"********PSO优化算法＋BP神经网络逼近二元分类问题程序********"<<endl;
		
		GetSet();
		ReadTrainData();
		char cmd[20];
		CPS *ps=new CPS(g_PopScale,g_ParamNumber,g_WL,g_WH);
		while(strcmp(cmd,"bye"))
		{
			cout<<"Input command:"<<endl;
			cout<<"'c' to optimize the neural network"<<endl;
			cout<<"'t' to test the neural network"<<endl;
			cout<<"'i' to reinitiate the particle swarm"<<endl;
			cout<<"'bye' to quit"<<endl;
			cout<<"command:";
			cin>>cmd;
			if(strcmp(cmd,"c")==NULL)
			{
				int start_time=timeGetTime();
				ps->Iterate(g_MaxIter);
				int end_time=timeGetTime();
				cout<<"Use "<<end_time-start_time<<" millisecond!"<<endl;
				SaveWParam(ps->m_gbest);
			}
			if(strcmp(cmd,"t")==NULL)
				Test();
			if(strcmp(cmd,"i")==NULL)
				ps->InitPop();
		}

		delete ps;
	}

	return nRetCode;
}