Merge pull request nest#2960 from ackurth/poisson_multapses

Add connection rule for Poisson-distributed synapse numbers—first merge of 2024!

doc/htmldoc/networks/spatially_structured_networks.rst

@@ -469,14 +469,16 @@ Pool
    +--------------------------------+--------------+--------------+
    | Connection parameters          | Driver       | Pool         |
    +================================+==============+==============+
-   | ``rule='pairwise_bernoulli'``  | source layer | target layer |
+   | ``rule='fixed_indegree'``      | target layer | source layer |
    +--------------------------------+--------------+--------------+
    | ``rule='fixed_outdegree'``     | source layer | target layer |
    +--------------------------------+--------------+--------------+
+   | ``rule='pairwise_bernoulli'``  | source layer | target layer |
+   +--------------------------------+--------------+--------------+
    | ``rule='pairwise_bernoulli'``  | target layer | source layer |
    | and ``use_on_source=True``     |              |              |
    +--------------------------------+--------------+--------------+
-   | ``rule='fixed_indegree'``      | target layer | source layer |
+   | ``rule='pairwise_poisson'``    | source layer | target layer |
    +--------------------------------+--------------+--------------+
 
 Displacement
@@ -500,6 +502,13 @@ Connection probability or ``p``
    The default probability is :math:`1`, i.e., connections are created with
    certainty. See section :ref:`sec_conn_kernels` for details.
 
+Pairwise average number of connections or ``pairwise_avg_num_conns``
+   The *pairwise average number of connections* between a driver and a pool node.
+   It is used in the ``pairwise_poisson`` connection rule and determines
+   the mean value of the Poisson distribution from which the number of
+   connections between the nodes is sampled. See section
+   :ref:`sec_conn_kernels` for details.
+
 Autapse
    An *autapse* is a synapse (connection) from a node onto itself.
    Autapses are permitted by default, but can be disabled by adding
@@ -832,10 +841,14 @@ Probabilistic connection rules
 
 Many neuronal network models employ probabilistic connection rules.
 NEST supports probabilistic connections through the
-``pairwise_bernoulli`` connection rule. The probability can then be a constant,
-depend on the position of the source or the target neuron, or on the
-distance between a driver and a pool node to a connection probability. To
-create dependencies on neuron positions, NEST parameters objects are used.
+``pairwise_bernoulli`` and the ``pairwise_poisson`` connection rule.
+For the ``pairwise_bernoulli`` rule, the probability can then be a
+constant, depend on the position of the source or the target neuron,
+or on the distance between a driver and a pool node to a connection probability.
+For the ``pairwise_poisson`` rule, the pairwise average number of connections
+can depend on the position of the source or the target neuron,
+or on the distance between a driver and a pool node to a connection probability.
+To create dependencies on neuron positions, NEST :py:class:`.Parameter` objects are used.
 NEST then generates a connection according to this probability.
 
 Probabilistic connections between layers can be generated in two different
@@ -849,6 +862,15 @@ Free probabilistic connections using ``pairwise_bernoulli``
    driver-pool pair is inspected exactly once* and that there will be *at
    most one connection between each driver-pool pair*.
 
+Free probabilistic connections using ``pairwise_poisson``
+   In this case, :py:func:`.Connect` considers each driver node :math:`D` in turn.
+   For each :math:`D`, it evaluates the parameter value for each pool node
+   :math:`P` within the mask and creates a number of connections that are sampled
+   from a Poisson distribution with mean of the parameter value. This means in
+   particular that *each possible driver-pool pair is inspected exactly once*
+   and that *more than one connection between each driver-pool pair is possible*.
+   Additionally, the parameter may be larger than :math:`1`.
+
 Prescribed number of connections
    can be obtained by using ``fixed_indegree`` or ``fixed_outdegree``
    connection rule, and specifying the number of connections to create

doc/htmldoc/synapses/connection_management.rst

@@ -256,6 +256,26 @@ must be ``True``.
                       'allow_autapses': False, 'make_symmetric': True}
     nest.Connect(A, B, conn_spec_dict)
 
+pairwise poisson
+~~~~~~~~~~~~~~~~
+
+For each possible pair of nodes from ``A`` and ``B``, a number of
+connections is created following a Poisson distribution with mean
+``pairwise_avg_num_conns``. This means that even for a small
+average number of connections between single neurons in ``A`` and
+``B`` multiple connections are possible. Thus, for this rule
+``allow_multapses`` cannot be ``False``.
+The ``pairwise_avg_num_conns`` can be greater than one.
+
+.. code-block:: python
+
+    n, m, p_avg_num_conns = 10, 12, 0.2
+    A = nest.Create('iaf_psc_alpha', n)
+    B = nest.Create('iaf_psc_alpha', m)
+    conn_spec_dict = {'rule': 'pairwise_poisson',
+                      'pairwise_avg_num_conns': p_avg_num_conns}
+    nest.Connect(A, B, conn_spec_dict)
+
 .. _tripartite_connectivity:
 
 Tripartite connectivity

nestkernel/conn_builder.cpp

@@ -1589,6 +1589,124 @@ nest::BernoulliBuilder::inner_connect_( const int tid, RngPtr rng, Node* target,
 }
 
 
+nest::PoissonBuilder::PoissonBuilder( NodeCollectionPTR sources,
+  NodeCollectionPTR targets,
+  const DictionaryDatum& conn_spec,
+  const std::vector< DictionaryDatum >& syn_specs )
+  : ConnBuilder( sources, targets, conn_spec, syn_specs )
+{
+  ParameterDatum* pd = dynamic_cast< ParameterDatum* >( ( *conn_spec )[ names::pairwise_avg_num_conns ].datum() );
+  if ( pd )
+  {
+    pairwise_avg_num_conns_ = *pd;
+  }
+  else
+  {
+    // Assume pairwise_avg_num_conns is a scalar
+    const double value = ( *conn_spec )[ names::pairwise_avg_num_conns ];
+    if ( value < 0 )
+    {
+      throw BadProperty( "Connection parameter 0 ≤ pairwise_avg_num_conns required." );
+    }
+    if ( not allow_multapses_ )
+    {
+      throw BadProperty( "Multapses must be allowed for this connection rule." );
+    }
+    pairwise_avg_num_conns_ = std::shared_ptr< Parameter >( new ConstantParameter( value ) );
+  }
+}
+
+void
+nest::PoissonBuilder::connect_()
+{
+#pragma omp parallel
+  {
+    // get thread id
+    const size_t tid = kernel().vp_manager.get_thread_id();
+
+    try
+    {
+      RngPtr rng = get_vp_specific_rng( tid );
+
+      if ( loop_over_targets_() )
+      {
+        NodeCollection::const_iterator target_it = targets_->begin();
+        for ( ; target_it < targets_->end(); ++target_it )
+        {
+          const size_t tnode_id = ( *target_it ).node_id;
+          Node* const target = kernel().node_manager.get_node_or_proxy( tnode_id, tid );
+          if ( target->is_proxy() )
+          {
+            // skip parameters handled in other virtual processes
+            skip_conn_parameter_( tid );
+            continue;
+          }
+
+          inner_connect_( tid, rng, target, tnode_id );
+        }
+      }
+      else
+      {
+        const SparseNodeArray& local_nodes = kernel().node_manager.get_local_nodes( tid );
+        SparseNodeArray::const_iterator n;
+        for ( n = local_nodes.begin(); n != local_nodes.end(); ++n )
+        {
+          const size_t tnode_id = n->get_node_id();
+
+          // Is the local node in the targets list?
+          if ( targets_->get_lid( tnode_id ) < 0 )
+          {
+            continue;
+          }
+          inner_connect_( tid, rng, n->get_node(), tnode_id );
+        }
+      }
+    }
+    catch ( std::exception& err )
+    {
+      // We must create a new exception here, err's lifetime ends at
+      // the end of the catch block.
+      exceptions_raised_.at( tid ) = std::shared_ptr< WrappedThreadException >( new WrappedThreadException( err ) );
+    }
+  } // of omp parallel
+}
+
+void
+nest::PoissonBuilder::inner_connect_( const int tid, RngPtr rng, Node* target, size_t tnode_id )
+{
+  const size_t target_thread = target->get_thread();
+
+  // check whether the target is on our thread
+  if ( static_cast< size_t >( tid ) != target_thread )
+  {
+    return;
+  }
+
+  poisson_distribution poi_dist;
+
+  // It is not possible to disable multapses with the PoissonBuilder, already checked
+  NodeCollection::const_iterator source_it = sources_->begin();
+  for ( ; source_it < sources_->end(); ++source_it )
+  {
+    const size_t snode_id = ( *source_it ).node_id;
+
+    if ( not allow_autapses_ and snode_id == tnode_id )
+    {
+      continue;
+    }
+
+    // Sample to number of connections that are to be established
+    poisson_distribution::param_type param( pairwise_avg_num_conns_->value( rng, target ) );
+    const size_t num_conns = poi_dist( rng, param );
+
+    for ( size_t n = 0; n < num_conns; ++n )
+    {
+      single_connect_( snode_id, *target, target_thread, rng );
+    }
+  }
+}
+
+
 nest::AuxiliaryBuilder::AuxiliaryBuilder( NodeCollectionPTR sources,
   NodeCollectionPTR targets,
   const DictionaryDatum& conn_spec,

nestkernel/conn_builder.h

@@ -454,6 +454,18 @@ class BernoulliBuilder : public ConnBuilder
   ParameterDatum p_; //!< connection probability
 };
 
+class PoissonBuilder : public ConnBuilder
+{
+public:
+  PoissonBuilder( NodeCollectionPTR, NodeCollectionPTR, const DictionaryDatum&, const std::vector< DictionaryDatum >& );
+
+protected:
+  void connect_() override;
+
+private:
+  void inner_connect_( const int, RngPtr, Node*, size_t );
+  ParameterDatum pairwise_avg_num_conns_; //!< Mean number of connections
+};
 /**
  * Helper class to support parameter handling for tripartite builders.
  *

nestkernel/connection_creator.cpp

@@ -135,6 +135,10 @@ ConnectionCreator::ConnectionCreator( DictionaryDatum dict )
       type_ = Pairwise_bernoulli_on_source;
     }
   }
+  else if ( connection_type == names::pairwise_poisson )
+  {
+    type_ = Pairwise_poisson;
+  }
   else if ( connection_type == names::pairwise_bernoulli_on_target )
   {
 

nestkernel/connection_creator.h

@@ -68,6 +68,7 @@ class ConnectionCreator
   {
     Pairwise_bernoulli_on_source,
     Pairwise_bernoulli_on_target,
+    Pairwise_poisson,
     Fixed_indegree,
     Fixed_outdegree
   };
@@ -144,6 +145,14 @@ class ConnectionCreator
     size_t tgt_thread,
     const Layer< D >& source );
 
+  template < typename Iterator, int D >
+  void connect_to_target_poisson_( Iterator from,
+    Iterator to,
+    Node* tgt_ptr,
+    const Position< D >& tgt_pos,
+    size_t tgt_thread,
+    const Layer< D >& source );
+
   template < int D >
   void pairwise_bernoulli_on_source_( Layer< D >& source,
     NodeCollectionPTR source_nc,
@@ -156,6 +165,10 @@ class ConnectionCreator
     Layer< D >& target,
     NodeCollectionPTR target_nc );
 
+  template < int D >
+  void
+  pairwise_poisson_( Layer< D >& source, NodeCollectionPTR source_nc, Layer< D >& target, NodeCollectionPTR target_nc );
+
   template < int D >
   void
   fixed_indegree_( Layer< D >& source, NodeCollectionPTR source_nc, Layer< D >& target, NodeCollectionPTR target_nc );