Where communities thrive


  • Join over 1.5M+ people
  • Join over 100K+ communities
  • Free without limits
  • Create your own community
People
Repo info
Activity
  • Dec 28 2016 09:12
    coveralls commented #28
  • Dec 28 2016 09:12
    coveralls commented #28
  • Dec 28 2016 09:12
    coveralls commented #28
  • Dec 28 2016 09:12
    coveralls commented #28
  • Dec 28 2016 09:12
    codecov-io commented #28
  • Dec 28 2016 09:07
    codecov-io commented #28
  • Dec 28 2016 09:05
    coveralls commented #28
  • Dec 28 2016 08:59
    ChrisRackauckas commented #28
  • Dec 28 2016 08:58
    ChrisRackauckas synchronize #28
  • Dec 28 2016 08:55
    ChrisRackauckas synchronize #28
  • Dec 28 2016 08:54
    coveralls commented #28
  • Dec 28 2016 08:54
    coveralls commented #28
  • Dec 28 2016 08:54
    coveralls commented #28
  • Dec 28 2016 08:54
    coveralls commented #28
  • Dec 28 2016 08:54
    codecov-io commented #28
  • Dec 28 2016 08:43
    ChrisRackauckas synchronize #28
  • Dec 22 2016 05:11
    ChrisRackauckas commented #28
  • Dec 12 2016 06:35
    coveralls commented #28
  • Dec 12 2016 06:35
    coveralls commented #28
  • Dec 12 2016 06:35
    codecov-io commented #28
Paweł Biernat
@pwl
some of these solvers should be a part of ODE, would that be feasible for you to use them through our API? Are there any downsides/problems that you would expect with this approach? Any improvements that we should make to our API?
I've got to make myself a dinner, I might be unavailable for a while
Christopher Rackauckas
@ChrisRackauckas
Okay.
To be honest, I found your API kind of confusing so I don't really know what the downsides/problems would be. I know a lot of the ODE solvers could transfer directly. In fact, I told @mauro3 about the Dormand-Prince 7/8 tableau, and the tableaus pretty much tranfer directly (so when I implement the Verner tableaus ChrisRackuackas/DifferentialEquations.jl#28 , I'll let you know about these as well). The implicit methods aren't hard to do, you'd just have to accept a dependency (I chose to do so to NLsolve and ForwardDiff).
I mean, it didn't take more than 6 hours to implement the lot of them (except for the Feagin methods), so it surely wouldn't be hard to port them over.
Christopher Rackauckas
@ChrisRackauckas
One think I am not sure of is, is there a way to not save the values at each timestep in your API?
Also, when I was previously looking at your API, it seemed like it had restrictions on what kinds of numbers could be used. And it still seems to. I am trying to make sure everything I have can plug in with Arbs (ArbFloats/ArbReals), DecFP, etc.
Plus, although it's messier, I like to special case and unroll a bunch of things so that it's more efficient. It seems like a lot what you guys have going on uses different tableau types which is great for implementing, but can be limiting in terms of performance.
That said, there's nothing about the ODE methods that makes it so they can't work somewhere else.
Paweł Biernat
@pwl
I don't know what you mean by your last point, but I will address your other points
1) Not saving the values is one of the feats of iterators, ideally a solver wouldn't allocate any memory during iteration, after it has been initialized. There might still be some imperfections though.
2) The design principle behind our solvers is to have them work on generic data types, with some minimal requirements for supported operations (+, eps, etc.). I notoriously use BigFloat in my work and I tested most of the iterative solvers with them. We don't have tests for that in ODE.jl yet, but I tested the rk solvers with types that don't mix with floats to test if the solvers are not polluted with some sneaky floating point constants. The solvers with the old API (the calls via ode45 etc.) might not have this feature, but they will once we bring them up to standards and re-implement them as iterators (we will have tests for supporting arbitrary values).
3) Unrolling stuff is not an issue, you could have a custom integrator, called RKWhateverUnrolled and make it completely unrelated to our RKIntegrator integrator. You wouldn't even have to make it a part of ODE.jl package, it would be enough if you implement the right backend interface so other people could use your solvers via our API and compare them with standard ones.
4) sorry if the API is confusing (is that the backend that you are talking about?). We are still working on better docs (some are already available here).
anyway, you should probably wait with any decisions to after we merge with the official ODE.jl
I am just point out the direction in which we are heading.
Ad 1) when you are using collect you are explicitly saving the results (it does a copy of y and dy at each time step), if you only want a save some results just do a for (t,y) in sol loop and save whatever you like. We might also implement stuff like last(sol) to facilitate extracting the end result of the integration.
Paweł Biernat
@pwl
still, the iterator functionality will be limited due to the lack of length. @mauro3 actually had a great idea to implement length for dense output (it has prescribed number of output steps), we will be testing it soon.
sorry for the docs link confusion, I fixed it now
Ok, I've got to go for now
thanks for all the comments!
Christopher Rackauckas
@ChrisRackauckas
Aight peace out.
Mauro
@mauro3

@ChrisRackauckas: concerning (3), you could make it a RKIntegrator and then just overload the onstep! etc functions to use the unrolled stuff. For the user it would then be just like the other RK-solvers, just faster.

@pwl: I love the docs! Tnx!

Also, looking through both your and my PR (#19 & #17), I have the feeling we should make tstop part of Problem. It features in the options everywhere, so it should be in the most central place. Then tsteps would be part of the options for fixed-step solvers (with @assert tsteps[end]==tstop), adaptive solvers wouldn't need any time, and dense output has tout (with @assert tout[end]<=tstop).

Also, then tdir(ivp).
Mauro
@mauro3
Paweł Biernat
@pwl
@mauro3 do you mean to move tdir to Problem or IVP? On one hand, there are common options showing up in each solver/integrator, on the other hand keeping all options, without exceptions, in the solver/integrator makes our types isolated and keeps the communication between them to minimum.
but yeah, this means that tdir might be more painful to use and implement but I would like to see it as a part of an internal implementation of a solver. Not everyone will want to use tdir in exactly the same form, some solvers might even choose not implement the reverse time integration at all.
Mauro
@mauro3

Yes, sorry, I should have written tdir(prob), i.e. tdir(prob::Problem) = sign(prob.tstop-prob.ivp.t0).

Integrators/solvers not implementing reverse integration could check in their constructor that tdir>0.

Paweł Biernat
@pwl
this is a viable option, but tdir would be extremely lonely in Problem, and it would look odd
like you keep these abstract IVP and Solver types in there and a single int
I also think we don't need tdir outside of Solver
apart from dense output but I proposed a fix for that
are there any other use cases of tdir?
Mauro
@mauro3
I think it is reasonable that a user may want to query the direction of integration.
Paweł Biernat
@pwl
can you write down a simple example of that? The user is the one who decides which way the integration goes, so he has access to this information
Mauro
@mauro3
Bob gets sent a .jld with an Problem in there. Needs to find out which way it integrates.
I'm just saying that something as fundamental as which way the intergration goes should be queryable without having to dig into the internals of a solver.
Paweł Biernat
@pwl
If we decide on putting it somewhere, I would put it in IVP
you can say that this is still a mathematical property of the initial value problem
some equations might only be well posed forward or backward in time
but then IVP is immutable and it could get in conflict with tout and others.
but that would be ok to me
what I mean by that is that we won't be able to change tdir according to keyword arguments
Paweł Biernat
@pwl
This message was deleted
Paweł Biernat
@pwl

Regarding #17 :

I think having solve(prob) would be putting too much stuff into our interface, which exacerbates the problem with naming that we already have. I was imagining something like this (using the current names for clarity):

1) Have ode for users who "just came here to solve an equation".

2) Have solve as an iterator constructor and the recommended way to deal with differential equations.

We already have a function that "solves" the Problem, and it's called collect. It is clear what it does and it is consistent with how Julia deals with iterators. I don't think there is enough reason to introduce a special name for what you suggest as solve(prob), which is essentially a transpose of collect (tuple of arrays instead of array of tuples coming from collect), we could call it collect_transpose or something but calling it solve has several issues: it would be a waste of a good name (we are short on these), it would suggest that it does more then it actually does (it is collect in a disguise) and it would cause confusion in the long run (isn't iterating the same as solving? which one is the recommended way? What's the difference between ode and solve(prob)?).

And, although I agree that ode is not the best name I wouldn't rename it to solve either. It would suggest that it is more general then it is, in practice it can only solve explicit ODEs. If we want to call it solve we should think about generalizing it to other problems too, which is, in my opinion, beyond the scope of this PR.

I think that, independently of the name we choose for solve, we have to make it very clear in the docs that solve returns an iterable object and it can be, and should be, dealt with just as any other iterable object in Julia. This is way more important then a name we pick (let's face it, there is no perfect choice here). And if you don't know how iterators work, feel free to use ode, which is still powerful enough for most use cases. Later on we could add functions like dae or implicit to solve other types of equations but that's a different story.

I wouldn't mind renaming solve to problem that much (my only gripe with it is the problem(initial value problem) thing), but I would mind having solve(problem(ivp)) replace collect and similar utility functions. That would be bad because people would no longer think in terms of iterators and that would be just one too many different ways to do the same thing, which is bad because we would have to maintain all of them at the same time.

Mauro
@mauro3
I thought I hit enter on this a few hours ago. This was meant for your previous post: "Hmm, it's a tricky one... By your argument just now, I think it [tstop/tdir] should go into Problem. But I also agree that it is a bit awkward there too. Let's sleep over it some more..."
Mauro
@mauro3

I think it is fine and should be encouraged to not use the iterators if not needed. The mostly empty loops in our examples are testimony that usually iterators are not needed. Calling the function which solves an ode ode is fine, but solve would be good too.

The problem with collect is that its return is essentially unusable: a vector of tuples. collect_transpose's return is usable, but still suffers from y being a vector and the old hcat(y...) problem. Also collect_transpose is really obscure.
So, having something like ode, returning something useful, namely a vector for t and an array for y if applicable, will be important to make ODE.jl easily usable.

Paweł Biernat
@pwl
wait, so what kind of return type do you expect from your solve(prob)?
in particular, what do you mean by y being a vector? In contrast to what?
Mauro
@mauro3
y should be a matrix if y0 is a vector. See JuliaLang/ODE.jl#80
Paweł Biernat
@pwl
1) our empty loops are empty because they are simple examples. They are not actual use cases. On the inside of the loop you could have a plotting function to draw the solution at the given time, some write functions to save the results into a file etc. These are common use cases in my programs. It's kind of pointless to argue what the typical use case is, instead we should think on how to provide the most flexible interface.
2) regardless of the output format, this collect_transpose, or whatever you call it, is essentially an extremely simple wrapper around an iterator. It's still far from something that a name solve would suggest.
3) maybe we shouldn't fixate ourselves on a particular output type being the most useful, and instead provide an option to collect_transpose as to how to accumulate the output
4) I don't claim that collect_transpose is the right name for the function we need but it's closer to what it actually does then solve.
What about results(solve(...))?
Paweł Biernat
@pwl
Another idea: have iterate(ode;opts...) return an iterable object and have solve(ode;opts...) return a solution. This way you don't have to pass the return value of iterate to solve, and the names look pretty clear to me.
As for the output type I still would rather have the big matrix output as an option instead of the default (you could pass format options to solve).