Refactor nrn_mlh_gsort

src/ivoc/ivocvect.cpp

@@ -1,6 +1,7 @@
 #include <../../nrnconf.h>
 
 //#include <string.h>
+#include <algorithm>
 #include <cstdio>
 #include <cstdlib>
 #include <cmath>
@@ -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;
 }

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

test/unit_tests/iovec.cpp

@@ -0,0 +1,45 @@
+#include <algorithm>
+#include <vector>
+
+#include "oc_ansi.h"
+
+#include <catch2/catch.hpp>
+
+// 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<double> input{1.2, -2.5, 5.1};
+
+    {
+        std::vector<int> 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<int> indices{2, 1, 1};
+        std::vector<int> 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]);
+            }
+        }
+    }
+}