Correct errors in slicing of node collections

libnestutil/CMakeLists.txt

@@ -27,6 +27,7 @@ set( nestutil_sources
     lockptr.h
     logging_event.h logging_event.cpp
     logging.h
+    nest_types.h
     numerics.h numerics.cpp
     regula_falsi.h
     sort.h

libnestutil/numerics.cpp

@@ -20,6 +20,11 @@
  *
  */
 
+#include <cstdlib>
+#include <iostream>
+#include <numeric>
+
+#include "nest_types.h"
 #include "numerics.h"
 
 #ifndef HAVE_M_E
@@ -126,3 +131,142 @@ is_integer( double n )
   // factor 4 allows for two bits of rounding error
   return frac_part < 4 * n * std::numeric_limits< double >::epsilon();
 }
+
+long
+mod_inverse( long a, long m )
+{
+  /*
+   The implementation here is based on the extended Euclidean algorithm which
+   solves
+
+     a x + m y = gcd( a, m ) = 1 mod m
+
+   for x and y. Note that the m y term is zero mod m, so the equation is equivalent
+   to
+
+     a x = 1 mod m
+
+   Since we only need x, we can ignore y and use just half of the algorithm.
+
+   We can assume without loss of generality that a < m, because if a = a' + j m with
+   0 < a' < m, we have
+
+      a x mod m = (a' + j m) x mod m = a' x + j x m mod m = a' x.
+
+   This implies that m ≥ 2.
+
+   For details on the algorithm, see D. E. Knuth, The Art of Computer Programming,
+   ch 4.5.2, Algorithm X (vol 2), and ch 1.2.1, Algorithm E (vol 1).
+   */
+
+  assert( 0 < a );
+  assert( 2 <= m );
+
+  const long a_orig = a;
+  const long m_orig = m;
+
+  // If a ≥ m, the algorithm needs two extra rounds to transform this to
+  // a' < m, so we take care of this in a single step here.
+  a = a % m;
+
+  // Use half of extended Euclidean algorithm required to compute inverse
+  long s_0 = 1;
+  long s_1 = 0;
+
+  while ( a > 0 )
+  {
+    // get quotient and remainder in one go
+    const auto res = div( m, a );
+    m = a;
+    a = res.rem;
+
+    // line ordering matters here
+    const long s_0_new = -res.quot * s_0 + s_1;
+    s_1 = s_0;
+    s_0 = s_0_new;
+  }
+
+  // ensure positive result
+  s_1 = ( s_1 + m_orig ) % m_orig;
+
+  assert( m == 1 );                         // gcd() == 1 required
+  assert( ( a_orig * s_1 ) % m_orig == 1 ); // self-test
+
+  return s_1;
+}
+
+size_t
+first_index( long period, long phase0, long step, long phase )
+{
+  assert( period > 0 );
+  assert( step > 0 );
+  assert( 0 <= phase0 and phase0 < period );
+  assert( 0 <= phase and phase < period );
+
+  /*
+   The implementation here is based on
+   https://math.stackexchange.com/questions/25390/how-to-find-the-inverse-modulo-m
+
+   We first need to solve
+
+        phase0 + k step = phase mod period
+   <=>  k step = ( phase - phase0 ) = d_phase mod period
+
+   This has a solution iff d = gcd(step, period) divides d_phase.
+
+   Then, if d = 1, the solution is unique and given by
+
+        k' = mod_inv(step) * d_phase mod period
+
+   If d > 1, we need to divide the equation by it and solve
+
+        (step / d) k_0 = d_phase / d  mod (period / d)
+
+   The set of solutions is then given by
+
+        k_j = k_0 + j * period / d  for j = 0, 1, ..., d-1
+
+   and we need the smallest of these. Now we are interested in
+   an index given by k * step with a period of lcm(step, period)
+   (the outer_period below, marked by | in the illustration in
+   the doxygen comment), we immediately run over
+
+        k_j * step = k_0 * step + j * step * period / d mod lcm(step, period)
+
+   below and take the smallest. But since step * period / d = lcm(step, period),
+   the term in j above vanishes and k_0 * step mod lcm(step, period) is actually
+   the solution.
+
+   We do all calculations in signed long since we may encounter negative
+   values during the algorithm. The result will be non-negative and returned
+   as size_t. This is important because the "not found" case is signaled
+   by invalid_index, which is size_t.
+   */
+
+  // This check is not only a convenience: If step == k * period, we only match if
+  // phase == phase0 and the algorithm below will fail if we did not return here
+  // immediately, because we get d == period -> period_d = 1 and modular inverse
+  // for modulus 1 makes no sense.
+  if ( phase == phase0 )
+  {
+    return 0;
+  }
+
+  const long d_phase = ( phase - phase0 + period ) % period;
+  const long d = std::gcd( step, period );
+
+  if ( d_phase % d != 0 )
+  {
+    return nest::invalid_index; // no solution exists
+  }
+
+  // Scale by GCD, since modular inverse requires gcd==1
+  const long period_d = period / d;
+  const long step_d = step / d;
+  const long d_phase_d = d_phase / d;
+
+  // Compute k_0 and multiply by step, see explanation in introductory comment
+  const long idx = ( d_phase_d * mod_inverse( step_d, period_d ) * step ) % std::lcm( period, step );
+
+  return static_cast< size_t >( idx );
+}

libnestutil/numerics.h

@@ -131,4 +131,52 @@ double dtruncate( double );
  */
 bool is_integer( double );
 
+/**
+ * Returns inverse of integer a modulo m
+ *
+ * For integer a > 0, m ≥ 2, find x so that ( a * x ) mod m = 1.
+ */
+long mod_inverse( long a, long m );
+
+/**
+ * Return first matching index for entry in container with double periodicity.
+ *
+ * Consider
+ * - a container of arbitrary length containing elements (e.g. GIDs) which map to certain values,
+ *     e.g., VP(GID), with a given *period*, e.g., the number of virtual processes
+ * - that the phase (e.g., the VP number) of the first entry in the container is *phase0*
+ * - that we slice the container with a given *step*, causing a double periodicity
+ * - that we want to know the index of the first element in the container with a given *phase*,
+ *     e.g., the first element on a given VP
+ *
+ * If such an index x exists, it is given by
+ *
+ * x = phase0 + k' mod lcm(period, step)
+ * k' = min_k ( phase0 + k step = phase mod period )
+ *
+ * As an example, consider
+ *
+ * idx 0 1 2 3 4 5 6 7 8 9  10 11 | 12 13 14 15 16 17 18
+ * gid 1 2 3 4 5 6 7 8 9 10 11 12 | 13 14 15 16 17 18 19
+ * vp  1 2 3 0 1 2 3 0 1 2  3  0  | 1  2  3  0  1  2  3
+ *    *     *     *     *        | *        *        *
+ *    1     0     3     2        |
+ *
+ * Here, idx is a linear index into the container, gid neuron ids which map to the given vp numbers.
+ * The container is sliced with step=3, i.e., starting with the first element we take every third, as
+ * marked by the asterisks. The question is then at which index we find the first entry belonging to
+ * vp 0, 1, 2, or 3. The numbers in the bottom row show for clarity on which VP we find the respective
+ * asterisks. The | symbol marks where the pattern repeats itself.
+ *
+ * phase0 in this example is 1, i.e., the VP of the first element in the container.
+ *
+ * @note This function returns that index. It is the responsibility of the caller to check if the returned index
+ * is within the bounds of the actual container—the algorithm assumes an infinite container.
+ *
+ * @note The function returns *invalid_index* if no solution exists.
+ *
+ * See comments in the function definition for implementation details.
+ */
+size_t first_index( long period, long phase0, long step, long phase );
+
 #endif

nestkernel/CMakeLists.txt

@@ -44,7 +44,6 @@ set ( nestkernel_sources
       histentry.h histentry.cpp
       model.h model.cpp
       model_manager.h model_manager_impl.h model_manager.cpp
-      nest_types.h
       nest_datums.h nest_datums.cpp
       nest_names.cpp nest_names.h
       nestmodule.h nestmodule.cpp

nestkernel/conn_builder.cpp

@@ -623,7 +623,7 @@ nest::OneToOneBuilder::connect_()
           Node* target = n->get_node();
 
           const size_t tnode_id = n->get_node_id();
-          const long lid = targets_->get_lid( tnode_id );
+          const long lid = targets_->get_nc_index( tnode_id );
           if ( lid < 0 ) // Is local node in target list?
           {
             continue;
@@ -823,7 +823,7 @@ nest::AllToAllBuilder::connect_()
           const size_t tnode_id = n->get_node_id();
 
           // Is the local node in the targets list?
-          if ( targets_->get_lid( tnode_id ) < 0 )
+          if ( targets_->get_nc_index( tnode_id ) < 0 )
           {
             continue;
           }
@@ -1106,7 +1106,7 @@ nest::FixedInDegreeBuilder::connect_()
           const size_t tnode_id = n->get_node_id();
 
           // Is the local node in the targets list?
-          if ( targets_->get_lid( tnode_id ) < 0 )
+          if ( targets_->get_nc_index( tnode_id ) < 0 )
           {
             continue;
           }
@@ -1538,7 +1538,7 @@ nest::BernoulliBuilder::connect_()
           const size_t tnode_id = n->get_node_id();
 
           // Is the local node in the targets list?
-          if ( targets_->get_lid( tnode_id ) < 0 )
+          if ( targets_->get_nc_index( tnode_id ) < 0 )
           {
             continue;
           }
@@ -1654,7 +1654,7 @@ nest::PoissonBuilder::connect_()
           const size_t tnode_id = n->get_node_id();
 
           // Is the local node in the targets list?
-          if ( targets_->get_lid( tnode_id ) < 0 )
+          if ( targets_->get_nc_index( tnode_id ) < 0 )
           {
             continue;
           }
@@ -1847,7 +1847,8 @@ nest::TripartiteBernoulliWithPoolBuilder::connect_()
         }
         else
         {
-          std::copy_n( third_->begin() + get_first_pool_index_( target.lid ), pool_size_, std::back_inserter( pool ) );
+          std::copy_n(
+            third_->begin() + get_first_pool_index_( target.nc_index ), pool_size_, std::back_inserter( pool ) );
         }
 
         // step 3, iterate through indegree to make connections for this target

nestkernel/connection_creator.h

@@ -87,8 +87,6 @@ class ConnectionCreator
    * - "mask": Mask definition (dictionary or masktype).
    * - "kernel": Kernel definition (dictionary, parametertype, or double).
    * - "synapse_model": The synapse model to use.
-   * - "targets": Which targets (model or lid) to select (dictionary).
-   * - "sources": Which targets (model or lid) to select (dictionary).
    * - "weight": Synaptic weight (dictionary, parametertype, or double).
    * - "delay": Synaptic delays (dictionary, parametertype, or double).
    * - other parameters are interpreted as synapse parameters, and may

nestkernel/connection_creator_impl.h

@@ -264,7 +264,7 @@ ConnectionCreator::pairwise_bernoulli_on_source_( Layer< D >& source,
 
         if ( not tgt->is_proxy() )
         {
-          const Position< D > target_pos = target.get_position( ( *tgt_it ).lid );
+          const Position< D > target_pos = target.get_position( ( *tgt_it ).nc_index );
 
           if ( mask_.get() )
           {
@@ -339,7 +339,7 @@ ConnectionCreator::pairwise_bernoulli_on_target_( Layer< D >& source,
     const int thread_id = kernel().vp_manager.get_thread_id();
     try
     {
-      NodeCollection::const_iterator target_begin = target_nc->local_begin();
+      NodeCollection::const_iterator target_begin = target_nc->thread_local_begin();
       NodeCollection::const_iterator target_end = target_nc->end();
 
       for ( NodeCollection::const_iterator tgt_it = target_begin; tgt_it < target_end; ++tgt_it )
@@ -348,7 +348,7 @@ ConnectionCreator::pairwise_bernoulli_on_target_( Layer< D >& source,
 
         assert( not tgt->is_proxy() );
 
-        const Position< D > target_pos = target.get_position( ( *tgt_it ).lid );
+        const Position< D > target_pos = target.get_position( ( *tgt_it ).nc_index );
 
         if ( mask_.get() )
         {
@@ -423,7 +423,7 @@ ConnectionCreator::pairwise_poisson_( Layer< D >& source,
 
         if ( not tgt->is_proxy() )
         {
-          const Position< D > target_pos = target.get_position( ( *tgt_it ).lid );
+          const Position< D > target_pos = target.get_position( ( *tgt_it ).nc_index );
 
           if ( mask_.get() )
           {
@@ -477,7 +477,7 @@ ConnectionCreator::fixed_indegree_( Layer< D >& source,
     throw IllegalConnection( "Spatial Connect with fixed_indegree to devices is not possible." );
   }
 
-  NodeCollection::const_iterator target_begin = target_nc->MPI_local_begin();
+  NodeCollection::const_iterator target_begin = target_nc->rank_local_begin();
   NodeCollection::const_iterator target_end = target_nc->end();
 
   // protect against connecting to devices without proxies
@@ -504,7 +504,7 @@ ConnectionCreator::fixed_indegree_( Layer< D >& source,
 
       size_t target_thread = tgt->get_thread();
       RngPtr rng = get_vp_specific_rng( target_thread );
-      Position< D > target_pos = target.get_position( ( *tgt_it ).lid );
+      Position< D > target_pos = target.get_position( ( *tgt_it ).nc_index );
 
       // We create a source pos vector here that can be updated with the
       // source position. This is done to avoid creating and destroying
@@ -637,7 +637,7 @@ ConnectionCreator::fixed_indegree_( Layer< D >& source,
       Node* const tgt = kernel().node_manager.get_node_or_proxy( target_id );
       size_t target_thread = tgt->get_thread();
       RngPtr rng = get_vp_specific_rng( target_thread );
-      Position< D > target_pos = target.get_position( ( *tgt_it ).lid );
+      Position< D > target_pos = target.get_position( ( *tgt_it ).nc_index );
 
       unsigned long target_number_connections = std::round( number_of_connections_->value( rng, tgt ) );
 
@@ -769,7 +769,7 @@ ConnectionCreator::fixed_outdegree_( Layer< D >& source,
     throw IllegalConnection( "Spatial Connect with fixed_outdegree to devices is not possible." );
   }
 
-  NodeCollection::const_iterator target_begin = target_nc->MPI_local_begin();
+  NodeCollection::const_iterator target_begin = target_nc->rank_local_begin();
   NodeCollection::const_iterator target_end = target_nc->end();
 
   for ( NodeCollection::const_iterator tgt_it = target_begin; tgt_it < target_end; ++tgt_it )