Merge pull request #180 from JuliaDynamics/TestGenericEstimate

Test generic estimate

Author: orkolorko

examples/time_convergence_plot2.jl

@@ -137,7 +137,7 @@ function save_coarse_data(prefix, K)
     compute_coarse_grid_quantities(f, 4; m = 24)  # warm-up for precompilation
     for n in 2 .^ K
         print("Coarse+fine $n...")
-        F, C = compute_coarse_grid_quantities(f, n; m = 24)
+        C, F = compute_coarse_grid_quantities(f, n; m = 24)
         serialize("$prefix-$n-coarse.juliaserialize", C)
         serialize("$prefix-$n-fine.juliaserialize", F)
         println("done.")

src/Basis/IntervalHatBasis.jl

@@ -151,9 +151,9 @@ end
 
 Base.length(S::AverageZero{HatNP{T}}) where {T} = length(S.basis) - 1
 
-weak_projection_error(B::HatNP) = 0.5 ⊘₊ Float64(length(B)-1, RoundDown)
-aux_normalized_projection_error(B::HatNP) = 0.5 ⊘₊ Float64(length(B)-1, RoundDown)
-strong_weak_bound(B::HatNP) = 2.0 ⊗₊ Float64(length(B)-1, RoundDown)
+weak_projection_error(B::HatNP) = 0.5 ⊘₊ Float64(length(B) - 1, RoundDown)
+aux_normalized_projection_error(B::HatNP) = 0.5 ⊘₊ Float64(length(B) - 1, RoundDown)
+strong_weak_bound(B::HatNP) = 2.0 ⊗₊ Float64(length(B) - 1, RoundDown)
 aux_weak_bound(B::HatNP) = 1.0
 weak_by_strong_and_aux_bound(B::HatNP) = (1.0, 1.0)
 bound_weak_norm_from_linalg_norm(B::HatNP) = @error "TODO"
@@ -238,7 +238,8 @@ end
 
 function Dual(B::HatNP, D; ϵ = 0.0, max_iter = 100)
     @assert is_full_branch(D)
-    results = collect(HatNPDualBranch(B.p, b, 1:length(B.p)-1; ϵ, max_iter) for b in branches(D))
+    results =
+        collect(HatNPDualBranch(B.p, b, 1:length(B.p)-1; ϵ, max_iter) for b in branches(D))
     x = vcat((result[1] for result in results)...)
     xlabel = vcat((result[2] for result in results)...)
     x′ = vcat((result[3] for result in results)...)

src/GenericEstimate.jl

@@ -105,16 +105,18 @@ end
 #
 # @deprecate contractmatrix norms_of_powers
 
-"""
-This function returns the bound on the weak norm of the discretized operator
-"""
-function boundnorm(B::Basis, P::AbstractMatrix{Interval{T}}, m) where {T}
-    W₁, W₂ = bound_weak_norm_from_linalg_norm(B)
-    α₁ = bound_linalg_norm_L1_from_weak(B)
-    α₂ = bound_linalg_norm_L∞_from_weak(B)
-    C, tilde_C = contractmatrix(B, P, m)
-    return (W₁ / α₁) * C + (W₂ / α₂) * tilde_C
-end
+# """
+# This function returns the bound on the weak norm of the discretized operator
+# """
+# function boundnorm(B::Basis, P::AbstractMatrix{Interval{T}}, m) where {T}
+#     W₁, W₂ = bound_weak_norm_from_linalg_norm(B)
+#     α₁ = bound_linalg_norm_L1_from_weak(B)
+#     α₂ = bound_linalg_norm_L∞_from_weak(B)
+#     C, tilde_C = contractmatrix(B, P, m)
+#     return (W₁ / α₁) * C + (W₂ / α₂) * tilde_C
+# end
+# 
+# @deprecate boundnorm
 
 """
 Uses different strategies to compute power norm bounds.
@@ -128,7 +130,7 @@ A vector of length m_extend is returned, such that norms[k] ≥ ||Q_h^k|_{U_h^0}
 function powernormbounds(B, D, m, m_extend; Q = DiscretizedOperator(B, D))
     normQ = opnormbound(B, weak_norm(B), Q)
     trivial_norms = norms_of_powers_trivial(normQ, m)
-    computed_norms = norms_of_powers(weak_norm(B), m, Q, integral_covector(B))
+    computed_norms = norms_of_powers(B, weak_norm(B), m, Q, integral_covector(B))
 
     # not interesting at the moment
     #(dfly_strongs, dfly_norms) = norms_of_powers_dfly(B, D, m)
@@ -270,7 +272,16 @@ function compute_coarse_grid_quantities(f, n; m = 8)
         residualbound(B, weak_norm(B), Q, w), mag(integral_covector(B) * w - 1)
     time_norms = @elapsed norms = powernormbounds(B, D, Q = Q, m = m)
     time_dfly = @elapsed dfly_coefficients = dfly(strong_norm(B), aux_norm(B), D)
-    return FineGridQuantities(
+    return CoarseGridQuantities(
+        B,
+        D,
+        norms,
+        dfly_coefficients,
+        time_assembling1,
+        time_norms,
+        time_dfly,
+    ),
+    FineGridQuantities(
         B,
         D,
         normQ,
@@ -279,15 +290,6 @@ function compute_coarse_grid_quantities(f, n; m = 8)
         ε₂,
         time_assembling1 + time_assembling2,
         time_eigen1 + time_eigen2,
-    ),
-    CoarseGridQuantities(
-        B,
-        D,
-        norms,
-        dfly_coefficients,
-        time_assembling1,
-        time_norms,
-        time_dfly,
     )
 end
 

test/TestAssembleHat.jl test/TestBasis/TestAssembleHat.jltest/TestAssembleHat.jl renamed to test/TestBasis/TestAssembleHat.jl

@@ -1,4 +1,8 @@
 include("TestUlam.jl")
+include("TestAssembleHat.jl")
+include("TestHat.jl")
+include("TestHatNP.jl")
+include("TestUlam2DSP.jl")
 
 include("TestFourierCommon.jl")
 include("TestFourierAdjoint.jl")

test/TestBasis/TestBasisIndex.jl

@@ -60,13 +60,13 @@
     v = ones(4)
     @test RigorousInvariantMeasures.change_of_basis(BU, BH, v) == ones(5)
 
-    
+
     @test RigorousInvariantMeasures.evaluate_integral(B, 1, Float64) == 0.25
-    
+
     Q = DiscretizedOperator(B, D)
 
-    @test Q.L[:, 1] == [0.5; 0.25; 0;  0; 0]
-    @test Q.L[:, 2] == [0.0; 0.25; 0.5; 0.25;  0]
+    @test Q.L[:, 1] == [0.5; 0.25; 0; 0; 0]
+    @test Q.L[:, 2] == [0.0; 0.25; 0.5; 0.25; 0]
     @test Q.L[:, 3] == [0.5; 0.25; 0.0; 0.25; 0.5]
     @test Q.L[:, 4] == [0; 0.25; 0.5; 0.25; 0.0]
     @test Q.L[:, 5] == [0; 0; 0; 0.25; 0.5]

test/TestHat.jl test/TestBasis/TestHat.jltest/TestHat.jl renamed to test/TestBasis/TestHat.jl

@@ -1,9 +1,29 @@
 @testset "Estimator" begin
     using IntervalArithmetic
 
-    D = mod1_dynamic(x -> 2 * x)
-    B = Ulam(32)
-    P = RigorousInvariantMeasures.assemble(B, D; ϵ = 10^(-15), max_iter = 100, T = Float64)
+
+    function twoxexperiment(n)
+        D = mod1_dynamic(x -> 2 * x)
+        B = Ulam(32)
+        Q = DiscretizedOperator(B, D)
+        return B, D, Q
+    end
+
+    FGQ, CGQ = RigorousInvariantMeasures.compute_coarse_grid_quantities(twoxexperiment, 64)
+
+    Q = DiscretizedOperator(CGQ.B, CGQ.D)
+    norms = powernormbounds(CGQ.B, CGQ.D; Q = Q)
+
+    @test norms == CGQ.norms
+
+    norms_other = RigorousInvariantMeasures.powernormbounds(CGQ.B, CGQ.D, 8, 16; Q = Q)
+    @test norms_other == norms
+
+    error, times = one_grid_estimate(CGQ, FGQ)
+    @test error == 0.0
+
+    error, times = two_grid_estimate(CGQ, FGQ)
+    @test error == 0.0
 
     # C = RigorousInvariantMeasures.boundnorm(B, P, 10)