Improved code.

README.md

@@ -40,7 +40,7 @@ plot(g, t, l=(2,:black), lab="sin(x) + 0.5sin(3.5x) + 0.5sin(5.1x)")
 for ζ ∈ [0.1, 0.002]
     # Minimize ∫(f(x)-g(x))²dx + ζ⋅||d²f/dx²||₀
     # Will automatically integrate the function to compute the costs
-    I, Y, cost = fit_pwl_reguralized(g, t, ζ)
+    I, Y, cost = fit_pwl_regularized(g, t, ζ)
 
     plot!(t[I], Y, l=2, m=:circle, lab = "l2-norm=$(round(cost,3)), zeta=$ζ")
 end

src/DynamicApproximations.jl

@@ -1,7 +1,7 @@
 module DynamicApproximations
 
 export QuadraticPolynomial, PiecewiseQuadratic, QuadraticForm
-export fit_pwl_constrained, fit_pwl_reguralized
+export fit_pwl_constrained, fit_pwl_regularized
 export pwq_dp_regularized, pwq_dp_constrained
 #Should we export the following?
 export compute_transition_costs, compute_discrete_transition_costs
@@ -29,13 +29,12 @@ include(joinpath("types","PiecewiseQuadratic.jl"))
 include(joinpath("types","QuadraticForm.jl"))
 include(joinpath("types","transition_costs.jl"))
 
-include("brute_force_search.jl")
-
-
 AbstractTransitionCost{T} = Union{Array{QuadraticForm{T},2}, TransitionCostDiscrete{T}, TransitionCostContinuous{T}}
 
-include("solve.jl")
 include("transition_cost_computation.jl")
+include("brute_force_search.jl")
+include("solve.jl")
+
 
 snp500_data() = readdlm(joinpath(Pkg.dir("DynamicApproximations"),"examples","data","snp500.txt"))
 

src/brute_force_search.jl

@@ -4,25 +4,32 @@ index sets with m segemets. The costs are evaluated using least squares.
 
 (Intended for verifying the dynamic programming algorithms)
 """
-function brute_force_search(ℓ, V_0N::QuadraticPolynomial, m::Integer)
+function brute_force_search(ℓ::AbstractTransitionCost{T}, V_N::QuadraticPolynomial{T}, m::Integer) where {T}
     cost_best = Inf
 
-    I_best = []
-    Y_best = []
+    I_best = Vector{Int64}(m+1)
+    Y_best = Vector{T}(m+1)
 
     N = size(ℓ, 2)
 
-    for I=IterTools.subsets(2:N-1, m-1)
+    for I_inner=IterTools.subsets(2:N-1, m-1)
+        I = [1; I_inner; N]
+    # I = zeros(Int64, m+1)
+    # I[1] = 1
+    # I[end] = N
+    #
+    # for I_inner=IterTools.subsets(2:N-1, m-1)
+    #
+    #     I[2:m] = I_inner
 
-        I = [1; I; N]
         P = zeros(m+1, m+1)
         q = zeros(m+1)
         r = 0
 
         # Add cost at right endpoint
-        P[end,end] = V_0N.a
-        q[end]     = V_0N.b
-        r          = V_0N.c
+        P[end,end] = V_N.a
+        q[end]     = V_N.b
+        r          = V_N.c
 
         # Form quadratic cost function Y'*P*Y + q'*Y + r
         # corresponding to the y-values in the vector Y
@@ -38,8 +45,8 @@ function brute_force_search(ℓ, V_0N::QuadraticPolynomial, m::Integer)
 
         if cost < cost_best
             cost_best = cost
-            Y_best = Yopt
-            I_best = I
+            Y_best .= Yopt
+            I_best .= I
         end
     end
     return I_best, Y_best, cost_best

src/solve.jl

@@ -227,16 +227,16 @@ global counter1
 global counter2
 
 """
-    Λ = pwq_dp_constrained(ℓ::AbstractTransitionCost{T}, V_0N::QuadraticPolynomial{T}, M::Integer, upper_bound=Inf) where {T}
+    Λ = pwq_dp_constrained(ℓ::AbstractTransitionCost{T}, V_N::QuadraticPolynomial{T}, M::Integer, upper_bound=Inf) where {T}
 
-Given the transition costs `ℓ[i,j](y_i,y_j)` and the cost at the endpoint `V_0N(y_N)` find all solutions `f` with up to `M` segments for the problem
+Given the transition costs `ℓ[i,j](y_i,y_j)` and the cost at the endpoint `V_N(y_N)` find all solutions `f` with up to `M` segments for the problem
 
-V_i^m = minimize_f^M [ Σ_{k=1}^i { ℓ[k,k+1](f(k),f(k+1)) } + V_0N(f(N)) ]
+V_i^m = minimize_f^M [ Σ_{k=1}^i { ℓ[k,k+1](f(k),f(k+1)) } + V_N(f(N)) ]
 s.t          f(k) being continuous piecewise linear with `m` segements.
 
 i.e. `Λ[m,i]` contains the best (in `ℓ` cost) continuous piecewise linear function `f` with up to `M` segments over the interval `i` to `N`
 """
-function pwq_dp_constrained{T}(ℓ::AbstractTransitionCost{T}, V_0N::QuadraticPolynomial{T}, M::Integer, upper_bound=Inf)
+function pwq_dp_constrained{T}(ℓ::AbstractTransitionCost{T}, V_N::QuadraticPolynomial{T}, M::Integer, upper_bound=Inf)
     #global counter1
     #global counter2
     #counter1 = 0
@@ -249,7 +249,7 @@ function pwq_dp_constrained{T}(ℓ::AbstractTransitionCost{T}, V_0N::QuadraticPo
     Λ = Array{PiecewiseQuadratic{T}}(M, N)
 
     for i=1:N-1
-        p = minimize_wrt_x2(ℓ[i, N], V_0N)
+        p = minimize_wrt_x2(ℓ[i, N], V_N)
         p.time_index = N
         Λ[1, i] .= create_new_pwq(p)
     end
@@ -318,14 +318,14 @@ where ℓ are positive-definite quadratic forms.
 
 Actually it computes cost to go functions V_(i,m), with
 """
-function pwq_dp_regularized{T}(ℓ::AbstractTransitionCost{T}, V_0N::QuadraticPolynomial{T}, ζ::T)
+function pwq_dp_regularized{T}(ℓ::AbstractTransitionCost{T}, V_N::QuadraticPolynomial{T}, ζ::T)
     N = size(ℓ, 2)
 
     Λ = Vector{PiecewiseQuadratic{T}}(N)
 
-    V_0N = deepcopy(V_0N)
-    V_0N.time_index = -1
-    Λ[N] = create_new_pwq(V_0N)
+    V_N = deepcopy(V_N)
+    V_N.time_index = -1
+    Λ[N] = create_new_pwq(V_N)
 
     ρ = QuadraticPolynomial{T}()
 
@@ -396,21 +396,21 @@ function recover_optimal_index_set{T}(Λ::PiecewiseQuadratic{T}, first_index=1)
 end
 
 """
-    Y, f = find_optimal_y_values(ℓ, V_0N, I)
-Given transition costs `ℓ`, cost of right endpoint `V_0N`, and breakpoint indicies `I`
+    Y, f = find_optimal_y_values(ℓ, V_N, I)
+Given transition costs `ℓ`, cost of right endpoint `V_N`, and breakpoint indicies `I`
 the optimal y-values `Y` and the optimal cost `f` are computed.
 """
 
-function find_optimal_y_values(ℓ, V_0N::QuadraticPolynomial, I)
+function find_optimal_y_values(ℓ, V_N::QuadraticPolynomial, I)
     m = length(I) - 1
 
     P = zeros(m+1, m+1)
     q = zeros(m+1)
 
     # Add cost for the right endpoint
-    P[m+1, m+1] = V_0N.a
-    q[m+1] = V_0N.b
-    r = V_0N.c
+    P[m+1, m+1] = V_N.a
+    q[m+1] = V_N.b
+    r = V_N.c
 
     # Form quadratic cost function Y'*P*Y + q'*Y + r
     # corresponding to the y-values in the vector Y
@@ -429,16 +429,22 @@ end
 
 # TODO: Include mζ in the cost?!
 """
-    I, Y, f = recover_solution(Λ::PiecewiseQuadratic{T}, ℓ, V_0N::QuadraticPolynomial, first_index=1)
+    I, Y, f = recover_solution(Λ::PiecewiseQuadratic{T}, ℓ, V_N::QuadraticPolynomial, first_index=1)
 
 """
-function recover_solution(Λ::PiecewiseQuadratic, ℓ, V_0N::QuadraticPolynomial, ζ=0.0)
+function recover_solution(Λ::PiecewiseQuadratic, ℓ, V_N::QuadraticPolynomial, ζ=0.0)
     I, _, f_expected = recover_optimal_index_set(Λ, 1)
-    Y, f = find_optimal_y_values(ℓ, V_0N::QuadraticPolynomial, I)
+    Y, f = find_optimal_y_values(ℓ, V_N::QuadraticPolynomial, I)
 
     f_regularized = f + ζ*(length(I)-1) # Include regularization cost
 
-    !isapprox(f_regularized, f_expected, atol=1e-10) && warn("Recovered cost ($f_regularized) is not what was expected from value function ($f_expected). Solution might be incorrect.")
+    if !isapprox(f_regularized, f_expected, atol=1e-10)
+        warn("Recovered cost ($f_regularized) is not what was expected from value function ($f_expected). Solution might be incorrect.")
+    end
+
+    if f_regularized < 0
+        warn("Computed cost < 0($f_regularized), if ≈ 0, this is probably due to numerical errors and nothing to worry about.")
+    end
 
     return I, Y, f
 end
@@ -553,32 +559,32 @@ end
 
 
 # Time-series
-function  fit_pwl_reguralized(g::AbstractArray, ζ; lazy=true)
+function  fit_pwl_regularized(g::AbstractArray, ζ; lazy=true)
     ℓ = lazy ? TransitionCostDiscrete{Float64}(g) :
                compute_discrete_transition_costs(g)
     # Discrete case, so cost at endpoint is quadratic
     cost_last = QuadraticPolynomial(1.0, -2*g[end], g[end]^2)
-    fit_pwl_reguralized_internal(ℓ, cost_last, ζ)
+    fit_pwl_regularized_internal(ℓ, cost_last, ζ)
 end
 
 # Continuous function
-function  fit_pwl_reguralized(g, t, ζ, tol=1e-3; lazy=true)
+function  fit_pwl_regularized(g, t, ζ, tol=1e-3; lazy=true)
     ℓ = lazy ? TransitionCostContinuous{Float64}(g, t, tol) :
                compute_transition_costs(g, t, tol)
     # Continouous case, no cost at endpoint
     cost_last = zero(QuadraticPolynomial{Float64})
-    fit_pwl_reguralized_internal(ℓ, cost_last, ζ)
+    fit_pwl_regularized_internal(ℓ, cost_last, ζ)
 end
 
-function  fit_pwl_reguralized_internal(ℓ, cost_last, ζ)
+function  fit_pwl_regularized_internal(ℓ, cost_last, ζ)
     Λ_reg = pwq_dp_regularized(ℓ, cost_last, ζ)
     #Get solution that starts at first index
-    I, _, f_reg = recover_optimal_index_set(Λ_reg[1])
-    Y, f = find_optimal_y_values(ℓ, cost_last, I)
+    I, Y, f = recover_solution(Λ_reg[1], ℓ, cost_last, ζ)
     return I, Y, f
 end
 
 
+
 # Time-series
 function  fit_pwl_constrained(g::AbstractArray, M; lazy=false)
     ℓ = lazy ? TransitionCostDiscrete{Float64}(g) :

src/types/QuadraticPolynomial.jl

@@ -12,6 +12,7 @@ has_been_used::Bool
 time_index::Int32
 ancestor::QuadraticPolynomial{T}
 function QuadraticPolynomial{T}(a::T, b::T, c::T) where {T}
+    # @assert a ≥ 0, # Δ may have negative a ..., Δ is currently not used...
     new(a, b, c, false, -1)
 end
 

test/examples.jl

@@ -20,9 +20,9 @@ M = 10
                   7.258287507540115,7.316043144146352,7.215020359632374]
     @test fvec[M] ≈ 0.11431141407229006
 
-    #Test reguralize discrete
+    #Test regularize discrete
     # should generate solution with 10 segments
-    I, Y, f = fit_pwl_reguralized(data, 0.02, lazy=lazy)
+    I, Y, f = fit_pwl_regularized(data, 0.02, lazy=lazy)
     @test length(I) == 9
     @test I == Ivec[8]
     @test Y ≈ Yvec[8]     rtol=sqrt(eps())
@@ -34,7 +34,7 @@ M = 10
 
 
     ζ = 0.1
-    I, Y, cost = fit_pwl_reguralized(g_, t, ζ, lazy=lazy)
+    I, Y, cost = fit_pwl_regularized(g_, t, ζ, lazy=lazy)
 
     @test I == [1, 87, 107, 129, 147, 201]
     @test cost ≈ 0.35055983342102515
@@ -46,7 +46,7 @@ M = 10
     @test cost ≈ fvec[5]       rtol=sqrt(eps())
 
     ζ = 0.002
-    I, Y, cost = fit_pwl_reguralized(g_, t, ζ, lazy=lazy)
+    I, Y, cost = fit_pwl_regularized(g_, t, ζ, lazy=lazy)
 
     @test I == [1, 9, 15, 33, 52, 68, 87, 104, 111, 125, 132, 146, 153, 170, 188, 201]
     @test cost ≈ 0.006061712727833957

test/problems/exp_problem.jl

@@ -4,7 +4,7 @@ function exp_problem()
 
     M_bf = 10
 
-    ζvec = [100, 30, 10, 3, 1, 0.3, 0.1, 0.03, 0.01]
+    ζ_vec = [100, 30, 10, 3, 1, 0.3, 0.1, 0.03, 0.01]
 
     I_sols =
     [[1, 11],
@@ -30,5 +30,5 @@ function exp_problem()
     0.02951761758059246
     -8.098410830825742e-12]
 
-    return g, M_bf, ζvec, I_sols, f_sols
+    return g, ζ_vec, I_sols, f_sols
 end

test/problems/snp500_problem.jl

@@ -6,7 +6,7 @@ function snp500_problem()
 
     M_bf = 2
 
-    ζ = 0.005:0.005:0.05
+    ζ_vec = 0.005:0.005:0.05
 
     I_sols =
     [[1, 300],
@@ -40,5 +40,5 @@ function snp500_problem()
     0.06051435296103591,
     0.05690190505811188]
 
-    return g, M_bf, ζ, I_sols, f_sols
+    return g, ζ_vec, I_sols, f_sols
 end

test/problems/square_wave_problem.jl

@@ -4,7 +4,7 @@ function square_wave_problem()
 
     M_bf = 5
 
-    ζvec = [100, 30, 10, 3, 1, 0.3, 0.1, 0.03, 0.01]
+    ζ_vec = [100, 30, 10, 3, 1, 0.3, 0.1, 0.03, 0.01]
 
     # For each m there are many non-unique solution
     I_sols =
@@ -21,5 +21,5 @@ function square_wave_problem()
         2.1356521739130434];
         zeros(Float64, 16)]
 
-    return g, M_bf, ζvec, I_sols, f_sols
+    return g, ζ_vec, I_sols, f_sols
 end

test/problems/straight_line_problem.jl

@@ -4,7 +4,7 @@ function straight_line_problem()
 
     M_bf = 2
 
-    ζvec = [100, 30, 10, 3, 1, 0.3, 0.1, 0.03, 0.01]
+    ζ_vec = [100, 30, 10, 3, 1, 0.3, 0.1, 0.03, 0.01]
 
     # For each m there are many non-unique solution
     I_sols =
@@ -13,5 +13,5 @@ function straight_line_problem()
 
     f_sols = zeros(Float64, 10)
 
-    return g, M_bf, ζvec, I_sols, f_sols
+    return g, ζ_vec, I_sols, f_sols
 end

test/test_discrete_fit.jl

@@ -1,17 +1,12 @@
 using Base.Test
 using DynamicApproximations: find_optimal_y_values
 
-include(joinpath(Pkg.dir("DynamicApproximations"),"test","problems", "snp500_problem.jl"))
-include(joinpath(Pkg.dir("DynamicApproximations"),"test","problems", "exp_problem.jl"))
-include(joinpath(Pkg.dir("DynamicApproximations"),"test","problems", "square_wave_problem.jl"))
-include(joinpath(Pkg.dir("DynamicApproximations"),"test","problems", "straight_line_problem.jl"))
 
-##
-
-for problem_fcn in [straight_line_problem, square_wave_problem, exp_problem, snp500_problem]
+for problem_fcn in ["straight_line_problem", "square_wave_problem", "exp_problem", "snp500_problem"]
 
+    include(joinpath(Pkg.dir("DynamicApproximations"),"test","problems", problem_fcn * ".jl"))
 ##
-    g, M_bf, ζvec, I_sols, f_sols = problem_fcn()
+    g, ζvec, I_sols, f_sols = @eval $(Symbol(problem_fcn))()
 
     M = length(I_sols)
 
@@ -28,25 +23,6 @@ for problem_fcn in [straight_line_problem, square_wave_problem, exp_problem, snp
         @test f ≈ f_sols[m]   rtol=1e-10 atol=1e-10
     end
 
-
-    ##
-
-    @testset "Data set: $problem_fcn, bruteforce with m=$m" for m in 1:M_bf
-
-        I_bf, Y_bf, f_bf = brute_force_search(ℓ, cost_last, m)
-
-        @test isempty(I_sols[m]) || I_bf == I_sols[m]
-        @test f_bf ≈ f_sols[m]  rtol=1e-10 atol=1e-10
-
-        if !isempty(I_sols[m])
-            Y_recovered, f_recovered = find_optimal_y_values(ℓ, cost_last, I_sols[m])
-            @test Y_bf ≈ Y_recovered    rtol=1e-10
-        end
-    end
-
-    ##
-
-    #ζvec = [100, 30, 10, 3, 1, 0.3, 0.1, 0.03, 0.01]
     @testset "Data set: $problem_fcn, regularization with ζ=$ζ" for ζ in ζvec
 
         Λ_reg = pwq_dp_regularized(ℓ, cost_last, ζ)
@@ -61,8 +37,7 @@ for problem_fcn in [straight_line_problem, square_wave_problem, exp_problem, snp
         @test m_expected == length(I) - 1
         @test isempty(I_sols[m_expected]) || I == I_sols[m_expected]
         @test f_reg ≈ f_sols[m_expected] + ζ*m_expected
-
-        #println(problem_fcn, "  " , length(I))
+
     end
 
 end