/* *********************************************************** */
/* file real.c:  contains the network evaluation and weight    */
/* adjustment procedures for the floating point versions       */
/*                                                             */
/* Copyright (c) 1990-96 by Donald R. Tveter                   */
/*                                                             */
/* *********************************************************** */

#include "rbp.h"

extern char deriv, outstr[], zeroderiv;
extern REAL alpha, decay, eta, eta2, etamax, kappa, qpdecayh, qpdecayo;
extern REAL qpeta, mu;
extern REAL theta1, theta2, toler, totaldiff, poft;
extern char qpslope, outformat, update, debugoutput;
extern LAYER *last, *start;
extern short nlayers;

extern int pg();

void forward()       /* computes unit activations */
{
register REAL sum;
register WTNODE *w;
register UNIT *predu;
register REAL denom;
register REAL fract, x, D;
register char actfun;
#ifdef UNIX
REAL val; /* should be in a register, but UNIX PC C-compiler does */
          /* not handle it correctly */
#else
register REAL val;
#endif
register int intpart;
UNIT *u;
LAYER *layer;

layer = start->next;
while (layer != NULL)
 {
  actfun = layer->activation;
  D = layer->D;
  u = (UNIT *) layer->units;
  while (u != NULL)
   {
    if (!u->inuse) u->oj = 0.0;
    else
     {
      sum = 0.0;
      w = (WTNODE *) u->wtlist;
      while (w)
       {
        predu = (UNIT *) w->backunit;
#ifdef SYMMETRIC
        sum = sum + *(w->weight) * predu->oj;
#else
        sum = sum + w->weight * predu->oj;
#endif
        w = w->next;
       };
      sum = sum * D;
      if (actfun == 's') u->oj = 1.0 / (1.0 + exp((double) -sum));
      else if (actfun == 'l') u->oj = sum;
      else pg(stdout,"invalid activation function\n");
     } /* end of else part of: if (u->inuse == 0) */
    u = u->next;
   };
  layer = layer->next;
 };
}

short backoutput()  /* back propagate errors from the output units */
{                   /* send down errors for any previous layers    */
register REAL deltaj;
register WTNODE *w;
register UNIT *bunit;
register REAL diff, adiff;
register PATNODE *target;
register UNIT *u;
register short notclose;
char actfun;
PATLIST *pl;

notclose = last->unitcount;
if (update != 't' && update != 'T')
 {
  pl = last->currentpat[TRAIN];
  target = pl->pats;
 };
actfun = last->activation;
u = (UNIT *) last->units;
while (u != NULL)
 {
  if (update == 't' || update == 'T') diff = 1.0;
  else diff = target->val - u->oj;
  target++;
  if (diff > 0) adiff = diff; else adiff = -diff;
  if (adiff < toler) notclose = notclose - 1;
  totaldiff = totaldiff + adiff;
  if (deriv == 'd')
   {
    deltaj = diff;
   }
  else 
   {
    switch (actfun) {/* deltaj in the switch is ONLY the derivative */
    case 'l': deltaj = 1.0;
              break;
    case 's': deltaj = u->oj * (1.0 - u->oj);
              break;}
    if (deriv == 'f')
     {
      deltaj = diff * (deltaj + 0.1);
     }
    else deltaj = diff * deltaj;
   };
  w = (WTNODE *) u->wtlist;
#ifdef SYMMETRIC
   while (w->next != NULL)
#else
   while (w) /* while (w != NULL) */
#endif
    {
     bunit = (UNIT *) w->backunit;
#ifdef SYMMETRIC
     *(w->total) = *(w->total) + deltaj * bunit->oj;
#else
     w->total = w->total + deltaj * bunit->oj;
     if (bunit->layernumber > 1)  /* pass back the error */
        bunit->error = bunit->error + deltaj * w->weight;
#endif
     w = w->next;
    };
  u = u->next;
 }
return(notclose);
}

#ifndef SYMMETRIC

void backinner()  /* compute weight changes for hidden layers */
{                 /* send down errors for any previous layers */
register REAL deltaj;
register UNIT *bunit;
register WTNODE *w;
register UNIT *u;
char actfun;
LAYER *layer;

layer = last->backlayer;
while (layer->backlayer != NULL)
 {
  actfun = layer->activation;
  u = (UNIT *) layer->units;
  while (u != NULL)
   {
    switch (actfun) {
    case 'l': deltaj = 1.0;
              break;
    case 's': deltaj = (u->oj * (1.0 - u->oj));
              break;}
    deltaj = deltaj * u->error;
    w = (WTNODE *) u->wtlist;
     while (w) /* while (w != NULL) */
      {
       bunit = (UNIT *) w->backunit;
       w->total = w->total + deltaj * bunit->oj;
       if (bunit->layernumber > 1)
          bunit->error = bunit->error + deltaj * w->weight;
       w = w->next;
      }
    u = u->next;
   };
  layer = layer->backlayer;
 };
}

#endif

void periodic_update()  /* the original periodic method */
{
register REAL reta, ralpha;
register UNIT *u;
register WTNODE *w;
LAYER *layer;

ralpha = alpha;
layer = last;
while (layer->backlayer != NULL)
 {
  if (layer == last) reta = eta; else reta = eta2;
  u = (UNIT *) layer->units;
  while (u != NULL)
   {
    w = (WTNODE *) u->wtlist;
    while (w)
     {
      if (w->inuse > 0)
       {
#ifdef SYMMETRIC
        if (((UNIT *) w->backunit)->unitnumber > u->unitnumber)
         {
          *(w->olddw) = *(w->total) * reta + ralpha * *(w->olddw);
          *(w->weight) = *(w->weight) + *(w->olddw);
         };
#else
        w->olddw = w->total * reta + ralpha * w->olddw;
        w->weight = w->weight + w->olddw;
#endif
       };
      w = w->next;
     };
    u = u->next;
   };
  layer = layer->backlayer;
 };
}

#ifndef SYMMETRIC
void dbd_update() /* delta-bar-delta method for changing weights */
{
register short stotal,sdbarm1;
register UNIT *u;
register WTNODE *w;
LAYER *layer;

/* w->olddw is used for delta-bar minus 1 */

layer = last;
while (layer->backlayer != NULL)
 {
  u = (UNIT *) layer->units;
  while (u != NULL)
   {
    w = (WTNODE *) u->wtlist;
    while (w)
     {
      if (w->inuse > 0)
       {
        if (w->total > 0) stotal = 1;
          else if (w->total < 0) stotal = -1;
           else stotal = 0;
        if (w->olddw > 0) sdbarm1 = 1;
          else if (w->olddw < 0) sdbarm1 = -1;
           else sdbarm1 = 0;
        w->olddw = theta2 * w->total + theta1 * w->olddw;
        if ((stotal > 0) && (sdbarm1 > 0)) w->eta = w->eta + kappa;
        else if ((stotal < 0) && (sdbarm1 < 0)) w->eta = w->eta + kappa;
        else if ((stotal > 0) && (sdbarm1 < 0)) w->eta = w->eta * decay;
        else if ((stotal < 0) && (sdbarm1 > 0)) w->eta = w->eta * decay;
        if (w->eta > etamax) w->eta = etamax;
        w->weight = w->weight + w->eta * w->total;
       };
      w = w->next;
     };
    u = u->next;
   };
  layer = layer->backlayer;
 };
}

void qp_update()
{
register REAL reta, s, nextdw, shrinkfactor, rqpdecay;
register WTNODE *w;
register short addslope;
register REAL rmu;
register UNIT *u;
LAYER *layer;

rmu = mu;
shrinkfactor = rmu / (1.0 + rmu);
reta = qpeta;
if (qpslope == '+') addslope = 1; else addslope = 0;
layer = last;
while (layer->backlayer != NULL)
 {
  if (layer == last) rqpdecay = qpdecayo;
  else rqpdecay = qpdecayh;
  u = (UNIT *) layer->units;
  while (u != NULL)
   {
    w = (WTNODE *) u->wtlist;
    while (w)
     {
      if (w->inuse > 0) /* if 1 or 2 but not 0 or -1 */
       {
        s = rqpdecay * w->weight - w->total;
        if (w->olddw < 0.0)
         {
          if (s >= (shrinkfactor * w->eta)) nextdw = rmu * w->olddw;
          else nextdw = w->olddw * s / (w->eta - s);
          if (addslope && s > 0.0) nextdw = nextdw - reta * s;
         }
        else if (w->olddw > 0.0)
         {
          if (s <= (shrinkfactor * w->eta)) nextdw = rmu * w->olddw;
          else nextdw = w->olddw * s / (w->eta - s);
          if (addslope && s < 0.0) nextdw = nextdw - reta * s;
         }
        else nextdw = - reta * s;
        w->olddw = nextdw;
        w->weight = w->weight + nextdw;
        w->eta = s;
       };
      w = w->next;
     };
    u = u->next;
   };
  layer = layer->backlayer;
 };
}

#else
void dbd_update() {}
void qp_update() {}
void supersab() {}
#endif

short cbackoutput()  /* for the "wrong" continuous update (c) */
{
register REAL deltaj, etadeltaj, diff, adiff, reta, ralpha;
register UNIT *u, *bunit;
register WTNODE *w;
register short notclose;
register PATNODE *target;
char actfun;
PATLIST *pl;

actfun = last->activation;
reta = eta;
ralpha = alpha;