-
-
Notifications
You must be signed in to change notification settings - Fork 199
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Lazy W operator support #443
Conversation
Codecov Report
@@ Coverage Diff @@
## master #443 +/- ##
==========================================
- Coverage 95.38% 95.25% -0.13%
==========================================
Files 78 78
Lines 9379 9483 +104
==========================================
+ Hits 8946 9033 +87
- Misses 433 450 +17
Continue to review full report at Codecov.
|
I will be using ImplicitEuler as my test solver as opposed to KenCarp3 as suggested in #416 (comment) because apparently KenCarp3 doesn' support problems with custom mass matrices at the moment. (Speaking of which, I should also remind myself to give a similar warning for the ExpRK methods.) The out-of-place version turns out to be fairly simple: all we need is add an additional clause in I'm still using the
|
@ChrisRackauckas I have a question regarding the If I understand it correctly, a transformed W operator is right, basically dividing the original W by I can probably add a |
This might be nice to keep in OrdinaryDiffEqExtendedTests.jl if not the main repo.
Yes, some methods are written in one way and others in the other way. The transformed version is simply better, but I think Rosenbrock23 hasn't upgraded to it yet. |
Note we will want to upgrade StochasticDiffEq.jl to get similar implicit handling soon after. We want to keep this part of them in sync. |
I guess it doesn't need a cache anyone could just build a cache into the linsolve? |
I took an attempt at the second route in #416 (comment) (i.e. migrate the The compromise I came up with is as follows: we keep Also, instead of dropping Simple test script for In-place ImplicitEuler using lazy W can be found here. Output:
I think the plan going forward should be something like this:
|
I think it's more natural to keep |
The upgrade is pretty much done. Be bold, make the change at least to ODEs. The change then needs to be propagated to SDEs and DDEs though. |
Newest version intended for SciML/DiffEqBase.jl#126. |
src/derivative_utils.jl
Outdated
@@ -78,7 +78,7 @@ internal cache (can be specified in the constructor; default to regular `Vector` | |||
It supports all of `AbstractDiffEqLinearOperator`'s interface. | |||
""" | |||
mutable struct WOperator{T, | |||
MType, | |||
MType <: Union{UniformScaling,AbstractMatrix}, |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
diffeqoperators are <: AbstractMatrix?
The previously removed concrete W calculation section in `calc_W!` is needed for problems whose function does not have a jacobian but instead calculates it via forward diff/finite difference using `calc_J!`
if isa(f, Union{SplitFunction, DynamicalODEFunction}) | ||
error("WOperator does not support $(typeof(f)) yet") | ||
end | ||
# Convert mass matrix, if needed |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
That's fine for now. I am hoping to support time-dependent mass matrices but we can just add an issue that it's a feature to implement in the future.
Base.size(W::WOperator, args...) = size(W.J, args...) | ||
function Base.getindex(W::WOperator, i::Int) | ||
if W.transform | ||
W.mass_matrix[i] / W.gamma - W.J[i] |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Since this shows up everywhere and it's constant given by the algorithm, it might make sense to make this a type parameter so the code can simply dispatch off of it, or at least not have large performance cuts because of this conditioning in indexing. I'm not sure this indexing is really used all that much though.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I'm not sure this indexing is really used all that much though.
This is true both for factorization (it indexes the concrete matrix instead of W) and mul!
(which is lazy).
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Sounds good
else | ||
J = ForwardDiff.derivative(uf,uprev) | ||
end | ||
W = mass_matrix - dtgamma*J |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
out of place isn't using WOperator yet?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
It's used in the branch where DiffEqBase.has_jac(f) == true
(https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/pull/443/files/21b1a5b7fae562c8f514cc445002453fa916da6c#diff-35ab888388acf2e7239e17f33e5eefe7R257 and https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/pull/443/files/21b1a5b7fae562c8f514cc445002453fa916da6c#diff-35ab888388acf2e7239e17f33e5eefe7R272).
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Why not all of the time?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I guess it's not necessary all of the time, and if we don't allow a linsolve option on it then it's pointless for now.
1. Fix global `fun` being polluted by complex_tests.jl 2. Add Compat to REQUIRE (used by iterator_tests.jl)
Finally got everything to pass. @ChrisRackauckas If everything so far looks good to you, then we can close this after the CI tests finish. I'll be working on StochasticDiffEq in the meantime. |
test/differentiation_traits_tests.jl
Outdated
@@ -26,7 +26,7 @@ using OrdinaryDiffEq, Test | |||
nothing | |||
end | |||
|
|||
Lotka_f = ODEFunction(Lotka,jac=Lotka_jac,tgrad=Lotka_tgrad) | |||
Lotka_f = ODEFunction(Lotka,jac=Lotka_jac,jac_prototype=zeros(2,2),tgrad=Lotka_tgrad) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
but nothing
still works?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Unfortunately it doesn't now.
Basically, previously we're assuming that jac_prototype
is a dense matrix (allocated in alg_cache
). However, now that it can in theory be everything but can't be simply determined by the constructor of ODEFunction
, we need to explicitly provide one.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
we need to have this work to keep the user interface clean. There should just be one spot where there's a hook that nothing
means the same as dense.
test/REQUIRE
Outdated
@@ -4,3 +4,4 @@ ParameterizedFunctions 3.0.0 | |||
StaticArrays | |||
DiffEqCallbacks | |||
DiffEqOperators | |||
Compat |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
where is this used?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
iterator_test.jl
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
src/derivative_utils.jl
Outdated
@@ -240,6 +240,32 @@ function calc_W!(integrator, cache::OrdinaryDiffEqMutableCache, dtgamma, repeat_ | |||
else | |||
new_W = false | |||
end | |||
else # concrete W using jacobian from `calc_J!` |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This branch is taken when?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
It's when f neither has invW nor jac, in which case the jacobian is computed using ForwardDiff/finite difference via calc_J!.
The WOperator implementation doesn't yet include this type of jacobian, which is something we could improve on in the future (might be difficult though because jacobian via differencing uses u and uf, and is different than the (J,u,p,t) interface for diffeq operators.
An integration test where an ODE is integrated with the mass matrix being a sparse matrix is required to show that this all captures the new functionality. |
You're right. I can modify the temporary testing script I wrote a while ago for this purpose. |
test/utility_tests.jl
Outdated
fun2 = ODEFunction(_f; mass_matrix=mm, jac=(u,p,t) -> t*A) | ||
fun1_ip = ODEFunction(_f_ip; mass_matrix=mm) | ||
fun2_ip = ODEFunction(_f_ip; mass_matrix=mm, | ||
jac_prototype=DiffEqArrayOperator(similar(A); update_func=(J,u,p,t) -> (J .= t .* A; J))) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
extra tab
Looks great. Merge this when the tests pass. |
About the integration test: In the temporary test script, I compared solution using 1) lazy W (by providing a jac function/jac prototype) and 2) ForwardDiff (by providing an I have also used a time-dependent jacobian so that the update functionality can be checked. As per your suggestion, sparse mass matrix and jacobian are used. One caveat is that I have to manually set the algorithm's linsolve to |
https://github.com/JuliaDiffEq/DiffEqBase.jl/blob/master/src/linear_nonlinear.jl#L14-L17 We can make the default be |
Implements the lazy W operators as suggested in #416 (comment).
I decided to write
WOperator
in OrdinaryDiffEq.jl so that DiffEqOperators is not required for OrdinaryDiffEq. The fact that the interface forAbstractDiffEqLinearOperator
in DiffEqBase helps a lot.(By the way, it occurs to me that a lot of methods in common_defaults.jl for DiffEqOperators should, in fact, be written in DiffEqBase instead. I should do a separate PR for this in the future).
Currently the lazy part of
WOperator
is already in place and what remains is to modifyalg_cache
andcalc_W!
to use it. (Also I need to migratemass_matrix
to*DEFunction
in DiffEqBase). The suggestions here #416 (comment) is also worth looking into.I also intend to add more tests for the methods in
derivative_utils.jl
toutility_tests.jl
, starting withcalc_J!
andcalc_W!
.