// (C) 2002 Jasper Bedaux, ITFA, University of Amsterdam
// http://www.bedaux.net/jasper/
// C++ file for layered neural network with recurrent connections

#include "neural_net.h"
#include "mtrand.h"
using namespace std;

// connect source to destination with specified chance
void connect(Neuron* src_begin, Neuron* src_end,
    Neuron* dest_begin, Neuron* dest_end, float chance /* = 1. */) {
  MTRand drand;
  for (Neuron* n = src_begin; n < src_end; ++n)
    for (Neuron* m = dest_begin; m < dest_end; ++m)
      if ((m != n) && (drand() < chance)) // no self connections
        n->connect(m);
}

// chances[0] to chances[hiddenlayers] are downward connection chances
// chances[hiddenlayers + 1] is lateral connection chance
// chances[hiddenlayers + 2] to chances[2 * hiddenlayers] are upward connection chances
Network::Network(int ninput, int noutput, int nhidden, int hiddenlayers,
  const vector<float>& chances) throw (string) :
    neuron(new Neuron[ninput + noutput + nhidden * hiddenlayers]),
    n(ninput + noutput + nhidden * hiddenlayers),
    ni(ninput), nh(nhidden), no(noutput), nl(hiddenlayers),
    input_end(neuron + ni), hidden_end(neuron + n - no), output_end(neuron + n) {
// check arguments
  if ((ni < 0) || (nh < 0) || (no < 0) || (nl < 0)) throw string("negative layer size");
  if (static_cast<int>(chances.size()) != (2 * nl + 1)) throw string("wrong number of connection chances");
  for (int i = 0; i < static_cast<int>(chances.size()); ++i) if ((chances[i] < 0.) || (chances[i] > 1.))
    throw string("connection chances must be in the interval [0, 1]");
// connect input layer to hidden layers and output layer
  for (int i = 0; i < size(layers); ++i)
    connect(begin(input), end(input), begin(i), end(i), chances[i]);
  connect(begin(input), end(input), begin(output), end(output), chances[size(layers)]);
// connect hidden layers to self, other hidden layers and output
  for (int i = 0; i < size(layers); ++i) {
    connect(begin(i), end(i), begin(i), end(i), chances[size(layers) + 1]); // lateral
    for (int j = i + 1; j < size(layers); ++j) // down
      connect(begin(i), end(i), begin(j), end(j), chances[j - (i + 1)]);
    for (int j = i - 1; j >= 0; --j) // up (feed-back)
      connect(begin(i), end(i), begin(j), end(j), chances[size(layers) + 2 + (i - 1) - j]);
    connect(begin(i), end(i), begin(output), end(output), chances[size(layers) - (i + 1)]);
  }
}

// constructor for simple feed-forward net with one hidden layer
Network::Network(int ninput, int nhidden, int noutput, float ph, float po)
    throw (string) : neuron(new Neuron[ninput + nhidden + noutput]),
    n(ninput + nhidden + noutput), ni(ninput), nh(nhidden), no(noutput), nl(1),
    input_end(neuron + ni), hidden_end(neuron + n - no), output_end(neuron + n) {
// check arguments
  if ((ni < 0) || (nh < 0) || (no < 0)) throw string("negative layer size");
  if ((ph < 0.) || (ph > 1.) || (po < 0.) || (po > 1.))
    throw string("connection chances must be in the interval [0, 1]");
// connect layers
  connect(begin(input), end(input), begin(hidden), end(hidden), ph);
  connect(begin(hidden), end(hidden), begin(output), end(output), po);
}

void Network::setInput(const vector<bool>& in) {
  int length = in.size();
  if (length > size(input)) length = size(input);
  for (int i = 0; i < length; ++i)
    neuron[i].setActivity(in[i]);
// set other input to inactive if length < size(input)
  for (int i = length; i < size(input); ++i)
    neuron[i].setActivity(false);
}

const vector<bool>& Network::getOutput() const {
  static vector<bool> out(size(output));
  for (int i = 0; i < size(output); ++i)
    out[i] = neuron[offset(output) + i].isActive();
  return out;
}

// let the active neurons in the specified region fire: update potentials
void fire(Neuron* start, Neuron* end) {
  for (Neuron* n = start; n < end; ++n)
    if (n->isActive()) for (int i = 0; i < n->size(); ++i)
       (*n)[i].post()->raisePotential((*n)[i].weight());
}

// update states of neurons in region, if na != 0 use extremal dynamics
void update(Neuron* start, Neuron* end, int na /* = 0 */) {
  if (!na) for (Neuron* n = start; n < end; ++n) // default dynamics
    n->setActivity(n->potential() > n->threshold());
  else { // extremal dynamics
    for (Neuron* n = start; n < end; ++n) n->setActivity(false);
    for (int i = 0; i < na; ++i) {
      Neuron* hmax = start;
      while (hmax->isActive()) ++hmax; // find non-firing neuron
      for (Neuron* n = hmax + 1; n < end; ++n)
        if ((n->potential() > hmax->potential()) && !(n->isActive()))
          hmax = n;
      hmax->setActivity(true);
    }
  }
}

// simulate one timestep for feed-forward networks
// if nah and nao are non-zero, extremal dynamics is used
void Network::findFixed(const vector<bool>& in, int nah, int nao) {
  setInput(in);
  for (Neuron* n = begin(hidden); n < end(); ++n) n->resetPotential();
  fire(begin(input), end(input));
  for (int i = 0; i < size(layers); ++i) {
    update(begin(i), end(i), nah);
    fire(begin(i), end(i));
  }
  update(begin(output), end(output), nao);
}

// find fixedpoint for recurrent networks
// returns number of steps needed or 0 if not found in tmax
// if nah and nao are non-zero, extremal dynamics is used
int Network::findFixed(int tmax, const vector<bool>& in, int nah, int nao) {
  static vector<bool> state(size(hidden) * size(layers));
  static vector<bool> newstate(size(hidden) * size(layers));
  setInput(in);
// reset potentials in hidden layers
  for (Neuron* n = begin(hidden); n < end(hidden); ++n) n->resetPotential();
// effect of input layer
  fire(begin(input), end(input));
  for (int i = 0; i < size(layers); ++i)
    update(begin(i), end(i), nah); // update hidden layers
// determine state of hidden layers
  for (int j = 0; j < static_cast<int>(state.size()); ++j)
    state[j] = neuron[j + size(input)].isActive();
  for (int i = 0; i < tmax; ++i) {
// reset potentials in hidden and output layers
    for (Neuron* n = begin(hidden); n < end(); ++n) n->resetPotential();
    fire(begin(input), end(input)); // fire input and hidden layers
    fire(begin(hidden), end(hidden));
    for (int j = 0; j < size(layers); ++j)
      update(begin(j), end(j), nah); // update hidden layers
    for (int j = 0; j < static_cast<int>(newstate.size()); ++j)
      newstate[j] = neuron[j + size(input)].isActive();
    if (newstate == state) { // fixedpoint found
      update(begin(output), end(output), nao);
      return i + 1;
    }
    state = newstate;
  }
  return 0; // no fixed point found in tmax steps
}

// Hebbian: eta(kappa(2x_i - 1) - (h_i - theta_i))(2x_i - 1)x_j
void Hebb(Neuron* start, Neuron* end, float eta, float kappa, float delta) {
  MTRand_normal gaussian;
  for (Neuron* n = start; n < end; ++n) {
    if (n->isActive()) { // x_j = 1
      for (int i = 0; i < n->size(); ++i) {
        Neuron* post = (*n)[i].post();
        if (post->isActive()) // (2x_i - 1) = 1
          (*n)[i].grow(eta * ( kappa - (post->potential() - post->threshold()))
            * (1. + delta * gaussian()));
        else // (2x_i - 1) = -1
          (*n)[i].grow(eta * (-kappa - (post->potential() - post->threshold()))
            * (1. + delta * gaussian()));
      }
    }
  }
}

// anti-Hebbian rho(x_i - a)x_j +- random part
void antiHebb(Neuron* start, Neuron* end, float rho, float alpha, float delta) {
  MTRand_normal gaussian;
  for (Neuron* n = start; n < end; ++n)
    if (n->isActive()) // x_j = 1
      for (int i = 0; i < n->size(); ++i)
        if ((*n)[i].post()->isActive())
          (*n)[i].grow(rho * (alpha - 1.) * (1. + delta * gaussian()));
        else (*n)[i].grow(rho * alpha); // without noise for speed
}