/*  1 hidden layer only
if use elliott's function, define in the main code:
#define sig(x) x/(1.0 + fabs(x))
#define dsig(x) (1.0 - fabs(x))* (1.0 - fabs(x))

if use logistic function, define in the main code:
#define sig(x) 1/(1.0 + exp(-x))
#define dsig(x) (x - x*x)
*/

#include <math.h>
#include "nrutil.h"
#include <stdio.h>


#define dsig(x) (1.0 - fabs(x))* (1.0 - fabs(x))

void backprop1(x,x1,delta,dw0,dw1,w0,w1)
float *x,*x1,*delta,**dw0,**dw1,**w0,**w1;

/* backprop
w0: weight matrix from layer 0
w1: weight matrix from layer 1
x: input
x1: 1st hidden layer output
*/

{
  extern int NP1,P,PP1,R;
  int i,j;
  float *delta1,**wtemp1;
  extern float lrate0,alpha,lrate1;
  delta1 = vector(1,P);
  wtemp1 = matrix(1,R,1,P);

/*delta2 is (dJ/ds2) for the 2nd hidden layer
  wtemp2 is used to relate delta(dJ/ds3) of the output nodes to delta of
  the hidden nodes but not the bias input, therefore, wtemp2 is one column
  less than w2.  Note that the delta_w term represented by -lrate*.....
  has the oppsite sign of lrate*(dJ/dw).  If we want to compute the
  gradient dJ/dw, and not delta_w using this subroutine, make sure that
  the sign is consistent  */

/* output layer */
  for (i = 1;i <= R;i++)
    for (j = 1;j <= PP1;j++)
      dw1[i][j] = alpha*dw1[i][j]-lrate1*delta[i]*x1[j];

  for (i = 1;i <= R;i++)
    for (j = 1;j <= P;j++)
      wtemp1[i][j] = w1[i][j];

  for (i = 1;i <= R;i++)
    for (j = 1;j <= PP1;j++)
      w1[i][j] += dw1[i][j];


/* note that we have only P (dJ/ds1)'s because the bias node (node PP1)
is not connected back to layer 0*/

  vecmat(delta,wtemp1,delta1,R,P);

  for (i = 1;i <= P;i++)
     delta1[i] *= dsig(x1[i]);

  for (i = 1;i <= P;i++)
    for (j = 1;j <= NP1;j++)
      {
	dw0[i][j] = alpha*dw0[i][j]-lrate0*delta1[i]*x[j];
	w0[i][j] += dw0[i][j];
      }

  free_vector(delta1,1,P);
  free_matrix(wtemp1,1,R,1,P);
}