bpalgoritim.cpp

感知机程序c语言编写

/*******************************************
* 目的：实现BP算法，数据结构采用数组的方法
*******************************************/
#include   "iostream.h"
#include   "stdio.h"
#include   "math.h"
#include   "stdlib.h" 
#include   "time.h"
#define     threshold 0.0001 //判断收敛的阈值
/************************************
* 神经单元的sigmoid函数
*************************************/
double sigmoid(double s);
/************************************
* 能量误差函数
*************************************/
double errorEnergy(double *desiredOutput, double*y,int n);
/************************************
* 产生(0-1)分布随即数的函数
*************************************/
double uniform(int a,int b, long int &seed); 
double xor(int x,int y);
main()
{
  int i;    //指向输入的源
  int j;    //指向单元节点的神经元
  int l;    //指向神经网络的某一层
  
  double    ***weightValue;   //权值
  double    **midleOut;       //各个阴层单元的输出量
  double    **delta;          //修改权值使用的改变量
  double    ***modifyWeight;  //训练时，各个权值的改变量,它的最后一层就是实际的输出 
  double    *inputX  =NULL;          //输入向量 
  double    *outputY =NULL;         //实际的输出向量
  double    *desiredOutput;   //理想输出的向量
  int       *numberCell  ;   //每一层的单元数
  int       layer ;           //神经网络的总的层数
  double    learningStep =0.1;//学习步长
  long       s;
  
  printf("press numbers of the layer of NN\n");
  scanf("%d",&layer);
  
  numberCell   = new int[layer+1];
  //第一个是输入向量，最后一个是输出向量
  midleOut     = new double*[layer+1];
  weightValue  = new double**[layer];
  modifyWeight = new double**[layer];
  
  delta        = new double*[layer];
  /****************************************************
  *  为各个指针变量分配内存单元，注意：使用完后要释放掉
  *****************************************************/
  
  for (l=0;l<=layer;l++)
  {
      
	  printf("press %d layer's numbers of cells\n", l);
      scanf("%d",&numberCell[l]);     //需不需要加&号
	  midleOut[l] = new double[numberCell[l]+1]; //加一用来添加阈值
	
  }
  for (l=0;l<layer;l++)
  {
     delta[l]  = new double[numberCell[l]+1]; 
  }
     desiredOutput = new double[numberCell[layer]];
  //第0层是输入的单元
  for(l=0;l<layer;l++)
  {
      weightValue[l]  = new double*[numberCell[l+1]];
	  modifyWeight[l] = new double*[numberCell[l+1]];

	  for(j=0;j<numberCell[l+1];j++)
	  {
	     weightValue[l][j]  = new double[numberCell[l]+1];
		 modifyWeight[l][j] = new double[numberCell[l]+1];
	  }
  }
  
  /*******************************************************
  *为权值设置初始值 随机数
  *******************************************************/
  srand( (unsigned)time( NULL ) );
  s = rand();
  for (l=0;l<layer;l++)
  {
    for (j=0;j<numberCell[l+1];j++)
    {
       for (i=0;i<=numberCell[l];i++)
       {
          weightValue[l][j][i] = uniform(0,1, s);  
	   }
	 
	}
  }
  for (l=0;l<=layer;l++)
  {
    for (j=0;j<numberCell[l];j++)
    {
       midleOut[l][j] = 0.0 ;
	}
    midleOut[l][numberCell[l]]=1; //每层的最后一个输出都是1 当作阈值
  }
  
  
  int    loops =0 ;                //纪录do--while循环的次数
  int    open_close=0;
  double energyAtZero = 0.0;       //纪录第一次的误差能量   
  double energyAtZeroUniform =0.0; //纪录每一次的归一化的值
  FILE   *fp=NULL;
  double tep =0.0;
  do {
	     if (open_close==0)
	     {
	     midleOut[0][0]=0.0;
		 midleOut[0][1]=1.0;
		 }
         else if (open_close==1)
         {
           midleOut[0][0]=1.0;
		   midleOut[0][1]=0.0;
		 }
		 else if (open_close==2)
         {
           midleOut[0][0]=0.0;
		   midleOut[0][1]=0.0;
		 }
		 else
		 {
		   midleOut[0][0]=1.0;
		   midleOut[0][1]=1.0;
		 
		 }
		 
		 open_close++;
       /*****************************************************
        * 下面这部分代码是临时替代产生理想输出的
        *****************************************************/
        for (j=0;j<numberCell[layer];j++)
		{
           desiredOutput[j]=xor(midleOut[0][0],midleOut[0][1]);
		}
      /********************************************
	  * 计算实际的输出 ；采用了三重的循环，第一重计算出一个节点的输出
	  * 第二重计算出一层的输出，第三层计算出各层的输出
      ********************************************/
	  for(l=0;l<layer;l++)
	  {	  
		  for(j=0;j<numberCell[l+1];j++)   // 这一层不到最下面 
		  {	  
			  for(i=0;i<=numberCell[l];i++) // 这一层扫描到网络的最下面
			  {
			    midleOut[l+1][j]+=weightValue[l][j][i]*midleOut[l][i]   ;
			  }
		      midleOut[l+1][j]=sigmoid(midleOut[l+1][j]);
		  }
	  }
  	 /*********************************************
	 * 训练神经网络
	 *********************************************/
     
	 /*********************************************
	 * 1)修改最后一层的权值
	 *********************************************/
     for (j=0;j<numberCell[layer];j++)
     {
       delta[layer-1][j] = midleOut[layer][j]*(1-midleOut[layer][j])
		                 *(desiredOutput[j]-midleOut[layer][j]);
	   for (i=0;i<=numberCell[layer-1];i++)
	   {
	      modifyWeight[layer-1][j][i] = delta[layer-1][j]*midleOut[layer-1][i];
		  weightValue[layer-1][j][i] += learningStep*modifyWeight[layer-1][j][i];
	   }
	   
	 }
	 /********************************************
	 * 2)修改其它层的权值
	 *********************************************/
	 for (l=layer-2; l>=0;l--)
	 {
	   for(j=0;j<numberCell[l+1];j++)
	   {
	      double sum = 0.0;
		  for (int q = 0 ; q <numberCell[l+2] ; q++)    //此处不加到底，一层中的单元的个数
		  {
		     sum +=delta[l+1][q]*weightValue[l+1][q][j];
		  }
		  delta[l][j]=midleOut[l+1][j]*(1-midleOut[l+1][j])*sum;
		  for (i=0;i<=numberCell[l];i++)    //此处加到底，对应一个单元的权数
	      {
	        modifyWeight[l][j][i] = delta[l][j]*midleOut[l][i]; 
			weightValue[l][j][i] += learningStep*modifyWeight[l][j][i];
		  }
	   }
	 }
     energyAtZeroUniform += errorEnergy(desiredOutput,midleOut[layer],numberCell[layer]);
	 if(loops<4)
	 {   
		energyAtZero  += errorEnergy(desiredOutput,midleOut[layer],numberCell[layer]);
	    if(energyAtZero==0.0)
			break;
	 } 
	else
		energyAtZeroUniform = errorEnergy(desiredOutput,midleOut[layer],numberCell[layer]);
	 energyAtZeroUniform = energyAtZeroUniform/energyAtZero;
     
    loops++;
	if (open_close==4)
	{
	 open_close=0;
	}
	if(loops>10000000)
      break;
  } while( energyAtZeroUniform > threshold);
  
  /***************************************************
   * 网络训练好以后,要输入的数据
   **************************************************/
 // for(j=0;j<numberCell[0];j++)
 // {  
 //	  printf("press input data %d\n",j);
 //	  scanf("%f",midleOut[0][j] );
 // }
  printf("the loops is %d\n",loops);
  for(open_close=0 ; open_close<4; open_close++)
  {
     
	   if (open_close==0)
	     {
	     midleOut[0][0]=0.0;
		 midleOut[0][1]=1.0;
		 }
         else if (open_close==1)
         {
           midleOut[0][0]=1.0;
		   midleOut[0][1]=0.0;
		 }
		 else if (open_close==2)
         {
           midleOut[0][0]=0.0;
		   midleOut[0][1]=0.0;
		 }
		 else
		 {
		   midleOut[0][0]=1.0;
		   midleOut[0][1]=1.0;
		 }
	    
        for(l=0;l<layer;l++)
		{	  
		   for(j=0;j<numberCell[l+1];j++)
		   {	  
			  for(i=0;i<numberCell[l];i++)
			  {
			    midleOut[l+1][j]+=weightValue[l][j][i]*midleOut[l][i]   ; 
			  }
		      midleOut[l+1][j]= sigmoid(midleOut[l+1][j]);
		  }
		}
	    printf("the result is %f\n",midleOut[layer][0]);
  }
 

  for (l = 0;l<layer;l++)
  {
     for (j=0;j<numberCell[l+1];j++)
     {
       delete[]  weightValue[l][j]  ;
	   delete[]  modifyWeight[l][j] ; 
	 }
	 delete [] weightValue[l] ;
	 delete [] modifyWeight[l];
	 delete [] midleOut[l]    ; 
	 delete [] delta[l]       ;
  }

     delete [] weightValue ;
     delete [] modifyWeight;
     delete [] midleOut;

     delete [] numberCell;
     delete [] desiredOutput;
     delete [] delta   ;
     scanf("%d",&i);
}
/********************************************************
* double *desiredOutput  理想的输出
* double *y              实际的输出
  int n                  需要的循环的次数
*********************************************************/
double errorEnergy(double *desiredOutput, double*y, int n)
{
     double sum = 0.0;
	 for ( int i =0 ; i<n ; i++ )
	 {
	    sum +=(desiredOutput[i]-y[i])*(desiredOutput[i]-y[i]) ;
	 }
     sum = sum/2;
	 return sum;
}

double sigmoid(double s)
{
  double result=0.0;
  result = 1/(1+exp((-1)*s));
  return result;
}
double uniform(int a,int b, long int &seed)
{
   double t;
   seed = 2045*(seed)+1;
   seed = seed -(seed/1048576)*1048576;
   t = seed/1048576.0;
   t= a+(b-a)*t;
   return(t);
   
}
 double xor(int x,int y)
 {
    if(x==0.0 && y==0.0)
      return 0.0;
    if(x==1.0&&y==1.0)
      return 0.0;
    if(x*y==0.0)
      return 1.0 ;
 }