WIP: Overhaul DynamicalSystem #24

Datseris · 2018-02-12T20:48:46Z

Closes #20

This (currently WIP) is the overhauling of all DynamicalSystem types. I have managed to put everything under the type DynamicalSystem and use aliasing and type-parameters for dispatch.

At the moment everything is looking extremely beatiful. I am literally drooling over how beatiful julia is and how extremely transparent code you can write even for such a complicated situation.

The accuracy of methods like e.g. lyapunov increased. Also, a lyapunov computation with NOT autodifferentiated jacobian for a CONTINUOUS system takes 0.07 seconds. I think this is god damn impressive!!! (I am using out-of-place method of course)

To-do:

I only have to test out of place discrete. The others are tested already.

Still todo: Auxilary functions, performance boosts, time evolution functions (evolve / trajectory)

Datseris · 2018-02-12T20:53:29Z

Alright everybody! This is iiiiiiit. I am so liiiiiiiiiiiiit!

So the main file is at the moment dynamicalsystem.jl, a new one. All definitions are there!
The other file to look at is the /test/dynsys_tangent.jl which has at the moment everything. There I have tests that everything works. I have also a huge amounts of timings.

NOTICE: dynamicalsystem.jl is not in a module now. You just run it. Then you run /test/dynsys_tangent.jl and see the results.

The design is pretty straight-forward. DynamicalSystem has the first type parameter the is-in-place. Because most often you want to dispatch on that. CDS and DDS are aliases of DynamicalSystem for continuous and discrete systems.

tangent_integrator makes the tangent dynamics integrator!

@ChrisRackauckas I would truly, truly, truly appreciate it if you took a look at it as well and gave me some hints. I am really bad at performance optimization and maybe I am doing something terribly wrong at some point. The timings do not look very bad, but I also have no knowledge to judge if they are actually bad! (Note: everything is evolved/solved with Vern9() and abstol = reltol = 1e-9)

@RainerEngelken if you want you can also join the discussion! Long story short, we are reworking the fundamental way a dynamical system is defined.

Both "auxilary"

Datseris · 2018-02-12T21:32:54Z

Hm I think I am doing something super-wrong for discrete systems.

lyapunov_oop(hopjac, 2)
@time lyapunov_oop(hopjac, 2) # -> 0.014448 seconds
# There seems to be massive slowdown on lyapunov of discrete. I expected 2
# orders of magnitude speed up. Plus old chaos tools found the exponents in
# 1 milisecond if I recall (for discete). Maybe there is something
# wrong here.

ChrisRackauckas · 2018-02-13T02:43:20Z

Nothing jumps out to me.

tkf

I think this is a significant improvement from the last version. But it can be even better. My overall impression is that the code is still not very DRY. I'm not sure if this is the final version so I'll not point out all the details. But I suggest making it more concise in the final version, since otherwise it's hard to be sure that you are providing well-factorized API.

Also, I think there are too many exports, especially the auxiliary ones. I think you'd want to keep the APIs minimal (let alone the exported names), since it makes maintaining the library harder. Furthermore, it's better to keep exported names minimal in general since that's going to contaminate user's namespace. I think that's the reason why ForwardDiff does not export anything.

Another general comment: I'd suggest being more careful about the import and the semantics of the functions you extend.

tkf · 2018-02-13T03:34:47Z

src/dynamicalsystem.jl

+import OrdinaryDiffEq: isinplace, ODEIntegrator, step!
+import Base: eltype
+
+export dimension, state, DynamicalSystem, DS, integrator, tangent_integrator


I would be hesitant to export aliases like DS, CDS and DDS. Those are easily definable on the user side but they may want to use it differently. Like DS as type parameter of the DynamicalSystem:

strcut Something{DS <: DynamicalSystem} end const ContinuousSomething = Something{ContinuousDynamicalSystem} const DiscreteSomething = Something{DiscreteDynamicalSystem}

Of course, I'm not opposed to internal use.

Agreed. Removed these exports.

tkf · 2018-02-13T04:02:56Z

src/dynamicalsystem.jl

+"""
+set_parameter!(prob, index, value) = (prob.p[index] = value)
+set_parameter!(prob, values) = (prob.p .= values)
+set_parameter!(ds::DS, args...) = set_parameter!(ds.prob, args...)


What's the purpose of set_parameter!?

ChaosTools, obviously.

I mean, ChaosTools just can ds.prob.p[index] = value etc., right? Why would you define a function here and make it an API? I mean, it makes maintaining the package harder...

tkf · 2018-02-13T04:10:04Z

src/dynamicalsystem.jl

+    # Create tangent eom:
+    reducedeom = (du, u) -> ds.prob.f(du, u, ds.prob.p, ds.prob.t)
+    cfg = ForwardDiff.JacobianConfig(reducedeom, rand(3), rand(3))
+    tangenteom = TangentIIP(ds.prob.f, similar(ds.prob.u0, D, D), cfg)


Shouldn't similar(ds.prob.u0, D, D) be similar(ds.prob.u0, D, size(Q0, 2)) to allow solving for Lyapunov exponents smaller than the system's dimension?

tkf · 2018-02-13T04:23:21Z

src/dynamicalsystem.jl

+# or continuous from just the equations of motion!)
+function ContinuousDynamicalSystem(eom, state::AbstractVector, p, j = nothing)
+    IIP = isinplace(eom, 4)
+    u0 = IIP ? Vector(state) : SVector{length(state)}(state...)


I don't think the use of standard/static array has to be enforced. End-users should be able to choose whatever they like. For example, maybe some users want to use GPU arrays. I'd suggest to just use the array type they provided.

If you are desperate for handy static array initialization, how about using SVector when state is a tuple?

If you want to support something other than SVector you are welcome to do a PR both here and at ChaosTools that supports it. I won't do it.

You missed my point. DynamicalSystemsBase.jl has to be generic on array type. I'm not talking about ChaosTools.jl. Everything you rely on (DifferentialEquations.jl and ForwardDiff.jl) accept generic type so it's not difficult at all to define DynamicalSystemsBase.jl in a type-generic way. In fact, you are writing more code to make it less generic...

Why don't you do it then?

Don't you think you are over-using the wording "has to be" ? For somebody that has written 0 lines of code for this package, I don't think you can be in a position to say what "has" to be done.

tkf · 2018-02-13T04:29:55Z

src/dynamicalsystem.jl

+function DiscreteDynamicalSystem(eom, state::AbstractVector, p, j = nothing)
+    IIP = isinplace(eom, 4)
+    u0 = IIP ? Vector(state) : SVector{length(state)}(state...)
+    prob = DiscreteProblem(eom, u0, DDS_TSPAN, p)


DiscreteDynamicalSystem code is mostly identical to ContinuousDynamicalSystem... I'd suggest to define them by

ContinuousDynamicalSystem(f, u0, p, j = nothing) = DynamicalSystem(ODEProblem(f, u0, CDS_TSPAN, p), j) DiscreteDynamicalSystem(f, u0, p, j = nothing) = DynamicalSystem(DiscreteProblem(f, u0, DDS_TSPAN, p), j)

and do all the routine work in the DynamicalSystem constructor.

Already overhauled.

tkf · 2018-02-13T04:42:43Z

src/dynamicalsystem.jl

+    solver, newkw = extract_solver(diff_eq_kwargs)
+
+    prob = ODEProblem(ds.prob.f, U0, CDS_TSPAN, ds.prob.p)
+    integ = init(prob, solver; newkw..., save_everystep = false)


Put save_everystep = false in DEFAULT_DIFFEQ_KWARGS and let people override the setting. Same for DDS. Some end-users may want to draw attractors with it.

Use trajectory for that.

That exclude users from using DynamicalSystemsBase.jl ecosystem. For example, DynamicalSystemsBase.jl has a rich plotting pipeline. Besides, it's not at all difficult to support it (just one line).

No it doesnt. Use solve(ds.prob) with whatever arguments you want. This is less than 1 line, it is just 2 extra characters and you get all the capabilities you ask for.

tkf · 2018-02-13T04:44:06Z

src/dynamicalsystem.jl

+
+
+"""
+    integrator(DS, u0 = ds.prob.u0)


I suggest integrator(ds, tspan = infinite_tspan(ds), u0 = ds.prob.u0).

I also suggest integrator(ds, tspan = infinite_tspan(ds), u0 = ds.prob.u0; kwargs...) and just pass kwargs to init for simplicity but that's not a strong opinion.

tkf · 2018-02-13T05:01:49Z

src/dynamicalsystem.jl

+systemtype(::ODEProblem) = "continuous"
+systemtype(::DiscreteProblem) = "discrete"
+
+function extract_solver(diff_eq_kwargs)


Isn't it equivalent to this?

function extract_solver(diff_eq_kwargs) newkw = deepcopy(diff_eq_kwargs) solver = pop!(newkw, DEFAULT_SOLVER) return solver, newkw end

I mean, assuming this is a small dictionary so that deepcopy is super cheap.

tkf · 2018-02-13T05:16:32Z

src/dynamicalsystem.jl

+#######################################################################################
+
+dimension(prob::DEProblem) = length(prob.u0)
+eltype(prob::DEProblem) = eltype(prob.u0)


No... No... DEProblem is not the type defined in this library so you shouldn't change any of its API, especially for the Base-related APIs. (Monkey-patch is required in some cases but this is totally not the case.) Furthermore, eltype mostly is for container or iterator. DEProblem is neither.

The second reason applies to DS. Don't define eltype(::DS).

What? Why.

I am using dimension, not DiffEq. What you say doesn't make sense.

Ah, I now see why you are asking about dimension in discourse. But please read carefully. I'm talking about Base.eltype. Not dimension.

Yes I have already removed that, it was only for internal test.

tkf · 2018-02-13T05:22:49Z

src/dynamicalsystem.jl

@@ -0,0 +1,349 @@
+using OrdinaryDiffEq, ForwardDiff, StaticArrays
+import OrdinaryDiffEq: isinplace, ODEIntegrator, step!


ODEIntegrator is not used

taken care of.

tkf · 2018-02-13T05:41:00Z

As I mentioned in #20 (comment), I think the name/concept "tangent integrator" has to be reconsidered and so has to be in the TODOs. I think the most meaningful choice at the moment is "augmented integrator" but I still don't like it.

Datseris · 2018-02-13T22:25:55Z

Currently at the master branch of DynamicalSystemsBase I have my own implementation of DiscreteProblem, which I have written about 2 weeks ago. I was wrongly benchmarking versus DiffEq, because I was benchmarking versus solve!, which after its first run, does nothing, thus returning ridicoulously low times.

The results are in this gist: https://gist.github.com/Datseris/5123a27f5659d14344551a296b7e57d3

In general I always gets 1 order of magnitude speed, most of the time even more. I think we can say that because discrete problems are inherently fast, this speed can be sacrificed for the callback functionality.

@ChrisRackauckas ?

ChrisRackauckas · 2018-02-13T23:09:35Z

Yeah, pure discrete problems with cheap maps are tough because we need to make sure we skip as much as possible in the normal usage, but still allow the extra functionality. I'm not sure what the performance tradeoff is, but a straight loop would be faster, the question is how much faster. But what we should do is profile it and see if we can make this use case faster. My gut feeling is that we can make it significantly faster by making the internals of the integrator skip some parts via type checks, and by that it should be able to almost match a straight loop. I don't believe we are there yet though (but Julia v0.7 will be extremely helpful because of how it propagates constants).

On a sidenote, I've wanted to implement SimpleDiffEq.jl that has a few straight loop (no callbacks, etc) implementations of the most useful algorithms to be able to get a sense of what the maximum performance is, and use that as a benchmark target.

Datseris · 2018-02-13T23:12:11Z

On a sidenote, I've wanted to implement SimpleDiffEq.jl that has a few straight loop (no callbacks, etc) implementations of the most useful algorithms to be able to get a sense of what the maximum performance is, and use that as a benchmark target.

This is a great idea, not only for performance comparison but also for normal users that have use-cases that do not care about callbacks. I am willing to put my minimal DiscreteProblem implementation there, should you accept it.

Of course it has to be somehow refined/adapted to DiffEq styles, but I am guessing you would want to shape it yourself afterwards, since I have precisely 0 knowledge of DiffEq internals.

ChrisRackauckas · 2018-02-14T06:49:17Z

I created the library. Let's do it.

Todo left is to make tangent integrator and to use generic tangent function

and many other minor improvements

Datseris · 2018-02-15T11:15:01Z

https://discourse.julialang.org/t/even-more-clarification-on-type-piracy/9079

tkf · 2018-02-16T08:28:48Z

I commented on discourse but it looks like ~~it takes sometime for it to appear, since it's the first time I comment.~~

Edited: It's already approved.

tkf · 2018-02-16T08:45:19Z

Commented on:

These comment still apply:

I also finished tangent_integrator for discrete

tkf

I think I've commented too much. I'll hit the approve button for now.

Datseris · 2018-02-17T11:11:22Z

This is so god damn sad.

When I compute lyapunovs of discrete systems I get 1000x slowdown from the current stable version.

The problem is that the code has now become so complicated that I have no clue where the regressions are; probably in multiple places

I am separating discrete and continuous with this commit. The `DynamicalSystem` is now an abstract type.

Datseris · 2018-02-18T11:28:54Z

towel map, current version

@btime ($lyapunovs)($ds, 1000)

  234.013 μs (176 allocations: 11.28 KiB)
Float64[3]
0.429…
0.365…
-3.28…

Datseris added 4 commits February 11, 2018 17:13

[skip ci] init commit

53ed26f

all tangent integrators written

4e822b0

I only have to test out of place discrete. The others are tested already.

finalize fundamental interface

d82af2c

Still todo: Auxilary functions, performance boosts, time evolution functions (evolve / trajectory)

add lyapunov timing

1eb4d2d

Datseris self-assigned this Feb 12, 2018

Datseris requested a review from tkf February 12, 2018 20:48

Datseris added help wanted documentation testing This issue/PR needs tests precision_improvement breaking performance labels Feb 12, 2018

add set_parameter! and safe_state_space

7741f0f

Both "auxilary"

Datseris requested a review from ChrisRackauckas February 12, 2018 21:29

example of slow discrete lyapunov

cb0db84

tkf reviewed Feb 13, 2018

View reviewed changes

Datseris mentioned this pull request Feb 14, 2018

Minimal discrete implementation SciML/SimpleDiffEq.jl#6

Closed

Datseris added 3 commits February 14, 2018 20:49

successfully split files

45bd2de

Todo left is to make tangent integrator and to use generic tangent function

[skip ci] update famous systems

936f938

and many other minor improvements

[skip ci] add trajectory

a66244f

Datseris added 3 commits February 16, 2018 16:37

[skip ci] add tangent_integrator for continuous

a362686

[skip ci] add parallel_integrator

2ac2749

I also finished tangent_integrator for discrete

[skip ci] minor cleanup

bdb7800

tkf approved these changes Feb 17, 2018

View reviewed changes

Datseris added 2 commits February 17, 2018 09:45

tests pass

7dc41d8

fix parallel integrator for CDS{iip}

c651f80

Datseris added 2 commits February 18, 2018 11:49

[skip ci] superb inference and timing improvements on Discrete

a89f356

I am separating discrete and continuous with this commit. The `DynamicalSystem` is now an abstract type.

[skip ci] Massive tangent boost for OOP

7ffbbc2

Datseris added 5 commits February 18, 2018 15:48

all things updated

f561c4f

reinit tweaks

02cb4a1

documentation improvements

bd2ec78

make reinit! same as DiffEq

d82d38f

discrete integrator extention

b982b05

Datseris merged commit 6a015f7 into master Feb 19, 2018

Datseris deleted the dynamicalsystem branch February 19, 2018 19:40

		@@ -0,0 +1,349 @@
		using OrdinaryDiffEq, ForwardDiff, StaticArrays
		import OrdinaryDiffEq: isinplace, ODEIntegrator, step!

WIP: Overhaul DynamicalSystem #24

WIP: Overhaul DynamicalSystem #24

Conversation

Datseris commented Feb 12, 2018 • edited Loading

Datseris commented Feb 12, 2018 • edited Loading

Datseris commented Feb 12, 2018 • edited Loading

ChrisRackauckas commented Feb 13, 2018

tkf left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Datseris Feb 16, 2018 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

tkf commented Feb 13, 2018

Datseris commented Feb 13, 2018

ChrisRackauckas commented Feb 13, 2018

Datseris commented Feb 13, 2018

ChrisRackauckas commented Feb 14, 2018

Datseris commented Feb 15, 2018

tkf commented Feb 16, 2018 • edited Loading

tkf commented Feb 16, 2018

tkf left a comment

Choose a reason for hiding this comment

Datseris commented Feb 17, 2018

Datseris commented Feb 18, 2018

Datseris commented Feb 12, 2018 •

edited

Loading

Datseris commented Feb 12, 2018 •

edited

Loading

Datseris commented Feb 12, 2018 •

edited

Loading

Datseris Feb 16, 2018 •

edited

Loading

tkf commented Feb 16, 2018 •

edited

Loading