Add a test to makesure that the previous bug can't happen anymore!

nathanaelbosch · Oct 10, 2023 · d3c7d45 · d3c7d45
1 parent d0e590f
commit d3c7d45
Showing 1 changed file with 81 additions and 55 deletions.
diff --git a/test/core/priors.jl b/test/core/priors.jl
@@ -1,10 +1,11 @@
+using ProbNumDiffEq: make_transition_matrices!
 using ProbNumDiffEq
 import ProbNumDiffEq as PNDE
 using Test
 using LinearAlgebra
 
-h = rand()
-σ = rand()
+h = 0.1
+σ = 0.1
 
 @testset "General prior API" begin
     for prior in (IWP(2, 3), IOUP(2, 3, 1), Matern(2, 3, 1))
@@ -26,24 +27,64 @@ end
 
     prior = PNDE.IWP(d, q)
 
+    # true sde parameters
+    F = [
+        0 1 0 0 0 0
+        0 0 1 0 0 0
+        0 0 0 0 0 0
+        0 0 0 0 1 0
+        0 0 0 0 0 1
+        0 0 0 0 0 0
+    ]
+    L = [
+        0 0
+        0 0
+        1 0
+        0 0
+        0 0
+        0 1
+    ]
+
+    # true transition matrices
+    AH_22_IBM = [
+        1 h h^2/2 0 0 0
+        0 1 h 0 0 0
+        0 0 1 0 0 0
+        0 0 0 1 h h^2/2
+        0 0 0 0 1 h
+        0 0 0 0 0 1
+    ]
+
+    QH_22_IBM = σ^2 .* [
+        h^5/20 h^4/8 h^3/6 0 0 0
+        h^4/8 h^3/3 h^2/2 0 0 0
+        h^3/6 h^2/2 h 0 0 0
+        0 0 0 h^5/20 h^4/8 h^3/6
+        0 0 0 h^4/8 h^3/3 h^2/2
+        0 0 0 h^3/6 h^2/2 h
+    ]
+
+    # true preconditioned transition matrices
+    AH_22_PRE = [
+        1 2 1 0 0 0
+        0 1 1 0 0 0
+        0 0 1 0 0 0
+        0 0 0 1 2 1
+        0 0 0 0 1 1
+        0 0 0 0 0 1
+    ]
+
+    QH_22_PRE = σ^2 * [
+        1/5 1/4 1/3 0 0 0
+        1/4 1/3 1/2 0 0 0
+        1/3 1/2 1/1 0 0 0
+        0 0 0 1/5 1/4 1/3
+        0 0 0 1/4 1/3 1/2
+        0 0 0 1/3 1/2 1/1
+    ]
+
     @testset "Test SDE" begin
         sde = PNDE.to_sde(prior)
-        F = [
-            0 1 0 0 0 0
-            0 0 1 0 0 0
-            0 0 0 0 0 0
-            0 0 0 0 1 0
-            0 0 0 0 0 1
-            0 0 0 0 0 0
-        ]
-        L = [
-            0 0
-            0 0
-            1 0
-            0 0
-            0 0
-            0 1
-        ]
         @test sde.F ≈ F
         @test sde.L ≈ L
     end
@@ -52,54 +93,39 @@ end
         Ah, Qh = PNDE.discretize(prior, h)
         Qh = PNDE.apply_diffusion(Qh, σ^2)
 
-        AH_22_IBM = [
-            1 h h^2/2 0 0 0
-            0 1 h 0 0 0
-            0 0 1 0 0 0
-            0 0 0 1 h h^2/2
-            0 0 0 0 1 h
-            0 0 0 0 0 1
-        ]
         @test AH_22_IBM ≈ Ah
-
-        QH_22_IBM =
-            σ^2 .* [
-                h^5/20 h^4/8 h^3/6 0 0 0
-                h^4/8 h^3/3 h^2/2 0 0 0
-                h^3/6 h^2/2 h 0 0 0
-                0 0 0 h^5/20 h^4/8 h^3/6
-                0 0 0 h^4/8 h^3/3 h^2/2
-                0 0 0 h^3/6 h^2/2 h
-            ]
         @test QH_22_IBM ≈ Matrix(Qh)
     end
 
     @testset "Test with preconditioning" begin
         A, Q = PNDE.preconditioned_discretize(prior)
         Qh = PNDE.apply_diffusion(Q, σ^2)
 
-        AH_22_PRE = [
-            1 2 1 0 0 0
-            0 1 1 0 0 0
-            0 0 1 0 0 0
-            0 0 0 1 2 1
-            0 0 0 0 1 1
-            0 0 0 0 0 1
-        ]
-
-        QH_22_PRE =
-            σ^2 * [
-                1/5 1/4 1/3 0 0 0
-                1/4 1/3 1/2 0 0 0
-                1/3 1/2 1/1 0 0 0
-                0 0 0 1/5 1/4 1/3
-                0 0 0 1/4 1/3 1/2
-                0 0 0 1/3 1/2 1/1
-            ]
-
         @test AH_22_PRE ≈ Matrix(A)
         @test QH_22_PRE ≈ Matrix(Qh)
     end
+
+    @testset "Test `make_transition_matrices!`" begin
+        struct DummyCache <: PNDE.AbstractODEFilterCache
+            d::Any
+            q::Any
+            A::Any
+            Q::Any
+            P::Any
+            PI::Any
+            Ah::Any
+            Qh::Any
+        end
+
+        A, Q, Ah, Qh, P, PI = PNDE.initialize_transition_matrices(prior, h)
+        @test AH_22_PRE ≈ A
+        @test QH_22_PRE ≈ Matrix(PNDE.apply_diffusion(Q, σ^2))
+
+        cache = DummyCache(d, q, A, Q, P, PI, Ah, Qh)
+        make_transition_matrices!(cache, prior, h)
+        @test AH_22_IBM ≈ Ah
+        @test QH_22_IBM ≈ Matrix(PNDE.apply_diffusion(cache.Qh, σ^2))
+    end
 end
 
 @testset "Test IOUP (d=2,q=2)" begin