Vectorisation API changed (broken?) for adaptive line integrals?

[burnpanck]

This is either a follow-up to #188 or a bug report: I have been using quadpy for some time now to compute line integrals of vectorised functions adaptively using quadpy.line_segment.integrate_adaptive. Previously working code seems now to be broken due to either a bug or an API change. If it was the latter, it would be helpful to document what the API is.

The wiki states that the function should expect point arrays of the shape (d, n, p) and should return values of shape (a, n, p) to compute a vectorised result of shape (a,). For line-integrals d would be equal to 1, and it seems that in this case, the input would actually be (n, p) (this is undocumented in the wiki) - at least that is what my code seems to expect. However, with quadpy 0.14.11 it actually seems receive inputs of shape (n*p,) (or at least, a one-dimensional array of inputs). Therefore, my code currently outputs (a, 1, n*p) shaped function values, and quadpy returns an output of shape (a,1) instead of the expected (a,).

Because the details of the vectorisation shape for line integrals are currently not documented, I can't tell if this was an unintended change (i.e. bug?) or an intended API change.

[nschloe]

This needs an MWE, too.

[burnpanck]

Well, half of it is a documentation issue. The documentation is not 100 % clear how the vectorised API is supposed to look, or at least was it slightly inconsistent with the actual behaviour observed (see final comment in #188). Below, you will find a highlight of the issue with vectorisation over the result types. There are three implementation of the integrand function: func_as_documented (strictly following the documentation with a leading dimensionality axis), func_v0_12 (works in v0.12) and func_v0_15 (accepts the arguments in v0.15). Feel free to try other combinations.

I also attempted to exercise vectorisation over the domain, though I didn't manage to get it done in the time available, as the documentation is again not clear how to do that for line integrals (maybe it's not supported for adaptive integrals at all?)

Also note that the very simple integrand we have here (it's a polynomial with a step at zero) never terminates in the given configuration on v0.15 in the third case (matrices). Maybe that is related to #319?

import numpy as np

import quadpy

print("QuadPy version: ",quadpy.__version__)

quadpy_version = tuple(int(v) for v in quadpy.__version__.split('.'))

if quadpy_version < (0,15,0):
    adaptive_line_segment = quadpy.line_segment.integrate_adaptive
else:
    adaptive_line_segment = quadpy.c1.integrate_adaptive

ndim = 1
kronrod_degree = 5

npoints = 2*kronrod_degree+1

def true_value(lo,hi,A0,A1):
    return np.sign(hi)*A0*hi-np.sign(lo)*A0*lo + A1/2 * (hi**2 - lo**2)

def func_as_documented(x,A0,A1):
    d, n, p = x.shape
    assert d == ndim, f"{d} == {ndim}"
    assert p == npoints, f"{p} == {npoints}"
    
    # append two dimensions to make room for n and p
    A0 = A0[...,None,None]
    A1 = A1[...,None,None]
    return A0*np.sign(x) + A1*x

def fun_v0_12(x,A0,A1):
    n, p = x.shape
#    assert d == ndim, f"{d} == {ndim}"
    assert p == npoints, f"{p} == {npoints}"
    
    # append two dimensions to make room for n and p
    A0 = A0[...,None,None]
    A1 = A1[...,None,None]
    return A0*np.sign(x) + A1*x

def fun_v0_15(x,A0,A1):
    n_times_p, = x.shape
    print(x.shape)
#    assert d == ndim, f"{d} == {ndim}"
#    assert p == npoints, f"{p} == {npoints}"
    
    # append two dimensions to make room for n and p
    A0 = A0[...,None]
    A1 = A1[...,None]
    return A0*np.sign(x) + A1*x

def test_line_adaptive(func,lo,hi,A0,A1,eps=1e-7):
    lo,hi = np.broadcast_arrays(lo,hi)
    A0,A1 = np.broadcast_arrays(A0,A1)
    
    domain = np.array([lo,hi])
    sh = np.broadcast_arrays(lo,A0)[0].shape
    print(
        f"Starting adaptive integration using \"{func.__name__}\""
        f" over domain of shape {domain.shape[1:]}"
        f" and values of shape {A0.shape},"
        f" expecting result of shape {sh}."
    )
    if quadpy_version < (0,15,0):
        eps_arg = dict(eps=eps)
    else:
        eps_arg = dict(eps_abs=eps,eps_rel=eps)
    val,error = adaptive_line_segment(
        lambda x: func(x,A0,A1),
        domain,
        kronrod_degree=kronrod_degree,
        **eps_arg,
    )
    assert val.shape == sh, f"{val.shape} == {sh}"
    exp = true_value(lo,hi,A0,A1)
#    print(val.flat[0],error.flat[0],exp.flat[0])
    assert np.all((val-error <= exp+eps) & (exp-eps <= val+error))

    

if quadpy_version < (0,13,0):
    good_func = fun_v0_12
else:
    good_func = fun_v0_15
    
# attempt a single integration
test_line_adaptive(good_func,-2,2,1,1)

# attempt a integration over vectors
test_line_adaptive(good_func,-2,2,1,np.r_[0,1])

# attempt a integration over matrices
value_shape = 3,4
test_line_adaptive(good_func,-2,2,1,np.arange(np.prod(value_shape)).reshape(value_shape))

# attempt a integration two domains
#test_line_adaptive(func_apparent_1,np.r_[-1,-2],np.r_[1,2],1,0)

[burnpanck]

I haven't run the code on 0.13 or 0.14, so the version checks to make it work might need small tweaking, if you want to test those too.

[nschloe]

I'm just plowing through some old bugs. Can you provide a shorter MWE?

[burnpanck]

This issue is about an undocumented API change between two versions of quadpy, both of them slightly mismatching what the documentation says. Therefore, it's up to you to decide if this is a bug or a feature. I am not actively using quadpy anymore, so I have no strong feelings toward either. The MWP I provided is to highlight the difference in the APIs. If you tell me which of the three APIs is correct (v0.12, v0.15 or "as documented"), then I can shorten the MWE to just that API. But apart from that, there is really not much in that API; a test function and a call to adaptive_line_integral.

[nschloe]

I see now. I'll look at the functions your provide and add a test such that accidental API changes like this won't occur anymore. Thanks for the examplary reports!