#include <math.h>

#include "synapse.h"

Synapse::Synapse() :
  src(-1), dst(-1), // seeing this value indicates an unused synapse
  weight(0.0),
  eligibility(0.0),
  delay(0.0),
  constant(false),
  firstSpike(- INFINITY),
  lastSpike(- INFINITY),
  evolvedUntil(0.0)
{}

// given an input time and current calculate their output values
void Synapse::computePostsynapticPulse(double &time, double &current) {
  time += delay;
  current = weight;
}

void Synapse::evolve(double time) {
  if (time < evolvedUntil)
    return;

  // get the dopamin history starting from a dopamin event earlier than the synpase's last event
  std::list<pair<double, double> >::iterator i = s.da_history.begin();
  while (i->first > evolvedUntil)
    i++;

  // evolve synapse state for each interval given by the dopamin history
  while (evolvedUntil < time) {
    double da_level = i->second;
    if (i->first < evolvedUntil) // the first dopamin event partially lies outside our time window
      da_level = g.decay_dopamin(da_level, evolvedUntil - i->first);

    double td = - evolvedUntil;
    if (i == s.da_history.begin()) {
      // this is the last (=newest) element of the history
      td += time;
      evolvedUntil = time;
    }else{
      i--;
      if (i->first > time) {
	// dopamin history goes furher (in time) than we have to simulate
	td += time;
	evolvedUntil = time;
      }else{
	td += i->first;
	evolvedUntil += td;
      }
    }

    /* below is the exact solution of the equation system
         dc/dt = -Cc          c - eligibility trace, C - cf. tau
         dd/dt = -Dd          d - dopamin level, D - cf. tau
         dw/dt = cd           w - weight
    */
    
    if (!constant) {
      double tau = g.dopamin_tau * g.stdp_et_tau / (g.dopamin_tau + g.stdp_et_tau);
      double dw = da_level * eligibility * tau * (1 - exp( -td / tau ));
      weight += dw;
      weight = fmax( msg.Wmin, fmin( msg.Wmax, weight ));
    }
    eligibility *= exp( -td / g.stdp_et_tau );
  }
}

// modify eligibility trace according to a single STDP event given by
// a pre-post time difference (= t_post - t_pre)
void Synapse::stdp(double td) {
  if (td >= 0.0) {// WARN: this doesn't work with negative zero
    eligibility = g.stdp_lambda_plus * exp( - td / g.stdp_tau_plus );
  }else{
    eligibility = - g.stdp_lambda_minus * exp( td / g.stdp_tau_minus );
  }
}

static void Synapse::updateTopology() {
  // clear the sin/sout lists of all neurons
  for (int i=0; i<numNeurons; i++) {
    n[i].sout.clear();
    n[i].sin.clear();
  }

  // iterate over all synapses wich have a valid source and destination (neuron) and add them to sin/sout lists
  int i;
  for (i=0; (i<maxSynapses) && (syn[i].src != -1) && (syn[i].dst != -1); i++) {
    n[syn[i].src].sout.push_front(syn[i]);
    n[syn[i].dst].sin.push_front(syn[i]);
  }

  // store number of (used) synapses
  numSynapses = i;
}

template<class Action>
bool Synapse::reflect<Action>(Action &a) {
  return
    a(src, "src")
    && a(dst, "dst")
    && a(weight, "weight")
    && a(eligibility, "eligibility")
    && a(delay, "delay")
    && a(constant, "constant")
    && a(firstSpike, "firstSpike")
    && a(lastSpike, "lastSpike")
    && a(evolvedUntil, "evolvedUntil");
}

static id_type Synapse::numElements() {
  return numSynapsesDec;
}

static Synapse * Synapse::singleton(id_type id) {
  return &(synapses[id]);
}