?? nnoe.c
字號:
/*
>>>>>>>>>>>>>>>>>>>> Linearize network <<<<<<<<<<<<<<<<<<<<<
*/
for(t=0;t<N;t++){
/*-- Derivative of output wrt. hidden outputs --*/
if(louts==1) for(k=0;k<hidden;k++) rvput(dy2dy1,k,rvget(W2,k));
else for(k=0;k<hidden;k++) rvput(dy2dy1,k,rvget(W2,k)*(1-\
rvget(y2,t)*rvget(y2,t)));
/*-- Partial deriv. of output from each hidden unit wrt. past net outp. --*/
for(j=0;j<lhids;j++){
i=(int)cvget(L_hidden,j);
for(k=0;k<na;k++) put_val(dy1dy,i,k,\
get_val(W1,i,k));
}
for(j=0;j<hhids;j++){
i=(int)cvget(H_hidden,j);
for(k=0;k<na;k++) put_val(dy1dy,i,k,\
get_val(W1,i,k)*(1-get_val(y1,i,t)*get_val(y1,i,t)));
}
/*--Partial derivative of net output w.r.t. past net outputs --*/
mmul(dy2dy_vec,dy2dy1,dy1dy);
for(i=0;i<na;i++) put_val(dy2dy,i,t,rvget(dy2dy_vec,i));
}
/*
>>>>>>>>>>>>>>>>>>> Filter partial derivatives <<<<<<<<<<<<<<<<<<<<
*/
for(t=0;t<N;t++){
j=na;
if(t<na) j=t;
for(k=1;k<=j;k++){
for(i=0;i<reduced;i++) {
ii =(int)cvget(theta_index,i);
PSI->mat[ii][t] += get_val(dy2dy,k-1,t)*get_val(PSI,ii,t-k);
}
}
}
dw = 0;
/*
>>>>>>>>>>>> Gradient (G = PSI_red*E_vector - D*theta_red) <<<<<<<<<<<<<
*/
for(i=0; i<reduced; i++){
ii = (int)cvget(theta_index,i);
for(sum=0.0,k=skip; k<N; k++) sum+=get_val(PSI,ii,k)*rvget(E,k);
cvput(G,i,sum - cvget(D,i)*cvget(theta_red,i));
}
/*
>>>>>>>>>> Mean square error part of Hessian (PSI_red*PSI_red') <<<<<<<<<<
*/
for(i=0; i<reduced; i++){
ii = (int)cvget(theta_index,i);
for(j=i; j<reduced; j++){
jj = (int)cvget(theta_index,j);
for(sum=0.0,k=skip; k<N; k++)
sum += get_val(PSI,ii,k)*get_val(PSI,jj,k);
put_val(R,i,j,sum);
put_val(R,j,i,sum);
}
}
}
/*
>>>>>>>>>>>>>>>>>>>>>>>>>>> COMPUTE h_k <<<<<<<<<<<<<<<<<<<<<<<<<<<
*/
/* -- Hessian (H = R + lambda*I + D) --*/
mset(H,R);
for(i=0;i<reduced;i++) /* Add diagonal matrix */
put_val(H,i,i,get_val(H,i,i)+lambda+cvget(D,i));
/* -- Search direction -- */
choldc(H, Htmp);
cholsl(Htmp,h,G);
/* -- Compute 'apriori' iterate -- */
madd(theta_red_new,theta_red,h); /* Update parameter vector */
mcopyi(theta,theta_index,index0,theta_red_new,index7,index0);
/* -- Put the parameters back into the weight matrices -- */
v2mreshape(W1_new,theta,parameters2);
v2mreshape(W2_new,theta,0);
/*
>>>>>>>>>>>>> Compute network output y2(theta+h) <<<<<<<<<<<<<<
*/
for(t=0;t<N;t++){
mvmul(h1,W1_new,PHI,t);
vtanh(y1,H_hidden,t,h1,H_hidden,0);
vcopyi(y1,L_hidden,t,h1,L_hidden,0);
mvmul(h2,W2_new,y1,t);
vtanh(y2,H_output,t,h2,H_output,0);
vcopyi(y2,L_output,t,h2,L_output,0);
j=na;
if(N-t-1<na) j=N-t-1;
for(i=1;i<=j;i++){
put_val(PHI,i-1,t+i,rvget(y2,t));
}
}
for(t=0;t<N;t++) /* Prediction error */
rvput(E_new,t,rvget(Y2,t)-rvget(y2,t));
for(SSE_new=0,t=skip;t<N;t++) /* Sum of squared errors */
SSE_new+=rvget(E_new,t)*rvget(E_new,t);
for(tmp1=0,i=0;i<reduced;i++) tmp1+=cvget(theta_red,i)*cvget(theta_red,i)*cvget(D,i);
NSSE_new = (SSE_new+tmp1)/(2*N2); /* Value of cost function*/
/*
>>>>>>>>>>>>>>>>>>>>>>>>>>> UPDATE lambda <<<<<<<<<<<<<<<<<<<<<<<<<<
*/
for(tmp1=0,i=0;i<reduced;i++) tmp1+=cvget(h,i)*cvget(h,i)*(cvget(D,i)+lambda);
L = sprod3(h,G) + tmp1;
/* Decrease lambda if SSE has fallen 'sufficiently' */
if(2*N2*(NSSE - NSSE_new) > (0.75*L)) lambda = lambda/2;
/* Increase lambda if SSE has grown 'sufficiently' */
else if(2*N2*(NSSE-NSSE_new) <= (0.25*L)) lambda = 2*lambda;
/*
>>>>>>>>>>>>>>>>>>> UPDATES FOR NEXT ITERATION <<<<<<<<<<<<<<<<<<<<
*/
/* Update only if criterion has decreased */
if(NSSE_new<NSSE)
{
tmp = W1; W1=W1_new; W1_new=tmp;
tmp = W2; W2=W2_new; W2_new=tmp;
tmp = theta_red; theta_red=theta_red_new; theta_red_new = tmp;
tmp = E; E = E_new; E_new = tmp;
dw = 1;
NSSE = NSSE_new;
cvput(NSSEvec,iteration-1,NSSE);
printf("iteration # %i PI = %4.3e\r",iteration,NSSE); /* On-line information */
++iteration;
}
}
/*
>>>>>>>>>> RETURN POINTERS TO RETURN ARGUMENTS & FREE MEMORY <<<<<<<<<<<<
*/
/* Swap pointers if they have been messed up */
if ((iteration&1) == 0) {
mset(W1_new,W1);
tmp = W1; W1=W1_new; W1_new=tmp;
mset(W2_new,W2);
tmp = W2; W2=W2_new; W2_new=tmp;
}
iteration=iteration-1;
if(iteration==0){
*NSSEvecpp = mmake(1,1);
(*NSSEvecpp)->row=0;
(*NSSEvecpp)->col=0;
}
else {
*NSSEvecpp = mmake(iteration,1);
subvec(*NSSEvecpp,NSSEvec,0,iteration-1);
}
*iter = iteration;
*lam = lambda;
mfree(L_hidden); mfree(H_hidden); mfree(L_output); mfree(H_output);
mfree(h1); mfree(h2); mfree(y1); mfree(y2); mfree(Y2);
mfree(E); mfree(E_new); mfree(dy2dy1); mfree(dy2dy); mfree(dy1dy); mfree(dy2dy_vec);
mfree(W1_new); mfree(W2_new); mfree(D);mfree(Dtmp); mfree(NSSEvec);mfree(Htmp);
mfree(theta); mfree(thtmp); mfree(theta_index); mfree(theta_red); mfree(theta_red_new);
mfree(PSI); mfree(G); mfree(H); mfree(R); mfree(h);
mfree(all); mfree(index0); mfree(index7);mfree(onesvec);mfree(tmp0); mfree(tmp2);
mfree(tmp3); mfree(index); mfree(index2); mfree(PHI); mfree(nb); mfree(nk);
lt = time(NULL);
c = localtime(<);
printf("\n\nNetwork training ended at %.8s\n\n\n",asctime(c)+11);
}
/*
--------------------------------------------------------------------------------
---------------- END OF NETWORK TRAINING --------------
--------------------------------------------------------------------------------
*/
/*********************************************************************************
* *
* G A T E W A Y R O U T I N E *
* *
*********************************************************************************/
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
/*
>>>>>>>>>>>>>>>>>> VARIABLE DECLARATIONS <<<<<<<<<<<<<<<<<<<
*/
matrix *NSSEvec, *NN, *nb;
matrix *NetDef, *W1, *W2, *U, *Y, *trparms;
double *M, lambda;
int iter, skip, m, N, na, nab, nu, hidden, k;
/*
>>>>>>>>>>>>>>>> CHECK FOR PROPER NUMBER OF ARGUMENTS <<<<<<<<<<<<<<<
*/
if (nrhs!=8) mexErrMsgTxt("Wrong number of input arguments");
else if (nlhs > 5) mexErrMsgTxt("Two many output arguments");
nu = mxGetM(prhs[7]); /* # of control signals */
m = mxGetM(prhs[6]); /* Rows of vector Y */
if(m!=1) mexErrMsgTxt("Wrong dimension of vector of observed outputs");
N = mxGetN(prhs[6]); /* Columns of vector Y */
if(N!=mxGetN(prhs[6])) mexErrMsgTxt("U and Y should have the same number of columns!");
hidden = mxGetN(prhs[0]); /* # of hidden units */
if(mxGetM(prhs[0])!=2) mexErrMsgTxt("Error in architecture definition!");
if(mxGetM(prhs[1])!=1) mexErrMsgTxt("NN should only have 1 row!");
if(mxGetN(prhs[1])!=2*nu+1) mexErrMsgTxt("Mismatch between U and NN");
/*
>>>>>>>>>>>>>>>>> CONVERT INPUT ARGUMENTS TO SM FORMAT <<<<<<<<<<<<<<<<
*/
NetDef = matstring2sm(prhs[0]); /* Network architecture */
NN = mat2sm(prhs[1]); /* Regressor structure */
if(mxGetM(prhs[4])!=0) /* Parameters associated with minimization */
trparms = mat2sm(prhs[4]);
else { /* Use defaults if passed as [] */
trparms = mmake(1,4);
mput(trparms,0,0,500);
mput(trparms,0,1,0.0);
mput(trparms,0,2,1);
mput(trparms,0,3,0.0);
}
if(mxGetM(prhs[5])!=0) /* Skip first gradients */
skip = (int)(*mxGetPr(prhs[5]));
else skip=0; /* Use skip=0 if passed as [] */
Y = mat2sm(prhs[6]); /* Vector of observed outputs */
U = mat2sm(prhs[7]); /* Vector of inputs */
/*
>>>>>>>>>>>>>>>> CHECK FOR PROPER NUMBER OF ARGUMENTS <<<<<<<<<<<<<<<
*/
na = vget(NN,0); /* Past predictions used as inputs */
nb = mmake(1,nu); /* Past controls used as input */
subvec(nb,NN,1,nu);
nab = na; /* na+nb */
for(k=0;k<nu;k++) nab=nab+rvget(nb,k);
/* Initialize weight matrices if passed as [] */
if(mxGetM(prhs[2])==0 || mxGetN(prhs[2])==0 || mxGetM(prhs[3])==0\
|| mxGetN(prhs[3])==0){
W1 = mmake(hidden,nab+1);
W2 = mmake(1,hidden+1);
W2->row = 0; /* Hack telling that the weights should be initialized */
}
else{
if(mxGetM(prhs[2])!=hidden) mexErrMsgTxt("W1 has the wrong dimension");
if(mxGetN(prhs[2])!=nab+1) mexErrMsgTxt("W1 has the wrong dimension");
if(mxGetM(prhs[3])!=1) mexErrMsgTxt("W2 has the wrong dimension");
if(mxGetN(prhs[3])!=hidden+1) mexErrMsgTxt("W2 has the wrong dimension");
W1 = mat2sm(prhs[2]); /* Input-to-hidden layer weights */
W2 = mat2sm(prhs[3]); /* Hidden-to-output layer weights */
}
/*
>>>>>>>>>>>>>>>>>>>>>> CALL THE C-ROUTINE <<<<<<<<<<<<<<<<<<<<<
*/
nnoe(&NSSEvec, &iter, &lambda, NetDef, NN, W1, W2, trparms, skip, Y, U);
/*
>>>>>>>>>>>>>>>>>>> CREATE OUTPUT MATICES <<<<<<<<<<<<<<<<<<
*/
plhs[0] = mxCreateDoubleMatrix(getrows(W1),getcols(W1),mxREAL);
plhs[1] = mxCreateDoubleMatrix(getrows(W2),getcols(W2),mxREAL);
plhs[2] = mxCreateDoubleMatrix(getrows(NSSEvec),getcols(NSSEvec),mxREAL);
plhs[3] = mxCreateDoubleMatrix(1,1,mxREAL);
plhs[4] = mxCreateDoubleMatrix(1,1,mxREAL);
sm2mat(plhs[0],W1);
sm2mat(plhs[1],W2);
sm2mat(plhs[2],NSSEvec);
M = mxGetPr(plhs[3]); M[0] = (double)iter;
M = mxGetPr(plhs[4]); M[0] = (double)lambda;
/*
>>>>>>>>>>>>>>>>>>>> FREE ARGUMENT MATRICES <<<<<<<<<<<<<<<<<<<<<
*/
mfree(NetDef);
mfree(NN);
mfree(U);
mfree(Y);
mfree(trparms);
}
?? 快捷鍵說明
復(fù)制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -