avoid out of bounds in _frequencies! (#19)

* avoid out of bounds in _frequencies!

* make trajectory use Vern9

* use absolutely identical keyword arguments

* relax the small lorenz delay test

* Be sure to also test the code of #18

src/estimate_delay.jl

@@ -266,7 +266,7 @@ function _frequencies!(f::AbstractVector{T}, edges::AbstractVector,
         end
         # when `s[i]` goes after the upper limit of the current bin,
         # move to the next bin and repeat until `s[i]` is found
-        while s[i] > edges[b+1]
+        while s[i] > edges[b+1] && b < nbins
             b += 1
         end
         f[b] += 1 # add point to the current bin

test/REQUIRE

@@ -1 +1,2 @@
 ChaosTools
+OrdinaryDiffEq

test/R_delay.jl

@@ -1,5 +1,5 @@
 using ChaosTools, DelayEmbeddings
-using Test
+using Test, OrdinaryDiffEq
 using DelayEmbeddings: localextrema, exponential_decay_fit
 
 println("\nTesting delay estimation...")
@@ -23,22 +23,31 @@ testval = (val, vmin, vmax) -> @test vmin ≤ val ≤ vmax
 
     ds = Systems.roessler(ones(3))
     dt = 0.02
-    data = trajectory(ds,500,dt=dt)
+    data = trajectory(ds,500.0;dt=dt,diffeq...)
     x = data[:,1]
     @test 1.3 ≤ estimate_delay(x,"ac_zero", 1:2:500)*dt ≤ 1.7
     @test 2.6 ≤ estimate_delay(x,"ac_min", 1:2:500)*dt  ≤ 3.4
-    @test 1.15 ≤ estimate_delay(x,"mi_min", 1:2:500)*dt ≤ 1.6
+    @test 1.0 ≤ estimate_delay(x,"mi_min", 1:2:500)*dt ≤ 1.6
 
     dt = 0.1
-    data = trajectory(ds,2000,dt=dt)
+    data = trajectory(ds,2000.0;dt=dt,diffeq...)
     x = data[:,1]
     @test 1.3 ≤ estimate_delay(x,"ac_zero", 1:1:50)*dt ≤ 1.7
     @test 2.6 ≤ estimate_delay(x,"ac_min", 1:1:50)*dt  ≤ 3.4
     @test 1.15 ≤ estimate_delay(x,"mi_min", 1:1:50)*dt ≤ 1.6
 
     ds = Systems.lorenz()
     dt = 0.1
-    data = trajectory(ds,5000;dt=dt)
+    data = trajectory(ds,5000;dt=dt,diffeq...)
     x = data[500:end,1]
     @test 0 < estimate_delay(x,"exp_extrema", 0:2:200)  < 200
+
+    # Issue #18
+    ds = Systems.gissinger(ones(3)) # 3D continuous chaotic system, also shown in orbit diagrams tutorial
+    dt = 0.05
+    data = trajectory(ds, 1000.0, dt = dt)
+    s = data[:, 1]
+
+    τ = estimate_delay(s, "mi_min", 0:1:400) # this was the erroring line
+    @test τ > 0
 end

test/R_dimension.jl

@@ -1,4 +1,4 @@
-using DelayEmbeddings, ChaosTools
+using DelayEmbeddings, ChaosTools, OrdinaryDiffEq
 using Test
 
 println("\nTesting R-param. estimation...")
@@ -30,7 +30,7 @@ test_value = (val, vmin, vmax) -> @test vmin <= val <= vmax
     # @test minimum(E2s) < 0.1 # THIS TEST FAILS
 
     ds = Systems.roessler();τ=15; dt=0.1
-    data = trajectory(ds,1000;dt=dt)
+    data = trajectory(ds,1000;dt=dt,diffeq...)
     s_roessler = data[:,1]
     Ds = 1:5
     E1s = DelayEmbeddings.estimate_dimension(s_roessler, τ, Ds)
@@ -42,7 +42,7 @@ test_value = (val, vmin, vmax) -> @test vmin <= val <= vmax
     @test saturation_point(Ds,E1s; threshold=0.1) ∈ [2, 3]
 
     ds = Systems.lorenz();τ=5; dt=0.01
-    data = trajectory(ds,500;dt=dt)
+    data = trajectory(ds,500;dt=dt,diffeq...)
     s_lorenz = data[:,1]
     Ds = 1:5
     E1s = DelayEmbeddings.estimate_dimension(s_lorenz, τ, Ds)

test/neighborhood_tests.jl

@@ -41,7 +41,7 @@ println("\nTesting neighborhoods...")
         end
 
         ds = Systems.lorenz()
-        data = trajectory(ds, 200, dt = 0.001)
+        data = trajectory(ds, 200, dt = 0.001, alg = Vern9())
         tree = KDTree(data)
 
         point1 = data[10000]

test/runtests.jl

@@ -1,7 +1,9 @@
-using DelayEmbeddings
+using DelayEmbeddings, OrdinaryDiffEq
 
 ti = time()
 
+const diffeq = (alg = Vern9(), atol = 1e-9, rtol = 1e-9, maxiters = typemax(Int))
+
 include("dataset_tests.jl")
 include("reconstruction_tests.jl")
 include("R_delay.jl")