From fe9ef106a3f9c86463ef107e9daabe54d34b419f Mon Sep 17 00:00:00 2001 From: Goran Jelic-Cizmek Date: Tue, 23 Jan 2024 17:17:16 +0100 Subject: [PATCH 1/4] Refactor `nrn_mlh_gsort` * use `std::sort` instead of a custom one * remove unnecessary macros and inline functions * add tests * fixes #2681 --- src/ivoc/ivocvect.cpp | 201 ++------------------------------------ src/oc/oc_ansi.h | 3 +- test/CMakeLists.txt | 1 + test/unit_tests/iovec.cpp | 26 +++++ 4 files changed, 36 insertions(+), 195 deletions(-) create mode 100644 test/unit_tests/iovec.cpp diff --git a/src/ivoc/ivocvect.cpp b/src/ivoc/ivocvect.cpp index cc6565bea3..e1d5bb53b0 100644 --- a/src/ivoc/ivocvect.cpp +++ b/src/ivoc/ivocvect.cpp @@ -1,6 +1,7 @@ #include <../../nrnconf.h> //#include +#include #include #include #include @@ -3875,198 +3876,10 @@ void Vector_reg() { #endif } -// hacked version of gsort from ../gnu/d_vec.cpp -// the transformation is that everything that used to be a double* becomes -// an int* and cmp(*arg1, *arg2) becomes cmp(vec[*arg1], vec[*arg2]) -// I am not sure what to do about the BYTES_PER_WORD - -// An adaptation of Schmidt's new quicksort - -static inline void SWAP(int* A, int* B) { - int tmp = *A; - *A = *B; - *B = tmp; -} - -/* This should be replaced by a standard ANSI macro. */ -#define BYTES_PER_WORD 8 -#define BYTES_PER_LONG 4 - -/* The next 4 #defines implement a very fast in-line stack abstraction. */ - -#define STACK_SIZE (BYTES_PER_WORD * BYTES_PER_LONG) -#define PUSH(LOW, HIGH) \ - do { \ - top->lo = LOW; \ - top++->hi = HIGH; \ - } while (0) -#define POP(LOW, HIGH) \ - do { \ - LOW = (--top)->lo; \ - HIGH = top->hi; \ - } while (0) -#define STACK_NOT_EMPTY (stack < top) - -/* Discontinue quicksort algorithm when partition gets below this size. - This particular magic number was chosen to work best on a Sun 4/260. */ -#define MAX_THRESH 4 - - -/* Order size using quicksort. This implementation incorporates - four optimizations discussed in Sedgewick: - - 1. Non-recursive, using an explicit stack of pointer that - store the next array partition to sort. To save time, this - maximum amount of space required to store an array of - MAX_INT is allocated on the stack. Assuming a 32-bit integer, - this needs only 32 * sizeof (stack_node) == 136 bits. Pretty - cheap, actually. - - 2. Chose the pivot element using a median-of-three decision tree. - This reduces the probability of selecting a bad pivot value and - eliminates certain extraneous comparisons. - - 3. Only quicksorts TOTAL_ELEMS / MAX_THRESH partitions, leaving - insertion sort to order the MAX_THRESH items within each partition. - This is a big win, since insertion sort is faster for small, mostly - sorted array segements. - - 4. The larger of the two sub-partitions is always pushed onto the - stack first, with the algorithm then concentrating on the - smaller partition. This *guarantees* no more than log (n) - stack size is needed! */ - -int nrn_mlh_gsort(double* vec, int* base_ptr, int total_elems, int (*cmp)(double, double)) { - /* Stack node declarations used to store unfulfilled partition obligations. */ - struct stack_node { - int* lo; - int* hi; - }; - int pivot_buffer; - int max_thresh = MAX_THRESH; - - if (total_elems > MAX_THRESH) { - int* lo = base_ptr; - int* hi = lo + (total_elems - 1); - int* left_ptr; - int* right_ptr; - stack_node stack[STACK_SIZE]; /* Largest size needed for 32-bit int!!! */ - stack_node* top = stack + 1; - - while (STACK_NOT_EMPTY) { - { - int* pivot = &pivot_buffer; - { - /* Select median value from among LO, MID, and HI. Rearrange - LO and HI so the three values are sorted. This lowers the - probability of picking a pathological pivot value and - skips a comparison for both the LEFT_PTR and RIGHT_PTR. */ - - int* mid = lo + ((hi - lo) >> 1); - - if (cmp(vec[*mid], vec[*lo]) < 0) - SWAP(mid, lo); - if (cmp(vec[*hi], vec[*mid]) < 0) { - SWAP(mid, hi); - if (cmp(vec[*mid], vec[*lo]) < 0) - SWAP(mid, lo); - } - *pivot = *mid; - pivot = &pivot_buffer; - } - left_ptr = lo + 1; - right_ptr = hi - 1; - - /* Here's the famous ``collapse the walls'' section of quicksort. - Gotta like those tight inner loops! They are the main reason - that this algorithm runs much faster than others. */ - do { - while (cmp(vec[*left_ptr], vec[*pivot]) < 0) - left_ptr += 1; - - while (cmp(vec[*pivot], vec[*right_ptr]) < 0) - right_ptr -= 1; - - if (left_ptr < right_ptr) { - SWAP(left_ptr, right_ptr); - left_ptr += 1; - right_ptr -= 1; - } else if (left_ptr == right_ptr) { - left_ptr += 1; - right_ptr -= 1; - break; - } - } while (left_ptr <= right_ptr); - } - - /* Set up pointers for next iteration. First determine whether - left and right partitions are below the threshold size. If so, - ignore one or both. Otherwise, push the larger partition's - bounds on the stack and continue sorting the smaller one. */ - - if ((right_ptr - lo) <= max_thresh) { - if ((hi - left_ptr) <= max_thresh) /* Ignore both small partitions. */ - POP(lo, hi); - else /* Ignore small left partition. */ - lo = left_ptr; - } else if ((hi - left_ptr) <= max_thresh) /* Ignore small right partition. */ - hi = right_ptr; - else if ((right_ptr - lo) > (hi - left_ptr)) /* Push larger left partition indices. */ - { - PUSH(lo, right_ptr); - lo = left_ptr; - } else /* Push larger right partition indices. */ - { - PUSH(left_ptr, hi); - hi = right_ptr; - } - } - } - - /* Once the BASE_PTR array is partially sorted by quicksort the rest - is completely sorted using insertion sort, since this is efficient - for partitions below MAX_THRESH size. BASE_PTR points to the beginning - of the array to sort, and END_PTR points at the very last element in - the array (*not* one beyond it!). */ - - - { - int* end_ptr = base_ptr + 1 * (total_elems - 1); - int* run_ptr; - int* tmp_ptr = base_ptr; - int* thresh = (end_ptr < (base_ptr + max_thresh)) ? end_ptr : (base_ptr + max_thresh); - - /* Find smallest element in first threshold and place it at the - array's beginning. This is the smallest array element, - and the operation speeds up insertion sort's inner loop. */ - - for (run_ptr = tmp_ptr + 1; run_ptr <= thresh; run_ptr += 1) - if (cmp(vec[*run_ptr], vec[*tmp_ptr]) < 0) - tmp_ptr = run_ptr; - - if (tmp_ptr != base_ptr) - SWAP(tmp_ptr, base_ptr); - - /* Insertion sort, running from left-hand-side up to `right-hand-side.' - Pretty much straight out of the original GNU qsort routine. */ - - for (run_ptr = base_ptr + 1; (tmp_ptr = run_ptr += 1) <= end_ptr;) { - while (cmp(vec[*run_ptr], vec[*(tmp_ptr -= 1)]) < 0) - ; - - if ((tmp_ptr += 1) != run_ptr) { - int* trav; - - for (trav = run_ptr + 1; --trav >= run_ptr;) { - int c = *trav; - int *hi, *lo; - - for (hi = lo = trav; (lo -= 1) >= tmp_ptr; hi = lo) - *hi = *lo; - *hi = c; - } - } - } - } - return 1; +template +int nrn_mlh_gsort(double *vec, int *base_ptr, int total_elems, F cmp) { + std::sort(base_ptr, base_ptr + total_elems, [&](const int &a, const int &b) { + return cmp(vec[a], vec[b]) < 0; + }); + return 1; } diff --git a/src/oc/oc_ansi.h b/src/oc/oc_ansi.h index 34cc37062f..2ef49ee3bc 100644 --- a/src/oc/oc_ansi.h +++ b/src/oc/oc_ansi.h @@ -104,7 +104,8 @@ std::FILE* hoc_obj_file_arg(int i); void hoc_reg_nmodl_text(int type, const char* txt); void hoc_reg_nmodl_filename(int type, const char* filename); std::size_t nrn_mallinfo(int item); -int nrn_mlh_gsort(double* vec, int* base_ptr, int total_elems, int (*cmp)(double, double)); +template +int nrn_mlh_gsort(double* vec, int* base_ptr, int total_elems, F cmp); void state_discontinuity(int i, double* pd, double d); IvocVect* vector_arg(int); diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index d5e9762cd7..e5731f55fb 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -17,6 +17,7 @@ add_executable( testneuron common/catch2_main.cpp unit_tests/basic.cpp + unit_tests/iovec.cpp unit_tests/container/container.cpp unit_tests/container/generic_data_handle.cpp unit_tests/container/mechanism.cpp diff --git a/test/unit_tests/iovec.cpp b/test/unit_tests/iovec.cpp new file mode 100644 index 0000000000..4d25e16b4b --- /dev/null +++ b/test/unit_tests/iovec.cpp @@ -0,0 +1,26 @@ +#include +#include + +#include "oc_ansi.h" + +#include + +TEST_CASE("Test nrn_mlh_gsort output", "[nrn_gsort]") { + std::vector input{1.2, -2.5, 5.1}; + + std::vector indices(input.size()); + // all values from 0 to size - 1 + std::iota(indices.begin(), indices.end(), 0); + + // for comparison + auto sorted_input = input; + std::sort(sorted_input.begin(), sorted_input.end()); + + SECTION("Test sorting") { + nrn_mlh_gsort(input.data(), indices.data(), input.size(), + [](double a, double b) { return a < b; }); + for (auto i = 0; i < input.size(); ++i) { + REQUIRE(sorted_input[i] == input[indices[i]]); + } + } +} From ad9bc90a225befecf29d6aa3cd7e5c48d9e4a6ab Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Mon, 29 Jan 2024 09:53:44 +0100 Subject: [PATCH 2/4] Format --- src/ivoc/ivocvect.cpp | 10 +++++----- test/unit_tests/iovec.cpp | 27 ++++++++++++++------------- 2 files changed, 19 insertions(+), 18 deletions(-) diff --git a/src/ivoc/ivocvect.cpp b/src/ivoc/ivocvect.cpp index e1d5bb53b0..efd6c914ad 100644 --- a/src/ivoc/ivocvect.cpp +++ b/src/ivoc/ivocvect.cpp @@ -3877,9 +3877,9 @@ void Vector_reg() { } template -int nrn_mlh_gsort(double *vec, int *base_ptr, int total_elems, F cmp) { - std::sort(base_ptr, base_ptr + total_elems, [&](const int &a, const int &b) { - return cmp(vec[a], vec[b]) < 0; - }); - return 1; +int nrn_mlh_gsort(double* vec, int* base_ptr, int total_elems, F cmp) { + std::sort(base_ptr, base_ptr + total_elems, [&](const int& a, const int& b) { + return cmp(vec[a], vec[b]) < 0; + }); + return 1; } diff --git a/test/unit_tests/iovec.cpp b/test/unit_tests/iovec.cpp index 4d25e16b4b..5a09b35379 100644 --- a/test/unit_tests/iovec.cpp +++ b/test/unit_tests/iovec.cpp @@ -6,21 +6,22 @@ #include TEST_CASE("Test nrn_mlh_gsort output", "[nrn_gsort]") { - std::vector input{1.2, -2.5, 5.1}; + std::vector input{1.2, -2.5, 5.1}; - std::vector indices(input.size()); - // all values from 0 to size - 1 - std::iota(indices.begin(), indices.end(), 0); + std::vector indices(input.size()); + // all values from 0 to size - 1 + std::iota(indices.begin(), indices.end(), 0); - // for comparison - auto sorted_input = input; - std::sort(sorted_input.begin(), sorted_input.end()); + // for comparison + auto sorted_input = input; + std::sort(sorted_input.begin(), sorted_input.end()); - SECTION("Test sorting") { - nrn_mlh_gsort(input.data(), indices.data(), input.size(), - [](double a, double b) { return a < b; }); - for (auto i = 0; i < input.size(); ++i) { - REQUIRE(sorted_input[i] == input[indices[i]]); + SECTION("Test sorting") { + nrn_mlh_gsort(input.data(), indices.data(), input.size(), [](double a, double b) { + return a < b; + }); + for (auto i = 0; i < input.size(); ++i) { + REQUIRE(sorted_input[i] == input[indices[i]]); + } } - } } From 4a00c66709656c1d5d618dad3251611b1ed5bc7d Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Mon, 29 Jan 2024 10:07:46 +0100 Subject: [PATCH 3/4] Remove template, add a test --- src/ivoc/ivocvect.cpp | 3 +-- src/oc/oc_ansi.h | 3 +-- test/unit_tests/iovec.cpp | 42 ++++++++++++++++++++++++++++----------- 3 files changed, 32 insertions(+), 16 deletions(-) diff --git a/src/ivoc/ivocvect.cpp b/src/ivoc/ivocvect.cpp index ae2d29ba6c..b124e516fe 100644 --- a/src/ivoc/ivocvect.cpp +++ b/src/ivoc/ivocvect.cpp @@ -3876,8 +3876,7 @@ void Vector_reg() { #endif } -template -int nrn_mlh_gsort(double* vec, int* base_ptr, int total_elems, F cmp) { +int nrn_mlh_gsort(double* vec, int* base_ptr, int total_elems, int (*cmp)(double, double)) { std::sort(base_ptr, base_ptr + total_elems, [&](const int& a, const int& b) { return cmp(vec[a], vec[b]) < 0; }); diff --git a/src/oc/oc_ansi.h b/src/oc/oc_ansi.h index 2ef49ee3bc..34cc37062f 100644 --- a/src/oc/oc_ansi.h +++ b/src/oc/oc_ansi.h @@ -104,8 +104,7 @@ std::FILE* hoc_obj_file_arg(int i); void hoc_reg_nmodl_text(int type, const char* txt); void hoc_reg_nmodl_filename(int type, const char* filename); std::size_t nrn_mallinfo(int item); -template -int nrn_mlh_gsort(double* vec, int* base_ptr, int total_elems, F cmp); +int nrn_mlh_gsort(double* vec, int* base_ptr, int total_elems, int (*cmp)(double, double)); void state_discontinuity(int i, double* pd, double d); IvocVect* vector_arg(int); diff --git a/test/unit_tests/iovec.cpp b/test/unit_tests/iovec.cpp index 5a09b35379..02e8558d2c 100644 --- a/test/unit_tests/iovec.cpp +++ b/test/unit_tests/iovec.cpp @@ -5,23 +5,41 @@ #include +// This function is the one that is used in all nrn-modeldb-ci +// Keep as is +int cmpdfn(double a, double b) { + return ((a) <= (b)) ? (((a) == (b)) ? 0 : -1) : 1; +} + TEST_CASE("Test nrn_mlh_gsort output", "[nrn_gsort]") { std::vector input{1.2, -2.5, 5.1}; - std::vector indices(input.size()); - // all values from 0 to size - 1 - std::iota(indices.begin(), indices.end(), 0); + { + std::vector indices(input.size()); + // all values from 0 to size - 1 + std::iota(indices.begin(), indices.end(), 0); + + // for comparison + auto sorted_input = input; + std::sort(sorted_input.begin(), sorted_input.end()); + + SECTION("Test sorting") { + nrn_mlh_gsort(input.data(), indices.data(), input.size(), cmpdfn); + for (auto i = 0; i < input.size(); ++i) { + REQUIRE(sorted_input[i] == input[indices[i]]); + } + } + } - // for comparison - auto sorted_input = input; - std::sort(sorted_input.begin(), sorted_input.end()); + { + std::vector indices{2, 1, 1}; + std::vector expected_result{1, 1, 2}; // as -2,5 < 5.1 - SECTION("Test sorting") { - nrn_mlh_gsort(input.data(), indices.data(), input.size(), [](double a, double b) { - return a < b; - }); - for (auto i = 0; i < input.size(); ++i) { - REQUIRE(sorted_input[i] == input[indices[i]]); + SECTION("Test sorting with repeted indices") { + nrn_mlh_gsort(input.data(), indices.data(), input.size(), cmpdfn); + for (auto i = 0; i < input.size(); ++i) { + REQUIRE(indices[i] == expected_result[i]); + } } } } From 2a7de92b48a27660329462da0b5302dbe1e8d537 Mon Sep 17 00:00:00 2001 From: Nicolas Cornu Date: Mon, 29 Jan 2024 16:52:25 +0100 Subject: [PATCH 4/4] fix review --- src/ivoc/ivocvect.cpp | 2 +- test/unit_tests/iovec.cpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/ivoc/ivocvect.cpp b/src/ivoc/ivocvect.cpp index b124e516fe..866fab6836 100644 --- a/src/ivoc/ivocvect.cpp +++ b/src/ivoc/ivocvect.cpp @@ -3877,7 +3877,7 @@ void Vector_reg() { } int nrn_mlh_gsort(double* vec, int* base_ptr, int total_elems, int (*cmp)(double, double)) { - std::sort(base_ptr, base_ptr + total_elems, [&](const int& a, const int& b) { + std::sort(base_ptr, base_ptr + total_elems, [&](int a, int b) { return cmp(vec[a], vec[b]) < 0; }); return 1; diff --git a/test/unit_tests/iovec.cpp b/test/unit_tests/iovec.cpp index 02e8558d2c..946ad7c7fa 100644 --- a/test/unit_tests/iovec.cpp +++ b/test/unit_tests/iovec.cpp @@ -35,7 +35,7 @@ TEST_CASE("Test nrn_mlh_gsort output", "[nrn_gsort]") { std::vector indices{2, 1, 1}; std::vector expected_result{1, 1, 2}; // as -2,5 < 5.1 - SECTION("Test sorting with repeted indices") { + SECTION("Test sorting with repeated indices") { nrn_mlh_gsort(input.data(), indices.data(), input.size(), cmpdfn); for (auto i = 0; i < input.size(); ++i) { REQUIRE(indices[i] == expected_result[i]);