The tag name "v0.6" is not of the appropriate SemVer form (vX.Y.Z).
cc: @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?

To allow identical interface between both system types, continuouys and disscrete

Goal

We request a minimal and self-contained integrator for systems of ODEs. This could be added here, but it makes much more sense for it to be in SimpleDiffEq repo. This issue is here for the bounty.

This minimal integrator, because it will not have most of the features of a standard DiffEq integrator, will be faster to compile and maybe faster to run as well. Who knows... But, besides any performance benefits, it will have educational value. A principal of DynamicalSystems.jl is that "you know whats going on", and a simple integrator with only a couple of .jl files will make this easier.

Our main goal is twofold:

1. Reduce the "time until first result", which is mostly compilation.
2. Make something that is easy to read and easy to understand.

Detailed Description

We want a continuous system ODE solver that has the following properties:

1. Extends init from DiffEqBase to create/initialize the integrator object given an ODEProblem.
2. Extends step! from DiffEqBase to "step" the integrator for one step. Adaptive stepping and error checking must be included.
3. It must not support callbacks, event handling, etc.
4. It must not support storing / creating a solution object. We will be exclusively using step!. Only the previous state and current state need to be stored.
5. It must support both in-place and out-of-place type of equations of motion (see diff eq for more).
6. It must support the following three types of states: Vector of numbers, Vector of Vectors (both Vector and SVector) and Matrix. Notice that for state-type VectorOfVectors only in-place version is necessary. With this I mean that it doesn't matter if the state is Vector{Vector} or Vector{SVector}, the outer vector is modified always in-place. For the other two cases both in-place (Vector/Matrix {Number}) and out-of-place (SVector/SMatrix {Number} ) are necessary.
7. It must support interpolation in between the previous and current state. The code from OrdinaryDiffEq does this already.
8. It must support "setting"/"changing" the current state. This allows not only reinit! but also allows one to change the state during integration which is crucial for JuliaDynamics.
9. You are also very welcomed to cut high-level features like allowing a distribution as initial state, etc.

You can choose whichever you want DOPRI5/Verner9 (I vote for Verner9 but who am I).
I am suggesting those two so you can actually copy the code that @ChrisRackauckas has pasted there, and have the least effort for this issue. Please try to copy as much from the actual step! implementation of OrdinaryDiffEq as possible, since it is highly optimized. This is open source in the end, just be sure to copy the license ;)

The Bounty

The 100$ bounty of this issue is posted by user @Datseris , a member of JuliaDynamics. @Datseris does not have merge or admin privileges to repositories of the DiffEq org. Once the integrator is merged you can do a one-liner PR here that changes the default integrator to the one you contributed. Even if the integrator is not accepted by DiffEq but it is accepted in JuliaDynamics, you can still move it here and you will get the bounty.

Did you help close this issue? Go claim the $100 bounty on Bountysource.

This is not that hard, but an extra step of solving this would be making a full conversion when creating a tree from this Dataset, as NEarestNeighbors does not allow views.

Want to back this issue? Post a bounty on it! We accept bounties via Bountysource.

Implement what is discussed here: JuliaDiff/ForwardDiff.jl#296

So that we can have datasets of vectors, matrices etc, instead of only svectors.

Also, define the default Dataset{D,T} = Dataset{D, T, SVector{D, T}.

Hey there,
I'm not sure how much you can do about this, but there is a serious issue with
compilation of the Dataset(::Matrix) function.

julia> a = rand(10,1000);

julia> @time Dataset(a);
 10.324998 seconds (315.09 k allocations: 15.435 MiB, 0.07% gc time)

julia> @btime Dataset(a);
  92.777 μs (23 allocations: 87.20 KiB)

julia> b = rand(10,2000);

julia> @time Dataset(b)
 39.327667 seconds (315.11 k allocations: 15.536 MiB, 0.01% gc time)

julia> @time Dataset(b)
  0.000229 seconds (27 allocations: 173.297 KiB)```

<bountysource-plugin>

---
Want to back this issue? **[Post a bounty on it!](https://www.bountysource.com/issues/61089189-throw-a-warning-when-creating-dataset-of-d-100?utm_campaign=plugin&utm_content=tracker%2F81663394&utm_medium=issues&utm_source=github)** We accept bounties via [Bountysource](https://www.bountysource.com/?utm_campaign=plugin&utm_content=tracker%2F81663394&utm_medium=issues&utm_source=github).
</bountysource-plugin>

This is a good change and an easy one, but it will require some amount of new code.

In principle, the following have to be modified

• Constructors of DynamicalSystem
• Tangent integrator to create MMatrix when the state is MVector
• Parallel integrator to create Vector[MArray]

I hope I am not missing something. Because DynamicalSystem is parameterized on the state type as well (see advanced documentation) expanding the required methods is very straight-forward!

Hey!

I'm trying to fix all my test cases to your D → D+1 change with the number of neighbors.
Why did you want to do this again? I find it super confusing that D is used as number
of neighbors D and as type parameter D=D+1 as well.

Anyway, there seems to be an inconsistency as
MTDelayEmbedding needs D = 4
for τ=[0 0; 15 15; 30 29; 45 45]
Even though there are only three actual neighbors.

Straight forward issue, easy to implement but unfortunately one has to write many lines of code! 🗡

This allows to passing matrix-valued ODEProblem as well.

Allowing to evolve the system for some transient amount of time Ttr and then recording would be very convenient.

the DynamicalSystem and DiscreteDynamicalSystem implementations are pure (only StaticArray dependence).

They can easily be updated! While this is done, be sure that:

• Discrete tests pass
• Stepping any discrete integrator is always non allocating

From @Datseris on August 3, 2017 10:58

Hybrid systems are a special type of dynamical system that switches between being continuous and discrete. An example could be the bouncing ball.

It will probably be a tiny bit tricky to implement them, but it has to be made as general as possible.

One need to have 2 equations of motion, one for continuous and one for discrete. Then a callback type function has to be created that will stop the ODE solving of DifferentialEquations and switch to discrete evolution (until some other kind of callback stops or for a specified amount of steps... implementation will show the best way...).

Copied from original issue: JuliaDynamics/DynamicalSystems.jl#14

Test the columns function using the gissinger model.

The theiler window is now part of neighborhood. However, it is still not perfectly correct. Because, it is still possible to get let's say 5 neighbors of a point, all of which are not temporal neighbors with the point itself, but they can still be temporal neighbors between themselfs.

If e.g. n=10 then [666, 667, 668] is still a valid neighborhood.

We should further make the neighborhood function to also not have neighbors that are temporal neighbors with each other.

This is a difficult decision though because you have to decide which point to keep from within each individual temporal neighborhood.

UPDATED ISSUE DESCRIPTION

Some kind of interface has to be constructed for using Hamiltonian continuous systems and integrating them while conserving energy.

DifferentialEquations.jl has some ways to already integrate while conserving energy:

Either as a callback or as splitting the system to 2 equations and using dynamical odes

dynamical odes:
http://docs.juliadiffeq.org/stable/types/dynamical_types.html
http://docs.juliadiffeq.org/stable/solvers/dynamical_solve.html

callbacks:
http://docs.juliadiffeq.org/stable/features/callback_library.html

What we need is a struct HamiltonianDynamicalSystem (shortcut HDS) that is the same as CDS but also has one extra field hamiltonian which is a function.

Questions

• How will the state be in HDS?
• How will the jacobian be expected? It must be able to take in a state.
• tangent and parallel integrators have to be created. Of course, both of these must also conserve energy!
• can we use the function hamiltonian to conserve energy without using these split ode types? This would allow HDS to have same behavior as CDS and we would write much less code.

Sad Answers

After talking with Chris, at the moment it is not possible to solve this issue. It is not possible to have a symplectic integrator that evolves states and tangent space.

However, it may be worth it to make HamiltonianDynamicalSystem never-the-less and just define

tangent_integrator(ds::HDS, args...;kwrags...) = error("Tangent space can't be evolved for 
Hamiltonian systems. Create a standard ContinuousDynamicalSystem instead."

Want to back this issue? Post a bounty on it! We accept bounties via Bountysource.

Makes 0 sense.

data = Dataset(rand(1000, 2))
write_dataset("data.txt", data, ',')

using BenchmarkTools

@btime read_dataset("data.txt", Dataset{2, Float64}, ',');
@btime Dataset(readdlm("data.txt", ',', Float64));

  3.279 ms (12996 allocations: 532.67 KiB)
  2.719 ms (4971 allocations: 227.22 KiB)

we should switch the implementation to use readdlm instead.

Allen, @halleysfifthinc this seems to be something we have to do, not for dispatch but for memory purposes. When I say memory I mean the actual "remembering the past" memory.

Algorithms that we want to implement actually want to know how many different timeseries are used to reconstruct. With the current type-information this is not possible.

(I am not assigning you to this issue or anything like that, just letting you know of an upcoming change)

The change will simply be
MDReconstruction{S2, D, τ, T} <: AbstractReconstruction{D*S2, τ, T} <: AbstractDataset{D*S2, T}

and it will also include parameterization on types τ::Vector{Int} in accordance with #32 .

And change default tolerances to 1e-9.

The only thing that has to be done is to change function get_solver!

So, dataset -> matrix conversion is superfast. But maybe it is faster to reinterpret to matrix, do svd, and then only transpose the U matrix, and return transposed U and S,, which would be the svd of the original matrix.

Ooookay.

I succumb to the peer pressure. Time to change this...

This was present in older versions and I removed it in humongous lapse of judgement.

Super cool conserevative 3D system (continuous)

Can be found here: https://books.google.de/books?id=xOnFCgAAQBAJ&pg=PA95&source=gbs_toc_r&cad=3#v=onepage&q&f=false

[book "Elegant Chaos" by Sprott]

This would make lyapunovs and gali also conserve energy.

Reading how ContinuousDS(prob::ODEProblem) w/o Jacobian:

# Constructors without Jacobian:
function ContinuousDS(prob::ODEProblem)
    state = prob.u0
    eom! = prob.f

    D = length(state); T = eltype(state)
    du = copy(state)
    J = zeros(T, D, D)

    jeom! = (du, u) -> eom!(0, u, du)
    jcf = ForwardDiff.JacobianConfig(jeom!, du, state)
    ForwardDiff_jacob! = (t, u, J) -> ForwardDiff.jacobian!(
    J, jeom!, du, u, jcf)
    ForwardDiff_jacob!(0, state, J)

    return ContinuousDS(prob, ForwardDiff_jacob!, J)
end

and variational_integrator:

function variational_integrator(ds::ContinuousDS, k::Int, T,
    S::AbstractMatrix; diff_eq_kwargs = DEFAULT_DIFFEQ_KWARGS)

    f! = ds.prob.f
    jac! = ds.jacob!
    J = ds.J
    # ... chop ....
    veom! = (t, u, du) -> begin
        us = view(u, :, 1)
        f!(t, us, view(du, :, 1))
        jac!(t, us, J)
        A_mul_B!(view(du, :, 2:k+1), J, view(u, :, 2:k+1))
    end
    # ... chop ....
    return vintegrator
end

are defined, I started worry that this implementation may be inefficient, since you call eom! (prob.f) twice. Discrete system also has the same problem. I haven't done any benchmark on this, though.

Here's an implementation for not calling eom! twice (as you know 😉 ):
https://gist.github.com/tkf/014787f386e99b7b8b42683cfbd8da01#file-lyapunovexponentswithforwarddiff-jl-L27

The implementation itself is not hard. The main problem is that you need to change the API. That is to say, you need to store the pair eom! and (say) "eom_and_jacob!" rather than eom! and jacob!. Of course, it is easy to make a backward-compatible layer by calling eom! and jacob! from the default "eom_and_jacob!". This also helps user to provide a highly efficient implementation.

(A bit tangential notes: I know that the name eom_and_jacob! is not cool. Not sure what's the best name. tangent_dynamics!? coeom! (co-eom)? teom! (tangent-eom)? veom! (variational-)? peom! (paired-)?)

I can do the fix but I thought I'd ask you first, since this is going to change the API (and there is a chance that you have done benchmarks and know that the performance benefit is negligible).

In a reconstruction interface, allow τ to also be a vector of integers, so that D == length(τ) + 1.

If this happens, each reconstruction component has delay time τ[i] instead of a constant τ::Int.

This reference: K. Judd and A.Mees, Physica D 120, 273–286 (1998) introduces this concept as well as its usefulness in multi-scale dynamics.

We have our own implementation of Discrete solvers, which is as fast as possible.

If you see the definition of TangentOOP, it does:

@inline function (tan::TangentOOP)(u, p, t)
    du = tan.f(u[:, 1], p, t)
    J = tan.jacobian(u[:, 1], p, t)
    dW = J*u[:, tan.ws]
    return hcat(du, dW)
end

this is necessary for the continuous case, because the state has to be one coherent sturcture.

For discrete it is not necessary. We can save the following computations:
accessing the state u[:, 1] since they will be separate fields of the new DiscreteTangentIntegrator. We also completely save the hcat operation.

Trying to test DynamicalSystems on Julia 0.7.beta2, I get this error after Pkg.resolve()

           _

_ _ ()_ | A fresh approach to technical computing
() | () () | Documentation: https://docs.julialang.org
_ _ | | __ _ | Type "?" for help, "]?" for Pkg help.
| | | | | | |/ ` | |
| | || | | | (| | | Version 0.7.0-beta2.0 (2018-07-13 19:54 UTC)
/ |_'|||_'_| | Official http://julialang.org/ release
|__/ | x86_64-w64-mingw32

(v0.7) pkg> add DynamicalSystemsBase
Updating registry at C:\Users\Denis\.julia\registries\General
Updating git-repo https://github.com/JuliaRegistries/General.git
Resolving package versions...
Updating C:\Users\Denis\.julia\environments\v0.7\Project.toml
[6e36e845] + DynamicalSystemsBase v0.10.3
Updating C:\Users\Denis\.julia\environments\v0.7\Manifest.toml
[no changes]

(v0.7) pkg> test DynamicalSystems
Testing DynamicalSystems
WARNING: Base.iteratorsize is deprecated, use IteratorSize instead.
likely near C:\Users\Denis.julia\packages\DynamicalSystemsBase\dHai\src\dataset.jl:14
WARNING: importing deprecated binding Base.Range into DynamicalSystemsBase.
WARNING: Base.Range is deprecated, use AbstractRange instead.
likely near C:\Users\Denis.julia\packages\DynamicalSystemsBase\dHai\src\dataset.jl:25
WARNING: Base.Range is deprecated, use AbstractRange instead.
likely near C:\Users\Denis.julia\packages\DynamicalSystemsBase\dHai\src\dataset.jl:27
WARNING: Base.Range is deprecated, use AbstractRange instead.
likely near C:\Users\Denis.julia\packages\DynamicalSystemsBase\dHai\src\dataset.jl:29
WARNING: Base.Range is deprecated, use AbstractRange instead.
likely near C:\Users\Denis.julia\packages\DynamicalSystemsBase\dHai\src\dataset.jl:31
┌ Warning: Requires now needs a UUID; please see the readme for changes in 0.7.
└ @ DynamicalSystemsBase C:\Users\Denis.julia\packages\DynamicalSystemsBase\dHai\src\dataset.jl:230
WARNING: Base.svd is deprecated: it has been moved to the standard library package LinearAlgebra.
Add using LinearAlgebra to your imports.
likely near C:\Users\Denis.julia\packages\DynamicalSystemsBase\dHai\src\dataset.jl:305
ERROR: LoadError: LoadError: ArgumentError: Package DynamicalSystemsBase does not have Distances in its dependencies:

• If you have DynamicalSystemsBase checked out for development and have
  added Distances as a dependency but haven't updated your primary
  environment's manifest file, try Pkg.resolve().
• Otherwise you may need to report an issue with DynamicalSystemsBase.

Stacktrace:
[1] require(::Module, ::Symbol) at .\loading.jl:821
[2] include at .\boot.jl:317 [inlined]
[3] include_relative(::Module, ::String) at .\loading.jl:1034
[4] include at .\sysimg.jl:29 [inlined]
[5] include(::String) at C:\Users\Denis.julia\packages\DynamicalSystemsBase\dHai\src\DynamicalSystemsBase.jl:7
[6] top-level scope at none:0
[7] include at .\boot.jl:317 [inlined]
[8] include_relative(::Module, ::String) at .\loading.jl:1034
[9] include(::Module, ::String) at .\sysimg.jl:29
[10] top-level scope at none:0
[11] eval at .\boot.jl:319 [inlined]
[12] eval(::Expr) at .\client.jl:394
[13] top-level scope at .\none:3 [inlined]
[14] top-level scope at .<missing>:0
in expression starting at C:\Users\Denis.julia\packages\DynamicalSystemsBase\dHai\src\neighborhoods.jl:2
in expression starting at C:\Users\Denis.julia\packages\DynamicalSystemsBase\dHai\src\DynamicalSystemsBase.jl:12
ERROR: LoadError: Failed to precompile DynamicalSystemsBase to C:\Users\Denis.julia\compiled\v0.7\DynamicalSystemsBase\oAQS.ji.
Stacktrace:
[1] error(::String) at .\error.jl:33
[2] macro expansion at .\logging.jl:298 [inlined]
[3] compilecache(::Base.PkgId) at .\loading.jl:1173
[4] _require(::Base.PkgId) at .\logging.jl:298
[5] require(::Base.PkgId) at .\loading.jl:838
[6] require(::Module, ::Symbol) at .\loading.jl:833
[7] include at .\boot.jl:317 [inlined]
[8] include_relative(::Module, ::String) at .\loading.jl:1034
[9] include(::Module, ::String) at .\sysimg.jl:29
[10] top-level scope at none:0
[11] eval at .\boot.jl:319 [inlined]
[12] eval(::Expr) at .\client.jl:394
[13] top-level scope at .\none:3 [inlined]
[14] top-level scope at .<missing>:0
in expression starting at C:\Users\Denis.julia\packages\DynamicalSystems\Ii9Y\src\DynamicalSystems.jl:14
ERROR: LoadError: Failed to precompile DynamicalSystems to C:\Users\Denis.julia\compiled\v0.7\DynamicalSystems\TY1v.ji.
Stacktrace:
[1] error(::String) at .\error.jl:33
[2] macro expansion at .\logging.jl:298 [inlined]
[3] compilecache(::Base.PkgId) at .\loading.jl:1173
[4] _require(::Base.PkgId) at .\loading.jl:942
[5] require(::Base.PkgId) at .\loading.jl:838
[6] require(::Module, ::Symbol) at .\loading.jl:833
[7] include at .\boot.jl:317 [inlined]
[8] include_relative(::Module, ::String) at .\loading.jl:1034
[9] include(::Module, ::String) at .\sysimg.jl:29
[10] include(::String) at .\client.jl:393
[11] top-level scope at none:0
in expression starting at C:\Users\Denis.julia\packages\DynamicalSystems\Ii9Y\test\runtests.jl:1
ERROR: Package DynamicalSystems errored during testing

(v0.7) pkg> resolve
Resolving package versions...
Updating C:\Users\Denis\.julia\environments\v0.7\Project.toml
[no changes]
Updating C:\Users\Denis\.julia\environments\v0.7\Manifest.toml
[no changes]

(v0.7) pkg> test DynamicalSystems
Testing DynamicalSystems
WARNING: Base.iteratorsize is deprecated, use IteratorSize instead.
likely near C:\Users\Denis.julia\packages\DynamicalSystemsBase\dHai\src\dataset.jl:14
WARNING: importing deprecated binding Base.Range into DynamicalSystemsBase.
WARNING: Base.Range is deprecated, use AbstractRange instead.
likely near C:\Users\Denis.julia\packages\DynamicalSystemsBase\dHai\src\dataset.jl:25
WARNING: Base.Range is deprecated, use AbstractRange instead.
likely near C:\Users\Denis.julia\packages\DynamicalSystemsBase\dHai\src\dataset.jl:27
WARNING: Base.Range is deprecated, use AbstractRange instead.
likely near C:\Users\Denis.julia\packages\DynamicalSystemsBase\dHai\src\dataset.jl:29
WARNING: Base.Range is deprecated, use AbstractRange instead.
likely near C:\Users\Denis.julia\packages\DynamicalSystemsBase\dHai\src\dataset.jl:31
┌ Warning: Requires now needs a UUID; please see the readme for changes in 0.7.
└ @ DynamicalSystemsBase C:\Users\Denis.julia\packages\DynamicalSystemsBase\dHai\src\dataset.jl:230
WARNING: Base.svd is deprecated: it has been moved to the standard library package LinearAlgebra.
Add using LinearAlgebra to your imports.
likely near C:\Users\Denis.julia\packages\DynamicalSystemsBase\dHai\src\dataset.jl:305
ERROR: LoadError: LoadError: ArgumentError: Package DynamicalSystemsBase does not have Distances in its dependencies:

• If you have DynamicalSystemsBase checked out for development and have
  added Distances as a dependency but haven't updated your primary
  environment's manifest file, try Pkg.resolve().
• Otherwise you may need to report an issue with DynamicalSystemsBase.

Stacktrace:
[1] require(::Module, ::Symbol) at .\loading.jl:821
[2] include at .\boot.jl:317 [inlined]
[3] include_relative(::Module, ::String) at .\loading.jl:1034
[4] include at .\sysimg.jl:29 [inlined]
[5] include(::String) at C:\Users\Denis.julia\packages\DynamicalSystemsBase\dHai\src\DynamicalSystemsBase.jl:7
[6] top-level scope at none:0
[7] include at .\boot.jl:317 [inlined]
[8] include_relative(::Module, ::String) at .\loading.jl:1034
[9] include(::Module, ::String) at .\sysimg.jl:29
[10] top-level scope at none:0
[11] eval at .\boot.jl:319 [inlined]
[12] eval(::Expr) at .\client.jl:394
[13] top-level scope at .\none:3 [inlined]
[14] top-level scope at .<missing>:0
in expression starting at C:\Users\Denis.julia\packages\DynamicalSystemsBase\dHai\src\neighborhoods.jl:2
in expression starting at C:\Users\Denis.julia\packages\DynamicalSystemsBase\dHai\src\DynamicalSystemsBase.jl:12
ERROR: LoadError: Failed to precompile DynamicalSystemsBase to C:\Users\Denis.julia\compiled\v0.7\DynamicalSystemsBase\oAQS.ji.
Stacktrace:
[1] error(::String) at .\error.jl:33
[2] macro expansion at .\logging.jl:298 [inlined]
[3] compilecache(::Base.PkgId) at .\loading.jl:1173
[4] _require(::Base.PkgId) at .\logging.jl:298
[5] require(::Base.PkgId) at .\loading.jl:838
[6] require(::Module, ::Symbol) at .\loading.jl:833
[7] include at .\boot.jl:317 [inlined]
[8] include_relative(::Module, ::String) at .\loading.jl:1034
[9] include(::Module, ::String) at .\sysimg.jl:29
[10] top-level scope at none:0
[11] eval at .\boot.jl:319 [inlined]
[12] eval(::Expr) at .\client.jl:394
[13] top-level scope at .\none:3 [inlined]
[14] top-level scope at .<missing>:0
in expression starting at C:\Users\Denis.julia\packages\DynamicalSystems\Ii9Y\src\DynamicalSystems.jl:14
ERROR: LoadError: Failed to precompile DynamicalSystems to C:\Users\Denis.julia\compiled\v0.7\DynamicalSystems\TY1v.ji.
Stacktrace:
[1] error(::String) at .\error.jl:33
[2] macro expansion at .\logging.jl:298 [inlined]
[3] compilecache(::Base.PkgId) at .\loading.jl:1173
[4] _require(::Base.PkgId) at .\loading.jl:942
[5] require(::Base.PkgId) at .\loading.jl:838
[6] require(::Module, ::Symbol) at .\loading.jl:833
[7] include at .\boot.jl:317 [inlined]
[8] include_relative(::Module, ::String) at .\loading.jl:1034
[9] include(::Module, ::String) at .\sysimg.jl:29
[10] include(::String) at .\client.jl:393
[11] top-level scope at none:0
in expression starting at C:\Users\Denis.julia\packages\DynamicalSystems\Ii9Y\test\runtests.jl:1
ERROR: Package DynamicalSystems errored during testing

(v0.7) pkg>

Easy follow to #28

A dataset has the dimension as type-information, so no need for SizedArray here

The methods and equations of motion of our package will change greatly, following the discussion of SciML/DifferentialEquations.jl#235

In short, the API for all equations of motion will become something like f(du, u, t, p) or f(du, u, p, t) (for out of place methods remove du).

This will also remove all kind of Functor creation of systems since the parameters p are always passed implicitly to the equations of motion. The same syntax will be used for the Jacobian. Using ForwardDiff with this syntax will be very easy nontheless.

The systems now will have a field that also has the parameters. The ODEProblem for continuous systems will have this field, and for discrete they will each have a dedicated field params or something like that.

This will break everything by the way.

For prediction of spatio-temporal timeseries a special type of delay embedding is necessary as descibed in [1].
The embedding should work like this:

The call signature could look like this:
Reconstruction(s,D,τ,I,K, c)
where D and τ delay and (delay-) dimension as usual.. I is the number of spatial neighbors and K the
spatial shift.
c is used for boundary conditions as visualized in the picture. If a proper default value can be found this could be made an optional kwarg.

[1] U. Parlitz, Prediction of Spatiotemporal Time Series Bases on Reconstructed Local States (2000)

If we do this, it will also be easy to use callbacks etc with diffeq.

Will also solve JuliaDynamics/DynamicalSystems.jl#14

https://en.wikipedia.org/wiki/Chua%27s_circuit#cite_note-5

Write generated code for the case of accessing a dataset like Dataset[anything, SVector{Int}].

Want to back this issue? Post a bounty on it! We accept bounties via Bountysource.

This makes more sense in our neighborhood implementation, which anyways requires n to be an argument.

In addition, this is a good opportunity to make the neighborhood methods return neighbors that are 1 theiler window apart. Optionally.

This seems to be really important for continuous systems that are sampled very densely.

Maybe its good to define substraction, addition etc. between Datasets. Maybe norm as well.

While we are at it, let's also define a Dataset{D,T}(N::Int) constructor that starts a dummy dataset.

Prof. Ulrich Parlitz knows this.

DiffEq already allows one to pass saveidxs to solve to save only specific variables for the state.

It is also very easy to modify trajectory for discrete systems to use the same keyword argument and have the same behavior!

We will do this when moving to 0.7!

Alright, its about time for fundamentally improving the function estimate_delay.

For the Gissinger system at default parameters it gives delay 500 while very good values are in the range of 30-40 for timestep dt = 0.05.

Notice however that the delay time suggested by the function multiplied with dt gives the perfect delay of around 30.

Hmmmmmm

EDIT: THIS IS THE NEW AND CURRENT PLAN, IN PROGRESS BY ME IN A BRANCH.

Example at the end!

• Rename ContinuousDS to ContinuousDynamicalSystem with alias CDS.

• CDS will have 2 fields: prob::ODEProblem and jacobian (which is a function).

  • At a later point we might consider having tangent dynamics function as part of it as well.

• Define two functions: makeintegrator(CDS) and maketangentintegrator(CDS) (names to be decided), that return integrators about the system or about the system+tangent space.

• Internally, all functions of DynamicalSystemsBase.jl and ChaosTools.jl (and any other packages of the ecosystem) will use step!(integrator, Δt). An argument like dt will be given to functions that actually care about dt.

• Define a function step!(integrator, Δt) that steps until a span of Δt has passed.

• All high-level functions and internals for both DynamicalSystemsBase.jl and ChaosTools.jl will not care about the tspan property of ODEProblem.

• Be sure that current time and parameters is passed properly into the Jacobian function from ForwardDiff.

• Define type aliases for CDS and DDS, e.g. CDS{ODE, A,B,C} = DS{ODE, A, B, C} where {ODE<:ODEProblem, A, B, C}. This will reduce integrator constructors by half!

• Provide high level constructors CDS(eom, state, p [, jacobian]) and CDS(prob [, jacobian]). The default tspan will be made (0.0, Inf).

• Parameterize CDS with IIP (already part of the ODEProblem).

• All of the above are meant to work for 2 cases: both in-place (iip) and out-of-place versions. This is a lot of work but the performance benefits of SVector for small systems are worth it.

• Create an oop Lorenz63 model and implement tests, similar to the tests here: https://github.com/JuliaDynamics/DynamicalSystemsBase.jl/blob/master/test/discrete_types.jl

• Create a ParallelIntegrator (this has already started in the source code) that evolves pairs of vectors of states.

• All high level functions that want to use a CDS, will first call makeintegrator(CDS) and immediatelly pass that to a lower-level function that uses the integrator and step!. This is very important for parallelization stuff.

Here is the struct:

struct ContinuousDynamicalSystem{IIP, IAD, ODE<:ODEProblem, JAC}
    prob::ODEProblem
    jacobian::JAC
end

where IIP is "isinplace" and IAD is "is auto differentiated".

If the equations for the system are out of place, the Jacobian must be as well. The same for in-place. Supporting all 4 possible variants is a lot of work that I cannot do.

EDIT MORE: Try to combine DISCRETE + Continuous into DynamicalSystem, parameterized by the ODEProblem type. Meaning:

struct DynamicalSystem{
        DEP<:DEProblem, # problem type also gives dispatch purposes
        IIP, # is in place , for dispatch purposes
        JAC, # jacobian function (either user-provided or FD)
        IAD} # is auto differentiated? Only for constructing TangentEvolver
    prob::DEP
    jacobian::JAC
end

This needs OrdinaryDiffEq first. Also, it is a good idea to do #64 in parallel with this!

As this is a repository for basic type definitions, it shouldn't have methods like these.

Should we change (for performance benefits) the parallel integrator of continuous systems to always use a matrix with hcat states?

Maybe the internals of DifferentialEquations.jl do not optimize very well vectors of vectors as state...???

# Initialize the lorenz system with random initial condition
lor = Systems.lorenz()

Q0 = eye(3)

btime lyapunovs($lor, 2000, $Q0; Ttr = 10.0);

  121.728 ms (270 allocations: 30.73 KiB)

@btime lyapunov($lor, 2000; Ttr = 10.0);

  161.619 ms (408 allocations: 33.38 KiB)

@ChrisRackauckas maybe this is interesting for you, so I'll give you some context. lyapunovs evolves a 3x4 matrix (12 elements) and does a QR decomposition 2000 times. The equations of motion is 1 application of the lorenz system and 1 application of the jacobian (matrix multiplication). The state of the system in this case is a 3x4 SMatrix.

lyapunov evolves a Vector of 2 SVectors. The equations of motion are 2 applications of the lorenz system, which I assume to be more slow than 1 lorenz and 1 jacobian. However, no QR decomposition is done, only a trivial distance computation. The state of the system in this case is a Vector{SVector}

I know that QR from StaticArrays is super fast:

s = rand(SMatrix{3,3})
@btime qr($s)

  51.059 ns (0 allocations: 0 bytes)

but I would still assume lyapunovs to be slower than lyapunov!!!
When doing

@profiler lyapunovs(lor, 2000, Q0; Ttr = 10.0);

@profiler lyapunov(lor, 2000.0; Ttr = 10.0);

in Juno, I can see that in both cases "all" time seems to be spent on perform_step! and loopfooter! so I can't really tell if this is an issue from DifferentialEquations.jl or from me. If you need more info let me know. But don't be alarmed... In the end of the day it may just be that the QR of static arrays for 3x3 is so damn fast that it is faster to evolve the tangent space than evolving 2 "normal" trajectories.