From c1b49bf7dfc2eee1b04c9bc04208bd39f520a90f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mattias=20F=C3=A4lt?= Date: Thu, 31 Aug 2017 20:32:07 +0200 Subject: [PATCH] 19x times faster, 45 times less memory --- REQUIRE | 1 + src/dev.jl | 64 ++++++++++++++++++++------------ src/types/QuadraticForm2.jl | 2 +- src/types/QuadraticPolynomial.jl | 7 ++++ test/tst.jl | 30 ++++++++------- 5 files changed, 65 insertions(+), 39 deletions(-) diff --git a/REQUIRE b/REQUIRE index 406ecaf..4f61adf 100644 --- a/REQUIRE +++ b/REQUIRE @@ -1,2 +1,3 @@ julia 0.6 StaticArrays +QuadGK diff --git a/src/dev.jl b/src/dev.jl index bcb3735..8bddac8 100644 --- a/src/dev.jl +++ b/src/dev.jl @@ -63,6 +63,7 @@ function add_quadratic{T}(Λ::PiecewiseQuadratic{T}, ρ::QuadraticPolynomial{T}) λ_prev = Λ λ_curr = Λ.next + Δ = QuadraticPolynomial(0., 0., 0.) while ~isnull(λ_curr) left_endpoint = get(λ_curr).left_endpoint @@ -70,17 +71,17 @@ function add_quadratic{T}(Λ::PiecewiseQuadratic{T}, ρ::QuadraticPolynomial{T}) # TODO: This is probably not needed now.. if left_endpoint == -1e9 #println("minus-inf") - left_endpoint = -10000 + left_endpoint = -10000.0 end right_endpoint = get_right_endpoint(get(λ_curr)) if right_endpoint == 1e9 #println("inf") - right_endpoint = left_endpoint + 20000 + right_endpoint = left_endpoint + 20000.0 end - Δ = ρ - get(λ_curr).p + Δ .= ρ .- get(λ_curr).p root1,root2 = roots(Δ) @@ -108,7 +109,6 @@ function add_quadratic{T}(Λ::PiecewiseQuadratic{T}, ρ::QuadraticPolynomial{T}) end #println("case 3:") - λ_prev, λ_curr = impose_quadratic_to_endpoint(λ_prev, get(λ_curr), (get(λ_curr).left_endpoint + right_endpoint)/2, ρ, Δ) #println(Λ) @@ -123,13 +123,12 @@ end # TODO Possibly just insert all roots, and then do the comparisons? # Leads to some unnecessary nodes being created but is perhaps simplier? - -function impose_quadratic_to_root(λ_prev, λ_curr, root, midpoint, ρ, Δ) +function impose_quadratic_to_root{T}(λ_prev::PiecewiseQuadratic{T}, λ_curr::PiecewiseQuadratic{T}, root::T, midpoint::T, ρ::QuadraticPolynomial{T}, Δ::QuadraticPolynomial{T}) if Δ((λ_curr.left_endpoint + root)/2) < 0 # ρ is smallest, i.e., should be inserted if λ_prev.p === ρ #println("1") λ_curr.left_endpoint = root - return λ_prev, Nullable{PiecewiseQuadratic}(λ_curr) + return λ_prev, Nullable{PiecewiseQuadratic{Float64}}(λ_curr) else #println("2") @@ -148,28 +147,26 @@ function impose_quadratic_to_root(λ_prev, λ_curr, root, midpoint, ρ, Δ) end end -function impose_quadratic_to_endpoint(λ_prev, λ_curr, midpoint, ρ, Δ) - if Δ(midpoint) < 0 # ρ is smallest, i.e., should be inserted +@inline function impose_quadratic_to_endpoint{T}(λ_prev::PiecewiseQuadratic{T}, λ_curr::PiecewiseQuadratic{T}, midpoint, ρ, Δ) + local v1, v2 + if Δ(midpoint) < 0 # ρ is smallest, i.e., should be inserted if λ_prev.p === ρ #println("1") λ_new = delete_next(λ_prev) - return λ_prev, λ_new + v1, v2 = λ_prev, λ_new else #println("2") λ_curr.p = ρ - return λ_curr, λ_curr.next + v1, v2 = λ_curr, λ_curr.next end else # Do nothing #println("3") - return λ_curr, λ_curr.next + v1, v2 = λ_curr, λ_curr.next end + return v1, v2 #Single point of exit end - - - - function minimize_wrt_x2(qf::QuadraticForm2) P = qf.P q = qf.q @@ -183,7 +180,24 @@ function minimize_wrt_x2(qf::QuadraticForm2) QuadraticPolynomial(P[1,1], q[1], r) else # There are some special cases, but disregards these - QuadraticPolynomial(0,0,-Inf) + QuadraticPolynomial(0.,0.,-Inf) + end +end + +function minimize_wrt_x2_fast{T}(qf::QuadraticForm2{T},a,b,c) + P = qf.P + q = qf.q + r = qf.r + + if P[2,2] > 0 + QuadraticPolynomial(P[1,1] - P[1,2]^2 / (P[2,2]+a), + q[1] - P[1,2]*(q[2]+b) / (P[2,2]+a), + (r+c) - (q[2]+b)^2 / (P[2,2]+a)/ 4) + elseif (P[2,2]+a) == 0 || P[1,2] == 0 || (q[2]+b) == 0 + QuadraticPolynomial(P[1,1], q[1], r+c) + else + # There are some special cases, but disregards these + QuadraticPolynomial(0.,0.,-Inf) end end @@ -191,11 +205,11 @@ end """ Find optimal fit """ -function find_optimal_fit(Λ_0, ℓ, M) +function find_optimal_fit{T}(Λ_0::Array{PiecewiseQuadratic{T},1}, ℓ::Array{QuadraticForm2{T},2}, M) N = size(ℓ, 2) - Λ = Array{PiecewiseQuadratic}(M, N-1) + Λ = Array{PiecewiseQuadratic{T}}(M, N-1) Λ[1, :] .= Λ_0 @@ -207,8 +221,10 @@ function find_optimal_fit(Λ_0, ℓ, M) for λ in Λ[m-1, ip] p = λ.p - ρ = dev.minimize_wrt_x2( - ℓ[i,ip] + dev.QuadraticForm2(@SMatrix([0 0; 0 1])*p.a, @SVector([0, 1])*p.b, p.c)) + # ρ = dev.minimize_wrt_x2( + # ℓ[i,ip] + dev.QuadraticForm2{T}(@SMatrix([0. 0; 0 1])*p.a, @SVector([0., 1])*p.b, p.c)) + # # Avoid ceting two extra QuadraticForm2 + ρ = dev.minimize_wrt_x2_fast(ℓ[i,ip], p.a, p.b, p.c) ρ.time_index = ip ρ.ancestor = p @@ -257,7 +273,7 @@ function find_minimum(Λ::PiecewiseQuadratic) end -function recover_solution(Λ::PiecewiseQuadratic, first_index=1, last_index=-1) +function recover_solution{T}(Λ::PiecewiseQuadratic{T}, first_index=1, last_index=-1) p, y, f = find_minimum(Λ) @@ -288,7 +304,7 @@ Computes the transition costs ℓ given a polynomial and a sequence t """ function compute_transition_costs(g, t::AbstractArray) - + T = Float64 # Find primitive functions to g, t*g, and g^2 # and evaluate them at the break points #I_g = polyint(g).(t) @@ -309,7 +325,7 @@ function compute_transition_costs(g, t::AbstractArray) end - ℓ = Array{QuadraticForm2}(N-1,N) + ℓ = Array{QuadraticForm2{T}}(N-1,N) for i=1:N-1 for ip=i+1:N diff --git a/src/types/QuadraticForm2.jl b/src/types/QuadraticForm2.jl index 29ac94d..439c123 100644 --- a/src/types/QuadraticForm2.jl +++ b/src/types/QuadraticForm2.jl @@ -1,6 +1,6 @@ # Quadratic form of 2 variables, used for representing the transition costs ℓ type QuadraticForm2{T} - P::SMatrix{2,2,T} + P::SMatrix{2,2,T,4} q::SVector{2,T} r::T end diff --git a/src/types/QuadraticPolynomial.jl b/src/types/QuadraticPolynomial.jl index bad6d30..a0d3573 100644 --- a/src/types/QuadraticPolynomial.jl +++ b/src/types/QuadraticPolynomial.jl @@ -33,6 +33,13 @@ end -(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) diff --git a/test/tst.jl b/test/tst.jl index 5931824..a81fd58 100644 --- a/test/tst.jl +++ b/test/tst.jl @@ -5,7 +5,7 @@ using Base.Test -include("../src/dev.jl") +include(joinpath(Pkg.dir("DynamicApproximations"),"src","dev.jl")) @@ -22,18 +22,18 @@ t = linspace(0,2π,N) #g = Poly([0,0,1,-1]) g = sin -@time ℓ = dev.compute_transition_costs(g, t) +@time ℓ = dev.compute_transition_costs(g, t); K = 4 -@time I, Y, f = dev.brute_force_optimization(ℓ, K-1) +@time I, Y, f = dev.brute_force_optimization(ℓ, K-1); -Λ_0 = [dev.create_new_pwq(dev.minimize_wrt_x2(ℓ[i, N])) for i in 1:N-1] +Λ_0 = [dev.create_new_pwq(dev.minimize_wrt_x2(ℓ[i, N])) for i in 1:N-1]; -@time Λ = dev.find_optimal_fit(Λ_0, ℓ, 7) +@time Λ = dev.find_optimal_fit(Λ_0, ℓ, 7); -@time I2, y2, f2 = dev.recover_solution(Λ[3, 1], 1, N) +@time I2, y2, f2 = dev.recover_solution(Λ[K, 1], 1, N) Y2, _ = dev.find_optimal_y_values(ℓ, I2) println("Comparison: ", sqrt(f), " --- ", sqrt(f2)) @@ -55,11 +55,13 @@ end #x = linspace(0,1,100) #closefig() - -using PyPlot -#close() -figure(1) -plot(t, g.(t), "g", t[I2], Y2, "bo-") - -figure(2) -plot(1:length(t)-1, l') +using Plots +plot(t, g.(t), lab="g(t)") +plot!(t[I2],Y2, m=:circle, lab="approximation") +#using PyPlot +# #close() +# figure(1) +# plot(t, g.(t), "g", t[I2], Y2, "bo-") +# +# figure(2) +# plot(1:length(t)-1, l')