smooth_util.pyx 5.77 KB
Newer Older
csongor's avatar
csongor committed
1 2
import numpy as np
cimport numpy as np
csongor's avatar
csongor committed
3
#defining NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
csongor's avatar
csongor committed
4

csongor's avatar
csongor committed
5
cpdef np.float32_t gaussianKernel_f(np.ndarray[np.float32_t, ndim=1] mpower, np.ndarray[np.float32_t, ndim=1] mk, np.float32_t mu, np.float32_t smooth_length):
csongor's avatar
csongor committed
6 7
    cdef np.ndarray[np.float32_t, ndim=1] C = np.exp(-(mk - mu) ** 2 / (2. * smooth_length ** 2))
    return np.sum(C * mpower) / np.sum(C)
csongor's avatar
csongor committed
8

csongor's avatar
csongor committed
9
cpdef np.float64_t gaussianKernel(np.ndarray[np.float64_t, ndim=1] mpower, np.ndarray[np.float64_t, ndim=1] mk, np.float64_t mu, np.float64_t smooth_length):
csongor's avatar
csongor committed
10
    cdef np.ndarray[np.float64_t, ndim=1] C = np.exp(-(mk - mu) ** 2 / (2. * smooth_length ** 2))
csongor's avatar
csongor committed
11 12 13
    return np.sum(C * mpower) / np.sum(C)


csongor's avatar
csongor committed
14
cpdef np.ndarray[np.float64_t, ndim=1] apply_kernel_along_array(np.ndarray[np.float64_t, ndim=1] power, np.int_t startindex, np.int_t endindex, np.ndarray[np.float64_t, ndim=1] distances, np.float64_t smooth_length):
csongor's avatar
csongor committed
15

csongor's avatar
csongor committed
16 17
    if smooth_length == 0.0:
        return power[startindex:endindex]
csongor's avatar
csongor committed
18 19

    if (smooth_length is None) or (smooth_length < 0):
csongor's avatar
csongor committed
20
        smooth_length = distances[1]-distances[0]
csongor's avatar
csongor committed
21

csongor's avatar
csongor committed
22
    cdef np.ndarray[np.float64_t, ndim=1] p_smooth = np.empty(endindex-startindex, dtype=np.float64)
csongor's avatar
csongor committed
23 24 25 26
    cdef np.int i
    for i in xrange(startindex, endindex):
        l = max(i-int(2*smooth_length)-1,0)
        u = min(i+int(2*smooth_length)+2,len(p_smooth))
csongor's avatar
csongor committed
27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45
        p_smooth[i-startindex] = gaussianKernel(power[l:u], distances[l:u], distances[i], smooth_length)

    return p_smooth

cpdef np.ndarray[np.float32_t, ndim=1] apply_kernel_along_array_f(np.ndarray[np.float32_t, ndim=1] power, np.int_t startindex, np.int_t endindex, np.ndarray[np.float32_t, ndim=1] distances, np.float32_t smooth_length):

    if smooth_length == 0.0:
        return power[startindex:endindex]

    if (smooth_length is None) or (smooth_length < 0):
        smooth_length = distances[1]-distances[0]

    cdef np.ndarray[np.float32_t, ndim=1] p_smooth = np.empty(endindex-startindex, dtype=np.float32_t)
    cdef np.int i
    for i in xrange(startindex, endindex):
        l = max(i-int(2*smooth_length)-1,0)
        u = min(i+int(2*smooth_length)+2,len(p_smooth))
        p_smooth[i-startindex] = gaussianKernel_f(power[l:u], distances[l:u], distances[i], smooth_length)

csongor's avatar
csongor committed
46 47
    return p_smooth

csongor's avatar
csongor committed
48

csongor's avatar
csongor committed
49 50
def getShape(a):
    return tuple(a.shape)
csongor's avatar
csongor committed
51

csongor's avatar
csongor committed
52
cpdef np.ndarray[np.float64_t, ndim=1] apply_along_axis(np.int_t axis, np.ndarray arr, np.int_t startindex,  np.int_t endindex, np.ndarray[np.float64_t, ndim=1] distances, np.float64_t smooth_length):
csongor's avatar
csongor committed
53
    cdef np.int_t nd = arr.ndim
csongor's avatar
csongor committed
54 55 56 57 58
    if axis < 0:
        axis += nd
    if (axis >= nd):
        raise ValueError("axis must be less than arr.ndim; axis=%d, rank=%d."
            % (axis, nd))
csongor's avatar
csongor committed
59 60
    cdef np.ndarray[np.int_t, ndim=1] ind = np.zeros(nd-1, dtype=np.int)
    cdef np.ndarray i = np.zeros(nd, dtype=object)
csongor's avatar
csongor committed
61
    cdef tuple shape = getShape(arr)
csongor's avatar
csongor committed
62 63
    cdef np.ndarray[np.int_t, ndim=1] indlist = np.asarray(range(nd))
    indlist = np.delete(indlist,axis)
csongor's avatar
csongor committed
64
    i[axis] = slice(None, None)
csongor's avatar
csongor committed
65 66
    cdef np.ndarray[np.int_t, ndim=1] outshape = np.asarray(shape).take(indlist)

csongor's avatar
csongor committed
67 68
    i.put(indlist, ind)

csongor's avatar
csongor committed
69 70 71 72 73 74
    cdef np.ndarray[np.float64_t, ndim=1] slicedArr
    cdef np.ndarray[np.float64_t, ndim=1] res

    cdef np.int_t Ntot = np.product(outshape)
    cdef np.ndarray[np.int_t, ndim=1] holdshape = outshape
    slicedArr = arr[tuple(i.tolist())]
csongor's avatar
csongor committed
75 76 77

    res = apply_kernel_along_array(slicedArr, startindex, endindex, distances, smooth_length)

csongor's avatar
csongor committed
78 79 80
    outshape = np.asarray(getShape(arr))
    outshape[axis] = endindex - startindex
    outarr = np.zeros(outshape, dtype=np.float64)
csongor's avatar
csongor committed
81
    outarr[tuple(i.tolist())] = res
csongor's avatar
csongor committed
82
    cdef np.int_t n, k = 1
csongor's avatar
csongor committed
83 84 85 86 87 88 89 90 91
    while k < Ntot:
        # increment the index
        ind[-1] += 1
        n = -1
        while (ind[n] >= holdshape[n]) and (n > (1-nd)):
            ind[n-1] += 1
            ind[n] = 0
            n -= 1
        i.put(indlist, ind)
csongor's avatar
csongor committed
92
        slicedArr = arr[tuple(i.tolist())]
csongor's avatar
csongor committed
93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141
        res = apply_kernel_along_array(slicedArr, startindex, endindex, distances, smooth_length)
        outarr[tuple(i.tolist())] = res
        k += 1

    return outarr


cpdef np.ndarray[np.float32_t, ndim=1] apply_along_axis_f(np.int_t axis, np.ndarray arr, np.int_t startindex,  np.int_t endindex, np.ndarray[np.float32_t, ndim=1] distances, np.float32_t smooth_length):
    cdef np.int_t nd = arr.ndim
    if axis < 0:
        axis += nd
    if (axis >= nd):
        raise ValueError("axis must be less than arr.ndim; axis=%d, rank=%d."
            % (axis, nd))
    cdef np.ndarray[np.int_t, ndim=1] ind = np.zeros(nd-1, dtype=np.int)
    cdef np.ndarray i = np.zeros(nd, dtype=object)
    cdef tuple shape = getShape(arr)
    cdef np.ndarray[np.int_t, ndim=1] indlist = np.asarray(range(nd))
    indlist = np.delete(indlist,axis)
    i[axis] = slice(None, None)
    cdef np.ndarray[np.int_t, ndim=1] outshape = np.asarray(shape).take(indlist)

    i.put(indlist, ind)

    cdef np.ndarray[np.float32_t, ndim=1] slicedArr
    cdef np.ndarray[np.float32_t, ndim=1] res

    cdef np.int_t Ntot = np.product(outshape)
    cdef np.ndarray[np.int_t, ndim=1] holdshape = outshape
    slicedArr = arr[tuple(i.tolist())]

    res = apply_kernel_along_array_f(slicedArr, startindex, endindex, distances, smooth_length)

    outshape = np.asarray(getShape(arr))
    outshape[axis] = endindex - startindex
    outarr = np.zeros(outshape, dtype=np.float32_t)
    outarr[tuple(i.tolist())] = res
    cdef np.int_t n, k = 1
    while k < Ntot:
        # increment the index
        ind[-1] += 1
        n = -1
        while (ind[n] >= holdshape[n]) and (n > (1-nd)):
            ind[n-1] += 1
            ind[n] = 0
            n -= 1
        i.put(indlist, ind)
        slicedArr = arr[tuple(i.tolist())]
        res = apply_kernel_along_array(slicedArr, startindex, endindex, distances, smooth_length)
csongor's avatar
csongor committed
142 143
        outarr[tuple(i.tolist())] = res
        k += 1
csongor's avatar
csongor committed
144 145

    return outarr
csongor's avatar
csongor committed
146