Code used to compare looping performance of R/C++ and Julia.

gistfile1.cpp

/**
 * This code is subject to the unlicense (http://unlicense.org/).
 *
 * I used it to compare the performance of C++ code called from R against native Julia.
 * It is based on the blog at:
 *
 * http://feedly.com/e/B5JIqzW8?goback=%2Egde_77616_member_5811511767066509316#%21
 * (accessed on 29 Nov, 2013)
 *
 * The Julia equivalent is given below (just incase that link disappears).
 *
 * On my (windows) machine the Julia code ran apprimately twice as fast as the C++
 * linked from R.
 */


//----------------------------------------------------------------------------------------------------------
// R function
//----------------------------------------------------------------------------------------------------------
/*
cont.mod <- function(burn.in, reps, n, d, l, s) {
    r <- .Call('cont_model', as.integer(reps), as.integer(n), d, l, s)
    return(mean(r[burn.in:reps]))
}
*/


//----------------------------------------------------------------------------------------------------------
// C++ code
//----------------------------------------------------------------------------------------------------------
#include <vector>
#include <algorithm>
#include <R.h>
#include <Rdefines.h>

extern "C" SEXP cont_model(SEXP rreps, SEXP rn, SEXP rd, SEXP rl, SEXP rs) {
    typedef std::vector<double> DVect;
    typedef DVect::iterator DVI;

    SEXP
        res;

    int
        reps(INTEGER(rreps)[0]),
        n(INTEGER(rn)[0]),
        i;

    double
        s(REAL(rs)[0]),
        d(REAL(rs)[0]),
        lxn(n * REAL(rl)[0]),
        sig, mul, tot;

    DVect
        tr(n);

    PROTECT(res = allocVector(REALSXP, reps));
    double *r(REAL(res));

    std::fill(tr.begin(), tr.end(), 0.0);

    GetRNGstate();
    for(i=0; i<reps; ++i, ++r) {
        sig = norm_rand() * d;
        if(sig < 0.0) {
            sig = -sig;
            mul = -1.0;
        } else {
            mul = 1.0;
        }
        tot = 0.0;
        for(DVI tri(tr.begin()), tre(tr.end()); tri!=tre; ++tri) if(sig > *tri) ++tot;
        tot /= lxn;
        *r = mul * tot;
        for(DVI tri(tr.begin()), tre(tr.end()); tri!=tre; ++tri) if(s > unif_rand()) *tri = tot;
    }
    PutRNGstate();
    UNPROTECT(1);
    return res;
}





//----------------------------------------------------------------------------------------------------------
// Julia code:
//----------------------------------------------------------------------------------------------------------

/*
using Distributions

function cont_run(burn_in, reps, n, d, l, s)
    tr = Array(Float64, n)
    r = Array(Float64, reps)
    for i in 1:reps
        aris = 0
        sig = randn() * d
        mul = 1
        if sig < 0
            sig = -sig
            mul = -1
        end
        for k in 1:n
            if sig > tr[k]
                aris += 1
            end
        end
        ari = aris / (l * n)
        r[i] = mul * ari
        for j in 1:n
            if rand() < s
                tr[j] = ari
            end
        end
    end
    mean(r[burn_in:reps])
end

n = 10
t_start = time()
k = Array(Float64, n)
for i in 1:n
    k[i] = cont_run(1000, 10000, 1000, 0.005, 10.0, 0.01)
end
*/