# Load the packages 
using MosekTools, Mosek, Random, Convex, JuMP, Suppressor, Images, DelimitedFiles
# using COSMO, Random, Convex, JuMP, Suppressor, Images, DelimitedFiles

# Load the data
include("img//portfolio_data.jl");

x = Variable(n)

objective = γ*quadform(x,Σ) - dot(x,μ) 

constraint_linear = ( sum(x) == d + sum(x0) )

constraint_pos = ( x >= 0 )

problem = minimize(objective, [constraint_linear, constraint_pos])

solve!(problem, Mosek.Optimizer)

x_sol = Convex.evaluate(x)

p_star = Convex.evaluate(objective)

## Define the problem
x = Variable(n);
objective = γ*quadform(x,Σ) - dot(x,μ);
constraint_linear = ( sum(x) == d + sum(x0) );
constraint_pos = ( x >= 0 );
problem = minimize(objective, [constraint_linear, constraint_pos]);

## Solve the problem
solve!(problem, Mosek.Optimizer);

## Extract the solution
x_sol = Convex.evaluate(x);
p_star = Convex.evaluate(objective);

## Define the problem
using JuMP
model = Model(Mosek.Optimizer)
@variable(model, xs[1:n] >= 0)
@objective(model, Min,  γ * xs'*Σ*xs  - μ' * xs);
@constraint(model, sum(xs) .== d + sum(x0))

## Solve the problem
optimize!(model)

## Get the optimal solution
x_sol_JuMP = value.(xs)
p_star_JuMP = objective_value(model)

norm(x_sol - x_sol_JuMP)

## Define the problem (student)
δ = 0.5
x_sparse = Variable(n)
objective_sparse = - dot(x_sparse,μ) + γ*norm(x_sparse,1)
constraint_linear = ( sum(x_sparse) == d + sum(x0) )
constraint_pos = ( x_sparse >= 0 )
constraint_additional = ( quadform(x_sparse,Σ) <= δ)
problem_sparse = minimize(objective_sparse, [constraint_linear, constraint_pos, constraint_additional])

## Solve the problem
solve!(problem_sparse, Mosek.Optimizer)

## Extract the solution
x_sol_sparse = Convex.evaluate(x_sparse)
p_star_sparse = Convex.evaluate(objective_sparse)

using Plots, Suppressor
gaston()
# inspectdr()
plot([1:n],[x_sol_JuMP x_sol_sparse], xlabel = "k", ylabel = "x[k]", label = [ "x original" "x sparse"])

using DelimitedFiles
τ = readdlm("img//melbourne_temps.txt", ',')
n = size(τ, 1)
plot(1:n, τ[1:n], ylabel="Temperature (°C)", label="data", xlabel = "Time (days)", xticks=0:365:n)

## Solution using Convex.jl
x = Variable(n)

eq_constraints = [ x[i] == x[i - 365] for i in 365 + 1 : n ];

λ = 100 # smoothing parameter
# Define the matrix A
A = zeros(n-1,n)
for i in 1:n-1
    A[i,i] = -1
    A[i,i+1] = 1
end

smooth_objective = norm(x-τ,1) + λ*norm(A*x,2)

smooth_problem = minimize(smooth_objective, eq_constraints);

solve!(smooth_problem, Mosek.Optimizer)

# Plot smooth fit
plot(1:n, τ[1:n], label="data")
plot!(1:n, Convex.evaluate(x)[1:n], linewidth=2, label="smooth fit",  ylabel="Temperature (°C)", xticks=0:365:n, xlabel="Time (days)")

using Images
lenna = load("img//lena128missing.png")

# convert to real matrices
Y = Float64.(lenna)

observed_entries_Y = findall(x->x!=0.0, Y)

N, N = size(Y)

X = Variable(N,N)

obj_SDP = nuclearnorm(X)

constraints_SDP = Convex.Constraint[]

for index_i_j in observed_entries_Y
    i = index_i_j[1]
    j = index_i_j[2]
    push!(constraints_SDP,  X[i,j] == Y[i,j])
end

problem_SDP = minimize(obj_SDP, constraints_SDP)

# [💀] Do not run on a laptop, I am solving this on MIT Supercloud
# This is a problem requiring lot of memory
# Constraints            : 73665           
# Scalar variables       : 49153 
solve!(problem_SDP, Mosek.Optimizer)

X_sdp_sol = Convex.evaluate(X);
lenna_original = load("img//Lenna_(test_image).png")
# [lenna colorview(Gray, X_sdp_sol) lenna_original]
[lenna colorview(Gray, X_sdp_sol)]

lenna_original

using JuMP, Mosek, MosekTools, LinearAlgebra

model_SDP_JuMP = Model(Mosek.Optimizer)

@variable(model_SDP_JuMP, X2[1:N, 1:N])

@variable(model_SDP_JuMP, U[1:N, 1:N]);

@variable(model_SDP_JuMP, V[1:N, 1:N]);

@variable(model_SDP_JuMP, t);

for index_i_j in observed_entries_Y
    i = index_i_j[1]
    j = index_i_j[2]
    @constraint(model_SDP_JuMP, X2[i,j] == Y[i,j])
end

@constraint(model_SDP_JuMP, Symmetric([U X2; X2' V]) in PSDCone());

@constraint(model_SDP_JuMP, tr(U)+tr(V) <= 2*t);

@objective(model_SDP_JuMP, Min, t)

optimize!(model_SDP_JuMP) 

X_Lenna_JuMP = value.(X2)

# [lenna colorview(Gray, X_Lenna_JuMP) lenna_original]
[lenna colorview(Gray, X_Lenna_JuMP)]