diff --git a/REQUIRE b/REQUIRE index 25c707c..c402e37 100644 --- a/REQUIRE +++ b/REQUIRE @@ -2,4 +2,3 @@ julia 0.6 StaticArrays QuadGK IterTools -Interpolations \ No newline at end of file diff --git a/test/discrete_transition_costs.jl b/test/discrete_transition_costs.jl index 0a1f9ac..9d69d9e 100644 --- a/test/discrete_transition_costs.jl +++ b/test/discrete_transition_costs.jl @@ -28,33 +28,34 @@ I1 = [1,10,30,40,50] Y1 = g[I1] y1 = interpolate((I1,), Y1, Gridded(Linear()))[1:length(g)] +@testset "Discrete Transition Costs" begin + @test compute_cost(ℓ,I1,Y1) ≈ sum((y1[1:end-1]-g[1:end-1]).^2) -@test compute_cost(ℓ,I1,Y1) ≈ sum((y1[1:end-1]-g[1:end-1]).^2) + ## Test 1b + I2 = [1,15,30,36,50] + Y2 = randn(5) -## Test 1b -I2 = [1,15,30,36,50] -Y2 = randn(5) + y2 = interpolate((t[I2],), Y2, Gridded(Linear()))[t] -y2 = interpolate((t[I2],), Y2, Gridded(Linear()))[t] + @test compute_cost(ℓ,I2,Y2) ≈ sum((y2[1:end-1]-g[1:end-1]).^2) -@test compute_cost(ℓ,I2,Y2) ≈ sum((y2[1:end-1]-g[1:end-1]).^2) + #plot(t, sin.(t)) + #plot!(t[I1], Y1) + #plot!(t, y1, color="red") -#plot(t, sin.(t)) -#plot!(t[I1], Y1) -#plot!(t, y1, color="red") + ## Test 2 + # Check the transition costs for a simple example + ℓ = dev.compute_discrete_transition_costs([1.0, 2.0, 3.0]) -## Test 2 -# Check the transition costs for a simple example -ℓ = dev.compute_discrete_transition_costs([1.0, 2.0, 3.0]) + # (y - 1)^2 + @test ℓ[1,2].P == [1.0 0.0; 0.0 0.0] + @test ℓ[1,2].q == [-2.0, 0] + @test ℓ[1,2].r == 1 -# (y - 1)^2 -@test ℓ[1,2].P == [1.0 0.0; 0.0 0.0] -@test ℓ[1,2].q == [-2.0, 0] -@test ℓ[1,2].r == 1 - -# (y - 1)^2 + (y/2 + y'/2 - 2)^2 -@test ℓ[1,3].P == [1.25 0.25; 0.25 0.25] -@test ℓ[1,3].q == [-4.0, -2.0] -@test ℓ[1,3].r == 5.0 + # (y - 1)^2 + (y/2 + y'/2 - 2)^2 + @test ℓ[1,3].P == [1.25 0.25; 0.25 0.25] + @test ℓ[1,3].q == [-4.0, -2.0] + @test ℓ[1,3].r == 5.0 +end diff --git a/test/find_optimal_fit.jl b/test/find_optimal_fit.jl index dc9edf2..f5ecb45 100644 --- a/test/find_optimal_fit.jl +++ b/test/find_optimal_fit.jl @@ -16,16 +16,19 @@ g = sin @time I_dp4, _, f_dp4 = dev.recover_solution(Λ[4, 1], 1, N) @time I_dp5, _, f_dp5 = dev.recover_solution(Λ[5, 1], 1, N) -# Comparisons to solutions found by brute force optimization -@test I_dp3 == [1, 21, 30, 50] -@test I_dp4 ∈ ([1, 8, 19, 30, 50], [1, 21, 32, 43, 50]) # Two symmetric solutions -@test I_dp5 == [1, 8, 19, 32, 43, 50] - -@test f_dp3 ≈ 1.7402065022125042 -@test f_dp4 ≈ 0.9196880458290417 -@test f_dp5 ≈ 0.09174442455649423 - -## Test of brute force optimization -m = 3 -@time I_bf, _, f_bf = dev.brute_force_optimization(ℓ, m-1); -@test I_bf == [1, 21, 30, 50] +@testset "Optimal Fit" begin + # Comparisons to solutions found by brute force optimization + @test I_dp3 == [1, 21, 30, 50] + @test I_dp4 ∈ ([1, 8, 19, 30, 50], [1, 21, 32, 43, 50]) # Two symmetric solutions + @test I_dp5 == [1, 8, 19, 32, 43, 50] + + @test f_dp3 ≈ 1.7402065022125042 + @test f_dp4 ≈ 0.9196880458290417 + @test f_dp5 ≈ 0.09174442455649423 + + ## Test of brute force optimization + m = 3 + @time I_bf, _, f_bf = dev.brute_force_optimization(ℓ, m-1); + @test I_bf == [1, 21, 30, 50] + +end diff --git a/test/piecewise_quadratics.jl b/test/piecewise_quadratics.jl index 160d385..8032a95 100644 --- a/test/piecewise_quadratics.jl +++ b/test/piecewise_quadratics.jl @@ -2,39 +2,7 @@ using Base.Test include(joinpath(Pkg.dir("DynamicApproximations"),"src","dev.jl")) -pwq = dev.generate_PiecewiseQuadratic(([2.0, 2, 1], -Inf), ([2.0, -2, 1], 0.0)) - -@test length(pwq) == 2 -dev.add_quadratic(pwq, dev.QuadraticPolynomial([1, 0, 1.0])) - -@test length(pwq) == 3 - -#--- -# Tests with (2x^2 + 1) and (x^2 + 2) which have intersections in -1, +1 -# and 5/2*x(x-1) + 1 = 2.5x^2 - 2.5x + 1 - -p1 = dev.QuadraticPolynomial([2.0, 0, 1]) -p2 = dev.QuadraticPolynomial([1.0, 0, 2]) -p3 = dev.QuadraticPolynomial(2.5, -2.5, 1.0) - -pwq1 = dev.create_new_pwq(p1) -dev.add_quadratic(pwq1, p2) -@test length(pwq1) == 3 -dev.add_quadratic(pwq1, p3) -@test length(pwq1) == 4 - -pwq2 = dev.create_new_pwq(p2) -dev.add_quadratic(pwq2, p3) -@test length(pwq2) == 3 -dev.add_quadratic(pwq2, p1) -@test length(pwq2) == 4 - -@test pwq1[1].p === pwq2[1].p == p1 -@test pwq1[2].left_endpoint == pwq2[2].left_endpoint == -1 -@test pwq1[2].p === pwq2[2].p == p2 -@test pwq1[3].left_endpoint == pwq2[3].left_endpoint == 0 -@test pwq1[3].p === pwq2[3].p - +global const TEST_PRINT_LEVEL = 0 #--- # Given a matrix with coefficients for quadratic polynomials this function @@ -55,28 +23,28 @@ function verify_piecewise_quadratic(c_mat, show_plot=false) # Insert quadratics into piecewise quadfatic pwq = dev.create_new_pwq() for p in poly_list - println("Inserting: ", p) + TEST_PRINT_LEVEL > 0 && println("Inserting: ", p) dev.add_quadratic(pwq, p) - println("After Insertion: ") - println(pwq) + TEST_PRINT_LEVEL > 0 && println("After Insertion: ") + TEST_PRINT_LEVEL > 0 && println(pwq) end # Check that the intersections between the quadratics are in increasing order - println("Checking that intersections points are increasing: ") + TEST_PRINT_LEVEL > 0 && println("Checking that intersections points are increasing: ") for λ in pwq @test λ.left_endpoint < dev.get_right_endpoint(λ) end # Check for continuity of the piecewise quadratic - println("Checking continuity: ") + TEST_PRINT_LEVEL > 0 && println("Checking continuity: ") for λ in pwq x = dev.get_right_endpoint(λ) if x == Inf break end - println("Continutiy diff: ", λ.p(x) - λ.next.p(x)) + TEST_PRINT_LEVEL > 0 && println("Continutiy diff: ", λ.p(x) - λ.next.p(x)) @test λ.p(x) ≈ λ.next.p(x) end @@ -86,40 +54,77 @@ function verify_piecewise_quadratic(c_mat, show_plot=false) y_pwq = dev.evalPwq(pwq, x_grid) for p in poly_list - @printf("suboptimality diff: %2.3f\n", maximum(y_pwq - p.(x_grid))) + TEST_PRINT_LEVEL > 0 && @printf("suboptimality diff: %2.3f\n", maximum(y_pwq - p.(x_grid))) #println(" ", x_grid[ (y_pwq - p.(x_grid)) .> 0]) @test (y_pwq .<= p.(x_grid)) == trues(length(x_grid)) end end -#--- -c_mat1 = [ -2.53106 1.59097 1.60448 -2.51349 0.681205 1.60553 -2.09364 0.714858 1.08607 -] -verify_piecewise_quadratic(c_mat1) - -#--- -c_mat2 = [ -2.01532 1.59269 1.65765 -2.50071 1.56421 1.53899 -2.02767 0.706143 1.129 -] -verify_piecewise_quadratic(c_mat2) - -#--- +pwq = dev.generate_PiecewiseQuadratic(([2.0, 2, 1], -Inf), ([2.0, -2, 1], 0.0)) -c_mat3 = [ -2.42298 0.796953 1.02011 -2.45055 1.32235 1.4894 -2.1762 0.855014 1.0647 -] -verify_piecewise_quadratic(c_mat3) -#--- -# Random test case -N = 10 -c_mat_rand = [2+rand(N) 2*rand(N) 1+rand(N)] +@testset "Piecewise Quadratics" begin + @test length(pwq) == 2 + dev.add_quadratic(pwq, dev.QuadraticPolynomial([1, 0, 1.0])) + + @test length(pwq) == 3 + + #--- + # Tests with (2x^2 + 1) and (x^2 + 2) which have intersections in -1, +1 + # and 5/2*x(x-1) + 1 = 2.5x^2 - 2.5x + 1 + + p1 = dev.QuadraticPolynomial([2.0, 0, 1]) + p2 = dev.QuadraticPolynomial([1.0, 0, 2]) + p3 = dev.QuadraticPolynomial(2.5, -2.5, 1.0) + + pwq1 = dev.create_new_pwq(p1) + dev.add_quadratic(pwq1, p2) + @test length(pwq1) == 3 + dev.add_quadratic(pwq1, p3) + @test length(pwq1) == 4 + + pwq2 = dev.create_new_pwq(p2) + dev.add_quadratic(pwq2, p3) + @test length(pwq2) == 3 + dev.add_quadratic(pwq2, p1) + @test length(pwq2) == 4 + + @test pwq1[1].p === pwq2[1].p == p1 + @test pwq1[2].left_endpoint == pwq2[2].left_endpoint == -1 + @test pwq1[2].p === pwq2[2].p == p2 + @test pwq1[3].left_endpoint == pwq2[3].left_endpoint == 0 + @test pwq1[3].p === pwq2[3].p + + + #--- + c_mat1 = [ + 2.53106 1.59097 1.60448 + 2.51349 0.681205 1.60553 + 2.09364 0.714858 1.08607 + ] + verify_piecewise_quadratic(c_mat1) + + #--- + c_mat2 = [ + 2.01532 1.59269 1.65765 + 2.50071 1.56421 1.53899 + 2.02767 0.706143 1.129 + ] + verify_piecewise_quadratic(c_mat2) + + #--- + + c_mat3 = [ + 2.42298 0.796953 1.02011 + 2.45055 1.32235 1.4894 + 2.1762 0.855014 1.0647 + ] + verify_piecewise_quadratic(c_mat3) + #--- + # Random test case + N = 10 + c_mat_rand = [2+rand(N) 2*rand(N) 1+rand(N)] + + verify_piecewise_quadratic(c_mat_rand) -verify_piecewise_quadratic(c_mat_rand) +end diff --git a/test/quadratic_forms.jl b/test/quadratic_forms.jl index c9e4d78..4f07f69 100644 --- a/test/quadratic_forms.jl +++ b/test/quadratic_forms.jl @@ -8,34 +8,35 @@ include(joinpath(Pkg.dir("DynamicApproximations"),"src","dev.jl")) # Addition of quaratic forms Q = dev.QuadraticForm2(@SMatrix([2.0 0; 0 1]), @SVector([0.0,0]), 1.0) + dev.QuadraticForm2(@SMatrix([2.0 1; 1 2]), @SVector([1.0,0]), 0.0) +@testset "Quadratic Forms" begin + @test Q == dev.QuadraticForm2(@SMatrix([4.0 1; 1 3]), @SVector([1.0,0]), 1.0) -@test Q == dev.QuadraticForm2(@SMatrix([4.0 1; 1 3]), @SVector([1.0,0]), 1.0) + # Maybe somewhat surprising that the tests pass considering equalities between + # Floats, or is the aritmhetic exact -# Maybe somewhat surprising that the tests pass considering equalities between -# Floats, or is the aritmhetic exact + # x^2 + y^2 + 1 + Q1 = dev.QuadraticForm2(@SMatrix([1.0 0; 0 0]), @SVector([0.0, 0]), 0.5) + p1 = dev.QuadraticPolynomial(1.0, 0.0, 0.5) + @test dev.minimize_wrt_x2(Q1, p1) == dev.QuadraticPolynomial(1,0,1) + @test Q1(1, 2) == 1.5 -# x^2 + y^2 + 1 -Q1 = dev.QuadraticForm2(@SMatrix([1.0 0; 0 0]), @SVector([0.0, 0]), 0.5) -p1 = dev.QuadraticPolynomial(1.0, 0.0, 0.5) -@test dev.minimize_wrt_x2(Q1, p1) == dev.QuadraticPolynomial(1,0,1) -@test Q1(1, 2) == 1.5 + # (2x - y + 1)^2 + x^2 - 2x + 4 + # = 5x^2 + y^2 + 2 + 2x - 2y - 4xy + Q2 = dev.QuadraticForm2(@SMatrix([5.0 -2; -2 0.5]), @SVector([2.0, -2]), 1.0) + p2 = dev.QuadraticPolynomial(0.5, 0.0, 5.0) + @test dev.minimize_wrt_x2(Q2, p2) == dev.QuadraticPolynomial(1,-2,5) -# (2x - y + 1)^2 + x^2 - 2x + 4 -# = 5x^2 + y^2 + 2 + 2x - 2y - 4xy -Q2 = dev.QuadraticForm2(@SMatrix([5.0 -2; -2 0.5]), @SVector([2.0, -2]), 1.0) -p2 = dev.QuadraticPolynomial(0.5, 0.0, 5.0) -@test dev.minimize_wrt_x2(Q2, p2) == dev.QuadraticPolynomial(1,-2,5) + # x^2 + 2x + 2 + Q3 = dev.QuadraticForm2(@SMatrix([1.0 0; 0 0]), @SVector([2.0, 0]), 2.0) + p3 = dev.QuadraticPolynomial(0.0, 0.0, 0.0) + @test dev.minimize_wrt_x2(Q3, p3) == dev.QuadraticPolynomial(1.0,2.0,2.0) -# x^2 + 2x + 2 -Q3 = dev.QuadraticForm2(@SMatrix([1.0 0; 0 0]), @SVector([2.0, 0]), 2.0) -p3 = dev.QuadraticPolynomial(0.0, 0.0, 0.0) -@test dev.minimize_wrt_x2(Q3, p3) == dev.QuadraticPolynomial(1.0,2.0,2.0) - -# y^2 + 2y + 2 -Q4 = dev.QuadraticForm2(@SMatrix([0.0 0; 0 1]), @SVector([0.0, 2]),2.0) -p4 = dev.QuadraticPolynomial(0.0, 0.0, 0.0) -@test dev.minimize_wrt_x2(Q4, p4) == dev.QuadraticPolynomial(0,0,1) + # y^2 + 2y + 2 + Q4 = dev.QuadraticForm2(@SMatrix([0.0 0; 0 1]), @SVector([0.0, 2]),2.0) + p4 = dev.QuadraticPolynomial(0.0, 0.0, 0.0) + @test dev.minimize_wrt_x2(Q4, p4) == dev.QuadraticPolynomial(0,0,1) +end diff --git a/test/quadratic_polynomials.jl b/test/quadratic_polynomials.jl index b6341db..478de3d 100644 --- a/test/quadratic_polynomials.jl +++ b/test/quadratic_polynomials.jl @@ -2,25 +2,27 @@ using Base.Test include(joinpath(Pkg.dir("DynamicApproximations"),"src","dev.jl")) -# Subtraction of quadratic polynomials -@test dev.QuadraticPolynomial(4.0, -2.0, 3.0) - dev.QuadraticPolynomial(1.0, 1.0, 2.0) == dev.QuadraticPolynomial(3.0,-3.0,1.0) +@testset "Quadratic Polynomials" begin + # Subtraction of quadratic polynomials + @test dev.QuadraticPolynomial(4.0, -2.0, 3.0) - dev.QuadraticPolynomial(1.0, 1.0, 2.0) == dev.QuadraticPolynomial(3.0,-3.0,1.0) -### Minimization of quadratic polynomials ### -@test dev.find_minimum(dev.QuadraticPolynomial(1 ,0 ,0)) == (0, 0) -@test dev.find_minimum(dev.QuadraticPolynomial(1 ,0 ,1)) == (0, 1) + ### Minimization of quadratic polynomials ### + @test dev.find_minimum(dev.QuadraticPolynomial(1 ,0 ,0)) == (0, 0) + @test dev.find_minimum(dev.QuadraticPolynomial(1 ,0 ,1)) == (0, 1) -# 2*x*(x+1) + 1 -@test dev.find_minimum(dev.QuadraticPolynomial(2, 2, 1)) == (-0.5, 0.5) + # 2*x*(x+1) + 1 + @test dev.find_minimum(dev.QuadraticPolynomial(2, 2, 1)) == (-0.5, 0.5) -# (x+1)^2 + 1 -@test dev.find_minimum(dev.QuadraticPolynomial(1, 2, 2)) == (-1, 1) + # (x+1)^2 + 1 + @test dev.find_minimum(dev.QuadraticPolynomial(1, 2, 2)) == (-1, 1) -# 2*((x+1)^2 + 1) -@test dev.find_minimum(dev.QuadraticPolynomial(2, 4, 4)) == (-1, 2) + # 2*((x+1)^2 + 1) + @test dev.find_minimum(dev.QuadraticPolynomial(2, 4, 4)) == (-1, 2) -# (3x+1)^2 + 1 = 9x^2 + 6x + 2 -@test dev.find_minimum(dev.QuadraticPolynomial(9, 6, 2)) == (-1/3, 1) + # (3x+1)^2 + 1 = 9x^2 + 6x + 2 + @test dev.find_minimum(dev.QuadraticPolynomial(9, 6, 2)) == (-1/3, 1) -# (3x+1)^2 + 1 = 9x^2 + 6x + 2 -@test dev.unsafe_minimum(dev.QuadraticPolynomial(9.0, 6.0, 2.0)) == 1.0 + # (3x+1)^2 + 1 = 9x^2 + 6x + 2 + @test dev.unsafe_minimum(dev.QuadraticPolynomial(9.0, 6.0, 2.0)) == 1.0 +end diff --git a/test/runtests.jl b/test/runtests.jl index 5c70632..55cbb02 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,9 +1,13 @@ #using DynamicApproximations using Base.Test -include("quadratic_polynomials.jl") -include("quadratic_forms.jl") -include("piecewise_quadratics.jl") -include("transition_costs.jl") -include("discrete_transition_costs.jl") -include("find_optimal_fit.jl") +include(joinpath(Pkg.dir("DynamicApproximations"),"src","dev.jl")) + +@testset "All Tests" begin + include("quadratic_polynomials.jl") + include("quadratic_forms.jl") + include("piecewise_quadratics.jl") + include("transition_costs.jl") + include("discrete_transition_costs.jl") + include("find_optimal_fit.jl") +end diff --git a/test/transition_costs.jl b/test/transition_costs.jl index 627e7da..f259b19 100644 --- a/test/transition_costs.jl +++ b/test/transition_costs.jl @@ -3,27 +3,28 @@ using Base.Test include(joinpath(Pkg.dir("DynamicApproximations"),"src","dev.jl")) t = [-1; 0; 1] +@testset "Transition Costs" begin + g = x -> 0 + ℓ = dev.compute_transition_costs(g, t) + @test ℓ[1,2].P == 1/6*[2 1; 1 2] + @test ℓ[1,2].q == [0,0] + @test ℓ[1,2].r == 0 -g = x -> 0 -ℓ = dev.compute_transition_costs(g, t) -@test ℓ[1,2].P == 1/6*[2 1; 1 2] -@test ℓ[1,2].q == [0,0] -@test ℓ[1,2].r == 0 + g = x -> 1 + ℓ = dev.compute_transition_costs(g, t) + @test ℓ[1,2].P == 1/6*[2 1; 1 2] + @test ℓ[1,2].q == [-1,-1] + @test ℓ[1,2].r == 1 -g = x -> 1 -ℓ = dev.compute_transition_costs(g, t) -@test ℓ[1,2].P == 1/6*[2 1; 1 2] -@test ℓ[1,2].q == [-1,-1] -@test ℓ[1,2].r == 1 - -g = x -> x -ℓ = dev.compute_transition_costs(g, t) -@test ℓ[2,3].P == 1/6*[2 1; 1 2] -@test ℓ[2,3].q ≈ [-1/3,-2/3] -@test ℓ[2,3].r ≈ 1/3 -@test ℓ[2,3](0,0) == 1/3 -@test ℓ[2,3](1,1) == 1/3 -@test ℓ[2,3](1,2) ≈ 1 -@test ℓ[2,3](1,0) ≈ 1/3 + g = x -> x + ℓ = dev.compute_transition_costs(g, t) + @test ℓ[2,3].P == 1/6*[2 1; 1 2] + @test ℓ[2,3].q ≈ [-1/3,-2/3] + @test ℓ[2,3].r ≈ 1/3 + @test ℓ[2,3](0,0) == 1/3 + @test ℓ[2,3](1,1) == 1/3 + @test ℓ[2,3](1,2) ≈ 1 + @test ℓ[2,3](1,0) ≈ 1/3 +end