Feat/sorted spike output (#89)

- Spikes are now sorted across MPI ranks according to time and gid before being
written to file
- Update build system to use c++11 by default and fix wrong directives
- Use nrnmpi wrappers and extend wrappers by one function
- Removed post-process sorting from integration tests
- Updated endianness tests to conform with C++11
- Fix warnings with c++11/clang and update mod2c

Author: ohm314

CMakeLists.txt

@@ -60,6 +60,10 @@ option(ENABLE_CUDA_MODULES "Enable modules like Random123 written in CUDA" ON)
 option(ENABLE_NET_RECEIVE_BUFFERING "Enable event buffering in net_receive function" ON)
 option(ENABLE_OMP_RUNTIME_SCHEDULE "Use runtime schedule for OpenMP" OFF)
 
+## set C++11 standard to be default
+set(CMAKE_CXX_STANDARD 11)
+set(CMAKE_CXX_STANDARD_REQUIRED ON)
+set(CMAKE_CXX_EXTENSIONS OFF)
 
 set(FRONTEND_C_COMPILER gcc CACHE FILEPATH "C compiler for building mod2c [frontend executable]")
 set(FRONTEND_CXX_COMPILER g++ CACHE FILEPATH "C++ compiler for building mod2c [frontend executable]")

coreneuron/nrniv/nrn_stats.cpp

@@ -69,9 +69,11 @@ void report_cell_stats(void) {
     stat_array[5] = (long)spikevec_gid.size();  // number of spikes
 
     int spikevec_positive_gid_size = 0;
-    for (int i = 0; i < spikevec_gid.size(); ++i)
-        if (spikevec_gid[i] > -1)
+    for (std::size_t i = 0; i < spikevec_gid.size(); ++i) {
+        if (spikevec_gid[i] > -1) {
             spikevec_positive_gid_size++;
+        }
+    }
 
     stat_array[6] = (long)spikevec_positive_gid_size;  // number of non-negative gid spikes
 

coreneuron/nrniv/output_spikes.cpp

@@ -31,6 +31,8 @@ THE POSSIBILITY OF SUCH DAMAGE.
 #include <string.h>
 #include <stdexcept>  // std::lenght_error
 #include <vector>
+#include <algorithm>
+#include <numeric>
 #include "coreneuron/nrnconf.h"
 #include "coreneuron/nrniv/nrniv_decl.h"
 #include "coreneuron/nrniv/output_spikes.h"
@@ -65,6 +67,74 @@ void spikevec_unlock() {
 }
 
 #if NRNMPI
+
+void local_spikevec_sort(std::vector<double>& isvect,
+                         std::vector<int>& isvecg,
+                         std::vector<double>& osvect,
+                         std::vector<int>& osvecg) {
+    osvect.resize(isvect.size());
+    osvecg.resize(isvecg.size());
+    // first build a permutation vector
+    std::vector<std::size_t> perm(isvect.size());
+    std::iota(perm.begin(), perm.end(), 0);
+    // sort by gid (second predicate first)
+    std::stable_sort(perm.begin(), perm.end(),
+                     [&](std::size_t i, std::size_t j) { return isvecg[i] < isvecg[j]; });
+    // then sort by time
+    std::stable_sort(perm.begin(), perm.end(),
+                     [&](std::size_t i, std::size_t j) { return isvect[i] < isvect[j]; });
+    // now apply permutation to time and gid output vectors
+    std::transform(perm.begin(), perm.end(), osvect.begin(),
+                   [&](std::size_t i) { return isvect[i]; });
+    std::transform(perm.begin(), perm.end(), osvecg.begin(),
+                   [&](std::size_t i) { return isvecg[i]; });
+}
+
+void sort_spikes(std::vector<double>& spikevec_time, std::vector<int>& spikevec_gid) {
+    double lmin_time = *(std::min_element(spikevec_time.begin(), spikevec_time.end()));
+    double lmax_time = *(std::max_element(spikevec_time.begin(), spikevec_time.end()));
+    double min_time = nrnmpi_dbl_allmin(lmin_time);
+    double max_time = nrnmpi_dbl_allmax(lmax_time);
+
+    // allocate send and receive counts and displacements for MPI_Alltoallv
+    std::vector<int> snd_cnts(nrnmpi_numprocs);
+    std::vector<int> rcv_cnts(nrnmpi_numprocs);
+    std::vector<int> snd_dsps(nrnmpi_numprocs);
+    std::vector<int> rcv_dsps(nrnmpi_numprocs);
+
+    double bin_t = (max_time - min_time) / nrnmpi_numprocs;
+    // first find number of spikes in each time window
+    for (const auto& st : spikevec_time) {
+        int idx = (int)(st - min_time) / bin_t;
+        snd_cnts[idx]++;
+    }
+    for (int i = 1; i < nrnmpi_numprocs; i++) {
+        snd_dsps[i] = snd_dsps[i - 1] + snd_cnts[i - 1];
+    }
+
+    // now let each rank know how many spikes they will receive
+    // and get in turn all the buffer sizes to receive
+    nrnmpi_int_alltoall(&snd_cnts[0], &rcv_cnts[0], 1);
+    for (int i = 1; i < nrnmpi_numprocs; i++) {
+        rcv_dsps[i] = rcv_dsps[i - 1] + rcv_cnts[i - 1];
+    }
+    std::size_t new_sz = 0;
+    for (const auto& r : rcv_cnts) {
+        new_sz += r;
+    }
+    // prepare new sorted vectors
+    std::vector<double> svt_buf(new_sz, 0.0);
+    std::vector<int> svg_buf(new_sz, 0);
+
+    // now exchange data
+    nrnmpi_dbl_alltoallv(spikevec_time.data(), &snd_cnts[0], &snd_dsps[0], svt_buf.data(),
+                         &rcv_cnts[0], &rcv_dsps[0]);
+    nrnmpi_int_alltoallv(spikevec_gid.data(), &snd_cnts[0], &snd_dsps[0], svg_buf.data(),
+                         &rcv_cnts[0], &rcv_dsps[0]);
+
+    local_spikevec_sort(svt_buf, svg_buf, spikevec_time, spikevec_gid);
+}
+
 /** Write generated spikes to out.dat using mpi parallel i/o.
  *  \todo : MPI related code should be factored into nrnmpi.c
  *          Check spike record length which is set to 64 chars
@@ -78,6 +148,7 @@ void output_spikes_parallel(const char* outpath) {
     if (nrnmpi_myid == 0) {
         remove(fname.c_str());
     }
+    sort_spikes(spikevec_time, spikevec_gid);
     nrnmpi_barrier();
 
     // each spike record in the file is time + gid (64 chars sufficient)
@@ -136,6 +207,11 @@ void output_spikes_serial(const char* outpath) {
     ss << outpath << "/out.dat";
     std::string fname = ss.str();
 
+    // reserve some space for sorted spikevec buffers
+    std::vector<double> sorted_spikevec_time(spikevec_time.size());
+    std::vector<int> sorted_spikevec_gid(spikevec_gid.size());
+    local_spikevec_sort(spikevec_time, spikevec_gid, sorted_spikevec_time, sorted_spikevec_gid);
+
     // remove if file already exist
     remove(fname.c_str());
 
@@ -145,9 +221,9 @@ void output_spikes_serial(const char* outpath) {
         return;
     }
 
-    for (unsigned i = 0; i < spikevec_gid.size(); ++i)
-        if (spikevec_gid[i] > -1)
-            fprintf(f, "%.8g\t%d\n", spikevec_time[i], spikevec_gid[i]);
+    for (std::size_t i = 0; i < sorted_spikevec_gid.size(); ++i)
+        if (sorted_spikevec_gid[i] > -1)
+            fprintf(f, "%.8g\t%d\n", sorted_spikevec_time[i], sorted_spikevec_gid[i]);
 
     fclose(f);
 }

coreneuron/nrnmpi/mpispike.cpp

@@ -341,13 +341,22 @@ void nrnmpi_assert_opstep(int opstep, double tt) {
 
 double nrnmpi_dbl_allmin(double x) {
     double result;
-    if (nrnmpi_numprocs < 2) {
+    if (!nrnmpi_use || (nrnmpi_numprocs < 2)) {
         return x;
     }
     MPI_Allreduce(&x, &result, 1, MPI_DOUBLE, MPI_MIN, nrnmpi_comm);
     return result;
 }
 
+double nrnmpi_dbl_allmax(double x) {
+    double result;
+    if (!nrnmpi_use || (nrnmpi_numprocs < 2)) {
+        return x;
+    }
+    MPI_Allreduce(&x, &result, 1, MPI_DOUBLE, MPI_MAX, nrnmpi_comm);
+    return result;
+}
+
 static void pgvts_op(double* in, double* inout, int* len, MPI_Datatype* dptr) {
     int i, r = 0;
     if (*dptr != MPI_DOUBLE)

coreneuron/nrnmpi/nrnmpidec.h

@@ -110,6 +110,7 @@ extern void nrnmpi_char_broadcast(char* buf, int cnt, int root);
 extern int nrnmpi_int_sum_reduce(int in);
 extern void nrnmpi_assert_opstep(int opstep, double t);
 extern double nrnmpi_dbl_allmin(double x);
+extern double nrnmpi_dbl_allmax(double x);
 extern int nrnmpi_pgvts_least(double* t, int* op, int* init);
 extern void nrnmpi_send_doubles(double* pd, int cnt, int dest, int tag);
 extern void nrnmpi_recv_doubles(double* pd, int cnt, int src, int tag);

coreneuron/nrnoc/eion.cpp

@@ -71,7 +71,7 @@ THE POSSIBILITY OF SUCH DAMAGE.
 namespace coreneuron {
 
 #define nparm 5
-static char* mechanism[] = {/*just a template*/
+static const char* mechanism[] = {/*just a template*/
                             "0", "na_ion", "ena", "nao", "nai", 0, "ina", "dina_dv_", 0, 0};
 
 void nrn_init_ion(NrnThread*, Memb_list*, int);

coreneuron/scopmath_core/sparse_thread.cpp

@@ -224,7 +224,7 @@ int _cvode_sparse_thread(void** v, int n, int* x, SPFUN fun, _threadargsproto_)
 }
 
 static int matsol(SparseObj* so, int _iml) {
-    register Elm *pivot, *el;
+    Elm *pivot, *el;
     unsigned i;
 
     /* Upper triangularization */
@@ -244,7 +244,7 @@ static int matsol(SparseObj* so, int _iml) {
 
 static void subrow(SparseObj* so, Elm* pivot, Elm* rowsub, int _iml) {
     double r;
-    register Elm* el;
+    Elm* el;
 
     int _cntml_padded = so->_cntml_padded;
     r = rowsub->value[_iml] / pivot->value[_iml];
@@ -300,7 +300,7 @@ static void prmat(SparseObj* so) {
 
 static void initeqn(SparseObj* so, unsigned maxeqn) /* reallocate space for matrix */
 {
-    register unsigned i, nn;
+    unsigned i, nn;
 
     if (maxeqn == so->neqn)
         return;
@@ -361,7 +361,7 @@ saves much time allocating and freeing during the solve phase
 static Elm* getelm(SparseObj* so, unsigned row, unsigned col, Elm* new_elem)
 /* return pointer to row col element maintaining order in rows */
 {
-    register Elm *el, *elnext;
+    Elm *el, *elnext;
     unsigned vrow, vcol;
 
     vrow = so->varord[row];

external/mod2c

@@ -30,8 +30,7 @@ then
   exit 1
 fi
 
-sort -k 1n,1n -k 2n,2n out.dat > sort_out.dat
-diff -w sort_out.dat @CMAKE_CURRENT_SOURCE_DIR@/@SIM_NAME@/out.dat.ref > diff.dat 2>&1
+diff -w out.dat @CMAKE_CURRENT_SOURCE_DIR@/@SIM_NAME@/out.dat.ref > diff.dat 2>&1
 
 if [ -s diff.dat ]
 then

tests/integration/integration_test.sh.in

@@ -86,11 +86,11 @@ template <> const unsigned char check_data<T>::data[]
 
 CHECK_DATA(char,'x','x')={'x'};
 
-CHECK_DATA(float,(float)0x1.8a4782p+79,(float)-0x1.468acep+3)=
-    {0xc1,0x23,0x45,0x67};
+CHECK_DATA(float,(float)1.1,(float)-428967904.0)=
+    {0xcd,0xcc,0x8c,0x3f};
 
-CHECK_DATA(double,-0x1.dab89674523c1p+45,-0x1.3456789abcdc2p+19)=
-    {0xc1,0x23,0x45,0x67,0x89,0xab,0xcd,0xc2};
+CHECK_DATA(double,1.1,-1.5423487136706484e-180)=
+    {0x9a,0x99,0x99,0x99,0x99,0x99,0xf1,0x3f};
 
 CHECK_DATA(uint16_t,0xf12e,0x2ef1)=
     {0x2e,0xf1};