juliadynamics / dynamicalsystemsbase.jl Goto Github PK
View Code? Open in Web Editor NEWDefinition of dynamical systems and integrators for DynamicalSystems.jl
License: Other
Definition of dynamical systems and integrators for DynamicalSystems.jl
License: Other
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>
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:
Pkg.resolve()
.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:
Pkg.resolve()
.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>
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.
HDS
?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.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.
Allowing to evolve the system for some transient amount of time Ttr
and then recording would be very convenient.
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.
If we do this, it will also be easy to use callbacks etc with diffeq.
Will also solve JuliaDynamics/DynamicalSystems.jl#14
We will do this when moving to 0.7!
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
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).
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
And change default tolerances to 1e-9
.
The only thing that has to be done is to change function get_solver
!
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}
.
As this is a repository for basic type definitions, it shouldn't have methods like these.
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
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.
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.
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.
This would make lyapunovs
and gali
also conserve energy.
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!
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
Makes 0 sense.
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.
Easy follow to #28
A dataset has the dimension as type-information, so no need for SizedArray
here
This needs OrdinaryDiffEq first. Also, it is a good idea to do #64 in parallel with this!
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:
We want a continuous system ODE solver that has the following properties:
init
from DiffEqBase
to create/initialize the integrator object given an ODEProblem
.step!
from DiffEqBase
to "step" the integrator for one step. Adaptive stepping and error checking must be included.step!
. Only the previous state and current state need to be stored.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.OrdinaryDiffEq
does this already.reinit!
but also allows one to change the state during integration which is crucial for JuliaDynamics.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 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.
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?
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.
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
DynamicalSystem
MMatrix
when the state is MVector
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!
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.
Ooookay.
I succumb to the peer pressure. Time to change this...
Test the columns
function using the gissinger
model.
Maybe its good to define substraction, addition etc. between Dataset
s. 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.
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]
Prof. Ulrich Parlitz knows this.
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.
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.
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.
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.
the DynamicalSystem
and DiscreteDynamicalSystem
implementations are pure (only StaticArray dependence).
They can easily be updated! While this is done, be sure that:
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 .
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)
To allow identical interface between both system types, continuouys and disscrete
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.
Implement what is discussed here: JuliaDiff/ForwardDiff.jl#296
This was present in older versions and I removed it in humongous lapse of judgement.
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.