correlomatrix_detector.h

Go to the documentation of this file.

/*
 *  correlomatrix_detector.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 CORRELOMATRIX_DETECTOR_H
#define CORRELOMATRIX_DETECTOR_H
 
 
// C++ includes:
#include <deque>
#include <vector>
 
// Includes from nestkernel:
#include "event.h"
#include "nest_timeconverter.h"
#include "nest_types.h"
#include "node.h"
#include "pseudo_recording_device.h"
 
 
namespace nest
{
 
/* BeginUserDocs: device, detector
 
Short description
+++++++++++++++++
 
Device for measuring the covariance matrix from several inputs
 
Description
+++++++++++
 
The ``correlomatrix_detector`` is a device that receives spikes from several pools
of spike inputs and calculates the covariance matrix of inter-spike intervals
(raw auto and cross correlation) binned to bins of duration ``delta_tau``, which
defaults to 5 times the simulation resolution. The histogram is only recorded for
non-negative time lags. The negative part can be obtained by the symmetry of
the covariance matrix :math:` C(t) = C^T(-t)`.
 
The result can be obtained from the node's status dictionary under the key
``count_covariance``.
 
In parallel it records a weighted histogram, where the connection weight are
used to weight every count, which is available under the key ``/covariance``.
Both are matrices of size ``N_channels x N_channels``, with each entry :math:`C_{ij}` being
a vector of size :math:`\tau_{max}/\delta_\tau + 1` containing the (weighted) histogram
for non-negative time lags.
 
The bins are centered around the time difference they represent, and are
left-closed and right-open in the lower triangular part of the matrix. On the
diagonal and in the upper triangular part the intervals are left-open and
right-closed. This ensures proper counting of events at the border of bins,
allowing consistent integration of a histogram over negative and positive
time lags by stacking two parts of the histogram
 
.. math::
 
    (C(t)=[C[i][j][::-1],C[j][i][1:]]).
 
In this case one needs to exclude :math:`C[j][i][0]` to avoid counting the zero-lag
bin twice.
 
The ``correlomatrix_detector`` has a variable number of inputs which can be set
via ``SetStatus`` under the key ``N_channels``. All incoming connections to a
specified receptor will be pooled.
 
Correlomatrix detectors ignore any connection delays.
 
Parameters
++++++++++
 
================ ========= ====================================================
Tstart           real      Time when to start counting events. This time should
                           be set to at least start + tau_max in order to avoid
                           edge effects of the correlation counts.
Tstop            real      Time when to stop counting events. This time should
                           be set to at most Tsim - tau_max, where Tsim is the
                           duration of simulation, in order to avoid edge
                           effects of the correlation counts.
delta_tau        ms        Bin width. This has to be an odd multiple of
                           the resolution, to allow the symmetry between
                           positive and negative time-lags. Defaults to 5 times
                           the simulation resolution.
tau_max          ms        One-sided width. In the lower triagnular part
                           events with differences in [0, tau_max+delta_tau/2)
                           are counted. On the diagonal and in the upper
                           triangular part events with differences in
                           (0, tau_max+delta_tau/2]. Defaults to 10 times the
                           value of delta_tau.
N_channels       integer   The number of pools. This defines the range of
                           receptor_type. Default is 1.
                           Setting N_channels clears count_covariance,
                           covariance and n_events.
covariance       3D        matrix of read-only -raw, weighted, auto/cross
                 matrix of correlation
                 integers
count_covariance 3D        matrix of read-only -raw, auto/cross correlation
                 matrix of counts
                 integers
n_events         list of   number of events from all sources
                 integers
================ ========= ====================================================
 
Receives
++++++++
 
SpikeEvent
 
See also
++++++++
 
correlation_detector, spike_recorder
 
Examples using this model
+++++++++++++++++++++++++
 
.. listexamples:: correlomatrix_detector
 
EndUserDocs */
 
void register_correlomatrix_detector( const std::string& name );
 
class correlomatrix_detector : public Node
{
 
public:
  correlomatrix_detector();
  correlomatrix_detector( const correlomatrix_detector& );
 
  bool
  has_proxies() const override
  {
    return true;
  }
 
  std::string
  get_element_type() const override
  {
    return names::recorder;
  }
 
  using Node::handle;
  using Node::handles_test_event;
 
  void handle( SpikeEvent& ) override;
 
  size_t handles_test_event( SpikeEvent&, size_t ) override;
 
  void get_status( Dictionary& ) const override;
  void set_status( const Dictionary& ) override;
 
  void calibrate_time( const TimeConverter& tc ) override;
 
private:
  void init_state_() override;
  void init_buffers_() override;
  void pre_run_hook() override;
 
  void update( Time const&, const long, const long ) override;
 
  // ------------------------------------------------------------
 
  struct Spike_
  {
    long timestep_;
    double weight_;
    long receptor_channel_;
 
    Spike_( long timestep, double weight, long receptorchannel )
      : timestep_( timestep )
      , weight_( weight )
      , receptor_channel_( receptorchannel )
    {
    }
 
    inline bool
    operator>( const Spike_& second ) const
    {
      return timestep_ > second.timestep_;
    }
  };
 
  typedef std::deque< Spike_ > SpikelistType;
 
  // ------------------------------------------------------------
 
  struct State_;
 
  struct Parameters_
  {
    Time delta_tau_;     
    Time tau_max_;       
    Time Tstart_;        
    Time Tstop_;         
    size_t N_channels_;  
 
    Parameters_();                      
    Parameters_( const Parameters_& );  
 
    Parameters_& operator=( const Parameters_& );
 
    void get( Dictionary& ) const;  
 
    bool set( const Dictionary&, const correlomatrix_detector&, Node* node );
 
    Time get_default_delta_tau();
  };
 
  // ------------------------------------------------------------
 
  struct State_
  {
    std::vector< long > n_events_;  
    SpikelistType incoming_;        
    std::vector< std::vector< std::vector< double > > > covariance_;
 
    std::vector< std::vector< std::vector< long > > > count_covariance_;
 
    State_();  
 
    void get( Dictionary& ) const;
 
    void set( const Dictionary&, const Parameters_&, bool, Node* node );
 
    void reset( const Parameters_& );
  };
 
  // ------------------------------------------------------------
 
  PseudoRecordingDevice device_;
  Parameters_ P_;
  State_ S_;
};
 
inline size_t
correlomatrix_detector::handles_test_event( SpikeEvent&, size_t receptor_type )
{
  if ( receptor_type > P_.N_channels_ - 1 )
  {
    throw UnknownReceptorType( receptor_type, get_name() );
  }
  return receptor_type;
}
 
inline void
correlomatrix_detector::get_status( Dictionary& d ) const
{
  device_.get_status( d );
  P_.get( d );
  S_.get( d );
}
 
inline void
correlomatrix_detector::set_status( const Dictionary& d )
{
  Parameters_ ptmp = P_;
  const bool reset_required = ptmp.set( d, *this, this );
 
  device_.set_status( d );
  P_ = ptmp;
  if ( reset_required )
  {
    S_.reset( P_ );
  }
}
 
inline Time
correlomatrix_detector::Parameters_::get_default_delta_tau()
{
  return 5 * Time::get_resolution();
}
 
}  // namespace
 
#endif /* #ifndef CORRELOMATRIX_DETECTOR_H */