`_x_hat_using_finite_difference` investigation

[pavelkomarov]

In the finite difference code lives a function that takes a first order FD and then uses scipy's trapezoidal integration immediately after. These two operations seem to cancel each other, although there will be some numerical difference due to the trapezoids. We should find out what this is actually accomplishing and properly document it, or if we determine it's not helpful remove it. @florisvb and I spoke about it briefly in review of #77, but he couldn't recall what it was for.

[pavelkomarov]

This appears to indeed be a smoothing procedure.

If I change the code to

def _x_hat_using_finite_difference(x, dt):
    """Find a smoothed estimate of the true function by taking FD and then integrating with trapezoids
    """
    from matplotlib import pyplot
    print("x", x)
    pyplot.plot(x, label='x')
    x_hat, dxdt_hat = first_order(x, dt)
    print("x_hat1", x_hat)
    pyplot.plot(x_hat, '--', label='x_hat1')
    x_hat = utility.integrate_dxdt_hat(dxdt_hat, dt)
    print("x_hat2", x_hat)
    pyplot.plot(x_hat, label='x_hat2')
    x0 = utility.estimate_initial_condition(x, x_hat)
    print("x_hat3", x_hat + x0)
    pyplot.plot(x_hat + x0, label='x_hat3')
    pyplot.legend()
    pyplot.show()
    return x_hat + x0

then I see plots like the following at the first iteration:

And then after a few more iterations:

It's just an iterative smoothing procedure. I may rename the function to highlight this.

[pavelkomarov]

Is there a more efficient way to achieve the same smoothing, or something spiritually equivalent? Trapezoidal integration is effectively doing some kind of convolution with a 2-wide moving average filter, plus some scaling. I'd rather not solve a linear inverse problem to determine the optimal x0 at every iteration.

[pavelkomarov]

I need to address this one as part of #120, since I'll be changing how iteration works in the module anyway. That will complete its refresh.

[pavelkomarov]

The way this is implemented isn't ideal. Finding the integration constant on every iteration is wasted effort, because we just differentiate it away in the next iteration. The way this is implemented is strange. There is averaging of forward and backward passes, which is extra work and actually unnecessary because first_order was a symmetrical second order method. I've settled on a common pattern to replace it all with:

if num_iterations < 1: raise ValueError("num_iterations must be >0")
    
x_hat = x # preserve a reference to x, because if iterating we need it to find the final constant of integration
dxdt_hat = np.zeros(x.shape) # preallocate memory

for i in range(num_iterations):
    # calculate the derivative
    dxdt_hat = # finite differencing core, changes per the order

    if i < num_iterations - 1: # if not the last iteration
        x_hat = utility.integrate_dxdt_hat(dxdt_hat, dt) # estimate new x_hat by integrating derivative
        # no need to find integration constant until the very end, because we just differentiate again

if num_iterations > 1: # We've lost a constant of integration in the above
    x_hat += utility.estimate_integration_constant(x, x_hat) # uses least squares

return x_hat, dxdt_hat

[pavelkomarov]

Okay, I'm managing to get essentially identical results with my structure, so this function indeed was not necessary, and I'm going to take it out as part of #120. Although there is another wrinkle involving #104. Go there for discussion.