nnarmax2.c

类神经网路─MATLAB的应用(范例程式)

/*
 *     INCLUDE HEADERS
 */
#include <stdio.h>
#include <math.h>
#include <time.h>
#include "mex.h"
#include "matrix2.h"
#include "nnmisc.h"

void nnarmax2(matrix**, int*, double*, matrix*, matrix*, matrix*, matrix*, matrix*,\
 int, matrix*, matrix*);



/*********************************************************************************
 *                                                                               *
 *    NNARMAX2                                                                   *
 *    --------                                                                   *
 *                                                                               *
 *    This is a CMEX-version of the Matlab function nnarmax2.                    *
 *    Type 'help nnarmax2' from Matlab for information on                        *
 *    how to call this function.                                                 *
 *                                                                               *
 *                                                                               *
 *    Programmed by: Magnus Norgaard                                             *
 *    LastEditDate : sep. 04, 1995                                               *
 *                                                                               *
 *********************************************************************************/
void nnarmax2(matrix **NSSEvecpp, int *iter, double *lam,\
	matrix *NetDef, matrix *NN, matrix *W1, matrix *W2, matrix *trparms,\
	int skip, matrix *Y, matrix *U)
{
/*
----------------------------------------------------------------------------------- 
---------------              VARIABLE DECLARATIONS                    ------------- 
----------------------------------------------------------------------------------- 
*/ 
register i, j, k, t; 
int max_iter, outputs, N, Nout, layers, dummy, hidden, inputs, iteration; 
int parameters1, parameters2, parameters, reduced, index1, ii, jj; 
int lhids, hhids, louts, houts, index11;
int Ndat, N2, na, nc, nu, nab, nabc, nmax, index5, dummy2;
double stop_crit, lambda, SSE, SSE_new, NSSE, NSSE_new, L, tmp1, sum, dummy3; 
char dw; 
matrix *L_hidden, *H_hidden, *L_output, *H_output, *h1, *h2, *y1, *y2; 
matrix *E, *E_new, *W1_new, *W2_new, *PHI, *D, *Dtmp; 
matrix *NSSEvec, *miter, *tmp, *Htmp, *R, *W1tmp; 
matrix *theta, *thtmp, *theta_index, *theta_red, *theta_red_new, *PSI, *G, *H, *h; 
matrix *all, *index0, *index7, *onesvec, *tmp0, *tmp2, *tmp3, *index, *index2;
matrix *nb, *nk, *dy2dy1, *dy2de, *dy1de, *dy2de_vec, *Y2, *dummy1;
struct tm *c; 
time_t lt; 
 
 
/* 
----------------------------------------------------------------------------------- 
---------------             NETWORK INITIALIZATIONS                   ------------- 
----------------------------------------------------------------------------------- 
 */
Ndat      = getcols(Y);                  /* # of data                            */
na        = vget(NN,0);                  /* Past predictions used as inputs      */
nu	  = getrows(U);			 /* # of input signals                   */ 
nc        = vget(NN,nu+1);               /* Past prediction errors used as input */
if(nu!=0){
	nb = mmake(1,nu);                /* Past controls used as inputs	 */
            subvec(nb,NN,1,nu);
	nk = mmake(1,nu);                /* Time delays                          */
	    subvec(nk,NN,nu+2,2*nu+1);
}
nmax      = na;		                 /* Oldest signal used as input          */
if(nmax<nc) nmax=nc;
for(k=0;k<nu;k++){
  i=rvget(nb,k)+rvget(nk,k)-1;
  if(nmax<i) nmax=i;
}
N         = Ndat - nmax;                 /* Size of training set                 */
N2        = N-skip;
nab       = na; 			 /* na+nb                                */
for(k=0;k<nu;k++) nab=nab+rvget(nb,k);
nabc      = nab+nc;			 /* na+nb+nc                             */
Y2        = mmake(1,N);                  /* Observed outputs used for training   */
hidden    = getcols(NetDef);             /* # of hidden units                    */
inputs    = nabc;                 	 /* Number of inputs to network          */
outputs   = 1;		                 /* Always one outputs                   */
Nout      = N*outputs;                   /* N*outputs                            */
L_hidden  = neuvector(NetDef,1,'L');     /* Location of linear hidden units      */
H_hidden  = neuvector(NetDef,1,'H');     /* Location of tanh hidden units        */ 
L_output  = neuvector(NetDef,2,'L');     /* Location of linear output units      */ 
H_output  = neuvector(NetDef,2,'H');     /* Location of tanh output units        */ 
lhids     = getrows(L_hidden);           /* # of linear hidden units             */ 
hhids     = getrows(H_hidden);           /* # of tanh hidden units               */ 
louts     = getrows(L_output);           /* # of linear output units             */  
houts     = getrows(H_output);           /* # of tanh output units               */
miter     = mmake(1,1);                  /* Temp element                         */ 
h1        = mmake(hidden,1);             /* Argument to hidden layer act. fcts   */ 
h2        = mmake(outputs,1);            /* Argument to hidden layer act. fcts   */ 
onesvec   = mmake(1,N);                  /* Vector of all ones                   */
minitx(onesvec,1.0); 
y1        = mmake(hidden+1,N);           /* Hidden layer outputs                 */
minit(y1); 
mat2mat(y1,hidden,0,onesvec);            /* Add a row of ones (bias to outputs)  */ 
y2        = mmake(outputs,N);            /* Output layer output                  */ 
minit(y2);
E         = mmake(outputs,N);            /* Prediction error matrix              */
E_new     = mmake(outputs,N);            /* A priori E                           */ 
index     = mmake(hidden,1);             /* Index vector outputs*(hidden+1)+...  */
for(k=0;k<hidden;k++) cvput(index,k,(double)(outputs*(hidden+1)+k*(inputs+1))); 
index2    = mmake(N,1);                  /* Index vector (0:N-1)*outputs         */
for(k=0;k<N;k++) cvput(index2,k,(double)k*outputs); 
iteration = 1;                           /* Initialize iteration counter         */
dw        = 1;                           /* Flag telling that the weights are new*/ 
parameters1= hidden*(inputs+1);          /* # of input-to-hidden weights         */
parameters2= outputs*(hidden+1);         /* # of hidden-to-output weights        */ 
parameters = parameters1+parameters2;    /* Total # of weights                   */ 

/*
 >>>>>>>>>>>>>>>>>>>>  CONSTRUCT THE REGRESSION MATRIX PHI   <<<<<<<<<<<<<<<<<<<<<
 */
PHI = mmake(nabc+1,N);      	         /* Matrix of input vectors (incl. bias) */
minit(PHI);
mat2mat(PHI,nabc,0,onesvec);
for(k=0;k<na;k++){
	for(i=0;i<Ndat-nmax;i++) mput(PHI,k,i,vget(Y,i+nmax-k-1));
}
index5 = na;                             /* Insert controls in PHI               */
for(i=0;i<nu;i++){
	for(k=0;k<vget(nb,i);k++){
		for(j=0;j<Ndat-nmax;j++){
			mput(PHI,index5+k,j,mget(U,i,nmax+j-k-vget(nk,i)));
		}
	}
	index5=index5+vget(nb,i);
}
for(t=0;t<N;t++) rvput(Y2,t,rvget(Y,t+nmax));

/*
 >>>>>>>>>>>>>>>>>  INITIALIZE WEIGHTS WITH NNARX IF NECESSARY   <<<<<<<<<<<<<<<<<<
 */
if(getrows(W2)==0){
	W2->row=1;
	mrand(W1); smul(W1,W1,0.025);
   	mrand(W2); smul(W2,W2,0.5);
   	W1tmp = mmake(hidden,nab+1);
   	mrand(W1tmp); smul(W1tmp,W1tmp,0.5);
	PHI->row=nab;
	marqc(&dummy1, &dummy2, &dummy3, NetDef, W1tmp, W2, PHI, Y2, trparms);
	PHI->row=inputs+1;
   	mat2mat(W1,0,0,W1tmp);
	mfree(dummy1); mfree(W1tmp);
}

W1_new     = mmake(hidden,inputs+1);     /* A priori updated W1                  */
W2_new     = mmake(outputs,hidden+1);    /* A priori updated W2                  */ 
theta      = mmake(parameters,1);        /* Vector containing all weights        */ 
m2vreshape(theta,0,W2); 
m2vreshape(theta,parameters2,W1); 
thtmp      = mnofind(theta,0.0);         /* Find non-zero entries in theta       */ 
reduced    = getrows(thtmp);             /* # of non-zero elements               */ 
theta_index = mmake(reduced,1);          /* Indices to weights <> 0              */ 
submat(theta_index,thtmp,0,reduced-1,0,0); 
theta_red = mmake(reduced,1);            /* Reduced parameter vector             */ 
for(i=0;i<reduced;i++)                   /* theta_red = theta(theta_index)       */ 
  cvput(theta_red,i,cvget(theta,(int)cvget(theta_index,i))); 
theta_red_new = mmake(reduced,1);        /* A priori update of parameters        */ 
dy2de     = mmake(nc,N);                 /* Der. of output wrt. past pred. errors*/
dy1de     = mmake(hidden,nc);            /* Der.of hid. outp. wrt. past pred. err*/
dy2dy1    = mmake(1,hidden);             /* Der. of outp. wrt. hidden outp.      */
dy2de_vec = mmake(1,nc);		 /* For temp. results   		 */
PSI       = mmake(parameters,Nout);      /* Der. of each output wrt. each weight */
G         = mmake(reduced,1);            /* Gradient vector                      */ 
H         = mmake(reduced,reduced);      /* Hessian matrix                       */ 
R         = mmake(reduced,reduced);      /* Mean square error G-N Hessian        */
Htmp      = mmake(reduced,reduced);      /* Matrix used by the linear sys solver */ 
h         = mmake(reduced,1);            /* Update vector                        */ 
all       = mmake(N,1);                  /* Index vector (0:N-1)                 */ 
for(k=0;k<N;k++) cvput(all,k,(double)k); 
index0    = mmake(1,1);                  /* Index vector (0)                     */ 
put_val(index0,0,0,0); 
index7    = mmake(parameters,1);         /* Index vector (0:parameters-1)        */ 
for(k=0;k<parameters;k++) cvput(index7,k,(double)k); 
if (hhids>0) tmp0 = mmake(hhids,N);      /* Used to construct PSI                */
else tmp0 = mmake(1,1);
tmp2      = mmake(1,N);                  /* Used to construct PSI                */ 
tmp3      = mmake(1,N);                  /* Used to construct PSI                */ 
max_iter  = vget(trparms,0);             /* Max. no. iterations                  */
stop_crit = vget(trparms,1);             /* Error bound                          */ 
lambda    = vget(trparms,2);             /* Levenberg-Marquardt parameter        */ 
D         = mmake(reduced,1);            /* Initialize vector cont. weight decays*/ 
Dtmp      = mmake(parameters,1); 
if(length(trparms)==4)                   /* Scalar weight decay parameters       */ 
  for(i=0;i<reduced;i++) cvput(D,i,rvget(trparms,3)); 
else if(length(trparms)==5)              /* Two weight decay parameters          */ 
{ 
  for(i=0;i<parameters2;i++) cvput(Dtmp,i,rvget(trparms,3)); 
  for(i=parameters2;i<parameters;i++) cvput(Dtmp,i,rvget(trparms,4)); 
  mcopyi(D,theta_index,index0,Dtmp,index7,index0); 
} 
else{                                    /* Individual weight decays             */ 
  for(i=0;i<reduced;i++) cvput(D,i,rvget(trparms,3+i));    
}
NSSE      = stop_crit+1;                 /* Intialize cost function              */
NSSEvec = mmake(max_iter,1);             /* Vector containing normailzed SSEs    */ 
minit(NSSEvec); 


/*
----------------------------------------------------------------------------------- 
---------------                    TRAIN NETWORK                      ------------- 
----------------------------------------------------------------------------------- 
*/
lt = time(NULL); 
c  = localtime(&lt); 

/* Clear screen on HP systems.
Uncomment the following line and comment the subsequent one */
/*printf("\x1BH\x1BJNetwork training started at %.8s\n\n",asctime(c)+11);*/

printf("\nNetwork training started at %.8s\n\n",asctime(c)+11);


/* 
 >>>>>>>>>>>>>>       Compute network output y2(theta)          <<<<<<<<<<<<<<<  
*/
for(t=0;t<N;t++){
	mvmul(h1,W1,PHI,t);
	vtanh(y1,H_hidden,t,h1,H_hidden,0);
	vcopyi(y1,L_hidden,t,h1,L_hidden,0);

	mvmul(h2,W2,y1,t);
	vtanh(y2,H_output,t,h2,H_output,0);
	vcopyi(y2,L_output,t,h2,L_output,0);
	
	rvput(E,t,rvget(Y2,t)-rvget(y2,t));          /* Prediction error        */	
	j=nc;
	if(N-t-1<nc) j=N-t-1;
	for(i=1;i<=j;i++){
		put_val(PHI,nab+i-1,t+i,rvget(E,t));
	}
}
for(SSE=0,t=skip;t<N;t++) SSE+=rvget(E,t)*rvget(E,t);/* Sum of squared errors   */
for(tmp1=0,i=0;i<reduced;i++) tmp1+=cvget(theta_red,i)*cvget(theta_red,i)*cvget(D,i); 
NSSE = (SSE+tmp1)/(2*N2);                            /* Value of cost function  */

while (iteration<=max_iter && NSSE>stop_crit && lambda<1e7)
{
  if(dw==1)
  {
/*
 >>>>>>>>>>>>>>>>>>>>>>>>>>>   COMPUTE THE PSI MATRIX   <<<<<<<<<<<<<<<<<<<<<<<<<<
 (The derivative of each network output (y2) with respect to each weight)
*/
/* Some intermidiate computations */
    for(j=0;j<hhids;j++)
    {
      jj = (int)cvget(H_hidden,j);
      for(k=0;k<N;k++)
	put_val(tmp0,j,k,1-get_val(y1,jj,k)*get_val(y1,jj,k));
    }

/*   ==========   Elements corresponding to the linear output units   ===========*/
    for(i=0; i<louts; i++)
    {
      ii = (int)cvget(L_output,i);

      /***  The part of PSI corresponding to hidden-to-output layer weights ***/
      index1 = ii * (hidden+1);
      psi1(PSI, index1, index2, ii, y1);
      /************************************************************************/

      /**** The part of PSI corresponding to input-to-hidden layer weights ****/
      for(j=0; j<lhids; j++)
      {
	jj = (int)cvget(L_hidden,j);
        psi2(PSI, (int)cvget(index,jj), index2, ii, get_val(W2,ii,jj), PHI);
      }

      for(j=0; j<hhids;j++)
      {
        jj = (int)cvget(H_hidden,j);
	psi3(tmp3, tmp0, j, get_val(W2,ii,jj));
	psi4(PSI, (int)cvget(index,jj), index2, ii, tmp3, PHI);
      }
      /************************************************************************/    
    }


    /* ===========  Elements corresponding to the tanh output units   ===========*/
    for(i=0; i<houts; i++)
    {
      ii = (int)cvget(H_output,i);
      index1 = ii * (hidden + 1);
      for(k=0; k<N; k++)
	put_val(tmp2,0,k,1-get_val(y2,ii,k)*get_val(y2,ii,k));

      /* -- The part of PSI corresponding to hidden-to-output layer weights --*/
      psi4(PSI, index1, index2, ii, tmp2, y1);
      /* ---------------------------------------------------------------------*/
    
      /* -- The part of PSI corresponding to input-to-hidden layer weights ---*/
      for(j=0; j<lhids; j++)
      {
        jj = (int)cvget(L_hidden,j);
	smul(tmp3, tmp2, get_val(W2,ii,jj));
	psi4(PSI, (int)cvget(index,jj), index2, ii, tmp3, PHI);
      }
      
      for(j=0; j<hhids; j++)
      {
      	jj = (int)cvget(H_hidden,j);
	psi3(tmp3, tmp0, j, get_val(W2,ii,jj));
        psi5(PSI, (int)cvget(index,jj), index2, ii, tmp3, tmp2, PHI);
      }
      /* ---------------------------------------------------------------------*/
    }
      

    /*