From 28bae47fc2324fa77cceb905e66b6b6c4e0e5b37 Mon Sep 17 00:00:00 2001 From: oloft Date: Thu, 7 Sep 2017 16:32:35 +0200 Subject: [PATCH] Removed old code. --- examples/example2.jl | 33 ++++--- examples/example3.jl | 2 +- src/dev.jl | 156 +------------------------------ src/types/LList.jl | 58 ------------ src/types/QuadraticForm2.jl | 11 ++- src/types/QuadraticPolynomial.jl | 46 +-------- test/piecewise_quadratics.jl | 14 +-- 7 files changed, 42 insertions(+), 278 deletions(-) delete mode 100644 src/types/LList.jl diff --git a/examples/example2.jl b/examples/example2.jl index fed277a..eb2c712 100644 --- a/examples/example2.jl +++ b/examples/example2.jl @@ -1,11 +1,11 @@ using IterTools - +using Plots include(joinpath(Pkg.dir("DynamicApproximations"),"src","dev.jl")) - +## data = readdlm(joinpath(Pkg.dir("DynamicApproximations"),"examples","data","snp500.txt")) -data = data[1:800] +data = data[1:300] N = length(data) @@ -13,18 +13,20 @@ N = length(data) @time ℓ = dev.compute_discrete_transition_costs(data); - - -K = 4 #@time I, Y, f = dev.brute_force_optimization(ℓ, K-1); -### cost_last = dev.QuadraticPolynomial(1.0, -2*data[end], data[end]^2) -Λ_0 = [dev.create_new_pwq(dev.minimize_wrt_x2_fast(ℓ[i, N], cost_last)) for i in 1:N-1]; - +Λ_0 = [dev.create_new_pwq(dev.minimize_wrt_x2(ℓ[i, N], cost_last)) for i in 1:N-1]; start_time = time() -@profiler Λ = dev.find_optimal_fit(Λ_0, ℓ, 7); +gc() +@time Λ = dev.find_optimal_fit(Λ_0, ℓ, 10, Inf); println("Time: ", time()-start_time) + +#println(dev.counter1, " ", dev.counter2) + +for k=3:10 + println(k, " : ", dev.recover_solution(Λ[k, 1], 1, N)[3]) +end ### @time I2, y2, f2 = dev.recover_solution(Λ[7, 1], 1, N) @@ -36,7 +38,8 @@ Y2, _ = dev.find_optimal_y_values(ℓ, I2) #find_optimal_y_values - +plot(layout=(1,2)) +## l = zeros(size(Λ)) for i=1:size(l,1) for j=1:size(l,2) @@ -45,8 +48,13 @@ for i=1:size(l,1) end end end +sum(l) +## +plot!(l', subplot=2) +sum(l) +## I2, y2, f2 = dev.recover_solution(Λ[6, 1], 1, N) println(I2) @@ -54,6 +62,3 @@ Y2, _ = dev.find_optimal_y_values(ℓ, I2) plot(data) plot!(I2, Y2) - - -plot(l') diff --git a/examples/example3.jl b/examples/example3.jl index d4d2550..340e08f 100644 --- a/examples/example3.jl +++ b/examples/example3.jl @@ -126,4 +126,4 @@ plot!([1.5, 2.5], [1,1]*cost_l0[3:end]') 1.89786 1.76981 1.75176 1.59654] -plot(costs) +plot(3:10,costs[3:end,:]) diff --git a/src/dev.jl b/src/dev.jl index a531dee..ed3d964 100644 --- a/src/dev.jl +++ b/src/dev.jl @@ -34,142 +34,6 @@ Inserts a quadratic polynomial ρ into the linked list Λ which represents a pie """ function add_quadratic{T}(Λ::PiecewiseQuadratic{T}, ρ::QuadraticPolynomial{T}) - if isnull(Λ.next) # I.e. the piecewise quadratic object is empty, perhaps better to add dummy polynomial - insert(Λ, ρ, -1e9) - return - end - - - λ_prev = Λ - λ_curr = Λ.next - - Δ = QuadraticPolynomial(0., 0., 0.) - while λ_curr.left_endpoint != Inf - - 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 - - Δ .= ρ .- λ_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, λ_curr, root1, (left_endpoint+root1)/2, ρ, Δ) - #println(Λ) - end - - if isnull(λ_curr) - break - end - - if root2 > left_endpoint && root2 < right_endpoint - #println("case 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 - - if isnull(λ_curr) - break - end - - #println("case 3:") - λ_prev, λ_curr = impose_quadratic_to_endpoint(λ_prev, λ_curr, (λ_curr.left_endpoint + right_endpoint)/2, ρ, Δ) - #println(Λ) - - if isnull(λ_curr) - break - end - end - -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{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{Float64}}(λ_curr) - else - #println("2") - - #λ_new = insert(λ_prev, ρ, λ_curr.left_endpoint) - #λ_curr.left_endpoint = root - #return λ_new, λ_curr - - λ_new = insert(λ_curr, λ_curr.p, root) - λ_curr.p = ρ - return λ_curr, λ_new - end - else - #println("3") - λ_new = insert(λ_curr, λ_curr.p, root) # insert new interval, but defer deciding on wheather λ_curr.p or ρ is smallest - return λ_curr, λ_new - end -end - -@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) - v1, v2 = λ_prev, λ_new - else - #println("2") - λ_curr.p = ρ - v1, v2 = λ_curr, λ_curr.next - end - else - # Do nothing - #println("3") - 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 - r = qf.r - - if P[2,2] > 0 - QuadraticPolynomial(P[1,1] - P[1,2]^2 / P[2,2], - q[1] - P[1,2]*q[2] / P[2,2], - r - q[2]^2 / P[2,2]/ 4) - elseif P[2,2] == 0 || P[1,2] == 0 || q[2] == 0 - QuadraticPolynomial(P[1,1], q[1], r) - else - # There are some special cases, but disregards these - QuadraticPolynomial(0.,0.,-Inf) - end -end - - - -## - -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 @@ -351,19 +215,10 @@ end return second_old_pwq_segment, second_old_pwq_segment.next end -### - - - - - # Takes a quadratic form in [x1; x2] and a polynomial in x2 # and returns the minimum of the sum wrt to x2, # i.e. a polynomial of x1 -@inline function minimize_wrt_x2_fast{T}(qf::QuadraticForm2{T},p::QuadraticPolynomial{T},ρ=QuadraticPolynomial{T}()) - - # Create quadratic form representing the sum of qf and p - +@inline function minimize_wrt_x2{T}(qf::QuadraticForm2{T},p::QuadraticPolynomial{T},ρ=QuadraticPolynomial{T}()) P = qf.P q = qf.q r = qf.r @@ -421,13 +276,13 @@ function find_optimal_fit{T}(Λ_0::Array{PiecewiseQuadratic{T},1}, ℓ::Array{Qu #counter1 += 1 - dev.minimize_wrt_x2_fast(ℓ[i,ip], p, ρ) + minimize_wrt_x2(ℓ[i,ip], p, ρ) if unsafe_minimum(ρ) > upper_bound continue end - add_quadratic2(Λ_new, ρ) + add_quadratic(Λ_new, ρ) if ρ.has_been_used == true ρ.time_index = ip @@ -445,10 +300,6 @@ end - - - - function recover_solution{T}(Λ::PiecewiseQuadratic{T}, first_index=1, last_index=-1) p, y, f = find_minimum(Λ) @@ -501,7 +352,6 @@ function compute_transition_costs(g, t::AbstractArray) I_tg[i] = I_tg[i-1] + quadgk(t -> t*g(t), t[i-1], t[i], reltol=tol)[1] end - ℓ = Array{QuadraticForm2{T}}(N-1,N) for i=1:N-1 diff --git a/src/types/LList.jl b/src/types/LList.jl deleted file mode 100644 index f2a65b9..0000000 --- a/src/types/LList.jl +++ /dev/null @@ -1,58 +0,0 @@ -type LList{T} - val::T - next::LList{T} - LList{T}() where T = new() #Hack for empty list - LList{T}(val::T) where T = (l = new(val); l.next = l) - LList{T}(val::T, l::LList{T}) where T = new(val, l) -end -LList{T}(val::T) = LList{T}(val) -LList{T}(val::T, l::LList{T}) = LList{T}(val, l) - -llist{T}(val::T) = LList(val) -llist{T}(val::T, args...) = LList(val, llist(args...)) -#Hack for empty set -Base.isempty(l::LList) = !isdefined(l, :next) - -function Base.show{T}(io::IO, l::LList{T}) - print(io, "llist(") - if !isempty(l) - show(io, l.val) - if !islast(l) - for val in l.next - print(io, ",") - show(io, val) - end - end - end - print(io, ")") -end - -islast(l::LList) = l === l.next || isempty(l.next) - -""" -Insert `val` in `l=llist(v1,v2,...,vn)` so that -`l=llist(v1,val,v2,...,vn)`. -If `l` is empty then `l` becomes `llist(val)`. -""" -function insertafter!{T}(val::T, l::LList{T}) - if isempty(l) - l.val = val - l.next = l - else - l2 = islast(l) ? LList{T}(val) : LList{T}(val, l.next) - l.next = l2 - end - return l -end - -Base.start{T}(l::LList{T}) = (l, isempty(l)) -Base.done{T}(l::LList{T}, state::Tuple{LList{T},Bool}) = state[2] -Base.next{T}(l::LList{T}, state::Tuple{LList{T},Bool}) = (state[1].val, (state[1].next, islast(state[1]))) - -function Base.length(l::LList) - n = 0 - for i in l - n += 1 - end - n -end diff --git a/src/types/QuadraticForm2.jl b/src/types/QuadraticForm2.jl index 439c123..7ef13bf 100644 --- a/src/types/QuadraticForm2.jl +++ b/src/types/QuadraticForm2.jl @@ -1,10 +1,19 @@ # Quadratic form of 2 variables, used for representing the transition costs ℓ -type QuadraticForm2{T} +struct QuadraticForm2{T} P::SMatrix{2,2,T,4} q::SVector{2,T} r::T end +struct QuadraticFormTest{T} + P11::T + P12::T + P22::T + q1::T + q2::T + r::T +end + function (qf::QuadraticForm2)(y, yp) return [y,yp]'*qf.P*[y,yp] + qf.q'*[y, yp] + qf.r end diff --git a/src/types/QuadraticPolynomial.jl b/src/types/QuadraticPolynomial.jl index 616ac91..fea9233 100644 --- a/src/types/QuadraticPolynomial.jl +++ b/src/types/QuadraticPolynomial.jl @@ -1,7 +1,7 @@ # Quadratic Polynomial to represent cost-to-go function, # keeos track of its "ancestor" to simplify the recovery of the solution -# TODO: penality for making mutable, so ancestor can be set later on, -# or set ancestor at creation +# has_been_used keeps track of if the Quadratic polynimial has been inserted +# into a piecewise quadratic, otherwise it is reused mutable struct QuadraticPolynomial{T<:Real} a::T b::T @@ -88,45 +88,3 @@ function roots{T}(p::QuadraticPolynomial{T}) return term, Inf end end - - -""" -n, intersec = intersections{T}(p1::QuadraticPolynomial{T},p2::QuadraticPolynomial{T}) -returns number of intersections `n` and a 2-vector `intersec` containing the intersection points -in the first `n` elements. -""" -function intersections{T}(p1::QuadraticPolynomial{T},p2::QuadraticPolynomial{T}) - a = p1.a-p2.a - b = p1.b-p2.b - c = p1.c-p2.c - n = 0 - intersec = [T(NaN), T(NaN)] - if a == 0 - #Numerical problems if division, 0 or 1 intersections - if b == 0 - # 0 intersections - #return 0, [T(NaN), T(NaN)] - else - # 1 intersection - n = 1 - intersec[1] = -c/b - #return 1, [-c/b, T(NaN)] - end - else - # 0 or 2 intersections, possibly in same point - r = b^2 - 4c*a - if r < 0 - # 0 intersections - #return 0, [T(NaN), T(NaN)] - else - # 2 intersections - r = sqrt(r)/(2a) #sqrt(b^2-4ca)/2a - v = -b/(2a) - n = 2 - intersec[1] = v-r - intersec[2] = v+r - #return 2, [v-r, v+r] - end - end - return n, intersec -end diff --git a/test/piecewise_quadratics.jl b/test/piecewise_quadratics.jl index ad0e7db..e7c233d 100644 --- a/test/piecewise_quadratics.jl +++ b/test/piecewise_quadratics.jl @@ -5,7 +5,7 @@ include(joinpath(Pkg.dir("DynamicApproximations"),"src","dev.jl")) pwq = dev.generate_PiecewiseQuadratic([([2.0, 2, 1], -Inf), ([2.0, -2, 1], 0.0)]) @test length(pwq) == 2 -dev.add_quadratic2(pwq, dev.QuadraticPolynomial([1, 0, 1.0])) +dev.add_quadratic(pwq, dev.QuadraticPolynomial([1, 0, 1.0])) @test length(pwq) == 3 @@ -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_quadratic2(pwq1, p2) +dev.add_quadratic(pwq1, p2) @test length(pwq1) == 3 -dev.add_quadratic2(pwq1, p3) +dev.add_quadratic(pwq1, p3) @test length(pwq1) == 4 pwq2 = dev.create_new_pwq(p2) -dev.add_quadratic2(pwq2, p3) +dev.add_quadratic(pwq2, p3) @test length(pwq2) == 3 -dev.add_quadratic2(pwq2, p1) +dev.add_quadratic(pwq2, p1) @test length(pwq2) == 4 @test pwq1[1].p === pwq2[1].p == p1 @@ -52,12 +52,12 @@ function verify_piecewise_quadratic(c_mat, show_plot=false) plot!(xgrid, p.(xgrid)) end end - + # Insert quadratics into piecewise quadfatic pwq = dev.create_new_pwq() for p in poly_list println("Inserting: ", p) - dev.add_quadratic2(pwq, p) + dev.add_quadratic(pwq, p) println("After Insertion: ") println(pwq)