[SciPy-User] weave.inline version of scipy.stats.rankdata

[Keith Goodman]

scipy.stats.rankdata looks like an easy function to speed up with
weave.inline. Well, easy if I knew what I was doing. My attempt below
sometimes segfaults and sometimes gives the wrong result on the first
run (but gives the right result on subsequent runs). Can anyone spot
what I am doing wrong? (A note about the code: I used "if...else if"
since ||, or, wouldn't compile.)

import numpy as np
from scipy import weave
from scipy.weave import converters

def rankdata(a):

    a  = np.asarray(a, dtype=float)
    n = len(a)
    svec, ivec = fastsort(a)
    newarray = np.zeros(n, float)

    code = """
            int dupcount;
            float sumranks, averank;
            sumranks = 0.0;
            dupcount = 0;
            for (int i=0; i < n; ++i) {
                sumranks = sumranks + i;
                dupcount = dupcount + 1;
                if (i==n-1) {
                    averank = sumranks / float(dupcount) + 1;
                    for (int j=i-dupcount+1; j < i+2; ++j) {
                        newarray[ivec[j]] = averank;
                    }
                    sumranks = 0;
                    dupcount = 0;
                }
                else if (svec[i] != svec[i+1]) {
                    averank = sumranks / float(dupcount) + 1;
                    for (int j=i-dupcount+1; j < i+2; ++j) {
                        newarray[ivec[j]] = averank;
                    }
                    sumranks = 0;
                    dupcount = 0;
                }
            }
           """
    err = weave.inline(code,
                       ['svec', 'ivec', 'n', 'newarray'],
                       type_converters=converters.blitz,
                       compiler='gcc')


    return newarray


def fastsort(a):
    # fixme: the wording in the docstring is nonsense.
    """Sort an array and provide the argsort.

    Parameters
    ----------
    a : array

    Returns
    -------
    (sorted array,
     indices into the original array,
    )
    """
    it = np.argsort(a)
    as_ = a[it]
    return as_, it