From ed0b46126332e5f817e90368aff5ac14114a57dc Mon Sep 17 00:00:00 2001 From: JCGoran Date: Mon, 29 Jan 2024 19:12:45 +0100 Subject: [PATCH] Refactor `nrn_mlh_gsort` with `std::sort` (#2692) Co-authored-by: Nicolas Cornu --- src/ivoc/ivocvect.cpp | 196 +------------------------------------- test/CMakeLists.txt | 1 + test/unit_tests/iovec.cpp | 45 +++++++++ 3 files changed, 50 insertions(+), 192 deletions(-) create mode 100644 test/unit_tests/iovec.cpp diff --git a/src/ivoc/ivocvect.cpp b/src/ivoc/ivocvect.cpp index 7b9087d328..866fab6836 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,9 @@ 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; - } - } - } - } + 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/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..946ad7c7fa --- /dev/null +++ b/test/unit_tests/iovec.cpp @@ -0,0 +1,45 @@ +#include +#include + +#include "oc_ansi.h" + +#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); + + // 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]]); + } + } + } + + { + std::vector indices{2, 1, 1}; + std::vector expected_result{1, 1, 2}; // as -2,5 < 5.1 + + 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]); + } + } + } +}