bp3.c

《MATLAB实用指南》系列丛书DE源代码

//////////////////////头文件
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include "mex.h"
#define DN 30     // Number of Neural net cells
#define DM 500    // Number of Samples
#define DL 20    // Output frequency 
#define LP 0.125  // Learning velocity
#define NSS 1000   //

//////////////////////用户C程序-start /////////////////////////////////
//////////////////////////////// f1 is logsig function;  f3 is pureline function;
double f1(double x){return (double)(1/(1+exp(-x)));}
double df1(double x){double tp; tp=(double)exp(-x);return tp/((1+tp)*(1+tp));}
double f3(double x){return x;}
double df3(double x){return 1.0;}

////////////////////////////////////////////////////////////
Bp3_C()
{
int i,j,s0,s1,s2,s3,k,n,jj,n0,NStep; 
double E,arfa,gama,err,err0;
double pp[DN][DM],tt[DN][DM],tt0[DN],ee[DN];
double aa0[DN],aa1[DN],aa2[DN],aa3[DN],ss1[DN],ss2[DN],ss3[DN],nn1[DN],nn2[DN],nn3[DN];
double bb1[DN],bb2[DN],bb3[DN],dbb10[DN],dbb20[DN],dbb30[DN],dbb1[DN],dbb2[DN],dbb3[DN];
double ww1[DN][DN],ww2[DN][DN],ww3[DN][DN];
double dww1[DN][DN],dww2[DN][DN],dww3[DN][DN];
double dww10[DN][DN],dww20[DN][DN],dww30[DN][DN];
////////////////
double kk[NSS],Err[NSS]; int k_DL;
mxArray *rP[1],*lP[1],*rrP[1],*llP[1];
/////////////////////////////////////////////////////////////////

char *DataFile,*Fout,*Fss;FILE *fp0,*fp1,*fp2;
err0=5.0; E=(double)1e-10;  arfa=LP; gama=0.45; 
///////////////////////////////////////////////////

DataFile="J6.txt"; ///////////////////////////////////输入文件名
//scanf("%s",DataFile);
if((fp0=fopen(DataFile,"r"))==NULL){mexErrMsgTxt(" Can't open OutPutfile\n");}
Fout=strdup(DataFile);Fss="-r.txt";strcat(Fout,Fss);
if((fp1=fopen(Fout,"w"))==NULL){mexErrMsgTxt(" Can't open OutPutfile\n");}
Fout=strdup(DataFile);Fss="-rr.txt";strcat(Fout,Fss);
if((fp2=fopen(Fout,"w"))==NULL){mexErrMsgTxt(" Can't open OutPutfile\n");}

fscanf(fp0,"%d%d%d%d%d%d%lf%lf%d",&n,&n0,&s0,&s1,&s2,&s3,&arfa,&gama,&NStep); 
for(jj=0;jj<n;jj++){
	for(i=0;i<s0;i++){fscanf(fp0,"%lf",&pp[i][jj]);}  // pp[][] is Input  Sample
	for(i=0;i<s3;i++){fscanf(fp0,"%lf",&tt[i][jj]);}  // tt[][] is target Sample
}   n=n0; if(NStep>=NSS*DL)NStep=NSS*DL-1;

//////////////////  Initial WW & BB  within(0,1)
for(i=0;i<s1;i++){bb1[i]=rand()/32767.0;for(j=0;j<s0;j++)ww1[i][j]=rand()/32768.0;}
for(i=0;i<s2;i++){bb2[i]=rand()/32767.0;for(j=0;j<s1;j++)ww2[i][j]=rand()/32768.0;}
for(i=0;i<s3;i++){bb3[i]=rand()/32767.0;for(j=0;j<s2;j++)ww3[i][j]=rand()/32768.0;}
for(i=0;i<s1;i++){dbb10[i]=0.0;for(j=0;j<s0;j++)dww10[i][j]=0.0;}
for(i=0;i<s2;i++){dbb20[i]=0.0;for(j=0;j<s1;j++)dww20[i][j]=0.0;}
for(i=0;i<s3;i++){dbb30[i]=0.0;for(j=0;j<s2;j++)dww30[i][j]=0.0;}
/////////////////////////////////////////////////////
k_DL=k=0;
while(1){
for(jj=0;jj<n;jj++){
for(i=0;i<s0;i++)aa0[i]=pp[i][jj];
for(i=0;i<s3;i++)tt0[i]=tt[i][jj];
//////////////////  Prepropagating aa0[] to aa3[]
for(i=0;i<s1;i++){nn1[i]=0.0;for(j=0;j<s0;j++)nn1[i]=nn1[i]+ww1[i][j]*aa0[j];}
for(i=0;i<s1;i++){nn1[i]=nn1[i]+bb1[i];aa1[i]=f1(nn1[i]);}
for(i=0;i<s2;i++){nn2[i]=0.0;for(j=0;j<s1;j++)nn2[i]=nn2[i]+ww2[i][j]*aa1[j];}
for(i=0;i<s2;i++){nn2[i]=nn2[i]+bb2[i];aa2[i]=f1(nn2[i]);}
for(i=0;i<s3;i++){nn3[i]=0.0;for(j=0;j<s2;j++)nn3[i]=nn3[i]+ww3[i][j]*aa2[j];}
for(i=0;i<s3;i++){nn3[i]=nn3[i]+bb3[i];aa3[i]=f3(nn3[i]);ee[i]=tt0[i]-aa3[i];}
////////////////////////////////////////////
err=0.0; for(i=0;i<s3;i++){if(err<(double)fabs(ee[i]))err=(double)fabs(ee[i]);}
/////////////////////  Calculating sensitivity ss
for(i=0;i<s3;i++){ss3[i]=-2*df3(nn3[i])*ee[i];}
for(j=0;j<s2;j++){ss2[j]=0;for(i=0;i<s3;i++)ss2[j]=ss2[j]+df1(nn2[j])*ww3[i][j]*ss3[i];}
for(j=0;j<s1;j++){ss1[j]=0;for(i=0;i<s2;i++)ss1[j]=ss1[j]+df1(nn1[j])*ww2[i][j]*ss2[i];}
///////////////////// Calculating dww & dbb and and Renew the ww & bb
for(i=0;i<s1;i++)for(j=0;j<s0;j++)dww1[i][j]=gama*dww10[i][j]-(1-gama)*arfa*ss1[i]*aa0[j];
for(i=0;i<s2;i++)for(j=0;j<s1;j++)dww2[i][j]=gama*dww20[i][j]-(1-gama)*arfa*ss2[i]*aa1[j];
for(i=0;i<s3;i++)for(j=0;j<s2;j++)dww3[i][j]=gama*dww30[i][j]-(1-gama)*arfa*ss3[i]*aa2[j];
for(i=0;i<s1;i++)dbb1[i]=gama*dbb10[i]-(1-gama)*arfa*ss1[i];
for(i=0;i<s2;i++)dbb2[i]=gama*dbb20[i]-(1-gama)*arfa*ss2[i];
for(i=0;i<s3;i++)dbb3[i]=gama*dbb30[i]-(1-gama)*arfa*ss3[i];

for(i=0;i<s1;i++)for(j=0;j<s0;j++){ww1[i][j]=ww1[i][j]+dww1[i][j];dww10[i][j]=dww1[i][j];}
for(i=0;i<s2;i++)for(j=0;j<s1;j++){ww2[i][j]=ww2[i][j]+dww2[i][j];dww20[i][j]=dww2[i][j];}
for(i=0;i<s3;i++)for(j=0;j<s2;j++){ww3[i][j]=ww3[i][j]+dww3[i][j];dww30[i][j]=dww3[i][j];}
for(i=0;i<s1;i++){bb1[i]=bb1[i]+dbb1[i];dbb10[i]=dbb1[i];}
for(i=0;i<s2;i++){bb2[i]=bb2[i]+dbb2[i];dbb20[i]=dbb2[i];}
for(i=0;i<s3;i++){bb3[i]=bb3[i]+dbb3[i];dbb30[i]=dbb3[i];}
} k++; ///    end of for(jj=0;jj<n;jj++)

/////////////////////////////////////////////
if(k%DL==0){
kk[k_DL]=(double)(k); Err[k_DL]=err;	
fprintf(fp1,"\n%d  %lf ",k/DL,err);
//mexPrintf("\nstep=  %d  err= %10.7f  arfa= %6.4f  gama=  %6.4f  s1=%d s2=%d",k/DL,err,arfa,gama,s1,s2);
k_DL++;
}   /// end of if(k%DL==0)
if(err<E||k>NStep)break;
} ////  end of while(1)
//////////////////////////// 调用MATLAB命令作计算误差图
rP[0]=mxCreateDoubleMatrix(NSS,1,mxREAL);
memcpy(mxGetPr(rP[0]),kk,k_DL*sizeof(double));
lP[0]=mxCreateDoubleMatrix(NSS,1,mxREAL);
memcpy(mxGetPr(lP[0]),Err,k_DL*sizeof(double));
mexCallMATLAB(1,rP,1,lP,"plot");
mxDestroyArray(rP[0]);mxDestroyArray(lP[0]);
///////////////////////////////  Output The WW & BB 
fprintf(fp1,"\n\n %d %d %d %d  %d",s0,s1,s2,s3,k);
fprintf(fp1,"\n\n");for(j=0;j<s1;j++) {fprintf(fp1,"\n");
for(i=0;i<s0;i++)fprintf(fp1,"   %lf",ww1[j][i]);fprintf(fp1,"   %lf",bb1[j]);}
fprintf(fp1,"\n\n");for(j=0;j<s2;j++) {fprintf(fp1,"\n");
for(i=0;i<s1;i++)fprintf(fp1,"   %lf",ww2[j][i]);fprintf(fp1,"   %lf",bb2[j]);}
fprintf(fp1,"\n\n");for(j=0;j<s3;j++) {fprintf(fp1,"\n");
for(i=0;i<s2;i++)fprintf(fp1,"   %lf",ww3[j][i]);fprintf(fp1,"   %lf",bb3[j]);}
fprintf(fp1,"\n\n");

///////////   Forcasting   ///////////////////
/////////////////////////////////////////////////////////////////////
fscanf(fp0,"%d",&n);
for(jj=0;jj<n;jj++)for(i=0;i<s0;i++)fscanf(fp0,"%lf",&pp[i][jj]);
	for(jj=0;jj<n;jj++){
for(i=0;i<s0;i++)aa0[i]=pp[i][jj];
for(i=0;i<s1;i++){nn1[i]=0.0;for(j=0;j<s0;j++)nn1[i]=nn1[i]+ww1[i][j]*aa0[j];}
for(i=0;i<s1;i++){nn1[i]=nn1[i]+bb1[i];aa1[i]=f1(nn1[i]);}
for(i=0;i<s2;i++){nn2[i]=0.0;for(j=0;j<s1;j++)nn2[i]=nn2[i]+ww2[i][j]*aa1[j];}
for(i=0;i<s2;i++){nn2[i]=nn2[i]+bb2[i];aa2[i]=f1(nn2[i]);}
for(i=0;i<s3;i++){nn3[i]=0.0;for(j=0;j<s2;j++)nn3[i]=nn3[i]+ww3[i][j]*aa2[j];}
for(i=0;i<s3;i++){nn3[i]=nn3[i]+bb3[i];aa3[i]=f3(nn3[i]);
}/////////////////////////////////////////////////////////////////////
fprintf(fp1,"\n");for(i=0;i<s3;i++)fprintf(fp1," %lf",aa3[i]);
kk[jj]=pow(10,aa0[0]); Err[jj]=pow(10,aa3[0]);
} /// endof for(jj)
//////////////////////////// 调用MATLAB命令将预测结果以图形方式显出
rrP[0]=mxCreateDoubleMatrix(70,1,mxREAL);
memcpy(mxGetPr(rrP[0]),kk,n*sizeof(double));
llP[0]=mxCreateDoubleMatrix(70,1,mxREAL);
memcpy(mxGetPr(llP[0]),Err,n*sizeof(double));
mexCallMATLAB(1,rrP,1,llP,"plot");
mxDestroyArray(rrP[0]);mxDestroyArray(llP[0]);
///////////////////////////////////////////////////
fclose(fp0); fclose(fp1);
}

//////////////////////用户C程序-end /////////////////////////////////

//////////////////　接口程序mexFunction/////////////////////////////////////
void mexFunction(int nlhs,mxArray *plhs[],int nrhs, mxArray *prhs[])
{
int rP; 
rP=mxCalloc(1,sizeof(int));
rP=(int)mxGetScalar(prhs[0]);
if(rP==0) Bp3_C();
else mexErrMsgTxt(" To call User C_program, The paramter must be 0. \n");
}