# Computing the distribution of a Jonckheere-Terpstra statistic

I have the following C++ function which I'd like to speed up, if possible. This program will have many different users, so parallelization is not really an option.

The C++ function is wrapped by C code and then called from R, using:

R CMD SHLIB -o jonckterp.so wrapper.cpp inside.cpp


from the command line, then:

dyn.load("jonckterp.so")


from within R.

inside.cpp

#include <stdio.h>
#include <vector>
#include <iostream>
#include <R.h>
#include <Rmath.h>

extern "C" void jonckterp(int *x,int *c, int *ngrp, double *comb);

int *finalcountdown;
static int count;

static int *tabInit(int num){
int *w;
w = (int *) R_alloc(num, sizeof(int));
memset(w, '\0', num * sizeof(int));
return(w);
}

static int mannwhit(std::vector<int> c, std::vector<int> sizes){
int Ustat=0;
int start=0;
size_t Ngroups=sizes.size();
for(int i=0; i<Ngroups-1;i++){
for(int j=start;j<(start+sizes[i]);j++){
for(int k=start+sizes[i];k<c.size();k++){
if(c[j]<c[k])
Ustat+=1;
}
}
start+=sizes[i];
}
return(Ustat);
}

void allocate_to_groups_impl( std::vector<char>& usage, std::vector<int>& order, const           std::vector<int>& cumsum, int group, int place, int min )
{
std::vector<int> store(order.size());
std::vector<int> sizes(cumsum.size());
sizes[0]=cumsum[0];
for(int i=1;i<cumsum.size();i++)
sizes[i]=(cumsum[i]-cumsum[i-1]);
if (place == cumsum[group]) {
group++; min = 0;
if (group == cumsum.size()) {
int ij=0;
for( std::vector<int>::iterator it = order.begin(); it != order.end(); ++it ){
store[ij]=*it;
ij++;
}
finalcountdown[count]=mannwhit(store,sizes);
count++;
return;
}
}

for( int v = min, max = usage.size() + place - cumsum[group]; v <= max; ++v ) {
if (!usage[v]) {
order[place] = v;
usage[v] = 1;
allocate_to_groups_impl(usage, order, cumsum, group, place+1, v+1);
usage[v] = 0;
}
}
}

void allocate_to_groups( int *c, int Ngroups )
{
size_t sum_of_c = 0;
std::vector<int> cumsum_of_c;
for( int* it = c; it < c + Ngroups; ++it ){
cumsum_of_c.push_back(sum_of_c += *it);
}
std::vector<int> order(sum_of_c);
std::vector<char> usage(sum_of_c);
allocate_to_groups_impl(usage, order, cumsum_of_c, 0, 0, 0);
}

void jonckterp(int *x,int *c, int *ngrp, double *numcomb){
count=0;
float p=0;
int d[*ngrp];
for(int i=0; i< *ngrp; i++)
d[i]=c[i];
finalcountdown =tabInit(*numcomb);
allocate_to_groups(d,*ngrp);
for(int test=0;test<*numcomb;test++){
if(finalcountdown[test]>=*x){
p+=1;
}
}
*x=p;
}


wrapper.cpp

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
extern "C" {
void jonckterp(int *x,int *c, int *ngrp, double *numcomb);
}


jonckterp.R

early.fn<-function(stat,group.sizes){
numcomb<-factorial(sum(group.sizes))/prod(factorial(group.sizes))
.C("jonckterp",p=as.integer(stat),as.integer(group.sizes),as.integer(length(group.sizes)),as.double(numcomb))\$p/numcomb
}

• Use some profiling tools to find out where the most CPU is being used. Once you have found out that you can concentrate on that part of the code. Feb 25, 2012 at 20:04
• It would be easier if you described what it is supposed to do... Feb 25, 2012 at 20:23
• @ronag It's computing the distribution of a type of nonparametric statistic called the Jonckheere-Terpstra statistic. Without getting into too much detail, i have n1, n2,...nk groups of objects, and I need to compute every possible order that they could've appeared, then for each of those (*numcomb) orderings, I sum up the k(k-1)/2 mann whitney counts between every pair (i,j) groups such that 1<=i<j<=k.
– Gschneider
Feb 25, 2012 at 20:36
• After using gprof, it seems like mannwhit() could be an area to focus on. It looks like the program spends 15% of its time there, which makes sense, as it runs this on each of the rows of "finalcountdown", which is often a tens of millions. I have read that it's a good idea to have loops count down rather than up, so I've changed the outer two loops in mannwhit to do that, and have changed all of my "int"s to "unsigned int"s, which makes a small difference, but any other suggestions would be greatly appreciated. Feb 27, 2012 at 19:04
• @ronag and some spaces here and there might not go astray for readability. Feb 28, 2012 at 3:29

I didn't read in it's entirity, but change this line

static int mannwhit(std::vector<int> c, std::vector<int> sizes)


to

static int mannwhit(const std::vector<int>& c,
const std::vector<int>& sizes)


This won't speed up your code but you should use

#include <cstdio> //instead of stdio.h


and memset will need to be prefixed as std::memset.

Also ensure you compiler is compiling at 02 and with out debug flags.

• That is a micro optimization, meaning he will not notice it. Also the include is fine. As for the optimization, O3 would be better, but how do you know how is he compiling? Feb 25, 2012 at 20:22
• @VJovic: O3 isn't necessarily better than O2. It has more optimizations per se, but may increase the code size and can worsen performance by putting too much pressure on the instruction cache (with excessive inlining, etc.). It's better to test whether O3 is actually faster than O2 instead of assuming so.
– netcoder
Feb 25, 2012 at 20:26
• Yes, this code seems to perform at similar speeds to the original, and std::memset rather than memset gives me an error. Thanks for trying, though.
– Gschneider
Feb 25, 2012 at 20:32
• @VJovic coping a container is a heavy weight operation. stdio.h is C it ins't garunteed to work by the C++ standard., std c++ specify cstdio as the header. 03 is not better than 02 in all cases, 02 is recommended for most tasks. You should always compare 02 and 03 if you are going to use 03.
– 111111
Feb 25, 2012 at 23:37