In [1]:

    

import sys 
sys.path.append('..')


    

In [2]:

    

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline


    

In [3]:

    

%load_ext cython


    

In [4]:

    

def sparseconv(v, a):
    """
    Returns the discrete, linear convolution of two one-dimensional sequences.
    output is of length len(v) + len(a) -1 (same as the default for numpy.convolve)

    v is typically the kernel, a is the input to the system

    Can run faster than numpy.convolve if:
    (1) a is much longer than v
    (2) a is sparse (has lots of zeros)
    """
    v_len = v.shape[-1]
    a_len = a.shape[-1]
    out = np.zeros(a_len +  v_len - 1)

    pos = np.where(a != 0)[0]
    # add shifted and scaled copies of v only where a is nonzero
    for p in pos:
        out[p:p + v_len] = out[p:p + v_len] + v * a[p]

    return out


    

In [5]:

    

%%cython
import cython
import numpy as np
cimport numpy as np

@cython.boundscheck(False)
@cython.wraparound(False)
@cython.profile(False)
def sparseconvx(np.ndarray[double, ndim=1] v, np.ndarray[double, ndim=1] a):
    """
    Returns the discrete, linear convolution of two one-dimensional sequences.
    output is of length len(v) + len(a) -1 (same as the default for numpy.convolve)

    v is typically the kernel, a is the input to the system

    Can run faster than numpy.convolve if:
    (1) a is much longer than v
    (2) a is sparse (has lots of zeros)
    """
    cdef int v_len = v.shape[0]
    cdef int a_len = a.shape[0]
    cdef out = np.zeros(a_len +  v_len - 1)

    pos = np.where(a != 0)[0]
    # add shifted and scaled copies of v only where a is nonzero
    for p in pos:
        out[p:p + v_len] = out[p:p + v_len] + v * a[p]

    return out


    

In [6]:

    

from numba import jit
sparseconvj = jit(sparseconv)


    

In [7]:

    

import electrode2currentmap as e2cm
import effectivecurrent2brightness as ec2b


    

In [8]:

    

stimulus = e2cm.Stimulus()


    

In [9]:

    

tau1 = .42/1000


    

In [10]:

    

t = np.arange(0, 8 * tau1, stimulus.tsample)


    

In [11]:

    

g = ec2b.gamma(1, tau1, t)


    

In [12]:

    

%%timeit 
sparseconv(g, stimulus.data)


    

    




1000 loops, best of 3: 1.14 ms per loop

In [13]:

    

%%timeit 
sparseconvx(g, stimulus.data)


    

    




1000 loops, best of 3: 1.22 ms per loop

In [14]:

    

%%timeit
np.convolve(g, stimulus.data)


    

    




100 loops, best of 3: 10.2 ms per loop

In [15]:

    

from scipy.signal import fftconvolve


    

In [16]:

    

%%timeit
fftconvolve(g, stimulus.data)


    

    




100 loops, best of 3: 7.25 ms per loop

In [17]:

    

%%timeit
sparseconvj(g, stimulus.data)


    

    




The slowest run took 1561.47 times longer than the fastest. This could mean that an intermediate result is being cached 
1 loops, best of 3: 229 µs per loop

In [ ]: