Merge pull request #6 from mfalt/code-improvement

Code improvement

README.md

@@ -4,6 +4,12 @@
 
 [![codecov](https://codecov.io/gh/mfalt/DynamicApproximations.jl/branch/master/graph/badge.svg?token=nt4j2gNB2l)](https://codecov.io/gh/mfalt/DynamicApproximations.jl)
 
+
+This package solves the problem of piecewise linear, *continuous*, approximation subject to either a hard limit or a regularization penalty on the number of break points. An exact solution is obtained using dynamic program over piecewise quadratic function, which avoids the combinatorial complexity of a naive approach.
+
+(FIXME: Mathematical formulation)
+
+
 ### Example: Constrained approximation
 
 Find best continouous piecewise approximations with up to M segments

examples/example2.jl

@@ -1,33 +1,20 @@
 using IterTools
 using Plots
-
-# The problem Σerror^2 + card(I)
-# Heuristic
-
-include(joinpath(Pkg.dir("DynamicApproximations"),"src","DynamicApproximations.jl"))
-##
-DA = DynamicApproximations
+using DynamicApproximations
 
 data = readdlm(joinpath(Pkg.dir("DynamicApproximations"),"examples","data","snp500.txt"))
 data = data[1:300]
 
 N = length(data)
 
-
-
 @time ℓ = DA.compute_discrete_transition_costs(data);
-
-#@time I, Y, f = brute_force_search(ℓ, K-1);
-
 cost_last = DA.QuadraticPolynomial(1.0, -2*data[end], data[end]^2)
 
 start_time = time()
 gc()
-@time Λ = DA.pwq_dp_constrained(ℓ, cost_last, 10, 1.65);
+@time Λ = pwq_dp_constrained(ℓ, cost_last, 10, 1.65);
 println("Time: ", time()-start_time)
 
-#println(counter1, " ", counter2)
-
 for k=3:10
     println(k, " : ", DA.recover_optimal_index_set(Λ[k, 1])[3])
 end

src/brute_force_search.jl

@@ -6,24 +6,21 @@ index sets with m segemets. The costs are evaluated using least squares.
 """
 function brute_force_search(ℓ::AbstractTransitionCost{T}, V_N::QuadraticPolynomial{T}, m::Integer) where {T}
     cost_best = Inf
-
     I_best = Vector{Int64}(m+1)
     Y_best = Vector{T}(m+1)
 
     N = size(ℓ, 2)
 
+    I = zeros(Int64, m+1)
+    I[1] = 1
+    I[end] = N
+
     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
-
-        P = zeros(m+1, m+1)
-        q = zeros(m+1)
+
+        I[2:m] .= I_inner
+
+        P = SymTridiagonal(zeros(T, m+1), zeros(T, m+1))
+        q = zeros(T, m+1)
         r = 0
 
         # Add cost at right endpoint
@@ -34,14 +31,16 @@ function brute_force_search(ℓ::AbstractTransitionCost{T}, V_N::QuadraticPolyno
         # Form quadratic cost function Y'*P*Y + q'*Y + r
         # corresponding to the y-values in the vector Y
         for j=1:m
-            P[j:j+1,j:j+1] .+= ℓ[I[j], I[j+1]].P
+            P.dv[j] += ℓ[I[j], I[j+1]].P[1,1]
+            P.dv[j+1] += ℓ[I[j], I[j+1]].P[2,2]
+            P.ev[j] += ℓ[I[j], I[j+1]].P[1,2]
             q[j:j+1] .+= ℓ[I[j], I[j+1]].q
             r += ℓ[I[j], I[j+1]].r
         end
 
         # find the optimal Y-vector, and compute the correspinding error
         Yopt = -(P \ q) / 2
-        cost = Yopt' * P * Yopt + q' * Yopt + r
+        cost = dot(Yopt, P*Yopt) + q' * Yopt + r
 
         if cost < cost_best
             cost_best = cost

src/solve.jl

@@ -24,9 +24,9 @@ function add_quadratic!{T}(Λ::PiecewiseQuadratic{T}, ρ::QuadraticPolynomial{T}
         right_endpoint = get_right_endpoint(λ_curr)
 
 
-        Δa = ρ.a - λ_curr.p.a
-        Δb = ρ.b - λ_curr.p.b
-        Δc = ρ.c - λ_curr.p.c
+        Δa = ρ.a - λ_curr.π.a
+        Δb = ρ.b - λ_curr.π.b
+        Δc = ρ.c - λ_curr.π.c
 
         b2_minus_4ac =  Δb^2 - 4*Δa*Δc
 
@@ -142,11 +142,11 @@ end
 
 @inline function update_segment_new(λ_prev, λ_curr, ρ)
     ρ.has_been_used = true # TODO: Could be moved to the second clause
-    if λ_prev.p === ρ
+    if λ_prev.π === ρ
         λ_prev.next = λ_curr.next
         v1, v2 = λ_prev, λ_curr.next
     else
-        λ_curr.p = ρ
+        λ_curr.π = ρ
         v1, v2 = λ_curr, λ_curr.next
     end
     return v1, v2
@@ -162,7 +162,7 @@ end
 
 @inline function update_segment_new_old(λ_prev, λ_curr, ρ, root)
     # ... λ_prev | (new segment) | λ_curr ...
-    if λ_prev.p === ρ
+    if λ_prev.π === ρ
         λ_curr.left_endpoint = root
     else
         ρ.has_been_used = true
@@ -185,7 +185,7 @@ end
     # The second new piece contains a copy of the polynomial in λ_curr
 
     ρ.has_been_used = true
-    λ_2 =  PiecewiseQuadratic(λ_curr.p, root2, λ_curr.next)
+    λ_2 =  PiecewiseQuadratic(λ_curr.π, root2, λ_curr.next)
     λ_1 =  PiecewiseQuadratic(ρ, root1, λ_2)
     λ_curr.next = λ_1
     return λ_2, λ_2.next
@@ -268,11 +268,8 @@ function pwq_dp_constrained{T}(ℓ::AbstractTransitionCost{T}, V_N::QuadraticPol
             for ip=i+1:N-m+1
                 DEBUG && println("(m:$m, i:$i, ip:$ip)")
                 for λ in Λ[m-1, ip]
-                    p = λ.p
 
-                    #counter1 +=
-
-                    minimize_wrt_x2(ℓ[i,ip], p, ρ)
+                    minimize_wrt_x2(ℓ[i,ip], λ.π, ρ)
 
                     DEBUG && println("Obtained ρ = $ρ")
 
@@ -288,7 +285,7 @@ function pwq_dp_constrained{T}(ℓ::AbstractTransitionCost{T}, V_N::QuadraticPol
 
                     if ρ.has_been_used == true
                         ρ.time_index = ip
-                        ρ.ancestor = p
+                        ρ.ancestor = λ.π
                         ρ = QuadraticPolynomial{T}()
                         ρ.has_been_used = false
                     end
@@ -336,18 +333,17 @@ function pwq_dp_regularized{T}(ℓ::AbstractTransitionCost{T}, V_N::QuadraticPol
 
             ζ_level_insertion = false
             for λ in Λ[ip]
-                p = λ.p
                 #counter1 += 1
 
-                minimize_wrt_x2(ℓ[i,ip], p, ρ)
+                minimize_wrt_x2(ℓ[i,ip], λ.π, ρ)
                 ρ.c += ζ # add cost for break point
 
 
                 add_quadratic!(Λ_new, ρ)
 
                 if ρ.has_been_used
                     ρ.time_index = ip
-                    ρ.ancestor = p
+                    ρ.ancestor = λ.π
                     ρ = QuadraticPolynomial{T}()
                     ρ.has_been_used = false
 
@@ -469,9 +465,9 @@ function poly_minus_constant_is_greater{T}(Λ::PiecewiseQuadratic{T}, ρ::Quadra
         left_endpoint = λ_curr.left_endpoint
         right_endpoint = get_right_endpoint(λ_curr)
 
-        Δa = ρ.a - λ_curr.p.a
-        Δb = ρ.b - λ_curr.p.b
-        Δc = (ρ.c - ζ) - λ_curr.p.c
+        Δa = ρ.a - λ_curr.π.a
+        Δb = ρ.b - λ_curr.π.b
+        Δc = (ρ.c - ζ) - λ_curr.π.c
 
         b2_minus_4ac =  Δb^2 - 4*Δa*Δc
 

src/types/PiecewiseQuadratic.jl

@@ -11,9 +11,10 @@ Note that the first element of the list is sometimes assumed to be a list head w
 The list/element `x` is the last element in the list when `x.next === x` or `x.next.left_endpoint==Inf`.
 """
 mutable struct PiecewiseQuadratic{T}
-    p::QuadraticPolynomial{T}
+    π::QuadraticPolynomial{T}
     left_endpoint::Float64
     next::PiecewiseQuadratic{T}
+
     PiecewiseQuadratic{T}()  where T = (x=new(); x.left_endpoint=Inf; x.next = x)
     PiecewiseQuadratic{T}(p, left_endpoint, next) where T = new(p, left_endpoint, next)
 end
@@ -126,7 +127,7 @@ function show{T}(io::IO, Λ::PiecewiseQuadratic{T})
         end
 
         print(io, "\t  :   ")
-        show(io, λ.p)
+        show(io, λ.π)
         print(io, "\n")
     end
     return
@@ -135,7 +136,7 @@ end
 function (Λ::PiecewiseQuadratic)(x::Number)
     for λ in Λ
         if λ.left_endpoint <= x < get_right_endpoint(λ)
-            return λ.p(x)
+            return λ.π(x)
         end
     end
     return NaN
@@ -145,7 +146,7 @@ function (Λ::PiecewiseQuadratic)(x::AbstractArray)
     y = zeros(x)
     for λ in Λ
         inds = λ.left_endpoint .<= x .< get_right_endpoint(λ)
-        y[inds] .= λ.p.(x[inds])
+        y[inds] .= λ.π.(x[inds])
     end
     return y
 end
@@ -184,7 +185,7 @@ function get_vals(Λ::PiecewiseQuadratic)
         end
 
         y_grid = linspace(left_endpoint, right_endpoint)
-        vals =  λ.p.(y_grid)
+        vals =  λ.π.(y_grid)
         append!(x, y_grid)
         append!(y, vals)
         push!(x_all, y_grid)
@@ -216,25 +217,25 @@ function find_minimum(Λ::PiecewiseQuadratic)
     # TODO: True?
 
     f_opt = Inf
-    p_opt = Λ.p # The segment containing the smallest polynimal
+    π_opt = typeof(Λ.π)()
     x_opt = NaN
 
     for λ in Λ
-        x, f = find_minimum(λ.p)
+        x, f = find_minimum(λ.π)
         if f < f_opt
             f_opt = f
             x_opt = x
-            p_opt = λ.p
+            π_opt = λ.π
         end
     end
-    return p_opt, x_opt, f_opt
+    return π_opt, x_opt, f_opt
 end
 
 function find_minimum_value(Λ::PiecewiseQuadratic)
     f_opt = Inf
 
     for λ in Λ
-        _, f = find_minimum(λ.p)
+        _, f = find_minimum(λ.π)
         if f < f_opt
             f_opt = f
         end

src/types/QuadraticPolynomial.jl

@@ -37,22 +37,13 @@ Base.zero{T<:QuadraticPolynomial}(S::T) = zero(T)
 
 -(p1::QuadraticPolynomial, p2::QuadraticPolynomial) = QuadraticPolynomial(p1.a-p2.a, p1.b-p2.b, p1.c-p2.c)
 
-# # x .= y .- z
-# # x .= y .+ z
-# function Base.broadcast!{T}(op::Function, x::QuadraticPolynomial{T}, y::QuadraticPolynomial{T}, z::QuadraticPolynomial{T})
-#     x.a = op(y.a,z.a)
-#     x.b = op(y.b,z.b)
-#     x.c = op(y.c,z.c)
-# end
-
 ==(p1::QuadraticPolynomial, p2::QuadraticPolynomial) = (p1.a==p2.a) && (p1.b==p2.b) && (p1.c==p2.c)
 
 @inline function (p::QuadraticPolynomial)(x)
     return p.a*x^2 + p.b*x + p.c
 end
 
 
-
 # Finds the minimum of a positive definite quadratic one variable polynomial
 # the find_minimum fcn returns (opt_x, opt_val)
 function find_minimum(p::QuadraticPolynomial)

test/README.md

@@ -0,0 +1 @@
+Empty solution index sets indicates that the solution is non-unique. In these cases there is no test.

test/auxilliary_test_fcns.jl

@@ -0,0 +1,37 @@
+# Some help functions to generate interesting test cases
+circle_segment(N) = sin.(acos.(linspace(-1, 1, N)))
+linear_trend(N) = (0:N-1)
+
+# Finds the solution to the cardinality constrained problem for
+# m=1:M segments using brute force
+function brute_force_multi(g, t, M; tol=1e-3)
+
+    if length(g) <= M
+        warn("Specified M too large. Setting M=length(g)-1.")
+        M = length(g)-1
+    end
+
+    ℓ = []
+    V_N = []
+
+    if typeof(g) <: AbstractArray
+        ℓ = compute_discrete_transition_costs(g)
+        V_N = QuadraticPolynomial(1.0, -2*g[end], g[end]^2)
+    else
+        ℓ = TransitionCostContinuous{Float64}(g, t, tol=tol)
+        V_N = zero(QuadraticPolynomial{Float64})
+    end
+
+    I_vec = Vector{Vector{Int64}}(M)
+    Y_vec = Vector{Vector{Float64}}(M)
+    f_vec = Vector{Float64}(M)
+
+    for m=1:M
+        (I_bf, Y_bf, f_bf) = brute_force_search(ℓ, V_N, m)
+        I_vec[m] = I_bf
+        Y_vec[m] = Y_bf
+        f_vec[m] = f_bf
+    end
+
+    return I_vec, Y_vec, f_vec
+end

test/brute_force_search.jl

@@ -0,0 +1,27 @@
+using Base.Test
+using DynamicApproximations
+
+include("auxilliary_test_fcns.jl")
+
+# Consider a subset of tests that are relatively small
+for problem_fcn in ["discontinuous1",
+                    "discontinuous2",
+                    "super_exponential1",
+                    "super_exponential3",
+                    "white_noise",
+                    "exponential"]
+
+    include(joinpath(Pkg.dir("DynamicApproximations"),"test","problems", problem_fcn * ".jl"))
+
+    g, ζ_vec, I_sols, f_sols = @eval $(Symbol(problem_fcn))()
+
+    M = length(I_sols)
+    I_vec, _, f_vec = brute_force_multi(g, -1, M)
+
+    @testset "Data set: $problem_fcn, brute_force_search m=$m" for m in 1:length(f_sols)
+        @test f_vec[m] ≈ f_sols[m]  atol=1e-10
+        if !isempty(I_sols[m])
+            @test I_vec[m] == I_sols[m]
+        end
+    end
+end