Merge pull request #3132 from heplesser/fix_nc_slicing

Correct errors in slicing of node collections

.github/workflows/nestbuildmatrix.yml

@@ -707,6 +707,8 @@ jobs:
           echo "backend : svg" > $HOME/.matplotlib/matplotlibrc
 
       - name: "Run NEST testsuite"
+        env:
+          DO_TESTS_SKIP_TEST_REQUIRING_MANY_CORES: ${{ contains(matrix.use, 'mpi') && contains(matrix.use, 'openmp') }}
         run: |
           pwd
           cd "$NEST_VPATH"
@@ -718,10 +720,7 @@ jobs:
         if: always()
         with:
           name: "build-logs-${{ matrix.os }}-${{ matrix.cpp_compiler }}-${{ matrix.use }}"
-          path: |
-            install_manifest.txt
-            **.log
-            build/reports/**
+          path: /home/runner/work/nest-simulator/nest-simulator/build/test_report_*/
 
   build_macos:
     if: ${{ !contains(github.event.head_commit.message, 'ci skip') }}
@@ -846,8 +845,5 @@ jobs:
         uses: actions/upload-artifact@5d5d22a31266ced268874388b861e4b58bb5c2f3 # v4.3.1
         if: ${{ always() }}
         with:
-          name: "${{ matrix.NEST_BUILD_TYPE }}-build-logs-${{ matrix.os }}-${{ matrix.cpp_compiler }}"
-          path: |
-            install_manifest.txt
-            **.log
-            build/reports/**
+          name: "build_logs_macos"
+          path: /Users/runner/work/nest-simulator/nest-simulator/build/test_report_*/

lib/sli/nest-init.sli

@@ -488,9 +488,14 @@ def
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-% add conversion from NodeCollection
-/cva [/nodecollectiontype]
-  /cva_g load
+% add new functions to trie if it exists, else create new
+/cva dup lookup not
+{
+  trie
+} if
+  [/connectiontype] /cva_C load addtotrie
+  [/nodecollectiontype] { /all cva_g_l } addtotrie
+  [/nodecollectiontype /literaltype] /cva_g_l load addtotrie
 def
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -1426,8 +1431,6 @@ def
 [/integertype /integertype] /TimeCommunicationAlltoallv_i_i load addtotrie
 def
 
-/cva [/connectiontype] /cva_C load def
-
 /abort
 {
   statusdict /exitcodes get /userabort get

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

models/iaf_psc_delta.cpp

@@ -154,7 +154,7 @@ nest::iaf_psc_delta::Parameters_::set( const DictionaryDatum& d, Node* node )
   }
   if ( c_m_ <= 0 )
   {
-    throw BadProperty( "Capacitance must be >0." );
+    throw BadProperty( "Capacitance must be > 0." );
   }
   if ( t_ref_ < 0 )
   {

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