From 8c8c72503c8077ef44ea44934e463cf72ac98b17 Mon Sep 17 00:00:00 2001 From: Hans Ekkehard Plesser Date: Wed, 16 Apr 2025 00:52:38 +0200 Subject: [PATCH 1/3] Hike ubuntu version for isort to 22.04 --- .github/workflows/nestbuildmatrix.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/nestbuildmatrix.yml b/.github/workflows/nestbuildmatrix.yml index d78e75c069..240c30c743 100644 --- a/.github/workflows/nestbuildmatrix.yml +++ b/.github/workflows/nestbuildmatrix.yml @@ -358,7 +358,7 @@ jobs: pylint --jobs=$(nproc) pynest/ testsuite/pytests/*.py testsuite/regressiontests/*.py isort: - runs-on: "ubuntu-20.04" + runs-on: "ubuntu-22.04" steps: - name: "Checkout repository content" uses: actions/checkout@11bd71901bbe5b1630ceea73d27597364c9af683 # v4.2.2 From 4f46e34c36a93201a170b58770d49a99a725041f Mon Sep 17 00:00:00 2001 From: Hans Ekkehard Plesser Date: Thu, 26 Jun 2025 23:26:50 +0200 Subject: [PATCH 2/3] Add clustered fixed total number connection rule (experimental) --- nestkernel/conn_builder.cpp | 266 ++++++++++++++++++ nestkernel/conn_builder.h | 21 ++ nestkernel/connection_manager.cpp | 1 + .../2/test_clustered_fixed_total_number.py | 1 + .../test_clustered_fixed_total_number.py | 64 +++++ 5 files changed, 353 insertions(+) create mode 120000 testsuite/pytests/mpi/2/test_clustered_fixed_total_number.py create mode 100644 testsuite/pytests/test_clustered_fixed_total_number.py diff --git a/nestkernel/conn_builder.cpp b/nestkernel/conn_builder.cpp index f70f00f0b4..7aeda84fb0 100644 --- a/nestkernel/conn_builder.cpp +++ b/nestkernel/conn_builder.cpp @@ -445,6 +445,72 @@ nest::BipartiteConnBuilder::single_connect_( size_t snode_id, Node& target, size } } +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,206 @@ nest::FixedTotalNumberBuilder::connect_() } +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 ); + + 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(); + + 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(); + + const long num_conns = conn_prob * size_sources * size_targets; + + // drawing connection ids + + // 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++ ) + { + int 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 ] ); + } + } + + // 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::vector< long > num_conns_on_vp( M, 0 ); // corresponds to n[] + + // calculate exact multinomial distribution + // get global rng that is tested for synchronization for all threads + RngPtr grng = get_rank_synced_rng(); + + // 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 + + binomial_distribution bino_dist; + 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, diff --git a/nestkernel/conn_builder.h b/nestkernel/conn_builder.h index 801d4bb355..9b81709ac8 100644 --- a/nestkernel/conn_builder.h +++ b/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: diff --git a/nestkernel/connection_manager.cpp b/nestkernel/connection_manager.cpp index 7228915296..7a32b68702 100644 --- a/nestkernel/connection_manager.cpp +++ b/nestkernel/connection_manager.cpp @@ -111,6 +111,7 @@ nest::ConnectionManager::initialize( const bool adjust_number_of_threads_or_rng_ 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" ); diff --git a/testsuite/pytests/mpi/2/test_clustered_fixed_total_number.py b/testsuite/pytests/mpi/2/test_clustered_fixed_total_number.py new file mode 120000 index 0000000000..059dcd13d2 --- /dev/null +++ b/testsuite/pytests/mpi/2/test_clustered_fixed_total_number.py @@ -0,0 +1 @@ +../../test_clustered_fixed_total_number.py \ No newline at end of file diff --git a/testsuite/pytests/test_clustered_fixed_total_number.py b/testsuite/pytests/test_clustered_fixed_total_number.py new file mode 100644 index 0000000000..eadefb63c4 --- /dev/null +++ b/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 . + +""" +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) From 5ab05cfe349bf8e93100ae9cfcd039d5554c2b98 Mon Sep 17 00:00:00 2001 From: Hans Ekkehard Plesser Date: Fri, 27 Jun 2025 08:06:29 +0200 Subject: [PATCH 3/3] Re-organize loop structure --- nestkernel/conn_builder.cpp | 52 ++++++++++++++++++++----------------- 1 file changed, 28 insertions(+), 24 deletions(-) diff --git a/nestkernel/conn_builder.cpp b/nestkernel/conn_builder.cpp index 7aeda84fb0..65c202c018 100644 --- a/nestkernel/conn_builder.cpp +++ b/nestkernel/conn_builder.cpp @@ -1985,34 +1985,41 @@ nest::ClusteredFixedTotalNumberBuilder::connect_() const long total_size_targets = targets_->size(); const double conn_prob = static_cast< double >( N_ ) / ( total_size_sources * total_size_targets ); - for ( size_t sc = 0; sc < num_clusters_; ++sc ) + // 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 src_nrns = sources_->slice( sc, total_size_sources, num_clusters_ ); - const long size_sources = src_nrns->size(); + const auto tgt_nrns = targets_->slice( tc, total_size_targets, num_clusters_ ); + const long size_targets = tgt_nrns->size(); - for ( size_t tc = 0; tc < num_clusters_; ++tc ) + // 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++ ) { - const auto tgt_nrns = targets_->slice( tc, total_size_targets, num_clusters_ ); - const long size_targets = tgt_nrns->size(); + 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 - // 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++ ) - { - int 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 ] ); - } - } // We use the multinomial distribution to determine the number of // connections that will be made on one virtual process, i.e. we @@ -2026,18 +2033,15 @@ nest::ClusteredFixedTotalNumberBuilder::connect_() // K from gsl is equivalent to M = n_vps // N is already taken from stack // p[] is targets_on_vp - std::vector< long > num_conns_on_vp( M, 0 ); // corresponds to n[] + 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 - // get global rng that is tested for synchronization for all threads - RngPtr grng = get_rank_synced_rng(); // 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 - binomial_distribution bino_dist; 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