Overhaul DiscreteDS and BigDiscreteDS and make them one?

[Datseris]

Can we do something similar as with #20 and also type-parameterize DiscreteDS and BigDiscreteDS and make the method parameter be in-place/out of place, and only have a single struct?

[tkf]

I reply to the comment about dummystate https://github.com/JuliaDynamics/DynamicalSystemsBase.jl/issues/17#issuecomment-362818803 here since it is more specific place, I think. (Sorry it's a bit late)

I actually prefer how DifferentialEquations.jl treats the similar situation. It has clear distinction of problem, integrator, and solution types. The problem types are immutable and all computations involving mutations go to the integrator. I think dummystate can go into the integrator-like type for the discrete DS as well.

Actually, if https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/issues/251 is solved, you can use the integrator from OrdinaryDiffEq.jl to solve the discrete DS. The great thing is that plotting facility via Plots.jl is automatically usable.

Reading perform_step!(integrator,cache::FunctionMapConstantCache,repeat_step=false), it seems there is a bit of extra ifs for each iteration. So maybe it is not as fast as your DS types at the moment. But maybe such code go into OrdinaryDiffEq.jl? OrdinaryDiffEq.jl already has a very generic parameterized DiscreteProblem type so I don't think you need to replicate that.

[tkf]

Oops, you started to use DiscreteProblem in #23 already :)

[tkf]

Actually, I find that you define DiscreteProblem in discrete.jl in #23. I think it's bad since the name crashes with the one in DifferentialEquations.jl.

[Datseris]

Yes. DiffEq definition is flawed from my perspective.

It makes no sense to follow the continuous system syntax for discrete systems. Specifically it makes no sense to "use a solver" since discrete systems are solved by their definition.

It also makes no sense to have time in the equations of motion. I have never seen a discrete system like.

Finally , having tspan part of the problem is a huge issue that doesnt fit discrete systems (or continuous but that is another story).

I can rename DiscreteProblem to avoid conflictions. I will name it DiscreteLaw .

On Feb 8, 2018 5:08 AM, "Takafumi Arakaki" notifications@github.com wrote:

Actually, I find that you define DiscreteProblem in discrete.jl in #23 https://github.com/JuliaDynamics/DynamicalSystemsBase.jl/pull/23. I think it's bad since the name crashes with the one in DifferentialEquations.jl.

— You are receiving this because you modified the open/close state. Reply to this email directly, view it on GitHub https://github.com/JuliaDynamics/DynamicalSystemsBase.jl/issues/21#issuecomment-363996568, or mute the thread https://github.com/notifications/unsubscribe-auth/ASwgYaaLMTXdq4Za22SA1GfcOYYamMADks5tSnMlgaJpZM4R4nwZ .

[tkf]

I disagree. DiffEq is much better. For example, ParallelEvolver you defined is kind of solver/integrator. See also my comment on serialization overhead in https://github.com/JuliaDynamics/DynamicalSystemsBase.jl/issues/20#issuecomment-363998194

[ChrisRackauckas]

A discrete problem doesn't necessarily have integer time. The quintessential example of this are continuous-time Markov chain Gillespie processes, like are shown here:

http://docs.juliadiffeq.org/latest/tutorials/discrete_stochastic_example.html

Those are probably the most common use of the DiscreteProblem actually, with a trivial identity map.

And there's one big reason to make use of the same syntax for DiscreteProblems: then they are compatible with all of the callbacks, so any extensions you made work. All the same controls like saveat etc. just work. And with enough work it should be close to zero overhead, I'm not sure we're quite there with the discrete case but it's worth checking. If you don't want that event handling interface, you can just write an algorithm that's a quick loop with nothing else (I plan to make a repo of these just for benchmarking).

And yes @tkf's point about the immutables is very relevant to parallelization. The problem type is just a small bundle that describes what you want to solve. The algorithm is a small (pure) type that describes how you want to solve it. These are designed to be cheap enough to create as necessary and pass around. The integrator is much heavier.

[Datseris]

Ok.

But is it true that to evolve a different state you have to create a new problem ? What about evolving for different time ?

On Feb 8, 2018 6:20 AM, "Christopher Rackauckas" notifications@github.com wrote:

A discrete problem doesn't necessarily have integer time. The quintessential example of this are continuous-time Markov chain Gillespie processes, like are shown here:

http://docs.juliadiffeq.org/latest/tutorials/discrete_ stochastic_example.html

Those are probably the most common use of the DiscreteProblem actually, with a trivial identity map.

And there's one big reason to make use of the same syntax for DiscreteProblems: then they are compatible with all of the callbacks, so any extensions you made work. All the same controls like saveat etc. just work. And with enough work it should be close to zero overhead, I'm not sure we're quite there with the discrete case but it's worth checking. If you don't want that event handling interface, you can just write an algorithm that's a quick loop with nothing else (I plan to make a repo of these just for benchmarking).

And yes @tkf https://github.com/tkf's point about the immutables is very relevant to parallelization. The problem type is just a small bundle that describes what you want to solve. The algorithm is a small (pure) type that describes how you want to solve it. These are designed to be cheap enough to create as necessary and pass around. The integrator is much heavier.

— You are receiving this because you modified the open/close state. Reply to this email directly, view it on GitHub https://github.com/JuliaDynamics/DynamicalSystemsBase.jl/issues/21#issuecomment-364005720, or mute the thread https://github.com/notifications/unsubscribe-auth/ASwgYVJ4oM4ZzV-Tb4vjxuAETaDwTaWxks5tSoQhgaJpZM4R4nwZ .

[tkf]

Creating a new problem is cheap so why not. I think it can also be done when creating the integrator.

Related to this, I suggest to re-define the semantics of tspan: https://github.com/JuliaDynamics/DynamicalSystemsBase.jl/issues/20#issuecomment-364004097

[Datseris]

Maybe you haven't noticed , but creating an integrator takes an absurdly long time , at least for continuous systems. Many thousand times more time than actually solving the problem.

Creating a new problem means you have to create a new integrator as well, but this. Must be avoided. From my perspective the answer is to use the integrator directly , which is mutable.

On Feb 8, 2018 6:34 AM, "Takafumi Arakaki" notifications@github.com wrote:

Creating a new problem is cheap so why not. I think it can also be done when creating the integrator.

Related to this, I suggest to re-define the semantics of tspan:

20 (comment)

https://github.com/JuliaDynamics/DynamicalSystemsBase.jl/issues/20#issuecomment-364004097

— You are receiving this because you modified the open/close state. Reply to this email directly, view it on GitHub https://github.com/JuliaDynamics/DynamicalSystemsBase.jl/issues/21#issuecomment-364007647, or mute the thread https://github.com/notifications/unsubscribe-auth/ASwgYfJ4hG2uA1gsczdc59CSaKhVV31Fks5tSodZgaJpZM4R4nwZ .

[tkf]

I mean creating a new ODEProblem and DiscreteProblem is cheap. There is a set_u0(prob::ODEProblem, u0) for making a new problem with new inital condition as well.

I'm not exactly sure why you mention the overhead of creating integrator, since you need to do it at least once anyway. But my guess: you want to provide evolve! for the DS types. That is the bad decision we are talking about here. You should construct a mutable integrator-like object first then do the evolution on that object. What I talk about reinit! here https://github.com/JuliaDynamics/DynamicalSystemsBase.jl/issues/20#issuecomment-364004097 is also related.

[ChrisRackauckas]

Many thousand times more time than actually solving the problem.

Do you have an example of that? The cost isn't the integrator though, it's usually the integrator.cache. But for most small systems you should be using static arrays, in which case the cache is trivial since there's no cache vectors. However, I've never seen a case where it's thousands of times slower than solving the problem, so I'm curious what that looks like. 10x sounds reasonable for a quick enough problem with small arrays, 100x is not something I've seen, 1000x sounds absurd and needs proof.

[Datseris]

???

Just time the amount of time it takes to create an integrator and then time the amount of time it takes to 'solve!' It.

There is no way you haven't seen that right ? The latter takes microseconds or even less. The first can take even up to seconds.

On Feb 8, 2018 7:32 AM, "Christopher Rackauckas" notifications@github.com wrote:

Many thousand times more time than actually solving the problem.

Do you have an example of that? The cost isn't the integrator though, it's usually the integrator.cache. But for most small systems you should be using static arrays, in which case the cache is trivial since there's no cache vectors. However, I've never seen a case where it's thousands of times slower than solving the problem, so I'm curious what that looks like.

— You are receiving this because you modified the open/close state. Reply to this email directly, view it on GitHub https://github.com/JuliaDynamics/DynamicalSystemsBase.jl/issues/21#issuecomment-364016108, or mute the thread https://github.com/notifications/unsubscribe-auth/ASwgYUpbQjsr5uAptSKJRhLkczvlnOKgks5tSpUagaJpZM4R4nwZ .

[ChrisRackauckas]

I have never seen an integrator take seconds to build. Please show me an example.