Use long dtype for indexing and stdlib quicksort

ext/_core.c

@@ -16,126 +16,68 @@
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */
+#include "_core.h"
 
 #include <stdlib.h>
 #include <math.h>
 
-/***************************************************************************************************
-*   C99 compatibility for macros INFINITY and NAN
-***************************************************************************************************/
-
-#ifdef _MSC_VER
-    /* handle Microsofts C99 incompatibility */
-    #include <float.h>
-    #define INFINITY (DBL_MAX+DBL_MAX)
-    #define NAN (INFINITY-INFINITY)
-#else
-    /* if not available otherwise, define INFINITY/NAN in the GNU style */
-    #ifndef INFINITY
-        #define INFINITY (1.0/0.0)
-    #endif
-    #ifndef NAN
-        #define NAN (INFINITY-INFINITY)
-    #endif
-#endif
-
 /***************************************************************************************************
 *   static convenience functions (without cython wrappers)
 ***************************************************************************************************/
 
-static double sqr(double x)
-{
+static inline double sqr(double x) {
     return (x == 0.0) ? 0.0 : x * x;
 }
 
-static double distance(double *points, int n, int m, int ndim)
-{
-    int i, o = n * ndim, p = m * ndim;
+static double distance(double *points, size_t n, size_t m, size_t ndim) {
+    size_t i, o = n * ndim, p = m * ndim;
     double sum = 0.0;
-    for(i=0; i<ndim; ++i)
+    for (i = 0; i < ndim; ++i)
         sum += sqr(points[o + i] - points[p + i]);
     return sqrt(sum);
 }
 
-static void mixed_sort(double *array, int L, int R)
-/* mixed_sort() is based on examples from http://www.linux-related.de (2004) */
-{
-    int l, r;
-    double swap;
-    if(R - L > 25) /* use quicksort */
-    {
-        l = L - 1;
-        r = R;
-        for(;;)
-        {
-            while(array[++l] < array[R]);
-            while((array[--r] > array[R]) && (r > l));
-            if(l >= r) break;
-            swap = array[l];
-            array[l] = array[r];
-            array[r] = swap;
-        }
-        swap = array[l];
-        array[l] = array[R];
-        array[R] = swap;
-        mixed_sort(array, L, l - 1);
-        mixed_sort(array, l + 1, R);
-    }
-    else /* use insertion sort */
-    {
-        for(l=L+1; l<=R; ++l)
-        {
-            swap = array[l];
-            for(r=l-1; (r >= L) && (swap < array[r]); --r)
-                array[r + 1] = array[r];
-            array[r + 1] = swap;
-        }
-    }
+static int compare_double(const void * a, const void * b) {
+    double da = *(const double*) a;
+    double db = *(const double*) b;
+    return (da > db) - (da < db);
 }
 
 /***************************************************************************************************
 *   pydpc core functions (with cython wrappers)
 ***************************************************************************************************/
 
-extern void _get_distances(double *points, int npoints, int ndim, double *distances)
-{
-    int i, j, o;
-    for(i=0; i<npoints-1; ++i)
-    {
-        o = i * npoints;
-        for(j=i+1; j<npoints; ++j)
-        {
+extern void _get_distances(double *points, size_t npoints, size_t ndim, double *distances) {
+    size_t i, j;
+    for (i = 0; i < npoints - 1; ++i) {
+        size_t o = i * npoints;
+        for (j = i + 1; j < npoints; ++j) {
             distances[o + j] = distance(points, i, j, ndim);
             distances[j * npoints + i] = distances[o + j];
         }
     }
 }
 
-extern double _get_kernel_size(double *distances, int npoints, double fraction)
-{
-    int i, j, o, m = 0, n = (npoints * (npoints - 1)) / 2;
+extern double _get_kernel_size(double *distances, size_t npoints, double fraction) {
+    size_t i, j, m = 0, n = (npoints * (npoints - 1)) / 2;
     double kernel_size;
-    double *scratch = (double*) malloc(n * sizeof(double));
-    for(i=0; i<npoints-1; ++i)
-    {
-        o = i * npoints;
-        for(j=i+1; j<npoints; ++j)
+    double *scratch = (double *) malloc(n * sizeof(double));
+    for (i = 0; i < npoints - 1; ++i) {
+        size_t o = i * npoints;
+        for (j = i + 1; j < npoints; ++j)
             scratch[m++] = distances[o + j];
     }
-    mixed_sort(scratch, 0, n - 1);
-    kernel_size = scratch[(int) floor(0.5 + fraction * n)];
+    qsort(scratch, n, sizeof(double), compare_double);
+    kernel_size = scratch[(size_t) floor(0.5 + fraction * n)];
 
     /* kernel_size is used as a divisor, it can't be zero. Fall back to the first
      * nonzero value. Technically, they *could* all be zero (degenerate case
      * where all points are the same), in which case zero will still be returned.
      * Client code checks for this.
      */
-    if(kernel_size == 0.0)
-    {
-        for(int i = 0; i < n; i++)
-        {
-            if(scratch[i] != 0.0)
-            {
+    if (kernel_size == 0.0) {
+        for (size_t i = 0; i < n; i++) {
+            if (scratch[i] != 0.0) {
                 kernel_size = scratch[i];
                 break;
             }
@@ -146,15 +88,12 @@ extern double _get_kernel_size(double *distances, int npoints, double fraction)
     return kernel_size;
 }
 
-extern void _get_density(double kernel_size, double *distances, int npoints, double *density)
-{
-    int i, j, o;
+extern void _get_density(double kernel_size, double *distances, size_t npoints, double *density) {
+    size_t i, j;
     double rho;
-    for(i=0; i<npoints-1; ++i)
-    {
-        o = i * npoints;
-        for(j=i+1; j<npoints; ++j)
-        {
+    for (i = 0; i < npoints - 1; ++i) {
+        size_t o = i * npoints;
+        for (j = i + 1; j < npoints; ++j) {
             rho = exp(-sqr(distances[o + j] / kernel_size));
             density[i] += rho;
             density[j] += rho;
@@ -163,23 +102,18 @@ extern void _get_density(double kernel_size, double *distances, int npoints, dou
 }
 
 extern void _get_delta_and_neighbour(
-    double max_distance, double *distances, int *order, int npoints, double *delta, int *neighbour)
-{
-    int i, j, o;
+        double max_distance, double *distances, size_t *order, size_t npoints, double *delta, long *neighbour) {
+    size_t i, j, o;
     double max_delta = 0.0;
-    for(i=0; i<npoints; ++i)
-    {
+    for (i = 0; i < npoints; ++i) {
         delta[order[i]] = max_distance;
         neighbour[i] = -1;
     }
     delta[order[0]] = -1.0;
-    for(i=1; i<npoints; ++i)
-    {
+    for (i = 1; i < npoints; ++i) {
         o = order[i] * npoints;
-        for(j=0; j<i; ++j)
-        {
-            if(distances[o + order[j]] < delta[order[i]])
-            {
+        for (j = 0; j < i; ++j) {
+            if (distances[o + order[j]] < delta[order[i]]) {
                 delta[order[i]] = distances[o + order[j]];
                 neighbour[order[i]] = order[j];
             }
@@ -190,56 +124,46 @@ extern void _get_delta_and_neighbour(
 }
 
 extern void _get_membership(
-    int *clusters, int nclusters, int *order, int *neighbour, int npoints, int *membership)
-{
-    int i;
-    for(i=0; i<npoints; ++i)
+        size_t *clusters, size_t nclusters, size_t *order, long *neighbour, size_t npoints, long *membership) {
+    size_t i;
+    for (i = 0; i < npoints; ++i)
         membership[i] = -1;
-    for(i=0; i<nclusters; ++i)
+    for (i = 0; i < nclusters; ++i)
         membership[clusters[i]] = i;
-    for(i=0; i<npoints; ++i)
-    {
-        if(membership[order[i]] == -1)
+    for (i = 0; i < npoints; ++i) {
+        if (membership[order[i]] == - 1)
             membership[order[i]] = membership[neighbour[order[i]]];
     }
 }
 
 extern void _get_border(
-    double kernel_size, double *distances, double *density, int *membership, int npoints,
-    int *border_member, double *border_density)
-{
-    int i, j, o;
+        double kernel_size, double *distances, double *density, long *membership, size_t npoints,
+        bool *border_member, double *border_density) {
+    size_t i, j, o;
     double average_density;
-    for(i=0; i<npoints-1; ++i)
-    {
+    for (i = 0; i < npoints - 1; ++i) {
         o = i * npoints;
-        for(j=i+1; j<npoints; ++j)
-        {
-            if((membership[i] != membership[j]) && (distances[o + j] < kernel_size))
-            {
+        for (j = i + 1; j < npoints; ++j) {
+            if ((membership[i] != membership[j]) && (distances[o + j] < kernel_size)) {
                 average_density = 0.5 * (density[i] + density[j]);
-                if(border_density[membership[i]] < average_density)
+                if (border_density[membership[i]] < average_density)
                     border_density[membership[i]] = average_density;
-                if(border_density[membership[j]] < average_density)
+                if (border_density[membership[j]] < average_density)
                     border_density[membership[j]] = average_density;
-                border_member[i] = 1;
-                border_member[j] = 1;
+                border_member[i] = true;
+                border_member[j] = true;
             }
         }
     }
 }
 
 extern void _get_halo(
-    int border_only, double *border_density,
-    double *density, int *membership, int *border_member, int npoints, int *halo)
-{
-    int i;
-    for(i=0; i<npoints; ++i)
-    {
-        if(density[i] < border_density[membership[i]])
-        {
-            if((0 == border_only) || (1 == border_member[i]))
-                halo[i] = -1;
+        bool border_only, double *border_density,
+        double *density, long *membership, bool *border_member, size_t npoints, int *halo) {
+    size_t i;
+    for (i = 0; i < npoints; ++i) {
+        if(density[i] < border_density[membership[i]] && (! border_only || border_member[i])) {
+            halo[i] = -1;
         }
     }
 }

ext/_core.h

@@ -19,19 +19,21 @@
 
 #ifndef PYDPC_HEADER
 #define PYDPC_HEADER
+#include <stdbool.h>
+#include <stddef.h>
 
-extern void _get_distances(double *points, int npoints, int ndim, double *distances);
-extern double _get_kernel_size(double *distances, int npoints, double fraction);
-extern void _get_density(double kernel_size, double *distances, int npoints, double *density);
+extern void _get_distances(double *points, size_t npoints, size_t ndim, double *distances);
+extern double _get_kernel_size(double *distances, size_t npoints, double fraction);
+extern void _get_density(double kernel_size, double *distances, size_t npoints, double *density);
 extern void _get_delta_and_neighbour(
-    double max_distance, double *distances, int *order, int npoints, double *delta, int *neighbour);
+    double max_distance, double *distances, size_t *order, size_t npoints, double *delta, long *neighbour);
 extern void _get_membership(
-    int *clusters, int nclusters, int *order, int *neighbour, int npoints, int *membership);
+    size_t *clusters, size_t nclusters, size_t *order, long *neighbour, size_t npoints, long *membership);
 extern void _get_border(
-    double kernel_size, double *distances, double *density, int *membership, int npoints,
-    int *border_member, double *border_density);
+    double kernel_size, double *distances, double *density, long *membership, size_t npoints,
+    bool *border_member, double *border_density);
 extern void _get_halo(
-    int border_only, double *border_density,
-    double *density, int *membership, int *border_member, int npoints, int *halo);
+    bool border_only, double *border_density,
+    double *density, long *membership, bool *border_member, size_t npoints, int *halo);
 
 #endif