whichmin.cpp

#include <Rcpp.h>
using namespace Rcpp;

// For each row, find which column has the minimum value
// With ties, pick the first
// Result is indexed starting at 1
// [[Rcpp::export]]
IntegerVector whichmin_byrow(const NumericMatrix& x)
{
    const int n = x.rows();
    const int p = x.cols();
    IntegerVector result(n);
    NumericVector the_min(n);

    // initialize with the first column
    the_min = x(_,0);
    for(int i=0; i<n; i++)
        result[i] = 1;

    // best to have loop over rows on the inside
    for(int j=1; j<p; j++) {
        for(int i=0; i<n; i++) {
            if(x(i,j) < the_min[i]) {
                result[i] = j+1;
                the_min[i] = x(i,j);
            }
        }
    }

    return result;
}

// For each column, find which row has the minimum value
// With ties, pick the first
// Result is indexed starting at 1
// [[Rcpp::export]]
IntegerVector whichmin_bycol(const NumericMatrix& x)
{
    const int n = x.rows();
    const int p = x.cols();
    IntegerVector result(p);

    for(int j=0; j<p; j++) {
        double the_min=x(0,j);
        result[j] = 1;
        for(int i=1; i<n; i++) {
            if(x(i,j) < the_min) {
                result[j] = i+1;
                the_min = x(i,j);
            }
        }
    }

    return result;
}

/*** R
n <- 10000
p <- 5000
x <- matrix(runif(n*p), ncol=p)

system.time(wh1 <- apply(x, 1, which.min))
system.time(wh2 <- whichmin_byrow(x))

y <- t(x)
system.time(wh3 <- apply(y, 2, which.min))
system.time(wh4 <- whichmin_bycol(y))
*/