[BUG] Weak test function evaluation incorrect for odd polynomial order

[Jacob-Stevens-Haas]

Putting a placeholder here; @znicolaou, you might be able to provide some context or this might solve a problem for you. It looks like the alternating signs of monomials in the odd expansion of _phi (and _phi_int and _xphi_int) are reversed. This shouldn't affect weak SINDy in general, however, since -phi is a valid test function.

Reproducing code example:

This is a pure-function adaptation of the code in weak SINDy. It'll eventually be merged:

def _phi(x: Float1D, d: int, p: int) -> Float1D:
    """
    One-dimensional polynomial test function (1-x**2)**p,
    differentiated d times, calculated term-wise in the binomial
    expansion.
    """
    ks = np.arange(p + 1)
    ks = ks[np.where(2 * (p - ks) - d >= 0)][:, np.newaxis]
    poly_coef = binom(p, ks)
    sign = (-1) ** ks
    x_pows = 2 * (p - ks) - d
    deriv_coef = perm(2 * (p - ks), d)
    monomials = sign * poly_coef * deriv_coef * x[np.newaxis, :] ** x_pows
    return np.sum(monomials, axis=0)

In order to understand how to re-implement weak, I've added a test (the full test is a bit longer and tests additional derivative orders):

@pytest.mark.parametrize("p", [2, 3, 4])
def test_test_function_phi(p):
    # Convince yourself that these are the correct derivatives
    # by differentiating (1-x^2)^p by hand.
    # Built-in phi uses vectorized operations, generic by derivative,
    # So these are easy to read but slower and manually defined by 
    # derivative order.
    def d0(x):
        return (1 - x**2)**p

    rt2 = np.sqrt(2)
    x = np.array([0, rt2 / 2, 1])
    expected = np.array([d0(0), d0(rt2/2), d0(1)])
    result = _phi(x, 0, p)
    assert_allclose(result, expected)

Error message:

All the values for p=3 are -1 * what's expected