uncool PR for getMaxColIndex stuff

R/getMaxColIndex.R

@@ -31,7 +31,7 @@
 #' print(getMaxIndexOfRows(x))
 #' print(getMinIndexOfRows(x))
 getMaxIndexOfRows = function(x, weights = NULL, ties.method = "random", na.rm = FALSE) {
-  mode(x) = "numeric"
+  if (!is.double(x)) mode(x) = "numeric"
   ties.method = switch(ties.method, random = 1L, first = 2L, last = 3L,
                        stop("Unknown ties method"))
   assertFlag(na.rm)
@@ -48,7 +48,7 @@ getMinIndexOfRows = function(x, weights = NULL, ties.method = "random", na.rm =
 #' @export
 #' @rdname getMaxIndexOfRows
 getMaxIndexOfCols = function(x, weights = NULL, ties.method = "random", na.rm = FALSE) {
-  mode(x) = "numeric"
+  if (!is.double(x)) mode(x) = "numeric"
   ties.method = switch(ties.method, random = 1L, first = 2L, last = 3L,
                        stop("Unknown ties method"))
   assertFlag(na.rm)

src/getMaxColRowIndex.c

@@ -15,8 +15,10 @@ SEXP c_getMaxIndexOfRows(SEXP s_x, SEXP s_w, SEXP s_ties_method, SEXP s_na_rm) {
     for (R_len_t i = 0; i < nrow_x; i++)
         ret[i] = get_max_index(x + i, ncol_x, nrow_x, ties_method, na_rm);
   } else {
+    double *buf = (double *) malloc(ncol_x * sizeof(double));
     for (R_len_t i = 0; i < nrow_x; i++)
-      ret[i] = get_max_index_w(x + i, w, ncol_x, nrow_x, ties_method, na_rm);
+      ret[i] = get_max_index_w(x + i, w, buf, ncol_x, nrow_x, ties_method, na_rm);
+    free(buf);
   }
   PutRNGstate();
   UNPROTECT(1); /* s_ret */
@@ -36,8 +38,10 @@ SEXP c_getMaxIndexOfCols(SEXP s_x, SEXP s_w, SEXP s_ties_method, SEXP s_na_rm) {
     for (R_len_t i = 0; i < ncol_x; ++i)
       ret[i] = get_max_index(x + i*nrow_x, nrow_x, 1, ties_method, na_rm);
   } else {
+    double *buf = (double *) malloc(nrow_x * sizeof(double));
     for (R_len_t i = 0; i < ncol_x; ++i)
-      ret[i] = get_max_index_w(x + i*nrow_x, w, nrow_x, 1, ties_method, na_rm);
+      ret[i] = get_max_index_w(x + i*nrow_x, w, buf, nrow_x, 1, ties_method, na_rm);
+    free(buf);
   }
   PutRNGstate();
   UNPROTECT(1); /* s_ret */

src/getMaxIndex.c

@@ -14,14 +14,14 @@
   return  : Index of maximal element (1-based) or
             -1 if we did not find a maximal elemnt (empty vector or only removed NAs)
 */
-int get_max_index(double *x, R_len_t n, R_len_t step, int ties_method, Rboolean na_rm) {
+int get_max_index(double *x, R_len_t n, const R_len_t step, int ties_method, Rboolean na_rm) {
   R_len_t i;
   int max_index = -2;
   int number_of_ties = 0;
   double max_value = -DBL_MAX, current_value;
 
   for (i = 0; i < n; ++i) {
-    current_value = x[i*step];
+    current_value = x[i * step];
     if (!na_rm && ISNAN(current_value))
       return NA_INTEGER;
     if (current_value > max_value) {
@@ -45,14 +45,12 @@ int get_max_index(double *x, R_len_t n, R_len_t step, int ties_method, Rboolean
 // get_max_index with weights.
 // copy x vector (only elements indexed by steps) multiplied with w
 // then call get_max_index with step=1
-int get_max_index_w(double *x, double *w, R_len_t n, R_len_t step, int ties_method, Rboolean na_rm) {
-  double *xx = (double *) malloc(n * sizeof(double));
+int get_max_index_w(double *x, double *w, double *buf, R_len_t n, R_len_t step, int ties_method, Rboolean na_rm) {
   R_len_t i;
   for (i=0; i<n; i++) {
-    xx[i] = x[i*step] * w[i];
+    buf[i] = x[i*step] * w[i];
   }
-  int j = get_max_index(xx, n, 1, ties_method, na_rm);
-  free(xx);
+  int j = get_max_index(buf, n, 1, ties_method, na_rm);
   return j;
 }
 
@@ -71,7 +69,9 @@ SEXP c_getMaxIndex(SEXP s_x, SEXP s_w, SEXP s_ties_method, SEXP s_na_rm) {
     index = get_max_index(x, len_x, 1, ties_method, na_rm);
   } else {
     UNPACK_REAL_VECTOR(s_w, w, len_w);
-    index = get_max_index_w(x, w, len_x, 1, ties_method, na_rm);
+    double *buf = (double *) malloc(len_x * sizeof(double));
+    index = get_max_index_w(x, w, buf, len_x, 1, ties_method, na_rm);
+    free(buf);
   }
   PutRNGstate();
   if (index == -1)

src/getMaxIndex.h

@@ -6,7 +6,7 @@
 #include <Rdefines.h>
 
 int get_max_index(double *x, R_len_t n, R_len_t step, int ties_method, Rboolean na_rm);
-int get_max_index_w(double *x, double *w, R_len_t n, R_len_t step, int ties_method, Rboolean na_rm);
+/* int get_max_index_w(double *x, double *w, double *buf, R_len_t n, R_len_t step, int ties_method, Rboolean na_rm); */
 SEXP c_getMaxIndex(SEXP s_x, SEXP s_w, SEXP s_ties_method, SEXP s_na_rm);
 
 #endif

test3.R

@@ -0,0 +1,14 @@
+library(devtools)
+load_all()
+
+# x = 1:3
+# w = c(100, 1, 1)
+# j = getMaxIndex(x, weights = w)
+# print(j)
+
+
+x = matrix(c(1,1,2,1,1,3), ncol = 3)
+print(x)
+j = getMaxIndexOfRows(x, weights = c(2.5, 1, 1))
+print(j)
+

test5.R

@@ -0,0 +1,21 @@
+library(microbenchmark)
+
+load_all()
+
+
+n = 50000
+p = 50
+
+x = matrix(rnorm(n * p), nrow = n)
+w = rep(1, p)
+
+mb = microbenchmark(
+  getMaxIndexOfCols(x),
+  # getMaxIndexOfRows(x, weights = w),
+  colMaxs(x)
+)
+print(mb)
+
+
+
+