Allow d(x) and k(x, s) to be matrices

[JoshKarpel]

In #35 , support was added for multi-dimensional IDEs. However, it did not account for the case where d(x) and k(x,s) might mix the dimensions, only that c(y, x) might. That case mostly-silently worked because of how numpy handles addition broadcasting; multiplication broadcasting is more complex and it seems like the simplest solution is to hackily switch on the array size. Not really sure whether this is the "proper" way of doing things, but it passes the test case I cooked up for it. I suppose we could force d(x) and k(x, s) to be matrices if y_0(x) is multi-dimensional, but that would be annoying to write out in many cases.

[codecov]

Codecov Report

Merging #60 (14c8fe7) into master (b64590d) will increase coverage by 0.0%.
The diff coverage is 100.0%.

@@          Coverage Diff           @@
##           master     #60   +/-   ##
======================================
  Coverage    98.3%   98.4%           
======================================
  Files          11      11           
  Lines         308     316    +8     
  Branches       55      57    +2     
======================================
+ Hits          303     311    +8     
  Misses          2       2           
  Partials        3       3           

Impacted Files	Coverage Δ
tests/test_solver_against_analytic_solutions.py	100.0% <ø> (ø)
idesolver/idesolver.py	97.2% <100.0%> (+0.1%)	⬆️

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update b64590d...14c8fe7. Read the comment docs.

[JoshKarpel]

@nbrucy since you added the original multi-dimensional IDE code, would you mind taking a look at this and letting me know if anything looks off? Or if you know any tricks for handling the scalar-vs-matrix multiplication switch more smoothly...

[nbrucy]

I suggest to replace d.size by d.ndim

[nbrucy]

We can add the following test if d and k are vector valued.

 (  # equivalent to previous case, but with d(x) and k(x, s) as vectors
        IDESolver(
            x=np.linspace(0, 7, 100),
            y_0=[0, 1],
            c=lambda x, y: [0.5 * (y[1] + 1), -0.5 * y[0]],
            d=lambda x: [-0.5, -0.5],
            k=lambda x, s: [1, 1],
            f=lambda y: y,
            lower_bound=lambda x: 0,
            upper_bound=lambda x: x,
        ),
        lambda x: [np.sin(x), np.cos(x)],
    ),

[nbrucy]

Hello!

I don’t see better way of switching between scalar and matrix multiplications but have you considered the case when d and k are vector-valued?

I think this test case would work before the PR and won't now:

 (  # equivalent to previous case, but with d(x) and k(x, s) as vectors
        IDESolver(
            x=np.linspace(0, 7, 100),
            y_0=[0, 1],
            c=lambda x, y: [0.5 * (y[1] + 1), -0.5 * y[0]],
            d=lambda x: [-0.5,  -0.5],
            k=lambda x, s: [1, 1],
            f=lambda y: y,
            lower_bound=lambda x: 0,
            upper_bound=lambda x: x,
        ),
        lambda x: [np.sin(x), np.cos(x)],
    ),

@ will do the dot product while * will perform term by term multiplication (what we want in this case I think).

I guess a simple solution would be do the test on d and k on number of dimensions instead of the length (see my comments).

[nbrucy]

I suggest to replace k.size by k.ndim

[JoshKarpel]

Hello!

I don’t see better way of switching between scalar and matrix multiplications but have you considered the case when d and k are vector-valued?

I think this test case would work before the PR and won't now:

 (  # equivalent to previous case, but with d(x) and k(x, s) as vectors
        IDESolver(
            x=np.linspace(0, 7, 100),
            y_0=[0, 1],
            c=lambda x, y: [0.5 * (y[1] + 1), -0.5 * y[0]],
            d=lambda x: [-0.5,  -0.5],
            k=lambda x, s: [1, 1],
            f=lambda y: y,
            lower_bound=lambda x: 0,
            upper_bound=lambda x: x,
        ),
        lambda x: [np.sin(x), np.cos(x)],
    ),

@ will do the dot product while * will perform term by term multiplication (what we want in this case I think).

I guess a simple solution would be do the test on d and k on number of dimensions instead of the length (see my comments).

Thank you so much! I've implemented your suggestions.

[nbrucy]

Nice! Looks goods for me now.