Skip to content

Commit

Permalink
New version of linked list.
Browse files Browse the repository at this point in the history
  • Loading branch information
@olof3
olof3 committed Sep 7, 2017
1 parent e948bb6 commit 423c248
Show file tree
Hide file tree
Showing 5 changed files with 147 additions and 71 deletions.
82 changes: 48 additions & 34 deletions src/dev.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -65,32 +65,32 @@ function add_quadratic{T}(Λ::PiecewiseQuadratic{T}, ρ::QuadraticPolynomial{T})
λ_curr = Λ.next

Δ = QuadraticPolynomial(0., 0., 0.)
while ~isnull(λ_curr)
while λ_curr.left_endpoint != Inf

left_endpoint = unsafe_get(λ_curr).left_endpoint
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(unsafe_get(λ_curr))
right_endpoint = get_right_endpoint(λ_curr)

if right_endpoint == 1e9
#println("inf")
right_endpoint = left_endpoint + 20000.0
end

Δ .= ρ .- unsafe_get(λ_curr).p
Δ .= ρ .- λ_curr.p

root1,root2 = roots(Δ)

#println(root1, " : ", root2, " .... (", (left_endpoint, right_endpoint))

if root1 > left_endpoint && root1 < right_endpoint
#println("case 1:")
λ_prev, λ_curr = impose_quadratic_to_root(λ_prev, unsafe_get(λ_curr), root1, (left_endpoint+root1)/2, ρ, Δ)
λ_prev, λ_curr = impose_quadratic_to_root(λ_prev, λ_curr, root1, (left_endpoint+root1)/2, ρ, Δ)
#println(Λ)
end

Expand All @@ -100,7 +100,7 @@ function add_quadratic{T}(Λ::PiecewiseQuadratic{T}, ρ::QuadraticPolynomial{T})

if root2 > left_endpoint && root2 < right_endpoint
#println("case 2:")
λ_prev, λ_curr = impose_quadratic_to_root(λ_prev, unsafe_get(λ_curr), root2, (root1 + root2)/2, ρ, Δ)
λ_prev, λ_curr = impose_quadratic_to_root(λ_prev, λ_curr, root2, (root1 + root2)/2, ρ, Δ)
#println(Λ)
if Δ.a > 0; return; end # Saves perhaps 5% of computation time
end
Expand All @@ -110,7 +110,7 @@ function add_quadratic{T}(Λ::PiecewiseQuadratic{T}, ρ::QuadraticPolynomial{T})
end

#println("case 3:")
λ_prev, λ_curr = impose_quadratic_to_endpoint(λ_prev, unsafe_get(λ_curr), (unsafe_get(λ_curr).left_endpoint + right_endpoint)/2, ρ, Δ)
λ_prev, λ_curr = impose_quadratic_to_endpoint(λ_prev, λ_curr, (λ_curr.left_endpoint + right_endpoint)/2, ρ, Δ)
#println(Λ)

if isnull(λ_curr)
Expand Down Expand Up @@ -191,8 +191,9 @@ end

function add_quadratic2{T}(Λ::PiecewiseQuadratic{T}, ρ::QuadraticPolynomial{T})


DEBUG && println("Inserting: ", ρ)
if isnull(Λ.next) # I.e. the piecewise quadratic object is empty, perhaps better to add dummy polynomial
if Λ.next.left_endpoint == Inf # I.e. the piecewise quadratic object is empty, perhaps better to add dummy polynomial
insert(Λ, ρ, -1e9)
return
end
Expand All @@ -202,34 +203,35 @@ function add_quadratic2{T}(Λ::PiecewiseQuadratic{T}, ρ::QuadraticPolynomial{T}

#left_endpoint = NaN

while ~isnull(λ_curr)
while λ_curr.left_endpoint != Inf #???
#global counter2 += 1
DEBUG && println(Λ)

left_endpoint = unsafe_get(λ_curr).left_endpoint
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(unsafe_get(λ_curr))
right_endpoint = get_right_endpoint(λ_curr)

if right_endpoint == 1e9
#println("inf")
right_endpoint = left_endpoint + 20000.0
end

Δa = ρ.a - unsafe_get(λ_curr).p.a
Δb = ρ.b - unsafe_get(λ_curr).p.b
Δc = ρ.c - unsafe_get(λ_curr).p.c
Δa = ρ.a - λ_curr.p.a
Δb = ρ.b - λ_curr.p.b
Δc = ρ.c - λ_curr.p.c

b2_minus_4ac = Δb^2 - 4*Δa*Δc

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 = unsafe_get(λ_curr)
λ_prev = λ_curr
λ_curr = λ_prev.next
else

Expand All @@ -244,29 +246,29 @@ function add_quadratic2{T}(Λ::PiecewiseQuadratic{T}, ρ::QuadraticPolynomial{T}
if root1 >= right_endpoint || root2 <= left_endpoint
# No intersections, old quadratic is smallest, step forward
DEBUG && println("Two intersections to the side")
λ_prev = unsafe_get(λ_curr)
λ_prev = λ_curr
λ_curr = λ_prev.next
elseif root1 <= left_endpoint && root2 >= right_endpoint
# No intersections, new quadratic is smallest
DEBUG && println("One intersections on either side")
λ_prev, λ_curr = update_segment_new(λ_prev, unsafe_get(λ_curr), ρ)
λ_prev, λ_curr = update_segment_new(λ_prev, λ_curr, ρ)
elseif root1 > left_endpoint && root2 < right_endpoint
DEBUG && println("Two intersections within the interval")
λ_prev, λ_curr = update_segment_old_new_old(unsafe_get(λ_curr), ρ, root1, root2)
λ_prev, λ_curr = update_segment_old_new_old(λ_curr, ρ, root1, root2)
elseif root1 > left_endpoint
DEBUG && println("Root 1 within the interval")
λ_prev, λ_curr = update_segment_old_new(unsafe_get(λ_curr), ρ, root1)
λ_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, unsafe_get(λ_curr), ρ, root2)
λ_prev, λ_curr = update_segment_new_old(λ_prev, λ_curr, ρ, root2)
else
error("Shouldn't end up here")
end
end

elseif Δa < 0 # ρ has lower curvature, i.e., ρ is smallest on the sides
if b2_minus_4ac <= 0
λ_prev, λ_curr = update_segment_new(λ_prev, unsafe_get(λ_curr), ρ)
λ_prev, λ_curr = update_segment_new(λ_prev, λ_curr, ρ)
else
# Compute the intersections
term1 = -(Δb / 2 / Δa)
Expand All @@ -277,18 +279,18 @@ function add_quadratic2{T}(Λ::PiecewiseQuadratic{T}, ρ::QuadraticPolynomial{T}
# Check where the intersections are and act accordingly
if root1 >= right_endpoint || root2 <= left_endpoint
# No intersections, ρ is smallest
λ_prev, λ_curr = update_segment_new(λ_prev, unsafe_get(λ_curr), ρ)
λ_prev, λ_curr = update_segment_new(λ_prev, λ_curr, ρ)
elseif root1 <= left_endpoint && root2 >= right_endpoint
# No intersections, old quadratic is smallest, just step forward
λ_prev = unsafe_get(λ_curr)
λ_prev = λ_curr
λ_curr = λ_prev.next
elseif root1 > left_endpoint && root2 < right_endpoint
# Two intersections within the interval
λ_prev, λ_curr = update_segment_new_old_new(λ_prev, unsafe_get(λ_curr), ρ, root1, root2)
λ_prev, λ_curr = update_segment_new_old_new(λ_prev, λ_curr, ρ, root1, root2)
elseif root1 > left_endpoint
λ_prev, λ_curr = update_segment_new_old(λ_prev, unsafe_get(λ_curr), ρ, root1)
λ_prev, λ_curr = update_segment_new_old(λ_prev, λ_curr, ρ, root1)
elseif root2 < right_endpoint
λ_prev, λ_curr = update_segment_old_new(unsafe_get(λ_curr), ρ, root2)
λ_prev, λ_curr = update_segment_old_new(λ_curr, ρ, root2)
else
error("Shouldn't end up here")
end
Expand All @@ -298,29 +300,29 @@ function add_quadratic2{T}(Λ::PiecewiseQuadratic{T}, ρ::QuadraticPolynomial{T}

if Δb == 0
if Δc >= 0
λ_prev, λ_curr = update_segment_do_nothing(unsafe_get(λ_curr))
λ_prev, λ_curr = update_segment_do_nothing(λ_curr)
else
λ_prev, λ_curr = update_segment_new(λ_prev, unsafe_get(λ_curr), ρ)
λ_prev, λ_curr = update_segment_new(λ_prev, λ_curr, ρ)
end
continue
end

root = -Δc / Δb
if Δb > 0
if root < left_endpoint
λ_prev, λ_curr = update_segment_do_nothing(unsafe_get(λ_curr))
λ_prev, λ_curr = update_segment_do_nothing(λ_curr)
elseif root > right_endpoint
λ_prev, λ_curr = update_segment_new(λ_prev, unsafe_get(λ_curr), ρ)
λ_prev, λ_curr = update_segment_new(λ_prev, λ_curr, ρ)
else
λ_prev, λ_curr = update_segment_new_old(λ_prev, unsafe_get(λ_curr), ρ, root)
λ_prev, λ_curr = update_segment_new_old(λ_prev, λ_curr, ρ, root)
end
else
if root < left_endpoint
λ_prev, λ_curr = update_segment_new(λ_prev, unsafe_get(λ_curr), ρ)
λ_prev, λ_curr = update_segment_new(λ_prev, λ_curr, ρ)
elseif root > right_endpoint
λ_prev, λ_curr = update_segment_do_nothing(unsafe_get(λ_curr))
λ_prev, λ_curr = update_segment_do_nothing(λ_curr)
else
λ_prev, λ_curr = update_segment_old_new(unsafe_get(λ_curr), ρ, root)
λ_prev, λ_curr = update_segment_old_new(λ_curr, ρ, root)
end
end
end
Expand Down Expand Up @@ -405,10 +407,21 @@ end
return v
end




#global counter1
#global counter2

"""
Find optimal fit
"""
function find_optimal_fit{T}(Λ_0::Array{PiecewiseQuadratic{T},1}, ℓ::Array{QuadraticForm2{T},2}, M::Int, upper_bound=Inf)
#global counter1
#global counter2
#counter1 = 0
#counter2 = 0

N = size(ℓ, 2)

Λ = Array{PiecewiseQuadratic{T}}(M, N)
Expand All @@ -423,6 +436,7 @@ function find_optimal_fit{T}(Λ_0::Array{PiecewiseQuadratic{T},1}, ℓ::Array{Qu
for λ in Λ[m-1, ip]
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
Expand Down
65 changes: 36 additions & 29 deletions src/types/PiecewiseQuadratic.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,23 +7,55 @@
mutable struct PiecewiseQuadratic{T}
p::QuadraticPolynomial{T}
left_endpoint::Float64
next::Nullable{PiecewiseQuadratic{T}}
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

"""
Constructs piece of quadratic polynomial poiting to NULL, i.e.
a rightmost piece, or one that has not been inserted into a piecewise quadfatic
"""
function PiecewiseQuadratic{T}(p::QuadraticPolynomial{T}, left_endpoint)
PiecewiseQuadratic(p, left_endpoint, Nullable{PiecewiseQuadratic{T}}())
PiecewiseQuadratic{T}(p, left_endpoint, PiecewiseQuadratic{T}())
end
function PiecewiseQuadratic{T}(p::QuadraticPolynomial{T}, left_endpoint, next::PiecewiseQuadratic{T})
PiecewiseQuadratic{T}(p, left_endpoint, next)
end


"""
Creates empty piecewise-quadratic polynomial consisting only of the list head
"""
function create_new_pwq()
return PiecewiseQuadratic(QuadraticPolynomial([NaN,NaN,NaN]), NaN, PiecewiseQuadratic{Float64}())
end

"""
Creates piecewise-quadratic polynomial containing one element
"""
function create_new_pwq(p::QuadraticPolynomial)
pwq = PiecewiseQuadratic(p, -1e9, PiecewiseQuadratic{Float64}())
return PiecewiseQuadratic(QuadraticPolynomial([NaN,NaN,NaN]), NaN, pwq)
end

function generate_PiecewiseQuadratic(args)
return PiecewiseQuadratic(QuadraticPolynomial([NaN,NaN,NaN]), NaN, _generate_PiecewiseQuadratic_helper(args))
end

function _generate_PiecewiseQuadratic_helper(args)
if Base.length(args) > 1
return PiecewiseQuadratic(QuadraticPolynomial(args[1][1]), args[1][2], _generate_PiecewiseQuadratic_helper(args[2:end]))
else
return PiecewiseQuadratic(QuadraticPolynomial(args[1][1]), args[1][2])
end
end


start{T}(pwq::PiecewiseQuadratic{T}) = pwq.next
done{T}(pwq::PiecewiseQuadratic{T}, iterstate::Nullable{PiecewiseQuadratic{T}}) = isnull(iterstate)
next{T}(pwq::PiecewiseQuadratic{T}, iterstate::Nullable{PiecewiseQuadratic{T}}) = (unsafe_get(iterstate), unsafe_get(iterstate).next)
done{T}(pwq::PiecewiseQuadratic{T}, iterstate::PiecewiseQuadratic{T}) = (iterstate.left_endpoint == Inf)
next{T}(pwq::PiecewiseQuadratic{T}, iterstate::PiecewiseQuadratic{T}) = (iterstate, iterstate.next)


# For trouble-shooting etc.
function getindex(Λ::dev.PiecewiseQuadratic, n::Int64)
Expand Down Expand Up @@ -148,32 +180,7 @@ function plot_pwq(Λ::PiecewiseQuadratic)
return p
end

"""
Creates empty piecewise-quadratic polynomial consisting only of the list head
"""
function create_new_pwq()
return PiecewiseQuadratic(QuadraticPolynomial([NaN,NaN,NaN]), NaN, Nullable{PiecewiseQuadratic{Float64}}())
end

"""
Creates piecewise-quadratic polynomial containing one element
"""
function create_new_pwq(p::QuadraticPolynomial)
pwq = PiecewiseQuadratic(p, -1e9, Nullable{PiecewiseQuadratic{Float64}}())
return PiecewiseQuadratic(QuadraticPolynomial([NaN,NaN,NaN]), NaN, pwq)
end

function generate_PiecewiseQuadratic(args)
return PiecewiseQuadratic(QuadraticPolynomial([NaN,NaN,NaN]), NaN, _generate_PiecewiseQuadratic_helper(args))
end

function _generate_PiecewiseQuadratic_helper(args)
if Base.length(args) > 1
return PiecewiseQuadratic(QuadraticPolynomial(args[1][1]), args[1][2], _generate_PiecewiseQuadratic_helper(args[2:end]))
else
return PiecewiseQuadratic(QuadraticPolynomial(args[1][1]), args[1][2])
end
end



Expand Down
28 changes: 28 additions & 0 deletions test/find_optimal_fit.jl
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,28 @@
using Base.Test

include(joinpath(Pkg.dir("DynamicApproximations"),"src","dev.jl"))

N = 50
t = linspace(0,4π,N)
g = sin

@time ℓ = dev.compute_transition_costs(g, t);

Λ_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 I_dp3, _, f_dp3 = dev.recover_solution(Λ[3, 1], 1, N)
@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 of brute force optimization
m = 3
@time I_bf, _, f_bf = dev.brute_force_optimization(ℓ, m-1);
@test I_bf == [1, 21, 30, 50]
Loading

0 comments on commit 423c248

Please sign in to comment.