/*
 *  gamma_sup_generator.cpp
 *
 *  This file is part of NEST.
 *
 *  Copyright (C) 2004 The NEST Initiative
 *
 *  NEST is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 2 of the License, or
 *  (at your option) any later version.
 *
 *  NEST is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with NEST.  If not, see <http://www.gnu.org/licenses/>.
 *
 */

#include "gamma_sup_generator.h"
#include "network.h"
#include "dict.h"
#include "integerdatum.h"
#include "doubledatum.h"
#include "numerics.h"
#include "datum.h"
#include <algorithm>
#include <limits>


/* ---------------------------------------------------------------- 
 * Constructor of internal states class
 * ---------------------------------------------------------------- */
    
nest::gamma_sup_generator::Internal_states_::Internal_states_(size_t num_bins, ulong_t ini_occ_ref, ulong_t ini_occ_act)
{
    occ_.resize(num_bins, ini_occ_ref);
    occ_.back() += ini_occ_act;
}

/* ---------------------------------------------------------------- 
 * Propagate internal states one time step and generate spikes
 * ---------------------------------------------------------------- */

nest::ulong_t nest::gamma_sup_generator::Internal_states_::update(double_t transition_prob, librandom::RngPtr rng)
{
    std::vector<ulong_t> n_trans;
    n_trans.resize( occ_.size() );
    
    // go through all states and draw number of transitioning components
    for (ulong_t i=0; i<occ_.size(); i++)
        {
        if (occ_[i]>0) 
            {
            /*The binomial distribution converges towards the Poisson distribution as
             the number of trials goes to infinity while the product np remains fixed.
             Therefore the Poisson distribution with parameter λ = np can be used as 
             an approximation to B(n, p) of the binomial distribution if n is 
             sufficiently large and p is sufficiently small. According to two rules 
             of thumb, this approximation is good if n ≥ 20 and p ≤ 0.05, or if 
             n ≥ 100 and np ≤ 10. Source:
             http://en.wikipedia.org/wiki/Binomial_distribution#Poisson_approximation */
            if (( occ_[i] >= 100 && transition_prob <= 0.01 ) || \
                (  occ_[i] >= 500 && transition_prob * occ_[i] <= 0.1 ))
                {
                poisson_dev_.set_lambda( transition_prob * occ_[i] );
                n_trans[i] = poisson_dev_.uldev(rng);
                if ( n_trans[i] > occ_[i] )
                    {
                    n_trans[i] = occ_[i];
                    }
                else
                    {;}
                }
            else
                {
                bino_dev_.set_p_n( transition_prob, occ_[i]);
                n_trans[i] = bino_dev_.uldev(rng); 
                }
            }
        else
            {
            n_trans[i] = 0;
            }
        }
    
    // according to above numbers, change the occupation vector
    for (ulong_t i=0; i<occ_.size(); i++)
        {
        if (n_trans[i]>0) 
            {
            occ_[i] -= n_trans[i];
            if (i==occ_.size()-1)
                occ_.front() += n_trans[i];
            else
                occ_[i+1] += n_trans[i]; 
            }
        }
    return n_trans.back();
}


/* ---------------------------------------------------------------- 
 * Default constructors defining default parameter
 * ---------------------------------------------------------------- */
    
nest::gamma_sup_generator::Parameters_::Parameters_()
  : rate_(0.0),  // Hz
    gamma_shape_(1),
    n_proc_(1),
    num_targets_(0)
{}

/* ---------------------------------------------------------------- 
 * Parameter extraction and manipulation functions
 * ---------------------------------------------------------------- */

void nest::gamma_sup_generator::Parameters_::get(DictionaryDatum &d) const
{
  (*d)[names::rate] = rate_;
  (*d)[names::gamma_shape] = gamma_shape_;
  (*d)[names::n_proc] = n_proc_;
}  

void nest::gamma_sup_generator::Parameters_::set(const DictionaryDatum& d)
{
  updateValue<long_t>(d, names::gamma_shape, gamma_shape_);
  if ( gamma_shape_ < 1 )
    throw BadProperty("The shape must be larger or equal 1");
    
  updateValue<double_t>(d, names::rate, rate_);
  if ( rate_ < 0.0 )
    throw BadProperty("The rate must be larger than 0.");
  
  long n_proc_l = n_proc_;
  updateValue<long_t>(d, names::n_proc, n_proc_l);
  if ( n_proc_l < 1 )
    throw BadProperty("The number of component processes cannot be smaller than one");
  else
    n_proc_ = static_cast<ulong_t> (n_proc_l);
}


/* ---------------------------------------------------------------- 
 * Default and copy constructor for node
 * ---------------------------------------------------------------- */

nest::gamma_sup_generator::gamma_sup_generator()
  : Node(),
    device_(), 
    P_()
{}

nest::gamma_sup_generator::gamma_sup_generator(const gamma_sup_generator& n)
  : Node(n), 
    device_(n.device_),
    P_(n.P_)
{}


/* ---------------------------------------------------------------- 
 * Node initialization functions
 * ---------------------------------------------------------------- */

void nest::gamma_sup_generator::init_state_(const Node& proto)
{ 
  const gamma_sup_generator& pr = downcast<gamma_sup_generator>(proto);

  device_.init_state(pr.device_);
}

void nest::gamma_sup_generator::init_buffers_()
{ 
  device_.init_buffers();
}

void nest::gamma_sup_generator::calibrate()
{
  device_.calibrate();
  
  double_t h = Time::get_resolution().get_ms();
 
  // transition probability in each time step
  V_.transition_prob_ = P_.rate_ * P_.gamma_shape_ * h / 1000.0;
    
  // approximate equilibrium occupation to initialize to
  ulong_t ini_occ_0 = static_cast<ulong_t> (P_.n_proc_ / P_.gamma_shape_)  ;
  
  // If new targets have been added during a simulation break, we
  // initialize the new elements in Internal_states with the initial dist. The existing
  // elements are unchanged.
  Internal_states_ internal_states0 (P_.gamma_shape_, ini_occ_0, P_.n_proc_ - ini_occ_0 * P_.gamma_shape_);
  B_.internal_states_.resize( P_.num_targets_, internal_states0 );
}



/* ---------------------------------------------------------------- 
 * Update function and event hook
 * ---------------------------------------------------------------- */

void nest::gamma_sup_generator::update(Time const & T, const long_t from, const long_t to)
{
  assert(to >= 0 && (delay) from < Scheduler::get_min_delay());
  assert(from < to);

  if ( P_.rate_ <= 0 || P_.num_targets_ == 0 ) 
    return;

  for ( long_t lag = from ; lag < to ; ++lag )
  {
    Time t = T + Time::step(lag); 
    
    if ( !device_.is_active( t ) )
      continue;  // no spike at this lag
    
    DSSpikeEvent se;
    network()->send(*this, se, lag);
  }
}


void nest::gamma_sup_generator::event_hook(DSSpikeEvent& e)
{

  // get port number
  const port prt = e.get_port();

  // we handle only one port here, get reference to vector elem
  assert(0 <= prt && static_cast<size_t>(prt) < B_.internal_states_.size() );

  // age_distribution object propagates one time step and returns number of spikes
  ulong_t n_spikes = B_.internal_states_[prt].update( V_.transition_prob_, net_->get_rng(get_thread()) );
  
  if ( n_spikes > 0 ) // we must not send events with multiplicity 0
  {
    e.set_multiplicity(n_spikes);
    e.get_receiver().handle(e);
  }
}