Where is the bottleneck in my Cython code?

The profile tells me it took ~15s to run, but without telling me more.

Tue Aug 19 20:55:38 2014 Profile.prof

3 function calls in 15.623 seconds

Ordered by: internal time

ncalls tottime percall cumtime percall filename:lineno(function)
1 15.623 15.623 15.623 15.623 {singleLoan.genLoan}
1 0.000 0.000 15.623 15.623 :1()
1 0.000 0.000 0.000 0.000 {method 'disable' of '_lsprof.Profiler' objects}


`import numpy as np
cimport numpy as np
from libc.stdlib cimport malloc, free
from libc.stdlib cimport rand, srand, RAND_MAX
import cython
cimport cython
import StringIO

cdef extern from "math.h":
int floor(double x)
double pow(double x, double y)
double exp(double x)

cdef double[:] zeros = np.zeros(360)
cdef double[:] stepCoupons = np.array([2.0,60,3.0,12.0,4.0,12.0,5.0])
cdef double[:,:] zeros2 = np.empty(shape=(999,2))

paraC2P = StringIO.StringIO('''1 10 0 0 (Intercept) 0 -4.981792
1 10 0 0 lv 50 0.55139
1 10 0 0 lv 51 0.53667
1 10 0 0 lv 52 0.52194
1 10 0 0 lv 53 0.50722
1 10 0 0 lv 54 0.49249
1 10 0 0 lv 55 0.47776
1 10 0 0 lv 56 0.46301
1 10 0 0 lv 57 0.44825
1 10 0 0 lv 58 0.43347
1 10 0 0 lv 59 0.41867
1 10 0 0 lv 60 0.40384
1 10 0 0 lv 61 0.38897
1 10 0 0 lv 62 0.37405
1 10 0 0 lv 63 0.35908
1 10 0 0 lv 64 0.34406
1 10 0 0 lv 65 0.32897
1 10 0 0 lv 66 0.31381
1 10 0 0 lv 67 0.29856
1 10 0 0 lv 68 0.28322
1 10 0 0 lv 69 0.26778
1 10 0 0 lv 70 0.25224
1 10 0 0 lv 71 0.23657
1 10 0 0 lv 72 0.22078
1 10 0 0 lv 73 0.20486
1 10 0 0 lv 74 0.18879
1 10 0 0 lv 75 0.17258
1 10 0 0 lv 76 0.1562
1 10 0 0 lv 77 0.13966
1 10 0 0 lv 78 0.12294
1 10 0 0 lv 79 0.10604
1 10 0 0 lv 80 0.08896
1 10 0 0 lv 81 0.07167
1 10 0

There's a lot of code here, so I'm just going to review interp2d.

-
There's no docstring. What does this function do? How am I supposed to call it? Are there any constraints on the parameters (for example, does the x array need to be sorted)?.

-
The function seems to be misnamed: it interpolates the function that takes x to y, but this is a function of one argument, so surely interp1d or just interp would be better names. (Compare numpy.interp, which does something very similar to your function, but is documentated as "one-dimensional linear interpolation".)

-
The parameter ex is opaquely named. What does it mean? Reading the code, it seems that it controls whether or not to extrapolate for values of x outside the given range. So it should be named extrapolate and it should be a Boolean (True or False) not a number.

-
Python allows indexing from the end of an array, so you can write x[-1] for the last element of an array, instead of the clusmy x[x.shape[0] - 1].

-
You find the interval in which to do the interpolation by a linear search over the array x, which takes time proportional to the size of x and so will be slow when x is large. Since x needs to be sorted in order for this algorithm to make sense, you should use binary search (for example, numpy.searchsorted) so that the time taken is logarithmic in the size of x.

-
You have failed to vectorize this function. The whole point of NumPy is that it provides fast operations on arrays of fixed-size numbers. If you find yourself writing a function that operates on just one value at a time (here, the single value new_x) then you probably aren't getting much, if any, benefit from NumPy.

-
Putting all this together, I'd write something like this:

def interp1d(x, y, a, extrapolate=False):
    """Interpolate a 1-D function.

    x is a 1-dimensional array sorted into ascending order.
    y is an array whose first axis has the same length as x.
    a is an array of interpolants.
    If extrapolate is False, clamp all interpolants to the range of x.

    Let f be the piecewise-linear function such that f(x) = y.
    Then return f(a).

    >>> x = np.arange(0, 10, 0.5)
    >>> y = np.sin(x)
    >>> interp1d(x, y, np.pi)
    0.0018202391415163

    """
    if not extrapolate:
        a = np.clip(a, x[0], x[-1])
    i = np.clip(np.searchsorted(x, a), 1, len(x) - 1)
    j = i - 1
    xj, yj = x[j], y[j]
    return yj + (a - xj) * (y[i] - yj) / (x[i] - xj)

(The rest of the program may need to be reorganized to pass arrays of values instead of one value at a time.)

-
Finally, what was wrong with scipy.interpolate.interp1d? It doesn't provide quite the same handling of points outside the range, but you could call numpy.clip yourself before calling it.

StackExchange Code Review Q#60540, answer score: 4