/****************************************************************************************/
/*                                                                                      */
/*  Project     :  1D-BRNN                                                              */
/*  Release     :  3.3                                                                  */
/*                                                                                      */
/*  File        :  Layer.cpp                                                            */
/*  Description :  Single layer of a neural network                                     */
/*                                                                                      */
/*  Author(s)   :  Christophe Magnan  (2013) - New generic version                      */
/*              :  Jianlin Cheng      (2003) - New custom version for SCRATCH           */
/*              :  Gianluca Pollastri (2001) - Customized version for SCRATCH           */
/*              :  Paolo Frasconi     (1997) - Initial generic version                  */
/*                                                                                      */
/*  Copyright   :  Institute for Genomics and Bioinformatics                            */
/*                 University of California, Irvine                                     */
/*                                                                                      */
/*  Modified    :  2015/07/01                                                           */
/*                                                                                      */
/****************************************************************************************/

#include "Layer.h"

/*******************************************/
/*                Interface                */
/*******************************************/

// Constructor
Layer::Layer() :
  Num_Nodes(0),Num_Inputs(0),IsOutputLay(0),
  Bias(NULL),Weights(NULL),Direct_Outputs(NULL),
  SigApp_Outputs(NULL),Gradient(NULL),
  Back_Prop(NULL),dBias(NULL),dWeights(NULL){}

// Destructor
Layer::~Layer()
{
  for(int n=0;n<Num_Nodes;n++){ delete[] Weights[n]; }
  delete[] Bias; delete[] Weights;
  delete[] Direct_Outputs; delete[] SigApp_Outputs;
  if(Gradient!=NULL){ for(int n=0;n<Num_Nodes;n++){
  delete[] dWeights[n]; } delete[] Gradient;
  delete[] Back_Prop; delete[] dBias; delete[] dWeights; }
}

// Write layer
void Layer::write(ostream& os)
{
  os << Num_Nodes << " " << Num_Inputs << " " << IsOutputLay << "\n";
  for(int n=0;n<Num_Nodes;n++){ for(int i=0;i<Num_Inputs;i++){
  os << Weights[n][i] << " "; } os << Bias[n] << "\n"; }
}

// Load layer
void Layer::load(istream& is)
{
  is >> Num_Nodes >> Num_Inputs >> IsOutputLay;
  Bias=new float[Num_Nodes]; Weights=new float*[Num_Nodes];
  for(int n=0;n<Num_Nodes;n++){ Weights[n]=new float[Num_Inputs];
  for(int i=0;i<Num_Inputs;i++){ is >> Weights[n][i]; } is >> Bias[n]; }
  check_viability(); Direct_Outputs=new float[Num_Nodes];
  SigApp_Outputs=new float[Num_Nodes];
}

// New random layer
void Layer::initialize(int nodes,int inputs)
{
  Num_Nodes=nodes; Num_Inputs=inputs; IsOutputLay=0; drand48();
  Bias=new float[Num_Nodes]; Weights=new float*[Num_Nodes];
  for(int n=0;n<Num_Nodes;n++)
  { Bias[n]=(float)(0.5-drand48())/10.0;
    Weights[n]=new float[Num_Inputs];
    for(int i=0;i<Num_Inputs;i++){
    Weights[n][i]=(float)(0.5-drand48())/10.0; } }
  Direct_Outputs=new float[Num_Nodes];
  SigApp_Outputs=new float[Num_Nodes];
}

// Allocate memory for back-propagation
void Layer::alloc_backpropagation()
{
  Gradient=new float[Num_Nodes]; Back_Prop=new float[Num_Inputs];
  dBias=new float[Num_Nodes]; dWeights=new float*[Num_Nodes];
  for(int n=0;n<Num_Nodes;n++){ dWeights[n]=new float[Num_Inputs]; }
  reset_gradient();
}

// Reset the gradient matrices
void Layer::reset_gradient()
{
  for(int n=0;n<Num_Nodes;n++){ 
  for(int i=0;i<Num_Inputs;i++){
  dWeights[n][i]=0.0; } dBias[n]=0.0; }
}

// Propagate the input features into the layer
void Layer::propagate(float* inputs)
{
  for(int n=0;n<Num_Nodes;n++){ Direct_Outputs[n]=Bias[n];
  for(int i=0;i<Num_Inputs;i++){ Direct_Outputs[n]+=Weights[n][i]*inputs[i]; }
  SigApp_Outputs[n]=(float)tanh(Direct_Outputs[n]); } if(IsOutputLay){
  float norm=0,max=Direct_Outputs[0]; int overflow=0; int argmax=0;
  for(int n=0;n<Num_Nodes;n++){ float v=Direct_Outputs[n];
  if(v>max){ max=v; argmax=n; } if(v>85){ overflow=1; }
  else{ norm+=(float)exp(v); } } for(int n=0;n<Num_Nodes;n++){
  if(overflow){ SigApp_Outputs[n]=(n==argmax?1.0:0.0); }
  else{ SigApp_Outputs[n]=(float)exp(Direct_Outputs[n])/norm; } } }
}

// Back-propagate the target values into the layer
void Layer::back_propagate(float* target)
{
  if(IsOutputLay)
  { int overflow=0; for(int n=0;n<Num_Nodes;n++){
    if(Direct_Outputs[n]>85){ overflow=1; } }
    for(int n=0;n<Num_Nodes;n++){
    if(overflow){ Gradient[n]=0; }
    else{ float app=SigApp_Outputs[n];
    Gradient[n]=(app-target[n])*app*(1-app); } } }
  else
  { for(int n=0;n<Num_Nodes;n++){
    float app=SigApp_Outputs[n];
    Gradient[n]=target[n]*(1.0-app*app); } }
  for(int i=0;i<Num_Inputs;i++)
  { Back_Prop[i]=0.0; for(int n=0;n<Num_Nodes;n++){
    Back_Prop[i]+=Weights[n][i]*Gradient[n]; } }
}

// Update the gradient matrices after back-propagation
void Layer::update_gradient(float* inputs)
{
  for(int n=0;n<Num_Nodes;n++){
  for(int i=0;i<Num_Inputs;i++){
  dWeights[n][i]+=Gradient[n]*inputs[i]; }
  dBias[n]+=Gradient[n]; }
}

// Update the weights after training
void Layer::update_weights(float epsilon)
{
  for(int n=0;n<Num_Nodes;n++){
  for(int i=0;i<Num_Inputs;i++){
  Weights[n][i]-=epsilon*dWeights[n][i]; }
  Bias[n]-=epsilon*dBias[n]; }
}

// Check the viability of a loaded layer
void Layer::check_viability()
{
  int v=1; if((Num_Nodes<1)||(Num_Inputs<1)){ v=0; }
  if((IsOutputLay<0)||(IsOutputLay>1)){ v=0; }
  if(!v){ cerr << "Inconsistent model, cannot proceed\n"; exit(1); }
}