Announcing DynamicalSystems.jl & Inviting for JuliaDynamics

package
proposal
announcement

#1

Hi all. I have been working on a package on my free time that I have mentioned in a previous post in this category. Now it is at a decent state and has been accepted as a Julia package. It’s called DynamicalSystems.jl! With this post I would like to do the following:

  1. Introduce the package
  2. Hear constructive criticism about what could be done better regarding the internals
    (mainly, should all data be represented as a vector of static arrays?)
  3. Invite relevant groups to create JuliaDynamics organization with me!

Point 1

The main goal is for it to be a useful companion for students and scientists treading on the field of Chaos, nonlinear dynamics and dynamical systems in general. I am really hopping that it will expand to become something big, assuming enough people find it useful and in turn would like to contribute. Therefore, if you feel like it, don’t be shy! :smiley:

So far the package is very basic however it does offer a lot of fundamental stuff:

  • Interface for defining dynamical systems and smoothly interacting with DifferentialEquations (for continuous systems)
  • Computation of lyapunov spectra
  • Computation of entropies in a highly efficient manner
  • Automated computation of attractor dimension (by finding the linear region automatically)
  • Delay coordinates reconstruction
  • Very detailed documentation strings and in general documentations (I actually cite papers where the method is taken from because one of the goals for the package is to be as transparent as possible)

**If you want to test the package right now, please use checkout: Pkg.add("DynamicalSystems"); Pkg.checkout("DynamicalSystems") because there have been some major changes in automated attractor dimension estimation.

For example, to get the lyapunov spectrum of a map, and then the attractor dimension you would simply do:

using DynamicalSystems, StaticArrays
@inline @inbounds eom_henon(x) = SVector{2}(1.0 - a*x[1]^2 + x[2], b*x[1])
@inline @inbounds jacob_henon(x) = @SMatrix [-2*a*x[1] 1.0; b 0.0]
he = DiscreteDS(zeros(2), eom_henon, jacob_henon)

this defines a system, and if the Jacobian is missing, it is calculated automatically using ForwardDiff. Now,

ly = lyapunovs(ds, 100000)
# result: [0.420, -1.62]

gives the lyapunov spectrum by applying QR method 100000 times. Also, you can see that it is not very slow:

using BenchmarkTools
@btime lyapunovs(he, 100000);
# result: 298.302 ms (2100006 allocations: 119.02 MiB)

but can probably be even better. Moving on, it is also very simple to calculate generalized entropies and thus dimensions.

ts = timeseries(he, 100000) # a 100000 x 2 matrix that contains the data
α = 1; #generalized entropy exponent
ε = 0.001 # boxsize
g = genentropy(α, ε, ts)  #generalized entropy
D = generalized_dim(α, ts) #generalized dimension (aka Renyi dimension)
# Result: D = 1.2291586862831074
@btime generalized_dim(1, $ts);
# Result: 1.701 s (25138684 allocations: 522.62 MiB)

If you ever use the function generalized_dim() please give me feedback on how well it does, because it performs a lot of automated steps (see the documentation string). In essence it calculates the generalized entropy for many different boxsizes.

All these functions that I have used here are generic and will work for your systems as well, not only the Henon map!

Point 2

Data inside the package is represented as an N x D matrix with N the amount of data points and D the system dimension. There have been comments that the better off approach would be to completely change this and make all data be represented as a vector of StaticArrays. This would not have much of an impact on Discrete systems, since they already operate with static arrays, but it will be a significant change for Continuous systems that operate using in-place functions.

The “best” interface for data representation could be to use RecursiveArrayTools and have the main data structure be a recursive array of static arrays. (notice that for my package almost all data is in the form of timeseries)

@ChrisRackauckas and @kristoffer.carlsson had some comments on that. I would seriously appreciate it if you had a chance to only give a glance at the source code (see e.g. the function timeseries in the files discrete.jl and continuous.jl) and simply confirm if you feel this is best.

Besides this, if anybody that uses the package has ideas for improvement please go ahead and say them, or even better contribute them directly.

Point 3

I have made 2 packages so far; this one and DynamicalBilliards.jl. Both are about some kind of dynamical systems. I am sure there are other packages out there that fit this theme. Interest was already shown from my previous post in this topic. Some of the interested persons were @tkoolen with his package about rigid body dynamics.

However, there isn’t yet any Julia organization that unites these packages. I want to change that, and make this organization. However, I would really, really prefer to not start it by myself with only my 2 packages. So, if there is anybody that wants to participate in this, please let me know asap because I will probably do this organization thingy before the end of september.


#2

Looks great, thanks for the introduction to the package. Also, just want to say that looking at animations from DynamicalBilliards is really mesmerizing…


#3

haha! Thanks so much!

Edit: I always knew making a billiard that is like the Julia logo would be super-catchy!


#4

That looks great.

I see that DynamicalBilliards.jl can be used to visualize the trajectories of the system using pyplot. Do you intend to provide additional plotting abilities, e.g. through the plots package?
Does DynamicalSystems.jl also allow to visualize the trajectories?


#5

I don’t use Plots.jl since I have stuck with PyPlot and there is no reason for me to change currently. However, because Plots.jl seems to be the “standard” approach in Julia, if somebody was to do a Pull Request that changes all plotting to use Plots.jl I would immediately accept it.

What additional plotting abilities are you referring to? Can you give me an example?

Plotting in DynamicalSystems.jl is trivial:

ds = some_system
ts = timeseries(ds, 10000.0)
plot(ts[:, 1], ts[:, 2])

Because the output is not special in any way, you simply use the functions provided by your preferred package directly.
(also you can use the interface of DifferentialEquations to plot, by taking advantage of ODEProblem(ds::DynamicalSystem, t_final).


#6

I think I am just confused by the documentation of DynamicalBilliards. I thought that there were special functions to plot the billiard. I am looking specifically at https://datseris.github.io/DynamicalBilliards.jl/stable/basic/basic_usage/ the very last example that shows how to reproduce the logo billiard (or the ones in https://datseris.github.io/DynamicalBilliards.jl/stable/tutorials/visualizing/). Three functions are advertised
; plot_obstacle, plot_billiard and plot_particle.

But I get the following when I try to reproduce the example.

  julia> plot_billiard(bt)
  ERROR: UndefVarError: plot_billiard not defined

Other examples make reference to a library DynamicalBilliardPlotting but it seems to be deprecated. I have installed DynamicalBilliards 1.3.0 which I think is the stable version.


#7

If you find an issue with DynamicalBilliards.jl (which is not unlikely at all), please report it in the appropriate github issues page. The topic here is about DynamicalSystems.jl, which is unrelated to the billiards package.

To quickly answer, in the installation page you can see that you need to call DynamicalBilliards.enableplotting() to enable plotting. This allows me to not have PyPlot in the REQUIRE file. The same is said in the Vizualizing tutorial.

If you found somewhere reference to DynamicalBilliardsPlotting, please report it as an issue, because it is clearly a mistake (I could not find that somewhere in the docs).


#8

There is no need to drop PyPlot support to enable plotting with Plots.jl :slight_smile:


#9

Well, as long as https://github.com/JuliaPlots/Plots.jl/issues/778 is still open, I cannot do something :smiley:
It is pretty straight-forward to switch to plots, and I do agree that it is the better option.


#10

Well, I trust you know the best way for a talented programmer to close an issue :slight_smile:


#11

I feel really bad for bringing up DybamicBilliards here which is not at all the point of this post. Can we please go back to talking about DynamicalSystems.jl


#12

Good and exciting stuff. In a lot of places in the DifferentialEquations.jl manual we link over to modeling tools. While currently those are all our own, I don’t see why we shouldn’t show off this work. Feel free to PR to add a link in the appropriate spot so that way users who may want this kind of functionality “in the ODE solver” know where to go.