Updated to work well with project files, updated examples

mfalt · Apr 23, 2020 · 1b4d100 · 1b4d100
1 parent 6c77334
commit 1b4d100
Show file tree

Hide file tree

Showing 10 changed files with 54 additions and 29 deletions.
diff --git a/Project.toml b/Project.toml
@@ -0,0 +1,20 @@
+name = "EllZeroTrendFiltering"
+uuid = "21f7a131-1a4b-4cd0-abd1-bca5723008be"
+authors = ["Mattias Fält <faldt.mattias@gmail.com>, Olof Troeng"]
+version = "0.1.0"
+
+[deps]
+DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
+IterTools = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+
+[extras]
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
+
+[targets]
+test = ["Test", "Random", "Interpolations"]
diff --git a/README.md b/README.md
@@ -45,9 +45,9 @@ Ivec, Yvec, fvec = fit_pwl_constrained(data, M)
 #Plot original data
 plot(data, l=:black, lab = "SNP500")
 #Plot solution with 5 segments
-plot!(Ivec[5], Yvec[5], l=2, m=:circle, lab="m=5, cost = $(round(fvec[5],3))")
+plot!(Ivec[5], Yvec[5], l=2, m=:circle, lab="m=5, cost = $(round(fvec[5],digits=3))")
 #Plot solution with 10 segments
-plot!(Ivec[M], Yvec[M], l=2, m=:circle, lab="m=$M, cost = $(round(fvec[M],3))")
+plot!(Ivec[M], Yvec[M], l=2, m=:circle, lab="m=$M, cost = $(round(fvec[M],digits=3))")
 
 ```
 ![Example figure](figures/snp500.svg)
@@ -67,7 +67,7 @@ for ζ ∈ [0.1, 0.002]
     # Will automatically integrate the function to compute the costs
     I, Y, cost = fit_pwl_regularized(g, t, ζ)
 
-    plot!(t[I], Y, l=2, m=:circle, lab = "l2-norm=$(round(cost,3)), zeta=$ζ")
+    plot!(t[I], Y, l=2, m=:circle, lab = "l2-norm=$(round(cost,digits=3)), zeta=$ζ")
 end
 plot!() # Show plot
 ```

diff --git a/examples/example1.jl b/examples/example1.jl
@@ -1,17 +1,17 @@
 using Plots
 using Convex
-using Mosek
+using Mosek, MosekTools
 using SCS
 using DelimitedFiles
-
+using SparseArrays
 
 data = readdlm(joinpath(dirname(@__FILE__),"data","snp500.txt"))
 n = size(data)[1]
 
 #p = plot(1:length(data), data)
 
 N = length(data)
-H =  spdiagm((ones(N-2), -2*ones(N-2), ones(N-2)), (0,1,2))
+H =  spdiagm(0=>ones(N-2), 1=>-2*ones(N-2), 2=>ones(N-2))
 
 
 x = Variable(N)
@@ -20,7 +20,8 @@ x = Variable(N)
 problem = minimize(0.5*sumsquares(x - data) + λ*norm(H*x, 1))
 
 # Solve the problem by calling solve!
-solve!(problem, MosekSolver())
+solve!(problem, Mosek.Optimizer)
+solve!(problem, SCS.Optimizer)
 
 # Check the status of the problem
 problem.status # :Optimal, :Infeasible, :Unbounded etc.

diff --git a/examples/example2.jl b/examples/example2.jl
@@ -23,9 +23,9 @@ end
 
 
 
-@time I2, y2, f2 = recover_optimal_index_set(Λ[7, 1], 1, N)
+@time I2, y2, f2 = recover_optimal_index_set(Λ[7, 1])
 println(I2)
-Y2, f2_2 = find_optimal_y_values(ℓ, cost_last, I2)
+Y2, f2_2 = EllZeroTrendFiltering.find_optimal_y_values(ℓ, cost_last, I2)
 
 
 I2, Y2, f3 = recover_solution(Λ[7, 1], ℓ, cost_last, )
@@ -45,9 +45,9 @@ plot!(l', subplot=2)
 sum(l)
 ##
 
-I2, y2, f2 = recover_solution(Λ[6, 1], 1, N)
+I2, y2, f2 = EllZeroTrendFiltering.recover_solution(Λ[6, 1], ℓ, cost_last)
 println(I2)
-Y2, _ = find_optimal_y_values(ℓ, I2)
+Y2, _ = EllZeroTrendFiltering.find_optimal_y_values(ℓ, cost_last, I2)
 
 plot(data)
 plot!(I2, Y2)
diff --git a/examples/example3.jl b/examples/example3.jl
@@ -4,6 +4,7 @@ using Mosek
 using SCS
 using Interpolations
 using DelimitedFiles
+using SparseArrays
 
 # This is a julia implementaiton of an example
 # from "ℓ1 Trend Filtering" by Kim et al (2009)
@@ -21,7 +22,7 @@ I_sets = [
 [1, 365, 515, 545, 628, 769, 837, 987, 1190, 1853, 2000]
 ]
 
-cost_l0 = [find_optimal_y_values(ℓ, I_sets[k])[2] for k=1:10]
+#cost_l0 = [find_optimal_y_values(ℓ, I_sets[k])[2] for k=1:10]
 
 
 #---
@@ -35,18 +36,19 @@ function is_local_max(y)
 end
 
 #---
-data = readdlm(joinpath(joinpath(dirname(@__FILE__),"data",,"snp500.txt"))
+data = readdlm(joinpath(joinpath(dirname(@__FILE__),"data","snp500.txt")))
 
 N = length(data)
-H =  spdiagm((ones(N-2), -2*ones(N-2), ones(N-2)), (0,1,2))
+H =  spdiagm(0=>ones(N-2), 1=>-2*ones(N-2), 2=>ones(N-2))
 
 x = Variable(N)
 
-
+ℓ = compute_discrete_transition_costs(data);
+cost_last = QuadraticPolynomial(1.0, -2*data[end], data[end]^2)
 #---
 r = 6
 I = I_sets[r]
-Y, _ = find_optimal_y_values(ℓ, I)
+Y, _ = EllZeroTrendFiltering.find_optimal_y_values(ℓ, cost_last, I)
 y_l0 = interpolate((I,), Y, Gridded(Linear()))[1.0:N]
 
 
@@ -62,7 +64,7 @@ println("Opt. val: ",problem.optval, "   (problem status: ", problem.status, ")"
 
 y_l1 = evaluate(x)[:]
 
-
+### THE REST OF THIS FILE IS OUTDATED
 
 v = sortperm(abs.(H*y_l1) + 1*is_local_max(abs.(H*y_l1)))
 

diff --git a/examples/example4.jl b/examples/example4.jl
@@ -12,14 +12,15 @@ N = length(g)
 
 
 ## ℓ_1 trend filtering, see "ℓ_1 Trend Filtering", Kim et al. (2009)
-H =  spdiagm((ones(N-2), -2*ones(N-2), ones(N-2)), (0,1,2))
+H =  spdiagm(0=>ones(N-2), 1=>-2*ones(N-2), 2=>ones(N-2))
 
 x = Variable(N)
 
 λ = 2.0
 problem = minimize(sumsquares(x - g) + λ * norm(H*x, 1))
 
-solve!(problem, MosekSolver())
+solve!(problem, Mosek.Optimizer)
+#solve!(problem, SCS.Optimizer)
 println(problem.status)
 
 x_l1 = evaluate(x)[:]
@@ -34,8 +35,8 @@ I_l0, Y_l0, f_l0 = fit_pwl_regularized(g, ζ)
 
 plot(g, label="SP500", size=(1000,550));
 plot!(x_l1, lw=2.5, label="l1 reg. lambda = $λ");
-plot!(I_l0, Y_l0, lw=2.5, label="l0 reg. zeta = $ζ");
-gui()
+plot!(I_l0, Y_l0, lw=2.5, label="l0 reg. zeta = $ζ")
 
-plot(H*x_sol);
-plot!(I_vec[m][2:end], 0.1*diff(Y_vec[m]) ./ diff(I_vec[m]), m=:o, lw=0)
+# plot(H*x.value);
+# plot!(I_l0[2:end], 0.1*diff(Y_l0) ./ diff(I_l0), m=:o, lw=0)
+#
diff --git a/examples/example6.jl b/examples/example6.jl
@@ -1,3 +1,4 @@
+### Note this code is outdated
 using Plots
 using Convex
 using Mosek

diff --git a/examples/plotall.jl b/examples/plotall.jl
@@ -1,5 +1,5 @@
-
-using Base.Test
+# Note this file is outdates
+using Test, Printf, Random
 
 # TODO What is this?
 #@everywhere include(joinpath(dirname(@__FILE__),"..","src","jl"))
@@ -28,7 +28,7 @@ function plotall(Λ,ℓ,t,g, both=true)
                 plot!(p, y, x, l=:red, subplot=k, lab="$I, $ystr, $fstr", ylims=(-4,4), xlims=(0,6))
                 if both && length(I) > 1 && fI != NaN
                     Y, _ = find_optimal_y_values(ℓ, I)
-                    tcont = range(t[1], stop=t[end], length=100*length(t)))
+                    tcont = range(t[1], stop=t[end], length=100*length(t))
                     plot!(p, tcont/t[end]*(N-1)+1, g.(tcont), l=:blue, lab="", subplot=k+N-1)
                     plot!(p, I, Y, m=:circle, l=:red, lab="", subplot=k+N-1)
                 end
@@ -41,7 +41,7 @@ function plotall(Λ,ℓ,t,g, both=true)
 end
 
 ##
-srand(9)
+Random.seed!(9)
 N = 20
 ran = cumsum(0.1*randn(10*N))
 

diff --git a/src/solve.jl b/src/solve.jl
@@ -428,7 +428,7 @@ end
 
 # TODO: Include mζ in the cost?!
 """
-    I, Y, f = recover_solution(Λ::PiecewiseQuadratic{T}, ℓ, V_N::QuadraticPolynomial, first_index=1)
+    I, Y, f = recover_solution(Λ::PiecewiseQuadratic{T}, ℓ, V_N::QuadraticPolynomial, ζ=0.0)
 
 """
 function recover_solution(Λ::PiecewiseQuadratic, ℓ, V_N::QuadraticPolynomial, ζ=0.0)

diff --git a/test/examples.jl b/test/examples.jl
@@ -1,5 +1,5 @@
 using StaticArrays, EllZeroTrendFiltering
-using Test, Pkg
+using Test
 
 #SNP example
 N = 400