diff --git a/doc/htmldoc/examples/index.rst b/doc/htmldoc/examples/index.rst index 6d2356f7c2..a091a8fdf7 100644 --- a/doc/htmldoc/examples/index.rst +++ b/doc/htmldoc/examples/index.rst @@ -244,6 +244,8 @@ PyNEST examples * :doc:`../auto_examples/sensitivity_to_perturbation` * :doc:`../auto_examples/intrinsic_currents_spiking` * :doc:`../auto_examples/intrinsic_currents_subthreshold` + * :doc:`../auto_examples/ou_vs_white_noise` + * :doc:`../auto_examples/ou_noise_single_neuron` .. grid-item-card:: Handling output :img-top: ../static/img/pynest/store_restore.png @@ -343,3 +345,5 @@ PyNEST examples ../auto_examples/EI_clustered_network/index ../auto_examples/eprop_plasticity/index ../auto_examples/wang_decision_making + ../auto_examples/ou_vs_white_noise + ../auto_examples/ou_noise_single_neuron diff --git a/models/ou_noise_generator.cpp b/models/ou_noise_generator.cpp new file mode 100644 index 0000000000..ecd473a8b7 --- /dev/null +++ b/models/ou_noise_generator.cpp @@ -0,0 +1,372 @@ +/* + * ou_noise_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 . + * + */ + +#include "ou_noise_generator.h" + +// Includes from libnestutil: +#include "compose.hpp" +#include "dict_util.h" +#include "logging.h" +#include "numerics.h" + +// Includes from nestkernel: +#include "event_delivery_manager_impl.h" +#include "kernel_manager.h" +#include "nest_impl.h" +#include "universal_data_logger_impl.h" + +// Includes from sli: +#include "dict.h" +#include "dictutils.h" +#include "doubledatum.h" + +namespace nest +{ +void +register_ou_noise_generator( const std::string& name ) +{ + register_node_model< ou_noise_generator >( name ); +} + +RecordablesMap< ou_noise_generator > ou_noise_generator::recordablesMap_; + +template <> +void +RecordablesMap< ou_noise_generator >::create() +{ + insert_( Name( names::I ), &ou_noise_generator::get_I_avg_ ); +} +} + +/* ---------------------------------------------------------------- + * Default constructors defining default parameter + * ---------------------------------------------------------------- */ + +nest::ou_noise_generator::Parameters_::Parameters_() + : mean_( 0.0 ) // pA + , std_( 0.0 ) // pA + , tau_( 1.0 ) // ms + , dt_( get_default_dt() ) + , num_targets_( 0 ) +{ +} + +nest::ou_noise_generator::Parameters_::Parameters_( const Parameters_& p ) + : mean_( p.mean_ ) + , std_( p.std_ ) + , tau_( p.tau_ ) + , dt_( p.dt_ ) + , num_targets_( 0 ) // we do not copy connections +{ + if ( dt_.is_step() ) + { + dt_.calibrate(); + } + else + { + dt_ = get_default_dt(); + } +} + +nest::ou_noise_generator::Parameters_& +nest::ou_noise_generator::Parameters_::operator=( const Parameters_& p ) +{ + if ( this == &p ) + { + return *this; + } + + mean_ = p.mean_; + std_ = p.std_; + tau_ = p.tau_; + dt_ = p.dt_; + + return *this; +} + +nest::ou_noise_generator::State_::State_() + : I_avg_( 0.0 ) // pA +{ +} + +nest::ou_noise_generator::Buffers_::Buffers_( ou_noise_generator& n ) + : next_step_( 0 ) + , logger_( n ) +{ +} + +nest::ou_noise_generator::Buffers_::Buffers_( const Buffers_& b, ou_noise_generator& n ) + : next_step_( b.next_step_ ) + , logger_( n ) +{ +} + +/* ---------------------------------------------------------------- + * Parameter extraction and manipulation functions + * ---------------------------------------------------------------- */ + +void +nest::ou_noise_generator::Parameters_::get( DictionaryDatum& d ) const +{ + ( *d )[ names::mean ] = mean_; + ( *d )[ names::std ] = std_; + ( *d )[ names::dt ] = dt_.get_ms(); + ( *d )[ names::tau ] = tau_; +} + +void +nest::ou_noise_generator::State_::get( DictionaryDatum& d ) const +{ +} + +void +nest::ou_noise_generator::Parameters_::set( const DictionaryDatum& d, const ou_noise_generator& n, Node* node ) +{ + updateValueParam< double >( d, names::mean, mean_, node ); + updateValueParam< double >( d, names::std, std_, node ); + updateValueParam< double >( d, names::tau, tau_, node ); + if ( tau_ <= 0 ) + { + throw BadProperty( "tau > 0 required." ); + } + double dt; + if ( updateValueParam< double >( d, names::dt, dt, node ) ) + { + dt_ = Time::ms( dt ); + } + if ( std_ < 0 ) + { + throw BadProperty( "The standard deviation cannot be negative." ); + } + + if ( not dt_.is_step() ) + { + throw StepMultipleRequired( n.get_name(), names::dt, dt_ ); + } +} + + +/* ---------------------------------------------------------------- + * Default and copy constructor for node + * ---------------------------------------------------------------- */ + +nest::ou_noise_generator::ou_noise_generator() + : StimulationDevice() + , P_() + , S_() + , B_( *this ) +{ + recordablesMap_.create(); +} + +nest::ou_noise_generator::ou_noise_generator( const ou_noise_generator& n ) + : StimulationDevice( n ) + , P_( n.P_ ) + , S_( n.S_ ) + , B_( n.B_, *this ) +{ +} + + +/* ---------------------------------------------------------------- + * Node initialization functions + * ---------------------------------------------------------------- */ + +void +nest::ou_noise_generator::init_state_() +{ + StimulationDevice::init_state(); +} + +void +nest::ou_noise_generator::init_buffers_() +{ + StimulationDevice::init_buffers(); + B_.logger_.reset(); + + B_.next_step_ = 0; + B_.amps_.clear(); + B_.amps_.resize( P_.num_targets_, 0.0 ); +} + +void +nest::ou_noise_generator::pre_run_hook() +{ + B_.logger_.init(); + + StimulationDevice::pre_run_hook(); + if ( P_.num_targets_ != B_.amps_.size() ) + { + LOG( M_INFO, "ou_noise_generator::pre_run_hook()", "The number of targets has changed, drawing new amplitudes." ); + init_buffers_(); + } + + V_.dt_steps_ = P_.dt_.get_steps(); + + const double h = P_.dt_.get_ms(); + + V_.prop_ = std::exp( -h / P_.tau_ ); + V_.noise_amp_ = P_.std_ * std::sqrt( -std::expm1( -2.0 * h / P_.tau_ ) ); + V_.mean_weight_ = 1.0 - V_.prop_; + V_.mean_incr_ = P_.mean_ * V_.mean_weight_; +} + + +/* ---------------------------------------------------------------- + * Update function and event hook + * ---------------------------------------------------------------- */ + +size_t +nest::ou_noise_generator::send_test_event( Node& target, size_t receptor_type, synindex syn_id, bool dummy_target ) +{ + StimulationDevice::enforce_single_syn_type( syn_id ); + + if ( dummy_target ) + { + DSCurrentEvent e; + e.set_sender( *this ); + return target.handles_test_event( e, receptor_type ); + } + else + { + CurrentEvent e; + e.set_sender( *this ); + const size_t p = target.handles_test_event( e, receptor_type ); + if ( p != invalid_port and not is_model_prototype() ) + { + ++P_.num_targets_; + } + return p; + } +} + +// +// Time Evolution Operator +// +void +nest::ou_noise_generator::update( Time const& origin, const long from, const long to ) +{ + const long start = origin.get_steps(); + + for ( long offs = from; offs < to; ++offs ) + { + + const long now = start + offs; + + // Update OU process state if next update is due + if ( now >= B_.next_step_ ) + { + // compute new currents + for ( double& amp : B_.amps_ ) + { + amp = V_.mean_incr_ + amp * V_.prop_ + V_.noise_amp_ * V_.normal_dist_( get_vp_specific_rng( get_thread() ) ); + } + // use now as reference, in case we woke up from inactive period + B_.next_step_ = now + V_.dt_steps_; + } + + S_.I_avg_ = 0.0; + + // Skip sending events while generator is inactiv + if ( not StimulationDevice::is_active( Time::step( now ) ) ) + { + B_.logger_.record_data( origin.get_steps() + offs ); + continue; + } + + // record values + for ( double& amp : B_.amps_ ) + { + S_.I_avg_ += amp; + } + + S_.I_avg_ /= std::max( 1, int( B_.amps_.size() ) ); + B_.logger_.record_data( origin.get_steps() + offs ); + + DSCurrentEvent ce; + kernel().event_delivery_manager.send( *this, ce, offs ); + } +} + +void +nest::ou_noise_generator::event_hook( DSCurrentEvent& e ) +{ + // get port number + const size_t prt = e.get_port(); + + // we handle only one port here, get reference to vector elem + assert( prt < B_.amps_.size() ); + + e.set_current( B_.amps_[ prt ] ); + e.get_receiver().handle( e ); +} + +void +nest::ou_noise_generator::handle( DataLoggingRequest& e ) +{ + B_.logger_.handle( e ); +} + +/* ---------------------------------------------------------------- + * Other functions + * ---------------------------------------------------------------- */ + +void +nest::ou_noise_generator::set_data_from_stimulation_backend( std::vector< double >& input_param ) +{ + Parameters_ ptmp = P_; // temporary copy in case of errors + ptmp.num_targets_ = P_.num_targets_; + + // For the input backend + if ( not input_param.empty() ) + { + if ( input_param.size() != 3 ) + { + throw BadParameterValue( "The size of the data for the ou_noise_generator needs to be 3 [mean, std, tau]." ); + } + DictionaryDatum d = DictionaryDatum( new Dictionary ); + ( *d )[ names::mean ] = DoubleDatum( input_param[ 0 ] ); + ( *d )[ names::std ] = DoubleDatum( input_param[ 1 ] ); + ( *d )[ names::tau ] = DoubleDatum( input_param[ 2 ] ); + ptmp.set( d, *this, this ); + } + + // if we get here, temporary contains consistent set of properties + P_ = ptmp; + P_.num_targets_ = ptmp.num_targets_; +} + +void +nest::ou_noise_generator::calibrate_time( const TimeConverter& tc ) +{ + if ( P_.dt_.is_step() ) + { + P_.dt_ = tc.from_old_tics( P_.dt_.get_tics() ); + } + else + { + const double old = P_.dt_.get_ms(); + P_.dt_ = P_.get_default_dt(); + std::string msg = String::compose( "Default for dt changed from %1 to %2 ms", old, P_.dt_.get_ms() ); + LOG( M_INFO, get_name(), msg ); + } +} diff --git a/models/ou_noise_generator.h b/models/ou_noise_generator.h new file mode 100644 index 0000000000..13106c0556 --- /dev/null +++ b/models/ou_noise_generator.h @@ -0,0 +1,320 @@ +/* + * ou_noise_generator.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 . + * + */ + +#ifndef OU_NOISE_GENERATOR_H +#define OU_NOISE_GENERATOR_H + +// C++ includes: +#include + +// Includes from nestkernel: +#include "connection.h" +#include "device_node.h" +#include "event.h" +#include "nest_timeconverter.h" +#include "nest_types.h" +#include "random_generators.h" +#include "stimulation_device.h" +#include "universal_data_logger.h" + +namespace nest +{ + +/* BeginUserDocs: device, generator + +Short description ++++++++++++++++++ + +Generates a temporally correlated noise current based on an Ornstein-Uhlenbeck process. + +Description ++++++++++++ + +The `ou_noise_generator` can be used to inject a temporally correlated noise current into a node. +The current I(t) follows an Ornstein-Uhlenbeck (OU) process, which is described by the following stochastic differential +equation: + +.. math:: + + dI = \frac{1}{\tau}(\mu - I)dt + \sigma_{stat} \sqrt{\frac{2}{\tau}} dW + +where: + - :math:`\mu` is the long-term mean of the process (`mean` parameter). + - :math:`\tau` is the time constant of the correlation (`tau` parameter). + - :math:`\sigma_{stat}` is the stationary standard deviation of the process (`std` parameter). + - :math:`dW` is a Wiener process (Gaussian white noise). + +The generator integrates this process at a user-defined interval `dt` and delivers the resulting current to its targets. +A larger time constant :math:`\tau` results in a more slowly varying noise signal. + +All targets of a noise generator receive different, independent noise currents, but the currents for all targets are +updated at the same points in time. The interval `dt` between updates must be a multiple of the simulation time step. + +.. admonition:: Recording the generated current + + You can use a :doc:`multimeter ` to record the average current sent to all targets for each time step if + simulating on a single thread; multiple MPI processes with one thread each also work. In this case, the recording + interval of the multimeter should be equal to the `dt` of the generator to avoid aliasing effects. In multi-threaded + mode, recording of noise currents is prohibited for technical reasons. + + +.. include:: ../models/stimulation_device.rst + +mean + The mean value :math:`\mu` to which the process reverts (pA). + +std + The stationary standard deviation :math:`\sigma_{stat}` of the process (pA). + +tau + The correlation time constant :math:`\tau` of the process (ms). + +dt + The time interval between updates of the noise current (ms). + + +Setting parameters from a stimulation backend +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +The parameters in this stimulation device can be updated with input +coming from a stimulation backend. The data structure used for the +update holds one value for each of the parameters mentioned above. +The indexing is as follows: + + 0. mean + 1. std + 2. tau + + +EndUserDocs */ + +void register_ou_noise_generator( const std::string& name ); + +class ou_noise_generator : public StimulationDevice +{ + +public: + ou_noise_generator(); + ou_noise_generator( const ou_noise_generator& ); + + //! Allow multimeter to connect to local instances + bool local_receiver() const override; + + /** + * Import sets of overloaded virtual functions. + * @see Technical Issues / Virtual Functions: Overriding, Overloading, and + * Hiding + */ + using Node::event_hook; + using Node::handle; + using Node::handles_test_event; + using Node::sends_signal; + + size_t send_test_event( Node&, size_t, synindex, bool ) override; + + SignalType sends_signal() const override; + + void handle( DataLoggingRequest& ) override; + + size_t handles_test_event( DataLoggingRequest&, size_t ) override; + + void get_status( DictionaryDatum& ) const override; + void set_status( const DictionaryDatum& ) override; + + void calibrate_time( const TimeConverter& tc ) override; + + StimulationDevice::Type get_type() const override; + void set_data_from_stimulation_backend( std::vector< double >& input_param ) override; + +private: + void init_state_() override; + void init_buffers_() override; + + /** + * Recalculates parameters and forces reinitialization + * of amplitudes if number of targets has changed. + */ + void pre_run_hook() override; + + void update( Time const&, const long, const long ) override; + void event_hook( DSCurrentEvent& ) override; + + // ------------------------------------------------------------ + + typedef std::vector< double > AmpVec_; + + /** + * Store independent parameters of the model. + */ + struct Parameters_ + { + double mean_; //!< mean current, in pA + double std_; //!< standard deviation of current, in pA + double tau_; //!< OU time constant, in ms + Time dt_; //!< time interval between updates + + /** + * Number of targets. + * This is a hidden parameter; must be placed in parameters, + * even though it is an implementation detail, since it + * concerns the connections and must not be affected by resets. + */ + size_t num_targets_; + + Parameters_(); //!< Sets default parameter values + Parameters_( const Parameters_& ); + Parameters_& operator=( const Parameters_& p ); + + void get( DictionaryDatum& ) const; //!< Store current values in dictionary + //! Set values from dictionary + void set( const DictionaryDatum&, const ou_noise_generator&, Node* node ); + + Time get_default_dt(); + }; + + // ------------------------------------------------------------ + + struct State_ + { + double I_avg_; //!< Average of instantaneous currents computed + //!< Used for recording current + + State_(); //!< Sets default parameter values + + void get( DictionaryDatum& ) const; //!< Store current values in dictionary + }; + + // ------------------------------------------------------------ + + // The next two classes need to be friends to access the State_ class/member + friend class RecordablesMap< ou_noise_generator >; + friend class UniversalDataLogger< ou_noise_generator >; + + // ------------------------------------------------------------ + + struct Buffers_ + { + long next_step_; //!< time step of next change in current + AmpVec_ amps_; //!< amplitudes, one per target + explicit Buffers_( ou_noise_generator& ); + Buffers_( const Buffers_&, ou_noise_generator& ); + UniversalDataLogger< ou_noise_generator > logger_; + }; + + // ------------------------------------------------------------ + + struct Variables_ + { + normal_distribution normal_dist_; //!< normal distribution + + long dt_steps_; //!< update interval in steps + double prop_; //!< frequency [radian/s] + double noise_amp_; //!< noise amplitude for one update, in pA + double mean_weight_; //!< weight for mean contribution + double mean_incr_; //!< precomputed mean contribution per step + }; + + double + get_I_avg_() const + { + return S_.I_avg_; + } + + // ------------------------------------------------------------ + + static RecordablesMap< ou_noise_generator > recordablesMap_; + + Parameters_ P_; + State_ S_; + Buffers_ B_; + Variables_ V_; +}; + +inline size_t +ou_noise_generator::handles_test_event( DataLoggingRequest& dlr, size_t receptor_type ) +{ + if ( kernel().vp_manager.get_num_threads() > 1 ) + { + throw KernelException( "Recording from a ou_noise_generator is only possible in single-threaded mode." ); + } + + if ( receptor_type != 0 ) + { + throw UnknownReceptorType( receptor_type, get_name() ); + } + return B_.logger_.connect_logging_device( dlr, recordablesMap_ ); +} + +inline void +ou_noise_generator::get_status( DictionaryDatum& d ) const +{ + P_.get( d ); + S_.get( d ); + StimulationDevice::get_status( d ); + + ( *d )[ names::recordables ] = recordablesMap_.get_list(); +} + +inline void +ou_noise_generator::set_status( const DictionaryDatum& d ) +{ + Parameters_ ptmp = P_; // temporary copy in case of errors + ptmp.num_targets_ = P_.num_targets_; // Copy Constr. does not copy connections + ptmp.set( d, *this, this ); // throws if BadProperty + + // We now know that ptmp is consistent. We do not write it back + // to P_ before we are also sure that the properties to be set + // in the parent class are internally consistent. + StimulationDevice::set_status( d ); + + // if we get here, temporaries contain consistent set of properties + P_ = ptmp; + P_.num_targets_ = ptmp.num_targets_; +} + +inline SignalType +ou_noise_generator::sends_signal() const +{ + return ALL; +} + +inline bool +ou_noise_generator::local_receiver() const +{ + return true; +} + +inline StimulationDevice::Type +ou_noise_generator::get_type() const +{ + return StimulationDevice::Type::CURRENT_GENERATOR; +} + +inline Time +ou_noise_generator::Parameters_::get_default_dt() +{ + return 10 * Time::get_resolution(); +} + +} // namespace + +#endif // OU_NOISE_GENERATOR_H diff --git a/modelsets/full b/modelsets/full index 46d6c6ee5a..9d2821603d 100644 --- a/modelsets/full +++ b/modelsets/full @@ -89,6 +89,7 @@ music_rate_in_proxy music_rate_out_proxy music_message_in_proxy noise_generator +ou_noise_generator parrot_neuron parrot_neuron_ps inhomogeneous_poisson_generator diff --git a/pynest/examples/ou_noise_single_neuron.py b/pynest/examples/ou_noise_single_neuron.py new file mode 100644 index 0000000000..36bc23aa13 --- /dev/null +++ b/pynest/examples/ou_noise_single_neuron.py @@ -0,0 +1,78 @@ +# -*- coding: utf-8 -*- +# +# ou_noise_single_neuron.py +# +# 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 . + +""" +One neuron with Ornstein-Uhlenbeck noise +---------------------------------------- + +This script simulates a neuron with input from the ``ou_noise_generator`` and +records the neuron's membrane potential. +""" + +############################################################################### +# First, we import the modules needed to simulate, analyze, and plot the +# example. We also reset the kernel. + +import matplotlib.pyplot as plt +import nest + +nest.ResetKernel() +nest.resolution = 0.1 + + +############################################################################### +# Second, we create the neuron, noise generator and voltmeter. + +neuron = nest.Create("iaf_psc_alpha") +ou_noise = nest.Create("ou_noise_generator") +voltmeter = nest.Create("voltmeter") + +############################################################################### +# Third, we configure the OU-noise generator. The parameters are: +# ``mean`` and ``std`` in `pA`, ``tau`` (correlation time) and ``dt`` in `ms`. +# The update interval ``dt`` must be a multiple of the simulation resolution. + +ou_noise.set( + { + "mean": 250.0, # pA + "std": 50.0, # pA + "tau": 50.0, # ms + "dt": nest.resolution, # ms + } +) + +############################################################################### +# Fourth, we connect the generator to the neuron and the voltmeter to the +# neuron to record the membrane potential. + +nest.Connect(ou_noise, neuron) +nest.Connect(voltmeter, neuron) + +############################################################################### +# Now we simulate the network for 1000 ms. + +nest.Simulate(1000.0) + +############################################################################### +# Finally, we plot the neuron's membrane potential as a function of time. + +nest.voltage_trace.from_device(voltmeter) +plt.show() diff --git a/pynest/examples/ou_vs_white_noise.py b/pynest/examples/ou_vs_white_noise.py new file mode 100644 index 0000000000..b1a47dfb4d --- /dev/null +++ b/pynest/examples/ou_vs_white_noise.py @@ -0,0 +1,121 @@ +# -*- coding: utf-8 -*- +# +# ou_vs_white_noise.py +# +# 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 . + +""" +Comparing white and OU noise +---------------------------- + +This example compares two neurons: one driven by the ``noise_generator`` (white noise) +and one driven by the ``ou_noise_generator`` (temporally correlated noise). +We match the mean and variance of both inputs and report firing rate, +ISI coefficient of variation, and membrane-potential statistics. +""" + +import nest +import numpy as np + +############################################################################### +# First, we reset the kernel and define the simulation resolution. +# The update interval of the generators must be an +# integer multiple of the resolution. + +nest.ResetKernel() +nest.resolution = 0.1 +dt = nest.resolution + +############################################################################### +# Second, we create two identical LIF neurons, one multimeter and one spike +# recorder per neuron. + +simtime = 500_000.0 + +lif_params = dict(C_m=250.0, tau_m=20.0, t_ref=2.0, V_reset=-70.0, E_L=-70.0, V_th=-54.0) +n_white = nest.Create("iaf_psc_alpha", params=lif_params) +n_ou = nest.Create("iaf_psc_alpha", params=lif_params) + +mm_white = nest.Create("multimeter", params={"record_from": ["V_m"], "interval": dt}) +mm_ou = nest.Create("multimeter", params={"record_from": ["V_m"], "interval": dt}) +sr_white = nest.Create("spike_recorder") +sr_ou = nest.Create("spike_recorder") + +############################################################################### +# Third, we create and configure the generators. We use the same mean and +# stationary standard deviation for both inputs and set ``tau`` only for the +# OU generator. + +mean_I, std_I = 250.0, 50.0 +tau_ou = 50.0 +ng_white = nest.Create("noise_generator", {"mean": mean_I, "std": std_I, "dt": dt}) +ng_ou = nest.Create("ou_noise_generator", {"mean": mean_I, "std": std_I, "tau": tau_ou, "dt": dt}) + +############################################################################### +# Fourth, we connect generators to neurons and measurement devices to neurons. + +nest.Connect(ng_white, n_white) +nest.Connect(n_white, sr_white) +nest.Connect(mm_white, n_white) +nest.Connect(ng_ou, n_ou) +nest.Connect(n_ou, sr_ou) +nest.Connect(mm_ou, n_ou) + +############################################################################### +# Now we simulate the network. + +nest.Simulate(simtime) + +############################################################################### +# Finally, we compute the membrane-potential 'mean/std' (after a 1 s burn-in), +# firing rate, and ISI coefficient of variation. + + +def isi_cv(spike_times): + if len(spike_times) < 3: + return np.nan + isi = np.diff(np.asarray(spike_times)) + m = isi.mean() + return np.nan if m == 0 else isi.std(ddof=1) / m + + +def vm_stats(mm, burn_in_ms=1000.0): + ev = mm.get("events") + t = np.asarray(ev["times"]) + vm = np.asarray(ev["V_m"]) + keep = t >= burn_in_ms + vm = vm[keep] + return vm.mean(), vm.std(ddof=1) + + +def spike_stats(sr, simtime_ms): + ev = sr.get("events") + times = np.asarray(ev["times"]) + rate = 1000.0 * len(times) / simtime_ms + return rate, isi_cv(times) + + +vm_mean_w, vm_std_w = vm_stats(mm_white) +vm_mean_o, vm_std_o = vm_stats(mm_ou) +rate_w, cv_w = spike_stats(sr_white, simtime) +rate_o, cv_o = spike_stats(sr_ou, simtime) + +print("=== White noise (noise_generator) ===") +print(f"Rate: {rate_w:.2f} Hz, ISI CV: {cv_w:.2f}, Vm mean: {vm_mean_w:.2f} mV, Vm std: {vm_std_w:.2f} mV") +print("=== OU noise (ou_noise_generator) ===") +print(f"Rate: {rate_o:.2f} Hz, ISI CV: {cv_o:.2f}, Vm mean: {vm_mean_o:.2f} mV, Vm std: {vm_std_o:.2f} mV") diff --git a/testsuite/pytests/sli2py_recording/test_multimeter_stepping.py b/testsuite/pytests/sli2py_recording/test_multimeter_stepping.py index c800c5ae61..1cf7c80de6 100644 --- a/testsuite/pytests/sli2py_recording/test_multimeter_stepping.py +++ b/testsuite/pytests/sli2py_recording/test_multimeter_stepping.py @@ -52,6 +52,7 @@ "ac_generator", # generator device, does not support spike input "dc_generator", # generator device, does not support spike input "noise_generator", # generator device, does not support spike input + "ou_noise_generator", # generator device, does not support spike input "step_current_generator", # generator device, does not support spike input "step_rate_generator", # generator device, does not support spike input "sinusoidal_poisson_generator", # generator device, does not support spike input diff --git a/testsuite/pytests/test_changing_tic_base.py b/testsuite/pytests/test_changing_tic_base.py index 6caa62b353..535091785e 100644 --- a/testsuite/pytests/test_changing_tic_base.py +++ b/testsuite/pytests/test_changing_tic_base.py @@ -52,6 +52,7 @@ def test_models(self): "correlomatrix_detector": ["delta_tau"], "correlospinmatrix_detector": ["delta_tau"], "noise_generator": ["dt"], + "ou_noise_generator": ["dt"], } # Generate a dictionary of reference values for each model. diff --git a/testsuite/pytests/test_ou_noise_generator.py b/testsuite/pytests/test_ou_noise_generator.py new file mode 100644 index 0000000000..2ac0731ab8 --- /dev/null +++ b/testsuite/pytests/test_ou_noise_generator.py @@ -0,0 +1,155 @@ +# -*- coding: utf-8 -*- +# +# test_ou_noise_generator.py +# +# 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 . + +""" +Tests parameter setting and statistical correctness for one application. +""" + +import nest +import numpy as np +import pytest +from scipy.signal import fftconvolve + + +@pytest.fixture +def prepare_kernel(): + nest.ResetKernel() + nest.resolution = 0.1 + + +def burn_in_start(dt, tau, k=10): + """Return number of steps to discard for a k*tau burn-in.""" + return int(round((k * tau) / dt)) + + +def test_ou_noise_generator_set_parameters(prepare_kernel): + params = {"mean": 210.0, "std": 60.0, "dt": 0.1} + + oung1 = nest.Create("ou_noise_generator") + oung1.set(params) + + nest.SetDefaults("ou_noise_generator", params) + + oung2 = nest.Create("ou_noise_generator") + assert oung1.get(params) == oung2.get(params) + + +def test_ou_noise_generator_incorrect_noise_dt(prepare_kernel): + with pytest.raises(nest.kernel.NESTError, match="StepMultipleRequired"): + nest.Create("ou_noise_generator", {"dt": 0.25}) + + +def test_ou_noise_mean_and_variance(prepare_kernel): + # run for resolution dt=0.1 project to iaf_psc_alpha. + # create 100 repetitions of 1000ms simulations + + oung = nest.Create("ou_noise_generator", {"mean": 0.0, "std": 60.0, "tau": 1.0, "dt": 0.1}) + neuron = nest.Create("iaf_psc_alpha") + + # we need to connect to a neuron otherwise the generator does not generate + nest.Connect(oung, neuron) + mm = nest.Create("multimeter", 1, {"record_from": ["I"], "interval": 0.1}) + nest.Connect(mm, oung, syn_spec={"weight": 1}) + + # Simulate for 100 times + n_sims = 100 + ou_current = np.empty(n_sims) + for i in range(n_sims): + nest.Simulate(1000.0) + ou_current[i] = mm.get("events")["I"][-1] + + curr_mean = np.mean(ou_current) + curr_var = np.var(ou_current) + expected_curr_mean = oung.mean + expected_curr_var = oung.std**2 + + # require mean within 3 std dev, std dev within 0.2 std dev + assert np.abs(curr_mean - expected_curr_mean) < 3 * oung.std + assert np.abs(curr_var - expected_curr_var) < 0.2 * curr_var + + +def test_ou_noise_generator_autocorrelation(prepare_kernel): + # verify lag-1 autocorr = exp(-dt/tau) + dt = 0.1 + tau = 1.0 + + oung = nest.Create("ou_noise_generator", {"mean": 0.0, "std": 60.0, "tau": tau, "dt": nest.resolution}) + neuron = nest.Create("iaf_psc_alpha") + + # we need to connect to a neuron otherwise the generator does not generate + nest.Connect(oung, neuron) + mm = nest.Create("multimeter", 1, {"record_from": ["I"], "interval": dt}) + nest.Connect(mm, oung, syn_spec={"weight": 1}) + nest.Simulate(10000.0) + ou_current = mm.get("events")["I"] + + # drop first k*tau + cutoff = burn_in_start(dt=dt, tau=tau) + x = ou_current[cutoff:] + x -= x.mean() + + # empirical lag-1 autocorrelation + emp_ac1 = np.dot(x[:-1], x[1:]) / ((len(x) - 1) * x.var()) + # theoretical lag-1: exp(-dt/tau) + theor_ac1 = np.exp(-dt / tau) + + assert abs(emp_ac1 - theor_ac1) < 0.05, f"autocorr {emp_ac1:.3f} vs theoretical {theor_ac1:.3f}" + + +def test_cross_correlation(prepare_kernel): + # two neurons driven by the same OU noise should remain uncorrelated + dt, tau = 0.1, 1.0 + + ou = nest.Create("ou_noise_generator", {"mean": 0.0, "std": 50.0, "tau": tau, "dt": dt}) + neurons = nest.Create("iaf_psc_alpha", 2, {"E_L": 0.0, "V_th": 1e9, "C_m": 1.0, "tau_m": tau}) + nest.Connect(ou, neurons) + + mm = nest.Create("multimeter", params={"record_from": ["V_m"], "interval": dt}) + nest.Connect(mm, neurons) + + simtime = 50000.0 + nest.Simulate(simtime) + + ev = mm.get("events") + senders = np.asarray(ev["senders"], int) + times = np.asarray(ev["times"]) + vms = np.asarray(ev["V_m"]) + + # build time array + n_steps = int(round(simtime / dt)) + 1 + V = np.zeros((n_steps, 2)) + for t, s, v in zip(times, senders, vms): + idx = int(round(t / dt)) + col = neurons.index(s) + V[idx, col] = v + + # discard burn-in + start = burn_in_start(dt, tau) + X1 = V[start:, 0] - V[start:, 0].mean() + X2 = V[start:, 1] - V[start:, 1].mean() + + # cross-corr via FFT + corr = fftconvolve(X1, X2[::-1], mode="full") + corr /= np.sqrt(np.dot(X1, X1) * np.dot(X2, X2)) + mid = corr.size // 2 + + # maximum absolute cross-corr should be near zero + assert np.max(np.abs(corr[mid:])) < 0.05