ISISNeutronMuon / MDANSE

MDANSE: Molecular Dynamics Analysis for Neutron Scattering Experiments
https://www.isis.stfc.ac.uk/Pages/MDANSEproject.aspx
GNU General Public License v3.0
20 stars 4 forks source link

[BUG] Velocity interpolation order 1 and 2 are equivalent #369

Closed ChiCheng45 closed 1 month ago

ChiCheng45 commented 5 months ago

Description of the error Velocity interpolation order 1 and 2 are equivalent. This occurs because of the order=1 and order=2 numerical derivatives are equivalent in MDANSE.

see MDANSE/Mathematics/Signal.py

    if order == 1:
        ts[0] = np.add.reduce(coefs[0, :] * a[:3])
        ts[-1] = np.add.reduce(coefs[2, :] * a[-3:])

        gj = a[1:] - a[:-1]
        ts[1:-1] = gj[1:] + gj[:-1]
        # the equivalence is clear if i rewrite the above
        # gj = a[1:] - a[:-1]
        # ts[1:-1] = gj[1:] + gj[:-1]
        #          = (a[1:] - a[:-1])[1:] + (a[1:] - a[:-1])[:-1]
        #          = a[1:][1:] - a[:-1][1:] + a[1:][:-1] - a[:-1][:-1]
        #          = a[2:] - a[:-2]

        fact /= 2.0

    # Case of the order 2
    elif order == 2:
        ts[0] = np.add.reduce(coefs[0, :] * a[:3])
        ts[-1] = np.add.reduce(coefs[2, :] * a[-3:])

        gj = np.zeros((a.size - 2, 3), dtype=np.float64)
        gj[:, 0] = coefs[1, 0] * a[:-2]  # coefs[1, 1] = -1
        gj[:, 1] = coefs[1, 1] * a[1:-1] # coefs[1, 1] = 0
        gj[:, 2] = coefs[1, 2] * a[2:]   # coefs[1, 1] = 1
        ts[1:-1] = np.add.reduce(gj, -1)

        fact /= 2.0

Suggested fix Remove the order=1 numerical derivatives or add an actual first order derivative method either forward or backward.

Additional details Here is what you see if you run VACF with order=1, order=2 or order=3. order=1 and 2 are exactly the same while order=3 is slightly different.

order=1 image

order=2 image

order=3 image

ChiCheng45 commented 4 months ago

Also the documentation for the 1st order velocity is not the same as what is done in the code above. The code and/or documentation needs to be corrected.

image