using MonotoneDecomposition
using Plots
using Random
using StatsBase
using LaTeXStrings

Random.seed!(123)
x, y, x0, y0 = gen_data(100, 0.1, x -> exp(x));
μs = 10.0 .^ (-2:0.2:3)
nμ = length(μs)
res = Array{Any, 1}(undef, nμ)

26-element Vector{Any}:
 #undef
 #undef
 #undef
 #undef
 #undef
 #undef
 #undef
 #undef
 #undef
 #undef
   ⋮
 #undef
 #undef
 #undef
 #undef
 #undef
 #undef
 #undef
 #undef
 #undef

J = 8

for i in 1:nμ
    res[i] = mono_decomp_cs(x, y, J = 8, s = μs[i], s_is_μ=true)
end
γups = hcat([res[i].γup for i in 1:nμ]...)
γdowns = hcat([res[i].γdown for i in 1:nμ]...)

8×26 Matrix{Float64}:
 0.58437  0.58437  0.58437  0.58437  …  0.58437  0.58437  0.58437  0.58437
 0.58437  0.58437  0.58437  0.58437     0.58437  0.58437  0.58437  0.58437
 0.58437  0.58437  0.58437  0.58437     0.58437  0.58437  0.58437  0.58437
 0.58437  0.58437  0.58437  0.58437     0.58437  0.58437  0.58437  0.58437
 0.58437  0.58437  0.58437  0.58437     0.58437  0.58437  0.58437  0.58437
 0.58437  0.58437  0.58437  0.58437  …  0.58437  0.58437  0.58437  0.58437
 0.58437  0.58437  0.58437  0.58437     0.58437  0.58437  0.58437  0.58437
 0.58437  0.58437  0.58437  0.58437     0.58437  0.58437  0.58437  0.58437

save the results

serialize("../res/solution_path/expx-sigma0.1-J8.sil", [γups, γdowns, res, μs, x, y])

calculate the coefficients via the least square solution

γls = inv(res[1].workspace.B' * res[1].workspace.B) * res[1].workspace.B' * y
c = mean(res[1].workspace.B * γls) / 2
γups_calculated = [γls ./ (μ + 1) .+ (μ - 1) / (μ + 1) * c for μ in μs]
γups_ls = hcat(γups_calculated...)

8×26 Matrix{Float64}:
 -0.256847  -0.252003   -0.24444    …  0.582242  0.583026  0.583522
 -0.168348  -0.164015   -0.157247      0.582466  0.583167  0.583611
 -0.019505  -0.0160281  -0.0105988     0.582842  0.583405  0.583761
  0.212734   0.214874    0.218215      0.58343   0.583777  0.583996
  0.525671   0.526009    0.526537      0.584222  0.584277  0.584311
  1.33034    1.32605     1.31934    …  0.586258  0.585563  0.585123
  1.57868    1.57296     1.56402       0.586887  0.58596   0.585374
  2.23185    2.22237     2.20756       0.58854   0.587004  0.586033

plot the coefficients along the tuning parameter μ

default_colors = palette(:auto)
plot(log10.(μs), γups[1, :], label = L"\hat\gamma_1^u", xlab = L"\log_{10}\; \mu", color = default_colors[1])
plot!(log10.(μs), γups_ls[1, :], label = L"\gamma_1^{u}", ls = :dash, lw = 2, color = default_colors[1])
plot!(log10.(μs), γups[2, :], label = L"\hat\gamma_2^u", color = default_colors[2])
plot!(log10.(μs), γups_ls[2, :], label = L"\gamma_2^{u}", ls = :dash, lw = 2, color = default_colors[2])
plot!(log10.(μs), γups[3, :], label = L"\hat\gamma_3^u", color = default_colors[3])
plot!(log10.(μs), γups_ls[3, :], label = L"\gamma_3^{u}", ls = :dash, lw = 2, color = default_colors[3])
plot!(log10.(μs), γups[4, :], label = L"\hat\gamma_4^u", color = default_colors[4])
plot!(log10.(μs), γups_ls[4, :], label = L"\gamma_4^{u}", ls = :dash, lw = 2, color = default_colors[4])
plot!(log10.(μs), γups[5, :], label = L"\hat\gamma_5^u", color = default_colors[5])
plot!(log10.(μs), γups_ls[5, :], label = L"\gamma_5^{u}", ls = :dash, lw = 2, color = default_colors[5])
plot!(log10.(μs), γups[6, :], label = L"\hat\gamma_6^u", color = default_colors[6])
plot!(log10.(μs), γups_ls[6, :], label = L"\gamma_6^{u}", ls = :dash, lw = 2, color = default_colors[6])
plot!(log10.(μs), γups[7, :], label = L"\hat\gamma_7^u", color = default_colors[7])
plot!(log10.(μs), γups_ls[7, :], label = L"\gamma_7^{u}", ls = :dash, lw = 2, color = default_colors[7])
plot!(log10.(μs), γups[8, :], label = L"\hat\gamma_8^u", color = default_colors[8])
plot!(log10.(μs), γups_ls[8, :], label = L"\gamma_8^{u}", ls = :dash, lw = 2, color = default_colors[8])

save the results

savefig("../res/solution_path/expx-sigma0.1-J8-up.pdf")

plot the difference between the estimated coefficients and the theoretical constant vector

diff_downs = sum((γdowns .- c).^2, dims = 1)[:]
plot(log10.(μs), log10.(diff_downs), xlab = L"\log_{10}\; \mu", label = "", ylab = L"\log_{10}\Vert \hat\gamma^d - \gamma^d\Vert_2^2")

save the figure

savefig("../res/solution_path/expx-sigma0.1-J8-down.pdf")

J = 9

for i in 1:nμ
    res[i] = mono_decomp_cs(x, y, J = 9, s = μs[i], s_is_μ=true)
end
γups = hcat([res[i].γup for i in 1:nμ]...)
γdowns = hcat([res[i].γdown for i in 1:nμ]...)

9×26 Matrix{Float64}:
 0.58437  0.58437  0.58437  0.58437  …  0.58437  0.58437  0.58437  0.58437
 0.58437  0.58437  0.58437  0.58437     0.58437  0.58437  0.58437  0.58437
 0.58437  0.58437  0.58437  0.58437     0.58437  0.58437  0.58437  0.58437
 0.58437  0.58437  0.58437  0.58437     0.58437  0.58437  0.58437  0.58437
 0.58437  0.58437  0.58437  0.58437     0.58437  0.58437  0.58437  0.58437
 0.58437  0.58437  0.58437  0.58437  …  0.58437  0.58437  0.58437  0.58437
 0.58437  0.58437  0.58437  0.58437     0.58437  0.58437  0.58437  0.58437
 0.58437  0.58437  0.58437  0.58437     0.58437  0.58437  0.58437  0.58437
 0.58437  0.58437  0.58437  0.58437     0.58437  0.58437  0.58437  0.58437

save the results

serialize("../res/solution_path/expx-sigma0.1-J9.sil", [γups, γdowns, res, μs, x, y])

without tied information

calculate the coefficients via the least square solution

γls = inv(res[1].workspace.B' * res[1].workspace.B) * res[1].workspace.B' * y
c = mean(res[1].workspace.B * γls) / 2
γups_calculated = [γls ./ (μ + 1) .+ (μ - 1) / (μ + 1) * c for μ in μs]
γups_ls = hcat(γups_calculated...)

9×26 Matrix{Float64}:
 -0.223576    -0.218924   -0.21166    …  0.582326  0.583079  0.583555
 -0.244814    -0.24004    -0.232585      0.582272  0.583045  0.583534
  0.00942349   0.0127339   0.017903      0.582915  0.583452  0.58379
  0.0883209    0.091177    0.0956369     0.583115  0.583578  0.58387
  0.379559     0.380738    0.38258       0.583852  0.584043  0.584164
  0.740732     0.739832    0.738426   …  0.584766  0.58462   0.584528
  1.42454      1.41971     1.41215       0.586497  0.585713  0.585218
  1.65321      1.64706     1.63745       0.587075  0.586079  0.585449
  2.24483      2.23527     2.22034       0.588573  0.587024  0.586046

plot the coefficients along the tuning parameter μ

plot(log10.(μs), γups[1, :], label = L"\hat\gamma_1^u", xlab = L"\log_{10}\; \mu", color = default_colors[1], lw = 3, alpha = 0.5)
plot!(log10.(μs), γups_ls[1, :], label = L"\gamma_1^{u}", ls = :dash, lw = 2, color = default_colors[1])
plot!(log10.(μs), γups[2, :], label = L"\hat\gamma_2^u", color = default_colors[2])
plot!(log10.(μs), γups_ls[2, :], label = L"\gamma_2^{u}", ls = :dash, lw = 2, color = default_colors[2])
plot!(log10.(μs), γups[3, :], label = L"\hat\gamma_3^u", color = default_colors[3])
plot!(log10.(μs), γups_ls[3, :], label = L"\gamma_3^{u}", ls = :dash, lw = 2, color = default_colors[3])
plot!(log10.(μs), γups[4, :], label = L"\hat\gamma_4^u", color = default_colors[4])
plot!(log10.(μs), γups_ls[4, :], label = L"\gamma_4^{u}", ls = :dash, lw = 2, color = default_colors[4])
plot!(log10.(μs), γups[5, :], label = L"\hat\gamma_5^u", color = default_colors[5])
plot!(log10.(μs), γups_ls[5, :], label = L"\gamma_5^{u}", ls = :dash, lw = 2, color = default_colors[5])
plot!(log10.(μs), γups[6, :], label = L"\hat\gamma_6^u", color = default_colors[6])
plot!(log10.(μs), γups_ls[6, :], label = L"\gamma_6^{u}", ls = :dash, lw = 2, color = default_colors[6])
plot!(log10.(μs), γups[7, :], label = L"\hat\gamma_7^u", color = default_colors[7])
plot!(log10.(μs), γups_ls[7, :], label = L"\gamma_7^{u}", ls = :dash, lw = 2, color = default_colors[7])
plot!(log10.(μs), γups[8, :], label = L"\hat\gamma_8^u", color = default_colors[8])
plot!(log10.(μs), γups_ls[8, :], label = L"\gamma_8^{u}", ls = :dash, lw = 2, color = default_colors[8])
plot!(log10.(μs), γups[9, :], label = L"\hat\gamma_9^u", color = default_colors[9])
plot!(log10.(μs), γups_ls[9, :], label = L"\gamma_9^{u}", ls = :dash, lw = 2, color = default_colors[9])

save the figure

savefig("../res/solution_path/expx-sigma0.1-J9-up-ls.pdf")

diff_downs = sum((γdowns .- c).^2, dims = 1)[:]
plot(log10.(μs), log10.(diff_downs), xlab = L"\log_{10}\; \mu", label = "", ylab = L"\log_{10}\Vert \hat\gamma^d - \gamma^d\Vert_2^2")

save the figure

savefig("../res/solution_path/expx-sigma0.1-J9-down-ls.pdf")

with tied information

G = hcat(zeros(8), 1.0I(8))
G[1, 1] = 1
γGls = G' * inv(G * res[1].workspace.B' * res[1].workspace.B * G') * G * res[1].workspace.B' * y
c = mean(res[1].workspace.B * γGls) / 2
γups_calculated = [γGls ./ (μ + 1) .+ (μ - 1) / (μ + 1) * c for μ in μs]
γups_ls = hcat(γups_calculated...)

plot the results

plot(log10.(μs), γups[1, :], label = L"\hat\gamma_1^u", xlab = L"\log_{10}\; \mu", color = default_colors[1], lw = 3, alpha = 0.5)
plot!(log10.(μs), γups_ls[1, :], label = L"\gamma_1^{u}", ls = :dash, lw = 2, color = default_colors[1])
plot!(log10.(μs), γups[2, :], label = L"\hat\gamma_2^u", color = default_colors[2])
plot!(log10.(μs), γups_ls[2, :], label = L"\gamma_2^{u}", ls = :dash, lw = 2, color = default_colors[2])
plot!(log10.(μs), γups[3, :], label = L"\hat\gamma_3^u", color = default_colors[3])
plot!(log10.(μs), γups_ls[3, :], label = L"\gamma_3^{u}", ls = :dash, lw = 2, color = default_colors[3])
plot!(log10.(μs), γups[4, :], label = L"\hat\gamma_4^u", color = default_colors[4])
plot!(log10.(μs), γups_ls[4, :], label = L"\gamma_4^{u}", ls = :dash, lw = 2, color = default_colors[4])
plot!(log10.(μs), γups[5, :], label = L"\hat\gamma_5^u", color = default_colors[5])
plot!(log10.(μs), γups_ls[5, :], label = L"\gamma_5^{u}", ls = :dash, lw = 2, color = default_colors[5])
plot!(log10.(μs), γups[6, :], label = L"\hat\gamma_6^u", color = default_colors[6])
plot!(log10.(μs), γups_ls[6, :], label = L"\gamma_6^{u}", ls = :dash, lw = 2, color = default_colors[6])
plot!(log10.(μs), γups[7, :], label = L"\hat\gamma_7^u", color = default_colors[7])
plot!(log10.(μs), γups_ls[7, :], label = L"\gamma_7^{u}", ls = :dash, lw = 2, color = default_colors[7])
plot!(log10.(μs), γups[8, :], label = L"\hat\gamma_8^u", color = default_colors[8])
plot!(log10.(μs), γups_ls[8, :], label = L"\gamma_8^{u}", ls = :dash, lw = 2, color = default_colors[8])
plot!(log10.(μs), γups[9, :], label = L"\hat\gamma_9^u", color = default_colors[9])
plot!(log10.(μs), γups_ls[9, :], label = L"\gamma_9^{u}", ls = :dash, lw = 2, color = default_colors[9])

save the figure

savefig("../res/solution_path/expx-sigma0.1-J9-up-gls.pdf")