Skip to content
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

add common interface bindings #119

Merged
merged 11 commits into from Dec 22, 2016
Merged

add common interface bindings #119

merged 11 commits into from Dec 22, 2016

Conversation

ChrisRackauckas
Copy link
Member

This adds common interface bindings to ODE.jl. Currently it only includes the algorithms in the README:

  • ode23Alg
  • ode45Alg
  • ode78Alg
  • ode23sAlg

I had to put the Alg in the type name because otherwise it matches the function's name. If anyone has a better idea of how to handle this let me know.

@coveralls
Copy link

coveralls commented Dec 12, 2016

Coverage Status

Coverage increased (+0.7%) to 93.803% when pulling 14eca21 on common into a11e3fe on master.

@coveralls
Copy link

coveralls commented Dec 12, 2016

Coverage Status

Coverage increased (+0.7%) to 93.803% when pulling 14eca21 on common into a11e3fe on master.

@@ -0,0 +1,99 @@
abstract ODEJLAlgorithm <: AbstractODEAlgorithm
immutable ode23Alg <: ODEJLAlgorithm end
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'd make those Ode23Alg, etc.


u0 = prob.u0

Ts = sort(unique([tspan[1];saveat;tspan[2]]))
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why do you sort those?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Two reasons. One is that it allows you to check for errors when saveat goes past tspan (I actually didn't include that error check, oops!). Also, I've had a request for this in the chats (a long time ago) by someone who wanted to save at 100 random timepoints. When I profiled everything, this sort never even showed up on the radar. But what should probably happen is it should be a sort which is very fast when the array is already sorted, but I don't know how to choose sort algorithms in Julia.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hmm, I'm not a fan of this: it adds unnecessary complexity. It's not much complexity but feature creep needs to be avoided actively. That user could easily have done that himself with one or two lines of code. Now we need to support X lines of code for this instead.

Also, this is wrong when integration is backward in time. Ah, I see, that is not allowed. Why not? ODE.jl solver can do it. (I guess I have my feature-creep priorities elsewhere...)

@@ -0,0 +1,21 @@
using ODE, DiffEqBase, DiffEqProblemLibrary
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can you test something more here than just throws-no-error? Maybe shape of solution? Within certain tolerance?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Indeed. That's (and missing some algorithms) is why I had it as WIP.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Guess I never did name it WIP, my bad. Should be complete now though, sans the ode23s issue which I think is beyond this PR.

@mauro3
Copy link
Contributor

mauro3 commented Dec 12, 2016

Thanks!

I think you can just drop Julia 0.4 (to make the travis error go away).

@kmsquire
Copy link

I had to put the Alg in the type name because otherwise it matches the function's name. If anyone has a better idea of how to handle this let me know.

Aren't functions also types now, as of Julia 0.5?

@ChrisRackauckas
Copy link
Member Author

Aren't functions also types now, as of Julia 0.5?

Ahh yes. I actually see that the problem was that I was putting the type declaration after the function declaration. If you have that the actual ode23 function adds a dispatch to the already existing ode23 function (from the constructor), it works:

immutable ode23 end

ode23(fn, y0, tspan; kwargs...) = 

So then the question is just whether you prefer it as ode23() to match the function or Ode23 since it's a type. With LSODA.jl and ODEInterface.jl I chose to match the function names: dopri5() and lsoda().

@ChrisRackauckas
Copy link
Member Author

This was updated to have the ode23() syntax and the tests check a bit more now. Notice that the tests on ode23s() were disabled: I think that might be an error in the actual function itself? It's only returning the initial timepoint.

maxstep=dtmax,
minstep=dtmin,
initstep=dt,
points=points)
end

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Again, if you're willing to give up the type declaration, and since all of the algorithms accept a kwargs... parameter, you could simplify the code greatly:

function solve{uType,tType,isinplace,F}(prob::AbstractODEProblem{uType,tType,isinplace,F}, ode_alg,...* rest of args *)
    < snip >
    ts,timeseries_tmp = ode_alg(f,u0,Ts,
                           norm = norm,
                           abstol=abstol,
                           reltol=reltol,
                           maxstep=dtmax,
                           minstep=dtmin,
                           initstep=dt,
                           points=points)

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Again, if you're willing to give up the type declaration

No, that would be a huge change and would have some big downsides for usage with things like Sundials.jl.

I could also just do the call via typeof(alg) since now they are the same name. The problem is that not all of the algorithms do accept kwargs parameters (ode4 is an example which does not). Maybe the better solution to make them all accept it?

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm curious why removing the type declaration would be a problem for usage with Sundials.jl? Even if this function doesn't have a type on the parameter, those solvers could be passed in. Or are you trying to prevent, e.g., a DAE algorithm from being passed into a ODE solver?

Either way, if there are only a few ode* algorithms which don't accept kwargs, then yes, I think giving all of them kwargs (and ignoring thes parameters if necessary) is worthwhile, since it let's you simplify this code considerably.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Oh I thought you meant type constructor. My bad for the confusion. No, the type declaration ode_alg::T with T<ODEJLAlgorithm is necessary for dispatch to work correctly, otherwise it doesn't know which solve function to use.

And sure, I'll change the others to use kwargs.

@coveralls
Copy link

Coverage Status

Coverage increased (+0.03%) to 93.151% when pulling d328716 on common into a11e3fe on master.

@coveralls
Copy link

coveralls commented Dec 12, 2016

Coverage Status

Coverage increased (+0.03%) to 93.151% when pulling d328716 on common into a11e3fe on master.

@coveralls
Copy link

coveralls commented Dec 13, 2016

Coverage Status

Coverage increased (+0.03%) to 93.151% when pulling d328716 on common into a11e3fe on master.

@ChrisRackauckas
Copy link
Member Author

I did the suggested changes. Please take another look. Thanks!

@coveralls
Copy link

Coverage Status

Coverage increased (+0.06%) to 93.182% when pulling 092cb63 on common into a11e3fe on master.

1 similar comment
@coveralls
Copy link

Coverage Status

Coverage increased (+0.06%) to 93.182% when pulling 092cb63 on common into a11e3fe on master.

Copy link

@kmsquire kmsquire left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

LGTM!

@codecov-io
Copy link

codecov-io commented Dec 13, 2016

Current coverage is 93.39% (diff: 93.93%)

Merging #119 into master will increase coverage by 0.26%

@@             master       #119   diff @@
==========================================
  Files             2          4     +2   
  Lines           320        348    +28   
  Methods           0          0          
  Messages          0          0          
  Branches          0          0          
==========================================
+ Hits            298        325    +27   
- Misses           22         23     +1   
  Partials          0          0          

Powered by Codecov. Last update a11e3fe...8af616f

@coveralls
Copy link

coveralls commented Dec 13, 2016

Coverage Status

Coverage increased (+0.2%) to 93.277% when pulling 092cb63 on common into a11e3fe on master.

@coveralls
Copy link

coveralls commented Dec 13, 2016

Coverage Status

Coverage increased (+0.06%) to 93.182% when pulling 243fac6 on common into a11e3fe on master.

@ChrisRackauckas
Copy link
Member Author

I'll wait for @mauro3 's approval.

@coveralls
Copy link

coveralls commented Dec 14, 2016

Coverage Status

Coverage increased (+0.06%) to 93.182% when pulling 3d7a211 on common into a11e3fe on master.

@ChrisRackauckas
Copy link
Member Author

Bump.

@ChrisRackauckas
Copy link
Member Author

I'll give @mauro3 a few more days.

@@ -16,6 +21,9 @@ export ode23s
# non-adaptive stiff:
export ode4s

# Common Interface
export ODEJLAlgorithm
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

ODEjlAlgorithm would be easier to read.

@@ -152,7 +160,7 @@ include("runge_kutta.jl")

# ODE_MS Fixed-step, fixed-order multi-step numerical method
# with Adams-Bashforth-Moulton coefficients
function ode_ms(F, x0, tspan, order::Integer)
function ode_ms(F, x0, tspan, order::Integer; kwargs...)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Are these kwargs... really necessary? It just seems like the abstraction is a bit leaky if they now need to accept throw away kwargs.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

These kwargs are necessary for the call

AlgType(f,u0,Ts;
                      norm = norm,
                      abstol=abstol,
                      reltol=reltol,
                      maxstep=dtmax,
                      minstep=dtmin,
                      initstep=dt,
                      points=points)

to not throw an error.

@@ -0,0 +1,8 @@
abstract ODEJLAlgorithm <: AbstractODEAlgorithm
immutable ode23 <: ODEJLAlgorithm end
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is a bit magical. How about a comment like:

Making the ODE-solver functions into types lets us dispatch on them.  Used in the DiffEqBase interface.


u0 = prob.u0

Ts = sort(unique([tspan[1];saveat;tspan[2]]))
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hmm, I'm not a fan of this: it adds unnecessary complexity. It's not much complexity but feature creep needs to be avoided actively. That user could easily have done that himself with one or two lines of code. Now we need to support X lines of code for this instead.

Also, this is wrong when integration is backward in time. Ah, I see, that is not allowed. Why not? ODE.jl solver can do it. (I guess I have my feature-creep priorities elsewhere...)

u0 = vec(prob.u0)
else
u0 = prob.u0
end
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

To save a few LOCs u0 = uType <: AbstractArray ? vec(prob.u0) : prob.u0

u0 = prob.u0
end

ts,timeseries_tmp = AlgType(f,u0,Ts;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Add a comment along the lines: Calling the solver

@mauro3
Copy link
Contributor

mauro3 commented Dec 22, 2016

Thanks Chris! Have you run your bigger test suites using this?

@ChrisRackauckas
Copy link
Member Author

Hmm, I'm not a fan of this: it adds unnecessary complexity. It's not much complexity but feature creep needs to be avoided actively. That user could easily have done that himself with one or two lines of code. Now we need to support X lines of code for this instead.

Also, this is wrong when integration is backward in time. Ah, I see, that is not allowed. Why not? ODE.jl solver can do it. (I guess I have my feature-creep priorities elsewhere...)

Ahh yes, I forgot ODE.jl can integrate backwards in time. I got rid of the sort and the error for this.

@ChrisRackauckas
Copy link
Member Author

I pushed an update which addressed most of what you mentioned. Tests on the nightly will fail because of an unrelated thing (macro hygiene change, which has been fixed, but a tag needs to happen somewhere else in the ecosystem for it to pass).

I got rid of the sorting and time direction check because indeed, one of the nice things about ODE.jl is that it can integrate backwards in time (and I'm hoping to try and get that in other places as well).

I don't know what you want to do about the kwargs... though: they are necessary for this style of implementation, unless you want to if AlgType <: ... and pull out each different case. We can subtype the algorithms a bit more and do that, but this is mostly a code "styling" issue so I'll just go with whatever makes others happy here.

Thanks Chris! Have you run your bigger test suites using this?

Yes, I have already been using it some places. I haven't used it in the full benchmarking suite yet because I am waiting for the PR49 changes to go through so I can update the benchmarks to include all of them. I can't actually run them simultaneously until the package is renamed, so I've just been holding off and doing other things in the meantime.

@coveralls
Copy link

coveralls commented Dec 22, 2016

Coverage Status

Coverage increased (+0.03%) to 93.151% when pulling 8af616f on common into a11e3fe on master.

@ChrisRackauckas ChrisRackauckas merged commit 8d4827b into master Dec 22, 2016
@mauro3
Copy link
Contributor

mauro3 commented Dec 22, 2016

LGTM

@ChrisRackauckas ChrisRackauckas deleted the common branch December 22, 2016 11:29
@mauro3
Copy link
Contributor

mauro3 commented Dec 22, 2016

I renamed JuliaODE/ODE.jl to PR49.jl

@pwl
Copy link

pwl commented Dec 22, 2016

Don't you need to rename the module too?

@ChrisRackauckas
Copy link
Member Author

Yes, the module needs to be renamed for it to be referenced correctly.

@pwl
Copy link

pwl commented Dec 22, 2016

Hm, if we are going to put an effort and rename the whole module anyway, shouldn't we think about a better name then PR49? We can discuss it in SciML/Roadmap#7.

@mauro3
Copy link
Contributor

mauro3 commented Dec 22, 2016

Sorry, I'll revert it.

@tkelman
Copy link
Contributor

tkelman commented Jan 12, 2017

You could also dispatch on typeof(ode45) and leave them as functions instead of types, now these have zero-argument constructors that don't serve much purpose

@ChrisRackauckas
Copy link
Member Author

In other cases they are used as functions in order to specify more type information. Some examples include:

The ODE.jl algorithms are too simple to need this (at least for right now), but they are setup in a way that is standardized because it would be more difficult to know to use ode45 sometimes and Tsit5() others. It's easier to just always have it have the constructor, which is what JuMP, Optim, and JuliaML do as well (so users should see some cross-Julia continuity in how algorithm choice is done).

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

7 participants