From 7215bc38d84b1bc45b741a084ad28dee9c3a80a1 Mon Sep 17 00:00:00 2001 From: oloft Date: Thu, 7 Sep 2017 14:30:42 +0200 Subject: [PATCH] Exploiting breaking, fixed allocation for p.ancestor. using Da>0 for abs(Da). --- src/dev.jl | 153 +++++++++++-------------------- src/types/PiecewiseQuadratic.jl | 30 ++++-- src/types/QuadraticPolynomial.jl | 25 ++++- test/find_optimal_fit.jl | 5 +- test/piecewise_quadratics.jl | 21 +++-- 5 files changed, 112 insertions(+), 122 deletions(-) diff --git a/src/dev.jl b/src/dev.jl index 69a0b3c..f694016 100644 --- a/src/dev.jl +++ b/src/dev.jl @@ -26,28 +26,7 @@ include("types/PiecewiseQuadratic.jl") include("types/QuadraticForm2.jl") global const DEBUG = false - - -function roots{T}(p::QuadraticPolynomial{T}) - b2_minus_4ac = p.b^2 - 4*p.a*p.c - if b2_minus_4ac < 0 - return NaN, NaN - end - - if p.a != 0 - term1 = (p.b / 2 / p.a) - term2 = sqrt(b2_minus_4ac) / 2 / abs(p.a) - return -term1-term2, -term1+term2 - else - term = -p.c / p.b - return term, Inf - end -end - - - - - +global const COUNTER_TEST = false # TODO come up with better symbol for ρ """ @@ -194,33 +173,20 @@ function add_quadratic2{T}(Λ::PiecewiseQuadratic{T}, ρ::QuadraticPolynomial{T} DEBUG && println("Inserting: ", ρ) if Λ.next.left_endpoint == Inf # I.e. the piecewise quadratic object is empty, perhaps better to add dummy polynomial - insert(Λ, ρ, -1e9) + insert(Λ, ρ, -Inf) return end λ_prev = Λ λ_curr = Λ.next - #left_endpoint = NaN - while λ_curr.left_endpoint != Inf #??? #global counter2 += 1 DEBUG && println(Λ) left_endpoint = λ_curr.left_endpoint - - # TODO: This is probably not needed now.. - if left_endpoint == -1e9 - #println("minus-inf") - left_endpoint = -10000.0 - end - right_endpoint = get_right_endpoint(λ_curr) - if right_endpoint == 1e9 - #println("inf") - right_endpoint = left_endpoint + 20000.0 - end Δa = ρ.a - λ_curr.p.a Δb = ρ.b - λ_curr.p.b @@ -231,23 +197,25 @@ function add_quadratic2{T}(Λ::PiecewiseQuadratic{T}, ρ::QuadraticPolynomial{T} if Δa > 0 # ρ has greater curvature, i.e., ρ is smallest in the middle if b2_minus_4ac <= 0 # No intersections, old quadratic is smallest, just step forward - λ_prev = λ_curr - λ_curr = λ_prev.next + DEBUG && println("No intersections, old quadratic is smallest, Δa > 0, breaking.") + break else # Compute the intersections term1 = -(Δb / 2 / Δa) - term2 = sqrt(b2_minus_4ac) / 2 / abs(Δa) + term2 = sqrt(b2_minus_4ac) / 2 / Δa # Δa > 0 root1, root2 = term1-term2, term1+term2 DEBUG && println("Δa > 0 root1:", root1, " root2:", root2) # Check where the intersections are and act accordingly - if root1 >= right_endpoint || root2 <= left_endpoint + if root1 >= right_endpoint + DEBUG && println("Two intersections to the right") + λ_prev, λ_curr = update_segment_do_nothing(λ_curr) + elseif root2 <= left_endpoint # No intersections, old quadratic is smallest, step forward - DEBUG && println("Two intersections to the side") - λ_prev = λ_curr - λ_curr = λ_prev.next + DEBUG && println("Two intersections to the left") + break # There will be no more intersections since Δa > 0 elseif root1 <= left_endpoint && root2 >= right_endpoint # No intersections, new quadratic is smallest DEBUG && println("One intersections on either side") @@ -255,12 +223,14 @@ function add_quadratic2{T}(Λ::PiecewiseQuadratic{T}, ρ::QuadraticPolynomial{T} elseif root1 > left_endpoint && root2 < right_endpoint DEBUG && println("Two intersections within the interval") λ_prev, λ_curr = update_segment_old_new_old(λ_curr, ρ, root1, root2) + break # There will be no more intersections since Δa > 0 elseif root1 > left_endpoint DEBUG && println("Root 1 within the interval") λ_prev, λ_curr = update_segment_old_new(λ_curr, ρ, root1) elseif root2 < right_endpoint DEBUG && println("Root 2 within the interval") λ_prev, λ_curr = update_segment_new_old(λ_prev, λ_curr, ρ, root2) + break # There will be no more intersections since Δa > 0 else error("Shouldn't end up here") end @@ -272,8 +242,8 @@ function add_quadratic2{T}(Λ::PiecewiseQuadratic{T}, ρ::QuadraticPolynomial{T} else # Compute the intersections term1 = -(Δb / 2 / Δa) - term2 = sqrt(b2_minus_4ac) / 2 / abs(Δa) - root1, root2 = term1-term2, term1+term2 + term2 = sqrt(b2_minus_4ac) / 2 / Δa # < 0 + root1, root2 = term1+term2, term1-term2 DEBUG && println("Δa < 0 root1:", root1, " root2:", root2) # Check where the intersections are and act accordingly @@ -282,8 +252,7 @@ function add_quadratic2{T}(Λ::PiecewiseQuadratic{T}, ρ::QuadraticPolynomial{T} λ_prev, λ_curr = update_segment_new(λ_prev, λ_curr, ρ) elseif root1 <= left_endpoint && root2 >= right_endpoint # No intersections, old quadratic is smallest, just step forward - λ_prev = λ_curr - λ_curr = λ_prev.next + λ_prev, λ_curr = update_segment_do_nothing(λ_curr) elseif root1 > left_endpoint && root2 < right_endpoint # Two intersections within the interval λ_prev, λ_curr = update_segment_new_old_new(λ_prev, λ_curr, ρ, root1, root2) @@ -332,6 +301,28 @@ function add_quadratic2{T}(Λ::PiecewiseQuadratic{T}, ρ::QuadraticPolynomial{T} end + +@inline function update_segment_do_nothing(λ_curr) + return λ_curr, λ_curr.next +end + +@inline function update_segment_new(λ_prev, λ_curr, ρ) + if λ_prev.p === ρ + λ_prev.next = λ_curr.next + v1, v2 = λ_prev, λ_curr.next + else + λ_curr.p = ρ + v1, v2 = λ_curr, λ_curr.next + end + return v1, v2 +end + +@inline function update_segment_old_new(λ_curr, ρ, break1) + new_pwq_segment = PiecewiseQuadratic(ρ, break1, λ_curr.next) + λ_curr.next = new_pwq_segment + return new_pwq_segment, new_pwq_segment.next +end + @inline function update_segment_new_old(λ_prev, λ_curr, ρ, break1) if λ_prev.p === ρ λ_curr.left_endpoint = break1 @@ -354,26 +345,6 @@ end return second_old_pwq_segment, second_old_pwq_segment.next end -@inline function update_segment_old_new(λ_curr, ρ, break1) - new_pwq_segment = PiecewiseQuadratic(ρ, break1, λ_curr.next) - λ_curr.next = new_pwq_segment - return new_pwq_segment, new_pwq_segment.next -end - -@inline function update_segment_new(λ_prev, λ_curr, ρ) - if λ_prev.p === ρ - λ_prev.next = λ_curr.next - v1, v2 = λ_prev, λ_curr.next - else - λ_curr.p = ρ - v1, v2 = λ_curr, λ_curr.next - end - return v1, v2 #λ_curr, λ_curr.next -end - -@inline function update_segment_do_nothing(λ_curr) - return λ_curr, λ_curr.next #λ_curr, λ_curr.next -end ### @@ -386,23 +357,28 @@ end @inline function minimize_wrt_x2_fast{T}(qf::QuadraticForm2{T},p::QuadraticPolynomial{T}) # Create quadratic form representing the sum of qf and p + P = qf.P q = qf.q r = qf.r P22_new = P[2,2] + p.a - local v + v = QuadraticPolynomial{T}() if P22_new > 0 - v = QuadraticPolynomial(P[1,1] - P[1,2]^2 / P22_new, - q[1] - P[1,2]*(q[2]+p.b) / P22_new, - (r+p.c) - (q[2]+p.b)^2 / P22_new/ 4) - elseif P22_new == 0 #|| P[1,2] == 0 || (q[2]+p.b) == 0 #why are the two last conditions needed? - v = QuadraticPolynomial(P[1,1], q[1], r+p.c) + v.a = P[1,1] - P[1,2]^2 / P22_new + v.b = q[1] - P[1,2]*(q[2]+p.b) / P22_new + v.c = (r+p.c) - (q[2]+p.b)^2 / P22_new / 4 + elseif P22_new == 0 #|| qf.P11 == 0 || (qf.q2+p.b) == 0 #why are the two last conditions needed? + v.a = P[1,1] + v.b = q[1] + v.c = r+p.c else # FIXME: what are these condtions? # There are some special cases, but disregards these - v = QuadraticPolynomial(0.,0.,-Inf) + v.a = 0.0 + v.b = 0.0 + v.c = 0.0 end return v end @@ -410,8 +386,8 @@ end -#global counter1 -#global counter2 +global counter1 +global counter2 """ Find optimal fit @@ -437,9 +413,6 @@ function find_optimal_fit{T}(Λ_0::Array{PiecewiseQuadratic{T},1}, ℓ::Array{Qu p = λ.p #counter1 += 1 - # ρ = 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) if unsafe_minimum(ρ) > upper_bound @@ -462,24 +435,6 @@ end -function find_minimum(Λ::PiecewiseQuadratic) - # May assume that the minimum is at the staionary point of the polynomial - # TODO: True? - - f_opt = Inf - p_opt = Λ.p # The segment containing the smallest polynimal - x_opt = NaN - - for λ in Λ - x, f = find_minimum(λ.p) - if f < f_opt - f_opt = f - x_opt = x - p_opt = λ.p - end - end - return p_opt, x_opt, f_opt -end function recover_solution{T}(Λ::PiecewiseQuadratic{T}, first_index=1, last_index=-1) @@ -491,11 +446,11 @@ function recover_solution{T}(Λ::PiecewiseQuadratic{T}, first_index=1, last_inde while true push!(I, p.time_index) - if isnull(p.ancestor); + if !isdefined(p, :ancestor); break; end - p = unsafe_get(p.ancestor) + p = p.ancestor end if last_index != -1 diff --git a/src/types/PiecewiseQuadratic.jl b/src/types/PiecewiseQuadratic.jl index ede4450..e5e28d3 100644 --- a/src/types/PiecewiseQuadratic.jl +++ b/src/types/PiecewiseQuadratic.jl @@ -35,7 +35,7 @@ end Creates piecewise-quadratic polynomial containing one element """ function create_new_pwq(p::QuadraticPolynomial) - pwq = PiecewiseQuadratic(p, -1e9, PiecewiseQuadratic{Float64}()) + pwq = PiecewiseQuadratic(p, -Inf, PiecewiseQuadratic{Float64}()) return PiecewiseQuadratic(QuadraticPolynomial([NaN,NaN,NaN]), NaN, pwq) end @@ -85,11 +85,7 @@ function delete_next{T}(pwq::PiecewiseQuadratic{T}) end function get_right_endpoint(λ::PiecewiseQuadratic) - if isnull(λ.next) - return 1e9 - else - return unsafe_get(λ.next).left_endpoint - end + return unsafe_get(λ.next).left_endpoint end @@ -106,13 +102,13 @@ function show(io::IO, Λ::PiecewiseQuadratic) #@printf("PiecewiseQuadratic (%i elems):\n", length(Λ)) for λ in Λ - if λ.left_endpoint == -1e9 + if λ.left_endpoint == -Inf print(io, "[ -∞ ,") else @printf(io, "[%3.2f,", λ.left_endpoint) end - if get_right_endpoint(λ) == 1e9 + if get_right_endpoint(λ) == Inf print(io, " ∞ ]") else @printf(io, " %3.2f]", get_right_endpoint(λ)) @@ -180,6 +176,24 @@ function plot_pwq(Λ::PiecewiseQuadratic) return p end +function find_minimum(Λ::PiecewiseQuadratic) + # May assume that the minimum is at the staionary point of the polynomial + # TODO: True? + + f_opt = Inf + p_opt = Λ.p # The segment containing the smallest polynimal + x_opt = NaN + + for λ in Λ + x, f = find_minimum(λ.p) + if f < f_opt + f_opt = f + x_opt = x + p_opt = λ.p + end + end + return p_opt, x_opt, f_opt +end diff --git a/src/types/QuadraticPolynomial.jl b/src/types/QuadraticPolynomial.jl index d652d31..54f23e3 100644 --- a/src/types/QuadraticPolynomial.jl +++ b/src/types/QuadraticPolynomial.jl @@ -6,16 +6,17 @@ mutable struct QuadraticPolynomial{T<:Real} a::T b::T c::T -ancestor::Nullable{QuadraticPolynomial{T}} time_index::Int +ancestor::QuadraticPolynomial{T} function QuadraticPolynomial{T}(a::T, b::T, c::T) where {T} # @assert a ≥ 0, # Δ may have negative a ... - new(a, b, c, Nullable{QuadraticPolynomial{T}}(),-1) + new(a, b, c,-1) end -function QuadraticPolynomial{T}(a::T, b::T, c::T, ancestor, time_index) where {T} +function QuadraticPolynomial{T}(a::T, b::T, c::T, time_index, ancestor) where {T} @assert a ≥ 0 - new(a, b, c, ancestor, time_index) + new(a, b, c, time_index, ancestor) end +QuadraticPolynomial{T}() where {T} = new() end QuadraticPolynomial{T}(a::T,b::T,c::T) = QuadraticPolynomial{T}(a,b,c) @@ -71,6 +72,22 @@ No checks of this is done end +function roots{T}(p::QuadraticPolynomial{T}) + b2_minus_4ac = p.b^2 - 4*p.a*p.c + if b2_minus_4ac < 0 + return NaN, NaN + end + + if p.a != 0 + term1 = (p.b / 2 / p.a) + term2 = sqrt(b2_minus_4ac) / 2 / abs(p.a) + return -term1-term2, -term1+term2 + else + term = -p.c / p.b + return term, Inf + end +end + """ n, intersec = intersections{T}(p1::QuadraticPolynomial{T},p2::QuadraticPolynomial{T}) diff --git a/test/find_optimal_fit.jl b/test/find_optimal_fit.jl index 51391ae..b2532fb 100644 --- a/test/find_optimal_fit.jl +++ b/test/find_optimal_fit.jl @@ -21,8 +21,11 @@ g = sin @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 +m = 5 @time I_bf, _, f_bf = dev.brute_force_optimization(ℓ, m-1); @test I_bf == [1, 21, 30, 50] diff --git a/test/piecewise_quadratics.jl b/test/piecewise_quadratics.jl index 818c0ca..4224573 100644 --- a/test/piecewise_quadratics.jl +++ b/test/piecewise_quadratics.jl @@ -18,15 +18,15 @@ 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) +dev.add_quadratic2(pwq1, p2) @test length(pwq1) == 3 -dev.add_quadratic(pwq1, p3) +dev.add_quadratic2(pwq1, p3) @test length(pwq1) == 4 pwq2 = dev.create_new_pwq(p2) -dev.add_quadratic(pwq2, p3) +dev.add_quadratic2(pwq2, p3) @test length(pwq2) == 3 -dev.add_quadratic(pwq2, p1) +dev.add_quadratic2(pwq2, p1) @test length(pwq2) == 4 @test pwq1[1].p === pwq2[1].p == p1 @@ -59,18 +59,21 @@ function verify_piecewise_quadratic(c_mat) end # Check that the intersections between the quadratics are in increasing order + 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: ") for λ in pwq - if isnull(λ.next) + x = dev.get_right_endpoint(λ) + if x == Inf break end - x = get(λ.next).left_endpoint - println("Continutiy diff: ", λ.p(x) - get(λ.next).p(x)) - @test λ.p(x) ≈ get(λ.next).p(x) + + println("Continutiy diff: ", λ.p(x) - λ.next.p(x)) + @test λ.p(x) ≈ λ.next.p(x) end # Check that the piecewise quadratic is smaller than all the quadratics @@ -105,8 +108,6 @@ verify_piecewise_quadratic(c_mat2) #--- c_mat3 = [ -#2.69548 1.50706 1.95186 -#2.86858 0.990155 1.86534 2.42298 0.796953 1.02011 2.45055 1.32235 1.4894 2.1762 0.855014 1.0647