-
Notifications
You must be signed in to change notification settings - Fork 385
[BUG] Weak test function evaluation incorrect for odd polynomial order #651
Open
Description
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
Activity
Metadata
Metadata
Assignees
Labels
No labels