iaf_psc_alpha_canon.h

Go to the documentation of this file.

/*
 *  iaf_psc_alpha_canon.h
 *
 *  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/>.
 *
 */

#ifndef IAF_PSC_ALPHA_CANON_H
#define IAF_PSC_ALPHA_CANON_H

#include "config.h"

#include "nest.h"
#include "event.h"
#include "node.h"
#include "ring_buffer.h"
#include "slice_ring_buffer.h"
#include "connection.h"

#include "universal_data_logger.h"

#include <vector>

/*BeginDocumentation
Name: iaf_psc_alpha_canon - Leaky integrate-and-fire neuron
with alpha-shape postsynaptic currents; canoncial implementation.

Description:
iaf_psc_alpha_canon is the "canonical" implementatoin of the leaky
integrate-and-fire model neuron with alpha-shaped postsynaptic
currents in the sense of [1].  This is the most exact implementation
available.

PSCs are normalized to an amplitude of 1pA.

The canonical implementation handles neuronal dynamics in a locally
event-based manner with in coarse time grid defined by the minimum
delay in the network, see [1]. Incoming spikes are applied at the
precise moment of their arrival, while the precise time of outgoing
spikes is determined by interpolation once a threshold crossing has
been detected. Return from refractoriness occurs precisly at spike
time plus refractory period.

This implementation is more complex than the plain iaf_psc_alpha
neuron, but achieves much higher precision. In particular, it does not
suffer any binning of spike times to grid points. Depending on your
application, the canonical application may provide superior overall
performance given an accuracy goal; see [1] for details.  Subthreshold
dynamics are integrated using exact integration between events [2].

Remarks:
The iaf_psc_delta_canon neuron does not accept CurrentEvent connections.
This is because the present method for transmitting CurrentEvents in 
NEST (sending the current to be applied) is not compatible with off-grid
currents, if more than one CurrentEvent-connection exists. Once CurrentEvents
are changed to transmit change-of-current-strength, this problem will 
disappear and the canonical neuron will also be able to handle CurrentEvents.
For now, the only way to inject a current is the built-in current I_e.

Please note that this node is capable of sending precise spike times
to target nodes (on-grid spike time plus offset). If this node is
connected to a spike_detector, the property "precise_times" of the
spike_detector has to be set to true in order to record the offsets
in addition to the on-grid spike times.

A further improvement of precise simulation is implemented in iaf_psc_exp_ps
based on [3].

Parameters:
The following parameters can be set in the status dictionary.

  V_m          double - Membrane potential in mV 
  E_L          double - Resting membrane potential in mV. 
  V_min        double - Absolute lower value for the membrane potential.
  C_m          double - Capacity of the membrane in pF
  tau_m        double - Membrane time constant in ms.
  t_ref        double - Duration of refractory period in ms. 
  V_th         double - Spike threshold in mV.
  V_reset      double - Reset potential of the membrane in mV.
  tau_syn      double - Rise time of the synaptic alpha function in ms.
  I_e          double - Constant external input current in pA.
  Interpol_Order  int - Interpolation order for spike time: 
                        0-none, 1-linear, 2-quadratic, 3-cubic

Note:
  tau_m != tau_syn is required by the current implementation to avoid a
  degenerate case of the ODE describing the model [1]. For very similar values,
  numerics will be unstable.

References:
[1] Morrison A, Straube S, Plesser H E, & Diesmann M (2006) Exact Subthreshold 
    Integration with Continuous Spike Times in Discrete Time Neural Network 
    Simulations. To appear in Neural Computation.
[2] Rotter S & Diesmann M (1999) Exact simulation of time-invariant linear
    systems with applications to neuronal modeling. Biologial Cybernetics
    81:381-402.
[3] Hanuschkin A, Kunkel S, Helias M, Morrison A & Diesmann M (2010) 
    A general and efficient method for incorporating exact spike times in 
    globally time-driven simulations Front Neuroinformatics, 4:113
        
Author: Diesmann, Eppler, Morrison, Plesser, Straube

Sends: SpikeEvent

Receives: SpikeEvent, CurrentEvent, DataLoggingRequest

SeeAlso: iaf_psc_alpha, iaf_psc_alpha_presc, iaf_psc_exp_ps

*/ 

namespace nest{

  class iaf_psc_alpha_canon:
    public Node
  {

    class Network;
    
  public:        
    
    iaf_psc_alpha_canon();

    iaf_psc_alpha_canon(const iaf_psc_alpha_canon&);

   using Node::handle;
    using Node::handles_test_event;

    port send_test_event(Node &, rport, synindex, bool);

    port handles_test_event(SpikeEvent&, rport);
    port handles_test_event(CurrentEvent&, rport);
    port handles_test_event(DataLoggingRequest &, rport);
    
    void handle(SpikeEvent &);
    void handle(CurrentEvent&);
    void handle(DataLoggingRequest &);

    bool is_off_grid() const {return true;}  // uses off_grid events

    void get_status(DictionaryDatum &) const;
    void set_status(const DictionaryDatum &);

  private:

    void init_state_(const Node& proto);
    void init_buffers_();
    void calibrate();

    void update(Time const & origin, const long_t from, const long_t to);
    

    void set_spiketime(Time const &);

    void propagate_(const double_t dt);

    void emit_spike_(Time const& origin, const long_t lag, 
             const double_t t0, const double_t dt);

    void emit_instant_spike_(Time const& origin, const long_t lag, 
                 const double_t spike_offset);


    enum interpOrder { NO_INTERPOL, LINEAR, QUADRATIC, CUBIC, END_INTERP_ORDER };

    double_t thresh_find_(double_t const) const;
    double_t thresh_find1_(double_t const) const;
    double_t thresh_find2_(double_t const) const;
    double_t thresh_find3_(double_t const) const;


    // The next two classes need to be friends to access the State_ class/member
    friend class RecordablesMap<iaf_psc_alpha_canon>;
    friend class UniversalDataLogger<iaf_psc_alpha_canon>;

    // ---------------------------------------------------------------- 

    struct Parameters_ {

      double_t tau_m_; 

      double_t tau_syn_;

      double_t c_m_;
    
      double_t t_ref_;

      double_t E_L_;

      double_t I_e_;

      double_t U_th_;

      double_t U_min_;

      double_t U_reset_;

      interpOrder Interpol_;

      Parameters_();  

      void get(DictionaryDatum&) const;  

      double set(const DictionaryDatum&);
    };
    
    // ---------------------------------------------------------------- 

    struct State_ {
      double_t y0_; 
      double_t y1_; 
      double_t y2_; 
      double_t y3_; 
      bool     is_refractory_;    
      long_t   last_spike_step_;   
      double_t last_spike_offset_; 
      
      State_();  
      
      void get(DictionaryDatum&, const Parameters_&) const;

      void set(const DictionaryDatum&, const Parameters_&, double);
    };
    
    // ---------------------------------------------------------------- 

    struct Buffers_ {
      Buffers_(iaf_psc_alpha_canon&);
      Buffers_(const Buffers_&, iaf_psc_alpha_canon&);

      SliceRingBuffer events_;
      RingBuffer currents_; 

      UniversalDataLogger<iaf_psc_alpha_canon> logger_;
    };

    // ---------------------------------------------------------------- 

    struct Variables_ { 
      double_t h_ms_;             
      double_t PSCInitialValue_;  
      long_t   refractory_steps_; 
      double_t gamma_;            
      double_t gamma_sq_;         
      double_t expm1_tau_m_;      
      double_t expm1_tau_syn_;    
      double_t P30_;              
      double_t P31_;              
      double_t P32_;              
      double_t y0_before_;        
      double_t y2_before_;        
      double_t y3_before_;        
    };
    
    // Access functions for UniversalDataLogger -------------------------------

    double_t get_V_m_() const { return S_.y3_ + P_.E_L_; }

    double_t get_y1_() const { return S_.y1_; }

    double_t get_y2_() const { return S_.y2_; }

   // ---------------------------------------------------------------- 

   Parameters_ P_;
   State_      S_;
   Variables_  V_;
   Buffers_    B_;
   static RecordablesMap<iaf_psc_alpha_canon> recordablesMap_;
  };
  
inline
port nest::iaf_psc_alpha_canon::send_test_event(Node& target, rport receptor_type, synindex, bool)
{
  SpikeEvent e;
  e.set_sender(*this);
  return target.handles_test_event(e, receptor_type);
}
  
inline
port iaf_psc_alpha_canon::handles_test_event(SpikeEvent&, rport receptor_type)
{
  if (receptor_type != 0)
    throw UnknownReceptorType(receptor_type, get_name());
  return 0;
}
 
inline
port iaf_psc_alpha_canon::handles_test_event(CurrentEvent&, rport receptor_type)
{
  if (receptor_type != 0)
    throw UnknownReceptorType(receptor_type, get_name());
  return 0;
}
 
inline
port iaf_psc_alpha_canon::handles_test_event(DataLoggingRequest& dlr, 
                       rport receptor_type)
{
  if (receptor_type != 0)
    throw UnknownReceptorType(receptor_type, get_name());
  return B_.logger_.connect_logging_device(dlr, recordablesMap_);
}

inline
void iaf_psc_alpha_canon::get_status(DictionaryDatum &d) const
{
  P_.get(d);
  S_.get(d, P_);
  (*d)[names::recordables] = recordablesMap_.get_list();
}

inline
void iaf_psc_alpha_canon::set_status(const DictionaryDatum &d)
{
  Parameters_ ptmp = P_;  // temporary copy in case of errors
  const double delta_EL = ptmp.set(d);                       // throws if BadProperty
  State_      stmp = S_;  // temporary copy in case of errors
  stmp.set(d, ptmp, delta_EL);                 // throws if BadProperty

  // if we get here, temporaries contain consistent set of properties
  P_ = ptmp;
  S_ = stmp;
}
  
} // namespace

#endif //IAF_PSC_ALPHA_CANON_H