?? nnarmax2.c
字號:
/*
* 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(<);
/* 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);
}
/* ---------------------------------------------------------------------*/
}
/*
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -