tests pass, ready to merge

I will only check if documentation makes sense before merging

Author: Datseris

CHANGELOG.md

@@ -1,14 +1,14 @@
 # v0.6
 ## Massively Breaking
-* All systems now have the parameter container as the third argument, and
-  mandatory as well.
 * Complete overhaul of Discrete Systems: all systems are now one, there is
   no longer `DDS` and `BigDDS`. Everything is `DDS`!
 * `variational_integrator` renamed to `tangent_integrator`.
+* `DDS` is made up of 2 fundamental parts: a DiscreteProblem and a
+  TangentEvolver.
 ## Non-breaking
-* Massive performance boost of up to 8x in system evolution.
-* `ParallelEvolver` for `DDS`
-* `TangentEvolver` for `DDS`
+* Massive performance boost of up to 8x in system evolution of continuous systems.
+* `ParallelEvolver`
+* `TangentEvolver`
 
 # v0.5.0
 

src/discrete.jl

@@ -633,8 +633,9 @@ evolve(ds::DDS, u) = evolve(ds.prob, u)
 trajectory(ds::DDS, args...) = trajectory(ds.prob, args...)
 
 """
-    jacobian(ds::DynamicalSystem, u = state(ds))
-Return the Jacobian of the equations of motion at `u`.
+    jacobian([J, ] ds::DynamicalSystem, u = state(ds))
+Return the Jacobian of the equations of motion at `u`, optionally writting the
+result in-place in `J`.
 """
 function jacobian(ds::DDS{true}, u = state(ds))
     ds.tangent.jacobian(ds.tangent.J, u, ds.prob.p)
@@ -643,6 +644,10 @@ end
 
 jacobian(ds::DDS{false}, u = state(ds)) = ds.tangent.jacobian(u, ds.prob.p)
 
+function jacobian(J, ds::DDS, u = state(ds))
+    ds.tangent.jacobian(J, u, ds.prob.p)
+end
+
 #####################################################################################
 #                                Pretty-Printing                                    #
 #####################################################################################

src/famous_systems.jl

@@ -353,7 +353,7 @@ Default values are the ones used in the original paper.
 [1] : O. E. Rössler, Phys. Lett. **71A**, pp 155 (1979)
 """
 function towel(u0=[0.085, -0.121, 0.075])
-        return DDS(u0, eom_towel, nothing, jacob_towel)
+    return DDS(u0, eom_towel, jacob_towel)
 end# should result in lyapunovs: [0.432207,0.378834,-3.74638]
 @inline function eom_towel(x, p)
     @inbounds x1, x2, x3 = x[1], x[2], x[3]
@@ -402,7 +402,7 @@ Nuclear Physics, Novosibirsk (1969)
 [2] : J. M. Greene, J. Math. Phys. **20**, pp 1183 (1979)
 """
 function standardmap(u0=0.001rand(2); k = 0.971635)
-    return DDS(u0, standardmap_eom, [k], standardmap_jacob)
+    return DDS(u0, standardmap_eom, standardmap_jacob; p = [k])
 end
 @inline @inbounds function standardmap_eom(x, par)
     theta = x[1]; p = x[2]
@@ -457,7 +457,7 @@ function coupledstandardmaps(M::Int, u0 = 0.001rand(2M);
         J[i+M, i+M] = 1
     end
 
-    return DDS(u0, csm, [ks, Γ], csm, J)
+    return DDS(u0, csm, csm, J; p = [ks, Γ])
 end
 mutable struct CoupledStandardMaps{N}
     idxs::UnitRange{Int}
@@ -567,7 +567,7 @@ function's documentation string.
 [2] : M. J. Feigenbaum, J. Stat. Phys. **19**, pp 25 (1978)
 """
 function logistic(x0=rand(); r = 4.0)
-    return DDS(x0, logistic_eom, [r], logistic_jacob)
+    return DDS(x0, logistic_eom, logistic_jacob; p = [r])
 end
 @inline logistic_eom(x, p) = p[1]*x*(1-x)
 @inline logistic_jacob(x, p) = p[1]*(1-2x)
@@ -586,7 +586,7 @@ The parameter container has the parameters in the same order as stated in this
 function's documentation string.
 """
 function circlemap(x0=rand(); K = 0.99, Ω = 1.0, usemod::Bool = false)
-    return DDS(x0, circlemap_eom, [Ω, K], circlemap_jacob)
+    return DDS(x0, circlemap_eom, circlemap_jacob; p = [Ω, K])
 end
 @inline circlemap_eom(x, p) = x + twopi*p[1] - p[2]*sin(x)
 @inline circlemap_jacob(x, p) = -p[2]*cos(x)

test/discrete_systems.jl

@@ -7,8 +7,8 @@ using Base.Test, StaticArrays
 @testset "Logistic Map" begin
 
   d1 = Systems.logistic(0.1)
-  d2 = DDS(0.1, d1.prob.f, d1.prob.p)
-  d3 = DDS(big(0.1), d1.prob.f, d1.prob.p, d1.jacobian)
+  d2 = DDS(0.1, d1.prob.f; p = d1.prob.p)
+  d3 = DDS(big(0.1), d1.prob.f, d1.tangent.jacobian; p =  d1.prob.p)
 
   @testset "Evolution & trajectory" begin
     st1 = evolve(d1, 1)
@@ -42,8 +42,8 @@ end
     s1 = Systems.henon(0.1ones(N))
   end
 
-  s2 = DDS(0.1ones(N), s1.prob.f, s1.prob.p)
-  s4 = DDS(round.(big.(0.1ones(N)),3), s1.prob.f, s1.prob.p, s1.jacobian)
+  s2 = DDS(0.1ones(N), s1.prob.f; p = s1.prob.p)
+  s4 = DDS(big.(0.1ones(N)), s1.prob.f, s1.tangent.jacobian; p = s1.prob.p)
 
   @testset "Evolution & trajectory" begin
     st1 = evolve(s1, 1)
@@ -82,11 +82,10 @@ end
     st1 = evolve(ds, 100)
 
     @test st1 != u0
-    @test u0 == state(ds)
+    evolve!(ds, 100)
+    @test st1 == state(ds)
 
-    Jbef = deepcopy(ds.J)
-    ds.jacobian(ds.J, evolve(ds, 1), ds.prob.p)
-    @test Jbef != ds.J
-    ds.jacobian(Jbef, evolve(ds, 1), ds.prob.p)
-    @test Jbef == ds.J
+    Jbef = deepcopy(ds.tangent.J)
+    jacobian(Jbef, ds)
+    @test Jbef != ds.tangent.J
 end