diff --git a/nestkernel/nestmodule.cpp b/nestkernel/nestmodule.cpp index 77e8d3ae7d..63361bd8f6 100644 --- a/nestkernel/nestmodule.cpp +++ b/nestkernel/nestmodule.cpp @@ -1320,8 +1320,14 @@ NestModule::SetMaxBufferedFunction::execute( SLIInterpreter* i ) const void NestModule::EnableStructuralPlasticity_Function::execute( SLIInterpreter* i ) const { - kernel().sp_manager.enable_structural_plasticity(); + i->assert_stack_load( 3 ); + const double max_distance = getValue< double >( i->OStack.pick( 0 ) ); + const double sigma = getValue< double >( i->OStack.pick( 1 ) ); + const bool use_kernel = getValue< bool >( i->OStack.pick( 2 ) ); + + kernel().sp_manager.enable_structural_plasticity( use_kernel, sigma, max_distance ); + i->OStack.pop( 3 ); i->EStack.pop(); } void diff --git a/nestkernel/sp_manager.cpp b/nestkernel/sp_manager.cpp index 4705b7de64..d31f7e36fe 100644 --- a/nestkernel/sp_manager.cpp +++ b/nestkernel/sp_manager.cpp @@ -24,6 +24,7 @@ // C++ includes: #include +#include // Includes from nestkernel: #include "conn_builder.h" @@ -33,6 +34,7 @@ #include "kernel_manager.h" #include "nest_names.h" #include "sp_manager_impl.h" +#include "spatial.h" namespace nest { @@ -41,6 +43,9 @@ SPManager::SPManager() : ManagerInterface() , structural_plasticity_update_interval_( 10000. ) , structural_plasticity_enabled_( false ) + , structural_plasticity_use_gaussian_kernel_( false ) + , structural_plasticity_gaussian_kernel_sigma_( 1. ) + , structural_plasticity_max_distance_( std::numeric_limits< double >::infinity() ) , sp_conn_builders_() , growthcurve_factories_() , growthcurvedict_( new Dictionary() ) @@ -64,6 +69,8 @@ SPManager::initialize( const bool adjust_number_of_threads_or_rng_only ) structural_plasticity_update_interval_ = 10000.; structural_plasticity_enabled_ = false; + structural_plasticity_use_gaussian_kernel_ = false; + structural_plasticity_gaussian_kernel_sigma_ = 1.; } void @@ -101,7 +108,6 @@ SPManager::get_status( DictionaryDatum& d ) def< std::string >( sp_synapse, names::synapse_model, model ); def< bool >( sp_synapse, names::allow_autapses, ( *i )->allows_autapses() ); def< bool >( sp_synapse, names::allow_multapses, ( *i )->allows_multapses() ); - def< DictionaryDatum >( sp_synapses, ( *i )->get_name(), sp_synapse ); } @@ -168,6 +174,134 @@ SPManager::set_status( const DictionaryDatum& d ) } } +void +SPManager::gather_global_positions_and_ids() +{ + std::vector< double > local_positions; + std::vector< int > local_ids; + std::vector< int > displacements; + + // Collect local positions and IDs + for ( size_t tid = 0; tid < kernel().vp_manager.get_num_threads(); ++tid ) + { + const SparseNodeArray& local_nodes = kernel().node_manager.get_local_nodes( tid ); + + for ( auto node_it = local_nodes.begin(); node_it < local_nodes.end(); ++node_it ) + { + int node_id = node_it->get_node_id(); + if ( node_id < 1 ) + { + throw std::runtime_error( "Invalid neuron ID (must be >= 1)." ); + } + + std::vector< double > pos = get_position( node_id ); + + if ( std::none_of( pos.begin(), pos.end(), []( double v ) { return std::isnan( v ); } ) ) + { + local_ids.push_back( node_id ); + local_positions.insert( local_positions.end(), pos.begin(), pos.end() ); + } + } + } + + // Communicate positions and IDs + kernel().mpi_manager.communicate( local_positions, global_positions, displacements ); + kernel().mpi_manager.communicate( local_ids, global_ids, displacements ); + + // Validate global_positions size consistency with global_ids + size_t num_neurons = global_ids.size(); + size_t total_positions = global_positions.size(); + + if ( num_neurons == 0 ) + { + throw std::runtime_error( + "No neurons with valid positions found. Please provide valid positions, or disable distance dependency." ); + } + if ( total_positions == 0 ) + { + throw std::runtime_error( "No positions found. Please provide positions, or disable distance dependency." ); + } + + if ( total_positions % num_neurons != 0 ) + { + throw std::runtime_error( "Mismatch in global positions dimensionality." ); + } + + pos_dim = total_positions / num_neurons; + + // Pair global_ids with their positions + std::vector< std::pair< int, std::vector< double > > > id_pos_pairs; + id_pos_pairs.reserve( num_neurons ); + for ( size_t i = 0; i < num_neurons; ++i ) + { + int node_id = global_ids[ i ]; + std::vector< double > pos( global_positions.begin() + i * pos_dim, global_positions.begin() + ( i + 1 ) * pos_dim ); + id_pos_pairs.emplace_back( node_id, pos ); + } + + // Sort id_pos_pairs based on node_id to ensure ordering from 1 to num_neurons + std::sort( id_pos_pairs.begin(), + id_pos_pairs.end(), + []( const std::pair< int, std::vector< double > >& a, const std::pair< int, std::vector< double > >& b ) -> bool + { return a.first < b.first; } ); + + // Verify that IDs are sequential + for ( size_t i = 0; i < num_neurons; ++i ) + { + if ( id_pos_pairs[ i ].first != static_cast< int >( i + 1 ) ) + { + throw std::runtime_error( "Neuron IDs are not sequential after sorting." ); + } + } + + // Assign sorted positions to temp_positions + std::vector< double > temp_positions( num_neurons * pos_dim, 0.0 ); + for ( size_t i = 0; i < num_neurons; ++i ) + { + std::copy( id_pos_pairs[ i ].second.begin(), id_pos_pairs[ i ].second.end(), temp_positions.begin() + i * pos_dim ); + } + + // Update global_positions with sorted positions + global_positions = std::move( temp_positions ); +} + + +// Method to perform roulette wheel selection +int +SPManager::roulette_wheel_selection( const std::vector< double >& probabilities, double rnd ) +{ + if ( probabilities.empty() ) + { + throw std::runtime_error( "Probabilities vector is empty." ); + } + + std::vector< double > cumulative( probabilities.size() ); + std::partial_sum( probabilities.begin(), probabilities.end(), cumulative.begin() ); + + // Ensure the sum of probabilities is greater than zero + double sum = cumulative.back(); + if ( sum < 0.0 ) + { + throw std::runtime_error( "Sum of probabilities must be greater than zero." ); + } + + + // Generate a random number in the range [0, sum) + double randomValue = rnd * sum; + + // Perform binary search to find the selected index + auto it = std::lower_bound( cumulative.begin(), cumulative.end(), randomValue ); + return static_cast< int >( std::distance( cumulative.begin(), it ) ); +} + + +double +SPManager::gaussian_kernel( const std::vector< double >& pos1, const std::vector< double >& pos2, const double sigma ) +{ + const double d2 = squared_distance( pos1, pos2 ); + return std::exp( -d2 / ( sigma * sigma ) ); +} + long SPManager::builder_min_delay() const { @@ -409,26 +543,47 @@ SPManager::create_synapses( std::vector< size_t >& pre_id, serialize_id( pre_id, pre_n, pre_id_rnd ); serialize_id( post_id, post_n, post_id_rnd ); - // Shuffle only the largest vector - if ( pre_id_rnd.size() > post_id_rnd.size() ) + std::vector< size_t > pre_ids_results; + std::vector< size_t > post_ids_results; + + if ( !structural_plasticity_use_gaussian_kernel_ ) { - // we only shuffle the n first items, - // where n is the number of postsynaptic elements - global_shuffle( pre_id_rnd, post_id_rnd.size() ); - pre_id_rnd.resize( post_id_rnd.size() ); + // Shuffle only the largest vector + if ( pre_id_rnd.size() > post_id_rnd.size() ) + { + // we only shuffle the n first items, + // where n is the number of postsynaptic elements + global_shuffle( pre_id_rnd, post_id_rnd.size() ); + pre_id_rnd.resize( post_id_rnd.size() ); + } + else + { + // we only shuffle the n first items, + // where n is the number of pre synaptic elements + global_shuffle( post_id_rnd, pre_id_rnd.size() ); + post_id_rnd.resize( pre_id_rnd.size() ); + } + + pre_ids_results = pre_id_rnd; + post_ids_results = post_id_rnd; } else { - // we only shuffle the n first items, - // where n is the number of pre synaptic elements - global_shuffle( post_id_rnd, pre_id_rnd.size() ); - post_id_rnd.resize( pre_id_rnd.size() ); + + // Shuffle pre_ids (Fisher-Yates) to randomize postsynaptic assignment order + for ( size_t i = pre_id_rnd.size() - 1; i > 0; --i ) + { + size_t j = get_rank_synced_rng()->ulrand( i + 1 ); + std::swap( pre_id_rnd[ i ], pre_id_rnd[ j ] ); + } + global_shuffle_spatial( + pre_id_rnd, post_id_rnd, pre_ids_results, post_ids_results, sp_conn_builder->allows_autapses() ); } // create synapse - sp_conn_builder->sp_connect( pre_id_rnd, post_id_rnd ); + sp_conn_builder->sp_connect( pre_ids_results, post_ids_results ); - return not pre_id_rnd.empty(); + return not pre_ids_results.empty(); } void @@ -663,10 +818,114 @@ nest::SPManager::global_shuffle( std::vector< size_t >& v, size_t n ) } v = v2; } +void +SPManager::global_shuffle_spatial( std::vector< size_t >& pre_ids, + std::vector< size_t >& post_ids, + std::vector< size_t >& pre_ids_results, + std::vector< size_t >& post_ids_results, + bool allow_autapse ) +{ + size_t maxIterations = std::min( pre_ids.size(), post_ids.size() ); + + for ( size_t iteration = 0; iteration < maxIterations; ++iteration ) + { + if ( pre_ids.empty() || post_ids.empty() ) + { + break; // Stop if either vector is empty + } + + size_t pre_id = pre_ids.back(); + pre_ids.pop_back(); + + std::vector< double > pre_pos( + global_positions.begin() + ( pre_id - 1 ) * pos_dim, global_positions.begin() + pre_id * pos_dim ); + + std::vector< double > probabilities; + std::vector< size_t > valid_post_ids; + double rnd; + for ( size_t post_id : post_ids ) + { + if ( post_id == pre_id && !allow_autapse ) + { + continue; // Skip self-connections + } + + // fetch post position + std::vector< double > post_pos( + global_positions.begin() + ( post_id - 1 ) * pos_dim, global_positions.begin() + post_id * pos_dim ); + + if ( !within_max_distance( pre_pos, post_pos ) ) + { + continue; // distance > max_distance -> masked out entirely + } + + double prob; + prob = gaussian_kernel( pre_pos, post_pos, structural_plasticity_gaussian_kernel_sigma_ ); + + if ( prob > 0 ) + { + probabilities.push_back( prob ); + valid_post_ids.push_back( post_id ); + } + } + + if ( probabilities.empty() ) + { + continue; // Skip if no valid connections are found + } + + rnd = get_rank_synced_rng()->drand(); + + // Select a post-synaptic neuron using roulette wheel selection + int selected_post_idx = roulette_wheel_selection( probabilities, rnd ); + size_t selected_post_id = valid_post_ids[ selected_post_idx ]; + + // Remove the selected post-synaptic neuron from the list + auto post_it = std::find( post_ids.begin(), post_ids.end(), selected_post_id ); + if ( post_it != post_ids.end() ) + { + post_ids.erase( post_it ); + } + + pre_ids_results.push_back( pre_id ); + post_ids_results.push_back( selected_post_id ); + } +} + +double +SPManager::squared_distance( const std::vector< double >& a, const std::vector< double >& b ) const +{ + if ( a.size() != b.size() ) + { + throw std::runtime_error( "Position vectors must have the same dimension." ); + } + + double d2 = 0.0; + for ( std::size_t i = 0; i < a.size(); ++i ) + { + const double diff = b[ i ] - a[ i ]; + d2 += diff * diff; + } + return d2; +} + +bool +SPManager::within_max_distance( const std::vector< double >& a, const std::vector< double >& b ) const +{ + if ( !std::isfinite( structural_plasticity_max_distance_ ) ) + { + return true; + } + const double r2 = structural_plasticity_max_distance_ * structural_plasticity_max_distance_; + + return squared_distance( a, b ) <= r2; +} void -nest::SPManager::enable_structural_plasticity() +nest::SPManager::enable_structural_plasticity( bool use_gaussian_kernel, + double gaussian_kernel_sigma, + double max_distance ) { if ( kernel().vp_manager.get_num_threads() > 1 ) { @@ -684,7 +943,15 @@ nest::SPManager::enable_structural_plasticity() "Structural plasticity can not be enabled if use_compressed_spikes " "has been set to false." ); } + structural_plasticity_use_gaussian_kernel_ = use_gaussian_kernel; + structural_plasticity_gaussian_kernel_sigma_ = gaussian_kernel_sigma; structural_plasticity_enabled_ = true; + structural_plasticity_max_distance_ = max_distance; + + if ( use_gaussian_kernel ) + { + gather_global_positions_and_ids(); + } } void @@ -693,4 +960,4 @@ nest::SPManager::disable_structural_plasticity() structural_plasticity_enabled_ = false; } -} // namespace nest +} // namespace nest \ No newline at end of file diff --git a/nestkernel/sp_manager.h b/nestkernel/sp_manager.h index 79cc92a72a..be03d5e264 100644 --- a/nestkernel/sp_manager.h +++ b/nestkernel/sp_manager.h @@ -24,6 +24,7 @@ #define SP_MANAGER_H // C++ includes: +#include #include // Includes from libnestutil: @@ -126,7 +127,9 @@ class SPManager : public ManagerInterface /** * Enable structural plasticity */ - void enable_structural_plasticity(); + void enable_structural_plasticity( bool use_gaussian_kernel, + double gaussian_kernel_sigma, + double max_distance = std::numeric_limits< double >::infinity() ); /** * Disable structural plasticity @@ -135,8 +138,16 @@ class SPManager : public ManagerInterface bool is_structural_plasticity_enabled() const; + /** Squared distance. Uses current pos_dim. */ + double squared_distance( const std::vector< double >& pos1, const std::vector< double >& pos2 ) const; + + /** Hard cutoff check using current max distance and metric. */ + bool within_max_distance( const std::vector< double >& pos1, const std::vector< double >& pos2 ) const; + double get_structural_plasticity_update_interval() const; + double get_structural_plasticity_gaussian_kernel_sigma() const; + /** * Returns the minimum delay of all SP builders. * @@ -188,6 +199,55 @@ class SPManager : public ManagerInterface void global_shuffle( std::vector< size_t >& v ); void global_shuffle( std::vector< size_t >& v, size_t n ); + /** + * Compute the Gaussian kernel value between two positions. + * + * @param pos1 Position of the first neuron. + * @param pos2 Position of the second neuron. + * @param sigma Standard deviation for the Gaussian kernel. + * @return Gaussian kernel value. + */ + double gaussian_kernel( const std::vector< double >& pos1, const std::vector< double >& pos2, const double sigma ); + + /** + * Perform global shuffling of pre- and post-synaptic neurons based on spatial probabilities. + * + * @param pre_ids Vector of pre-synaptic neuron IDs. + * @param post_ids Vector of post-synaptic neuron IDs. + * @param pre_ids_results Vector to store shuffled pre-synaptic IDs. + * @param post_ids_results Vector to store shuffled post-synaptic IDs. + */ + void global_shuffle_spatial( std::vector< size_t >& pre_ids, + std::vector< size_t >& post_ids, + std::vector< size_t >& pre_ids_results, + std::vector< size_t >& post_ids_results, + bool allow_autapse ); + + /** + * Gather global neuron positions and IDs from all nodes. + */ + void gather_global_positions_and_ids(); + + /** + * Perform roulette wheel selection to randomly select an index based on probabilities. + * + * @param probabilities Vector of probabilities for selection. + * @param rnd Random number. + * @return Selected index. + */ + int roulette_wheel_selection( const std::vector< double >& probabilities, double rnd ); + + /** + * Global list of neuron IDs used for structural plasticity computations. + */ + std::vector< int > global_ids; + + /** + * Global list of neuron positions used for spatial computations in + * structural plasticity. + */ + std::vector< double > global_positions; + private: /** * Time interval for structural plasticity update (creation/deletion of @@ -200,6 +260,35 @@ class SPManager : public ManagerInterface * Off (False). */ bool structural_plasticity_enabled_; + + /** + * Flag indicating whether a Gaussian spatial kernel is used for connection + * probability computation. + */ + bool structural_plasticity_use_gaussian_kernel_; + + /** + * Standard deviation parameter for the Gaussian kernel used in + * spatial probability calculations. + */ + double structural_plasticity_gaussian_kernel_sigma_; + + /** + * List of precomputed probabilities for neuron connections, indexed + * by neuron pair indices for efficient lookup. + */ + std::vector< double > probability_list; + + /** + * Maximum allowed Euclidean distance between pre- and post-neurons. + */ + double structural_plasticity_max_distance_; + + /** + * Dimensionality of the neuron positions + */ + int pos_dim; + std::vector< SPBuilder* > sp_conn_builders_; /** @@ -228,7 +317,11 @@ SPManager::get_structural_plasticity_update_interval() const { return structural_plasticity_update_interval_; } - +inline double +SPManager::get_structural_plasticity_gaussian_kernel_sigma() const +{ + return structural_plasticity_gaussian_kernel_sigma_; +} } // namespace nest -#endif /* #ifndef SP_MANAGER_H */ +#endif /* #ifndef SP_MANAGER_H */ \ No newline at end of file diff --git a/pynest/nest/lib/hl_api_simulation.py b/pynest/nest/lib/hl_api_simulation.py index 9ce0d6f569..f3d30b2588 100644 --- a/pynest/nest/lib/hl_api_simulation.py +++ b/pynest/nest/lib/hl_api_simulation.py @@ -332,14 +332,22 @@ def Install(module_name): @check_stack -def EnableStructuralPlasticity(): +def EnableStructuralPlasticity( + use_gaussian_kernel=False, + gaussian_kernel_sigma=0.0, + max_distance=float("inf"), +): """Enable structural plasticity for the network simulation See Also -------- DisableStructuralPlasticity + """ + sps(bool(use_gaussian_kernel)) + sps(float(gaussian_kernel_sigma)) + sps(float(max_distance)) sr("EnableStructuralPlasticity") diff --git a/testsuite/cpptests/run_all.cpp b/testsuite/cpptests/run_all.cpp index fcc0868731..99b9390f2c 100644 --- a/testsuite/cpptests/run_all.cpp +++ b/testsuite/cpptests/run_all.cpp @@ -26,6 +26,7 @@ // Includes from cpptests #include "test_block_vector.h" +#include "test_distance_dependent_structural_plasticity.h" #include "test_enum_bitfield.h" #include "test_parameter.h" #include "test_sort.h" diff --git a/testsuite/cpptests/test_distance_dependent_structural_plasticity.h b/testsuite/cpptests/test_distance_dependent_structural_plasticity.h new file mode 100644 index 0000000000..67bcce7f1a --- /dev/null +++ b/testsuite/cpptests/test_distance_dependent_structural_plasticity.h @@ -0,0 +1,74 @@ +/* + * test_distance_dependent_structural_plasticity.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 TEST_DISTANCE_DEPENDENT_H +#define TEST_DISTANCE_DEPENDENT_H + +#define BOOST_TEST_DYN_LINK +#include + +// C++ includes: +#include +#include +#include + +// Includes from nestkernel: +#include "../nestkernel/random_manager.h" +#include "../nestkernel/sp_manager.h" + +namespace nest +{ + +/** + * Test cases: Distance-dependent connection methods in SPManager + */ +BOOST_AUTO_TEST_SUITE( test_distance_dependent ) + +BOOST_AUTO_TEST_CASE( test_gaussianKernel ) +{ + SPManager sp_manager; + + // Test for zero distance + std::vector< double > pos1 = { 0.0, 0.0 }; + std::vector< double > pos2 = { 0.0, 0.0 }; + double sigma = 1.0; + + double expected = 1.0; + BOOST_REQUIRE_CLOSE( sp_manager.gaussian_kernel( pos1, pos2, sigma ), expected, 1e-6 ); + + // Test for unit distance + pos2 = { 1.0, 0.0 }; + expected = std::exp( -1.0 ); + BOOST_REQUIRE_CLOSE( sp_manager.gaussian_kernel( pos1, pos2, sigma ), expected, 1e-6 ); + + // Test for negative sigma (will compute as if sigma were positive) + sigma = -1.0; + double result = sp_manager.gaussian_kernel( pos1, pos2, sigma ); + expected = std::exp( -1.0 ); // Same as sigma=1 since squared value is used + BOOST_REQUIRE_CLOSE( result, expected, 1e-6 ); +} + +BOOST_AUTO_TEST_SUITE_END() + +} // namespace nest + +#endif /* TEST_DISTANCE_DEPENDENT_H */ diff --git a/testsuite/pytests/test_sp/test_sp_manager.py b/testsuite/pytests/test_sp/test_sp_manager.py index 4631b613ab..271b088ae6 100644 --- a/testsuite/pytests/test_sp/test_sp_manager.py +++ b/testsuite/pytests/test_sp/test_sp_manager.py @@ -49,12 +49,18 @@ def setUp(self): "rate_connection_delayed_lbl", ] + HAVE_GSL = nest.ll_api.sli_func("statusdict/have_gsl ::") + def test_register_synapses(self): for syn_model in nest.synapse_models: if syn_model not in self.exclude_synapse_model: nest.ResetKernel() nest.SetDefaults(syn_model, {"delay": 0.5}) - syn_dict = {"synapse_model": syn_model, "pre_synaptic_element": "SE1", "post_synaptic_element": "SE2"} + syn_dict = { + "synapse_model": syn_model, + "pre_synaptic_element": "SE1", + "post_synaptic_element": "SE2", + } # For co-dependent properties, we use `set()` instead of kernel attributes nest.set(min_delay=0.1, max_delay=1.0) nest.structural_plasticity_synapses = {"syn1": syn_dict} @@ -121,7 +127,11 @@ def test_synapse_creation(self): for syn_model in nest.synapse_models: if syn_model not in self.exclude_synapse_model: nest.ResetKernel() - syn_dict = {"synapse_model": syn_model, "pre_synaptic_element": "SE1", "post_synaptic_element": "SE2"} + syn_dict = { + "synapse_model": syn_model, + "pre_synaptic_element": "SE1", + "post_synaptic_element": "SE2", + } nest.structural_plasticity_synapses = {"syn1": syn_dict} neurons = nest.Create( "iaf_psc_alpha", diff --git a/testsuite/pytests/test_sp/test_sp_manager_spatial.py b/testsuite/pytests/test_sp/test_sp_manager_spatial.py new file mode 100644 index 0000000000..c91c103499 --- /dev/null +++ b/testsuite/pytests/test_sp/test_sp_manager_spatial.py @@ -0,0 +1,172 @@ +# -*- coding: utf-8 -*- +# +# test_sp_manager_spatial.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 . + +import unittest + +import nest +import numpy as np + + +class TestStructuralPlasticityManagerSpatial(unittest.TestCase): + def setUp(self): + nest.ResetKernel() + nest.set_verbosity("M_INFO") + self.exclude_synapse_model = [ + "stdp_dopamine_synapse", + "stdp_dopamine_synapse_lbl", + "stdp_dopamine_synapse_hpc", + "stdp_dopamine_synapse_hpc_lbl", + "gap_junction", + "gap_junction_lbl", + "diffusion_connection", + "diffusion_connection_lbl", + "rate_connection_instantaneous", + "rate_connection_instantaneous_lbl", + "rate_connection_delayed", + "rate_connection_delayed_lbl", + ] + + HAVE_GSL = nest.ll_api.sli_func("statusdict/have_gsl ::") + + def test_synapse_creation_distance_dependent(self): + for syn_model in nest.synapse_models: + if syn_model not in self.exclude_synapse_model: + nest.ResetKernel() + # assign random 2d positions + positions = np.random.uniform([0, 0], [10, 10], size=(2, 2)) + + syn_dict = { + "synapse_model": syn_model, + "pre_synaptic_element": "SE1", + "post_synaptic_element": "SE2", + } + nest.structural_plasticity_synapses = {"syn1": syn_dict} + # create 2 neurons that with 10 presynaptic and 10 postsymaptic elements each. + neurons = nest.Create( + "iaf_psc_alpha", + 2, + { + "synaptic_elements": { + "SE1": {"z": 10.0, "growth_rate": 0.0}, + "SE2": {"z": 10.0, "growth_rate": 0.0}, + } + }, + positions=nest.spatial.free(pos=positions.tolist()), + ) + # set update interval to 10 seconds + nest.structural_plasticity_update_interval = 10000 + nest.EnableStructuralPlasticity(use_gaussian_kernel=True, gaussian_kernel_sigma=1.0) + nest.Simulate(10.0) + status = nest.GetStatus(neurons, "synaptic_elements") + for st_neuron in status: + assert st_neuron["SE1"]["z_connected"] == 10 + assert st_neuron["SE2"]["z_connected"] == 10 + + assert len(nest.GetConnections(neurons, neurons, syn_model)) == 20 + break + + def test_structural_plasticity_with_positions(self): + # Test structural plasticity when neurons have positions in multiple dimensions. + nest.ResetKernel() + nest.structural_plasticity_update_interval = 1.0 + + syn_model = "static_synapse" + nest.CopyModel(syn_model, "sp_synapse") + nest.SetDefaults("sp_synapse", {"weight": 1.0, "delay": 1.0}) + + nest.structural_plasticity_synapses = { + "sp_syn": { + "synapse_model": "sp_synapse", + "post_synaptic_element": "Den_ex", + "pre_synaptic_element": "Axon_ex", + } + } + + num_neurons = 30 + positions = np.random.uniform(0, 10, (num_neurons, 3)) # 3D positions + neuron_params = { + "synaptic_elements": { + "Axon_ex": {"z": 1.0, "growth_rate": 1.0}, + "Den_ex": {"z": 1.0, "growth_rate": 1.0}, + }, + } + + neurons = nest.Create( + "iaf_psc_alpha", + num_neurons, + neuron_params, + positions=nest.spatial.free(pos=positions.tolist()), + ) + + nest.EnableStructuralPlasticity(use_gaussian_kernel=True, gaussian_kernel_sigma=1.0) + nest.Simulate(10.0) + + connections = nest.GetConnections(neurons, neurons, synapse_model="sp_synapse") + self.assertGreater(len(connections), 0, "No connections were created") + + # Check that neurons closer together are more likely to be connected + distances = [] + for conn in connections: + pre_pos = positions[conn.source - 1] + post_pos = positions[conn.target - 1] + distance = np.linalg.norm(np.array(pre_pos) - np.array(post_pos)) + distances.append(distance) + + avg_distance = np.mean(distances) + self.assertLess(avg_distance, 10.0, "Average connection distance is too large") + + def test_distance_dependent_without_positions(self): + # Test if an error is raised when distance dependency is enabled without providing positions. + nest.ResetKernel() + + syn_model = "static_synapse" + nest.CopyModel(syn_model, "sp_synapse") + nest.SetDefaults("sp_synapse", {"weight": 1.0, "delay": 1.0}) + + nest.structural_plasticity_synapses = { + "sp_syn": { + "synapse_model": "sp_synapse", + "post_synaptic_element": "Den_ex", + "pre_synaptic_element": "Axon_ex", + } + } + + num_neurons = 10 + neuron_params = { + "synaptic_elements": { + "Axon_ex": {"z": 1.0, "growth_rate": 1.0}, + "Den_ex": {"z": 1.0, "growth_rate": 1.0}, + }, + } + + # Create neurons without positions + neurons = nest.Create("iaf_psc_alpha", num_neurons, params=neuron_params) + + # Expecting a NEST CppException due to missing positions + with self.assertRaises(nest.NESTErrors.CppException) as context: + nest.EnableStructuralPlasticity(use_gaussian_kernel=True, gaussian_kernel_sigma=1.0) + nest.Simulate(10.0) + + # Check the error message + self.assertIn( + "No neurons with valid positions found. Please provide valid positions, or disable distance dependency.", + str(context.exception), + )