Connection rule for direct support of clusterd connectivity

nestkernel/conn_builder.cpp

@@ -445,6 +445,72 @@
   }
 }
 
+void
+nest::BipartiteConnBuilder::clustered_single_connect_( size_t snode_id,
+  Node& target,
+  size_t target_thread,
+  RngPtr rng,
+  const bool is_intra )
+{
+  // TODO Almost verbatim copy of single_connect_, just if instead of loop.
+  if ( this->requires_proxies() and not target.has_proxies() )
+  {
+    throw IllegalConnection( "Cannot use this rule to connect to nodes without proxies (usually devices)." );
+  }
+
+  assert( synapse_model_id_.size() == 2 );
+
+  const size_t synapse_indx = is_intra ? 0 : 1;
+
+  update_param_dict_( snode_id, target, target_thread, rng, synapse_indx );
+
+  if ( default_weight_and_delay_[ synapse_indx ] )
+  {
+    kernel().connection_manager.connect( snode_id,
+      &target,
+      target_thread,
+      synapse_model_id_[ synapse_indx ],
+      param_dicts_[ synapse_indx ][ target_thread ] );
+  }
+  else if ( default_weight_[ synapse_indx ] )
+  {
+    kernel().connection_manager.connect( snode_id,
+      &target,
+      target_thread,
+      synapse_model_id_[ synapse_indx ],
+      param_dicts_[ synapse_indx ][ target_thread ],
+      delays_[ synapse_indx ]->value_double( target_thread, rng, snode_id, &target ) );
+  }
+  else if ( default_delay_[ synapse_indx ] )
+  {
+    kernel().connection_manager.connect( snode_id,
+      &target,
+      target_thread,
+      synapse_model_id_[ synapse_indx ],
+      param_dicts_[ synapse_indx ][ target_thread ],
+      numerics::nan,
+      weights_[ synapse_indx ]->value_double( target_thread, rng, snode_id, &target ) );
+  }
+  else
+  {
+    const double delay = delays_[ synapse_indx ]->value_double( target_thread, rng, snode_id, &target );
+    const double weight = weights_[ synapse_indx ]->value_double( target_thread, rng, snode_id, &target );
+    kernel().connection_manager.connect( snode_id,
+      &target,
+      target_thread,
+      synapse_model_id_[ synapse_indx ],
+      param_dicts_[ synapse_indx ][ target_thread ],
+      delay,
+      weight );
+  }
+
+  // We connect third-party only once per source-target pair, not per collocated synapse type
+  if ( third_out_ )
+  {
+    third_out_->third_connect( snode_id, target );
+  }
+}
+
 void
 nest::BipartiteConnBuilder::set_synaptic_element_names( const std::string& pre_name, const std::string& post_name )
 {
@@ -1865,6 +1931,210 @@
 }
 
 
+nest::ClusteredFixedTotalNumberBuilder::ClusteredFixedTotalNumberBuilder( NodeCollectionPTR sources,
+  NodeCollectionPTR targets,
+  ThirdOutBuilder* third_out,
+  const DictionaryDatum& conn_spec,
+  const std::vector< DictionaryDatum >& syn_specs )
+  : BipartiteConnBuilder( sources, targets, third_out, conn_spec, syn_specs )
+  , N_( ( *conn_spec )[ names::N ] )
+  , num_clusters_( ( *conn_spec )[ "num_clusters" ] )
+{
+  // check for potential errors
+
+  // verify that total number of connections is not larger than
+  // N_sources*N_targets
+  if ( not allow_multapses_ )
+  {
+    if ( ( N_ > static_cast< long >( sources_->size() * targets_->size() ) ) )
+    {
+      throw BadProperty( "Total number of connections cannot exceed product of source and target population sizes." );
+    }
+  }
+
+  if ( N_ < 0 )
+  {
+    throw BadProperty( "Total number of connections cannot be negative." );
+  }
+
+  if ( num_clusters_ < 1 )
+  {
+    throw BadProperty( "There must be at least one cluster." );
+  }
+
+  if ( synapse_model_id_.size() != 2 )
+  {
+    throw BadProperty( "For clustered connectivity, syn_specs must be CollocatedSynapse with two elements." );
+  }
+
+  // for now multapses cannot be forbidden
+  // TODO: Implement option for multapses_ = False, where already existing
+  // connections are stored in
+  // a bitmap
+  if ( not allow_multapses_ )
+  {
+    throw NotImplemented( "Connect doesn't support the suppression of multapses in the FixedTotalNumber connector." );
+  }
+}
+
+void
+nest::ClusteredFixedTotalNumberBuilder::connect_()
+{
+  const size_t M = kernel().vp_manager.get_num_virtual_processes();
+  const long total_size_sources = sources_->size();
+  const long total_size_targets = targets_->size();
+  const double conn_prob = static_cast< double >( N_ ) / ( total_size_sources * total_size_targets );
+
+  // get global rng that is tested for synchronization for all threads
+  RngPtr grng = get_rank_synced_rng();
+  std::vector< long > num_conns_on_vp( M, 0 ); // corresponds to n[] in GSL 1.8 binomial algo
+  binomial_distribution bino_dist;
+
+  for ( size_t tc = 0; tc < num_clusters_; ++tc )
+  {
+    const auto tgt_nrns = targets_->slice( tc, total_size_targets, num_clusters_ );
+    const long size_targets = tgt_nrns->size();
+
+    // Compute the distribution of targets over processes using the modulo
+    // function
+    std::vector< size_t > number_of_targets_on_vp( M, 0 );
+    std::vector< size_t > local_targets;
+    local_targets.reserve( size_targets / kernel().mpi_manager.get_num_processes() );
+    for ( size_t t = 0; t < tgt_nrns->size(); t++ )
+    {
+      size_t vp = kernel().vp_manager.node_id_to_vp( ( *tgt_nrns )[ t ] );
+      ++number_of_targets_on_vp[ vp ];
+      if ( kernel().vp_manager.is_local_vp( vp ) )
+      {
+        local_targets.push_back( ( *tgt_nrns )[ t ] );
+      }
+    }
+
+    for ( size_t sc = 0; sc < num_clusters_; ++sc )
+    {
+      const auto src_nrns = sources_->slice( sc, total_size_sources, num_clusters_ );
+      const long size_sources = src_nrns->size();
+
+
+      const long num_conns = conn_prob * size_sources * size_targets;
+
+      // drawing connection ids
+
+
+      // We use the multinomial distribution to determine the number of
+      // connections that will be made on one virtual process, i.e. we
+      // partition the set of edges into n_vps subsets. The number of
+      // edges on one virtual process is binomially distributed with
+      // the boundary condition that the sum of all edges over virtual
+      // processes is the total number of edges.
+      // To obtain the num_conns_on_vp we adapt the gsl
+      // implementation of the multinomial distribution.
+
+      // K from gsl is equivalent to M = n_vps
+      // N is already taken from stack
+      // p[] is targets_on_vp
+      std::fill( num_conns_on_vp.begin(), num_conns_on_vp.end(), 0 ); // corresponds to n[], reset to 0 for each use
+
+      // calculate exact multinomial distribution
+
+      // begin code adapted from gsl 1.8 //
+      double sum_dist = 0.0; // corresponds to sum_p
+      // norm is equivalent to size_targets
+      unsigned int sum_partitions = 0; // corresponds to sum_n
+
+      for ( int k = 0; k < M; k++ )
+      {
+        // If we have distributed all connections on the previous processes we exit the loop. It is important to
+        // have this check here, as N - sum_partition is set as n value for GSL, and this must be larger than 0.
+        if ( num_conns == sum_partitions )
+        {
+          break;
+        }
+        if ( number_of_targets_on_vp[ k ] > 0 )
+        {
+          double num_local_targets = static_cast< double >( number_of_targets_on_vp[ k ] );
+          double p_local = num_local_targets / ( size_targets - sum_dist );
+
+          binomial_distribution::param_type param( num_conns - sum_partitions, p_local );
+          num_conns_on_vp[ k ] = bino_dist( grng, param );
+        }
+
+        sum_dist += static_cast< double >( number_of_targets_on_vp[ k ] );
+        sum_partitions += static_cast< unsigned int >( num_conns_on_vp[ k ] );
+      }
+
+      // end code adapted from gsl 1.8
+
+#pragma omp parallel
+      {
+        // get thread id
+        const size_t tid = kernel().vp_manager.get_thread_id();
+
+        try
+        {
+          const size_t vp_id = kernel().vp_manager.thread_to_vp( tid );
+
+          if ( kernel().vp_manager.is_local_vp( vp_id ) )
+          {
+            RngPtr rng = get_vp_specific_rng( tid );
+
+            // gather local target node IDs
+            std::vector< size_t > thread_local_targets;
+            thread_local_targets.reserve( number_of_targets_on_vp[ vp_id ] );
+
+            std::vector< size_t >::const_iterator tnode_id_it = local_targets.begin();
+            for ( ; tnode_id_it != local_targets.end(); ++tnode_id_it )
+            {
+              if ( kernel().vp_manager.node_id_to_vp( *tnode_id_it ) == vp_id )
+              {
+                thread_local_targets.push_back( *tnode_id_it );
+              }
+            }
+
+            assert( thread_local_targets.size() == number_of_targets_on_vp[ vp_id ] );
+
+            while ( num_conns_on_vp[ vp_id ] > 0 )
+            {
+
+              // draw random numbers for source node from all source neurons
+              const long s_index = rng->ulrand( size_sources );
+              // draw random numbers for target node from
+              // targets_on_vp on this virtual process
+              const long t_index = rng->ulrand( thread_local_targets.size() );
+              // map random number of source node to node ID corresponding to
+              // the source_adr vector
+              const long snode_id = ( *src_nrns )[ s_index ];
+              // map random number of target node to node ID using the
+              // targets_on_vp vector
+              const long tnode_id = thread_local_targets[ t_index ];
+
+              Node* const target = kernel().node_manager.get_node_or_proxy( tnode_id, tid );
+              const size_t target_thread = target->get_thread();
+
+              if ( allow_autapses_ or snode_id != tnode_id )
+              {
+                clustered_single_connect_( snode_id,
+                  *target,
+                  target_thread,
+                  rng,
+                  /* is_intra */ sc == tc );
+                num_conns_on_vp[ vp_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 ) );
+        }
+      } // omp parallel
+    }   // for target
+  }     // for source
+}
+
+
 nest::BernoulliBuilder::BernoulliBuilder( NodeCollectionPTR sources,
   NodeCollectionPTR targets,
   ThirdOutBuilder* third_out,

nestkernel/conn_builder.h

@@ -189,6 +189,10 @@ class BipartiteConnBuilder
   void single_connect_( size_t, Node&, size_t, RngPtr );
   void single_disconnect_( size_t, Node&, size_t );
 
+  //! TODO Experimental
+  void
+  clustered_single_connect_( size_t snode_id, Node& target, size_t target_thread, RngPtr rng, const bool is_intra );
+
   /**
    * Moves pointer in parameter array.
    *
@@ -713,6 +717,23 @@ class FixedTotalNumberBuilder : public BipartiteConnBuilder
   long N_;
 };
 
+class ClusteredFixedTotalNumberBuilder : public BipartiteConnBuilder
+{
+public:
+  ClusteredFixedTotalNumberBuilder( NodeCollectionPTR,
+    NodeCollectionPTR,
+    ThirdOutBuilder* third_out,
+    const DictionaryDatum&,
+    const std::vector< DictionaryDatum >& );
+
+protected:
+  void connect_() override;
+
+private:
+  long N_;
+  long num_clusters_;
+};
+
 class BernoulliBuilder : public BipartiteConnBuilder
 {
 public:

nestkernel/connection_manager.cpp

@@ -111,6 +111,7 @@
     register_conn_builder< PoissonBuilder >( "pairwise_poisson" );
     register_conn_builder< SymmetricBernoulliBuilder >( "symmetric_pairwise_bernoulli" );
     register_conn_builder< FixedTotalNumberBuilder >( "fixed_total_number" );
+    register_conn_builder< ClusteredFixedTotalNumberBuilder >( "clustered_fixed_total_number" );
     register_third_conn_builder< ThirdBernoulliWithPoolBuilder >( "third_factor_bernoulli_with_pool" );
 #ifdef HAVE_LIBNEUROSIM
     register_conn_builder< ConnectionGeneratorBuilder >( "conngen" );
@@ -808,7 +809,7 @@
 }
 
 void
 nest::ConnectionManager::connect_sonata( const DictionaryDatum& graph_specs, const long hyberslab_size )
 {
 #ifdef HAVE_HDF5
   SonataConnector sonata_connector( graph_specs, hyberslab_size );

testsuite/pytests/mpi/2/test_clustered_fixed_total_number.py

@@ -0,0 +1 @@
+../../test_clustered_fixed_total_number.py

testsuite/pytests/test_clustered_fixed_total_number.py

@@ -0,0 +1,64 @@
+# -*- coding: utf-8 -*-
+#
+# test_clustered_fixed_total_number.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 <http://www.gnu.org/licenses/>.
+
+"""
+Test models with calcium concentration.
+
+This set of tests verify the behavior of the calcium concentration in models
+that inherit from the strutural plasticity node class in the kernel.
+"""
+
+import nest
+import numpy as np
+import pytest
+
+
+@pytest.mark.skipif_missing_threads
+@pytest.mark.parametrize("n_threads", [1, 2, 3])
+def test_clustered_fixed_total_number(n_threads):
+    """Minimal test for clustered connectivity."""
+
+    nest.ResetKernel()
+    nest.local_num_threads = n_threads
+
+    a = nest.Create("parrot_neuron", n=50)
+    b = nest.Create("parrot_neuron", n=50)
+
+    num_conns = 1000
+    w_intra = 100
+    w_cross = 3
+
+    nest.Connect(
+        a,
+        b,
+        {"rule": "clustered_fixed_total_number", "N": num_conns, "num_clusters": 2},
+        nest.CollocatedSynapses({"weight": w_intra}, {"weight": w_cross}),
+    )
+
+    assert nest.num_connections * nest.num_processes == num_conns
+
+    # All odd numbered neurons get assigned to one cluster and the even numbered
+    # ones to the other. Thus, the sum of source and target ids for will be even
+    # for intra-cluster connections and odd for cross-cluster connections.
+    c = nest.GetConnections().get(["source", "target", "weight"], output="pandas")
+
+    assert all(c.loc[(c["source"] + c["target"]) % 2 == 0, "weight"] == w_intra)
+    assert all(c.loc[(c["source"] + c["target"]) % 2 == 1, "weight"] == w_cross)