Epanechnikov kernel for use in Kernel Density Estimation. Note, for large arrays, performance may be slow. Originally designed for sparse grids.

epanechnikov.pyx

cimport cython
import numpy as np
cimport numpy as np

cdef extern from 'math.h':
    float exp(float x)
    float cos(float x)
    float sin(float x)
    float fabs(float x)

DTYPE64 = np.float64
ctypedef np.float64_t DTYPE64_t

@cython.boundscheck(False)
@cython.cdivision(True)
def epanechnikov(np.ndarray[DTYPE64_t, ndim=2] data, float bandwidth,
                 float dx, sig_as_grid_points=False):              
    '''
    The Epanechnikov kernel for use in Kernel Density Estimation.
    
    Note: Originally designed for sparse grids. Performance may be slow
          if used for large, dense grids.
    
    Parameters
    ----------
    data : np.float64 array
        The spatial data to "smooth"
    bandwidth : float
        The bandwidth to use in the Epanechnikov kernel.
    dx : float
        The distance between grid points in both the x and y direction
    sig_as_grid_points : bool
        A flag to tell the function if the bandwidth is given in grid points (True)
        or in physical distance (False). If bandwidth is given in physical space,
        the bandwidth is divided by dx to convert it to grid points.
        
    Returns
    -------
    A smoothed 2D array of type np.float64
    
    '''

    cdef unsigned int vlength = data.shape[0]
    cdef unsigned int ulength = data.shape[1]
    cdef unsigned int ng, nx, ny, nw
    cdef int iw, ie, js, jn, ngn
    cdef float sig_sq, dist_sq, ng_sq
    cdef Py_ssize_t i, j, ii, jj, nxx, nyy

    if not sig_as_grid_points:
        bandwidth = bandwidth / dx

    ng = int(bandwidth)
    ng_sq = float(ng * ng)
    ngn = -1 * ng
    nx = 2*ng+1
    ny = 2*ng+1

    cdef np.ndarray[DTYPE64_t, ndim=1] partweight = np.zeros([nx*ny], dtype=DTYPE64)
    cdef np.ndarray[DTYPE64_t, ndim=2] frc_data = np.zeros([vlength, ulength], dtype=DTYPE64)

    nw = -1
    for nyy in range(ngn, ng+1):
        for nxx in range(ngn, ng+1):
            nw = nw+1
            dist_sq = float(nxx*nxx) + float(nyy*nyy)
            if dist_sq <= ng_sq:
                partweight[nw] = 0.75 * (1 - (dist_sq / ng_sq)) * 1/bandwidth**2

    for j in range(0, vlength):
        for i in range(0, ulength):
            if data[j,i] > 0:
                iw=i-ng
                ie=i+ng
                js=j-ng
                jn=j+ng
                nw = -1
                for jj in range(js, jn+1):
                    for ii in range(iw, ie+1):
                        nw += 1
                        if jj < 0 or jj >= vlength or ii < 0 or ii >= ulength: continue
                        frc_data[jj,ii] = frc_data[jj,ii] + data[j,i]*partweight[nw]
    return frc_data