Skip to content

x_hat_b and x_hat_f are the same in convolutional_smoother #129

@pavelkomarov

Description

@pavelkomarov

Here is the current utility.convolutional_smoother code:

def convolutional_smoother(x, kernel, iterations=1):
    """Perform smoothing by convolving x with a kernel forward and backward and then combining

    :param np.array[float] x: 1D data
    :param np.array[float] kernel: kernel to use in convolution
    :param int iterations: number of iterations, >=1
    
    :return: **x_hat** (np.array[float]) -- smoothed x
    """
    x_hat = np.hstack((x[::-1], x, x[::-1])) # pad
    w = np.linspace(0, 1, len(x_hat)) # weights

    for _ in range(iterations):
        x_hat_f = np.convolve(x_hat, kernel, 'same')
        x_hat_b = np.convolve(x_hat[::-1], kernel, 'same')[::-1]
        
        x_hat = x_hat_f*w + x_hat_b*(1-w)

    return x_hat[len(x):len(x)*2]

But try this:

>>> import numpy as np
>>> x = np.arange(10)
>>> kernel = np.array([1/3, 1/3, 1/3])
>>> np.convolve(x, kernel, 'same')
array([0.33333333, 1.        , 2.        , 3.        , 4.        ,
       5.        , 6.        , 7.        , 8.        , 5.66666667])
>>> np.convolve(x[::-1], kernel, 'same')[::-1]
array([0.33333333, 1.        , 2.        , 3.        , 4.        ,
       5.        , 6.        , 7.        , 8.        , 5.66666667])
>>> x = np.random.randn(10)
>>> x
array([ 1.30231522,  1.28738889,  0.82759237, -0.41610496, -2.53010508,
        1.53911992,  0.1855655 ,  0.59859514, -0.48033393,  0.41225268])
>>> np.convolve(x, kernel, 'same')
array([ 0.8632347 ,  1.13909883,  0.5662921 , -0.70620589, -0.46903004,
       -0.26847322,  0.77442685,  0.10127557,  0.17683796, -0.02269375])
>>> np.convolve(x[::-1], kernel, 'same')[::-1]
array([ 0.8632347 ,  1.13909883,  0.5662921 , -0.70620589, -0.46903004,
       -0.26847322,  0.77442685,  0.10127557,  0.17683796, -0.02269375])
>>> from pynumdiff.utils.utility import friedrichs_kernel
>>> x = np.random.randn(100)
>>> a = np.convolve(x, kernel, 'same')
>>> b = np.convolve(x[::-1], kernel, 'same')[::-1]
>>> np.max(np.abs(a - b))
np.float64(3.1138286393783687e-16)

So essentially the code is doing more work than necessary. There is no need to do forward and backward and then weight them together.

Metadata

Metadata

Assignees

Labels

simplificationunifying, shortening, and cleaning tasks that make the modules and user interface more cohesive

Projects

No projects

Milestone

No milestone

Relationships

None yet

Development

No branches or pull requests

Issue actions