📄 nnssif.c
字号:
/*
* INCLUDE HEADERS
*/
#include <stdio.h>
#include <math.h>
#include <time.h>
#include "mex.h"
#include "matrix2.h"
#include "nnmisc.h"
void nnssif(matrix**, int*, double*, matrix*, int, matrix*, matrix*, matrix*,\
matrix*, int, matrix*, matrix*);
/*********************************************************************************
* *
* NNSSIF *
* ------ *
* *
* This is a CMEX-version of the Matlab function nnssif. *
* Type 'help nnssif' from Matlab for information on *
* how to call this function. *
* *
* *
* Programmed by: Magnus Norgaard *
* LastEditDate : sep. 04, 1995 *
* *
*********************************************************************************/
void nnssif(matrix **NSSEvecpp, int *iter, double *lam,\
matrix *NetDef, int nx, matrix *W1, matrix *W2, matrix *obsidx,\
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, hid1, hid2, Nny, nxu, skipstart;
int Ndat, N2, nu, ny, index5, index6, dummy2;
double stop_crit, lambda, SSE, SSE_new, NSSE, NSSE_new, L, tmp1, sum, dummy3;
char dw, stateflag;
matrix *L_hidden, *H_hidden, *L_output, *H_output, *h1, *h2, *y1, *y2;
matrix *E, *Evec, *Evec_new, *W1_new, *W2_new, *PHI, *D, *Dtmp;
matrix *NSSEvec, *miter, *tmp, *Htmp, *R, *PSIx, *Yhat;
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 *rowidx, *nrowidx, *dxdy1, *Ahat, *Khat, *dy1de, *dy1dx, *Y2, *dummy1, *C;
matrix *AKC;
struct tm *c;
time_t lt;
/*
-----------------------------------------------------------------------------------
--------------- NETWORK INITIALIZATIONS -------------
-----------------------------------------------------------------------------------
*/
Ndat = getcols(Y); /* # of data */
ny = getrows(Y); /* # of outputs */
nu = getrows(U); /* # of controls */
N = Ndat - 1; /* Size of training set */
Nny = N*ny;
nxu = nx+nu;
N2 = N-skip;
skipstart = ny*skip;
Y2 = mmake(ny,N); /* Observed outputs used for training */
Yhat = mmake(ny,1); /* Output prediction */
hidden = getcols(NetDef); /* # of hidden units */
inputs = nx+nu+ny; /* Number of inputs to network */
outputs = nx; /* Always one 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(ny,1); /* Prediction error */
Evec = mmake(1,Nny); /* Prediction error vector */
Evec_new = mmake(1,Nny); /* A priori Evec */
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 */
rowidx = mmake(ny,1); /* Row indices */
vput(rowidx,0,vget(obsidx,0));
for(k=1;k<ny;k++) vput(rowidx,k,vget(obsidx,k)+vget(rowidx,k-1));
for(k=0;k<ny;k++) vput(rowidx,k,vget(rowidx,k)-1.0);
if(nx-ny!=0){
nrowidx = mmake(nx-ny,1); /* Not row indices */
for(j=0,k=0;k<nx;k++){
stateflag=0;
for(i=0;i<ny;i++){
if((int)vget(rowidx,i)==k) stateflag=1;
}
if(stateflag==0) vput(nrowidx,j++,(double)k);
}
}
else{
nrowidx=mmake(1,1);
nrowidx->row = 0;
nrowidx->col = 0;
}
/* Initialize weights if necessary */
if(getrows(W2)==0){
W2->row=outputs;
mrand(W1); smul(W1,W1,0.025);
mrand(W2); smul(W2,W2,0.025);
}
/* Insert zeros to ensure observability */
hid1 = (int)floor(hidden/2+0.5);
hid2 = hidden - hid1;
for(i=hid1;i<hidden;i++){
for(j=0;j<nx;j++) mput(W1,i,j,0.0);
for(j=0;j<ny;j++) mput(W2,(int)vget(rowidx,j),i,0.0);
}
for(i=0;i<hid1;i++){
for(j=0;j<(nx-ny);j++) mput(W2,(int)vget(nrowidx,j),i,0.0);
}
/* Observation matrix */
C = mmake(ny,nx); minit(C);
mput(C,0,0,1.0);
for(i=1;i<ny;i++) mput(C,i,vget(rowidx,i-1)+1,1.0);
Ahat = mmake(nx,nx); /* Deriv. of states wrt. past states */
Khat = mmake(nx,ny); /* Deriv. of states wrt. past residuals */
AKC = mmake(nx,nx); /* Stores temp. results */
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 */
dy1dx = mmake(hidden,nx); /* Der. of hid. outp. wrt. past states */
dy1de = mmake(hidden,ny); /* Der.of hid. outp. wrt. past pred. err*/
dxdy1 = mmake(nx,hidden); /* Der. of state estim. wrt. hid. outp. */
PSI = mmake(parameters,N*ny); /* Der. of each output wrt. each weight */
PSIx = mmake(parameters,N*nx); /* Der. of est. states wrt. each weight */
minit(PSIx);
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);
/*
>>>>>>>>>>>>>>>>>>>> CONSTRUCT THE REGRESSION MATRIX PHI <<<<<<<<<<<<<<<<<<<<<
*/
PHI = mmake(inputs+1,N); /* Matrix of input vectors (incl. bias) */
minit(PHI);
mat2mat(PHI,inputs,0,onesvec); /* Insert biases in PHI */
for(i=0;i<nu;i++){ /* Insert controls in PHI */
for(t=0;t<N;t++)
mput(PHI,nx+i,t,mget(U,i,t));
}
for(i=0;i<ny;i++){
for(t=0;t<N;t++) mput(Y2,i,t,mget(Y,i,t+1));
}
/*
-----------------------------------------------------------------------------------
--------------- 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);
for(k=0;k<(nx-ny);k++){
i=(int)cvget(nrowidx,k);
y2->mat[i][t]+=get_val(PHI,i+1,t);
}
mvmul(Yhat,C,y2,t); /* Output prediction */
for(i=0;i<ny;i++){
cvput(E,i,get_val(Y2,i,t)-cvget(Yhat,i));/* Prediction error */
rvput(Evec,t*ny+i,cvget(E,i)); /*Store E in vector Evec*/
}
if(t<N-1){
for(i=0;i<nx;i++)
put_val(PHI,i,t+1,get_val(y2,i,t));
for(i=0;i<ny;i++)
put_val(PHI,nx+nu+i,t+1,cvget(E,i));
}
}
for(SSE=0,t=skipstart;t<Nny;t++)
SSE+=rvget(Evec,t)*rvget(Evec,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 PSIx corresponding to hidden-to-output layer weights ***/
index1 = ii * (hidden+1);
psi1(PSIx, index1, index2, ii, y1);
/************************************************************************/
/**** The part of PSIx corresponding to input-to-hidden layer weights ****/
for(j=0; j<lhids; j++)
{
jj = (int)cvget(L_hidden,j);
psi2(PSIx, (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(PSIx, (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 PSIx corresponding to hidden-to-output layer weights --*/
psi4(PSIx, index1, index2, ii, tmp2, y1);
/* ---------------------------------------------------------------------*/
/* -- The part of PSIx 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(PSIx, (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(PSIx, (int)cvget(index,jj), index2, ii, tmp3, tmp2, PHI);
}
/* ---------------------------------------------------------------------*/
}
/*
>>>>>>>>>>>>>>>>>>>> Linearize network <<<<<<<<<<<<<<<<<<<<<
*/
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -