transitioning_from_matlab.jl

# Copyright (c) 2024 Mateus Araújo and contributors                              #src
#                                                                                #src
# Permission is hereby granted, free of charge, to any person obtaining a copy   #src
# of this software and associated documentation files (the "Software"), to deal  #src
# in the Software without restriction, including without limitation the rights   #src
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell      #src
# copies of the Software, and to permit persons to whom the Software is          #src
# furnished to do so, subject to the following conditions:                       #src
#                                                                                #src
# The above copyright notice and this permission notice shall be included in all #src
# copies or substantial portions of the Software.                                #src
#                                                                                #src
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR     #src
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,       #src
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE    #src
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER         #src
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,  #src
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE  #src
# SOFTWARE.                                                                      #src

# # Transitioning from MATLAB

# [YALMIP](https://yalmip.github.io/) and [CVX](https://cvxr.com/cvx/) are two
# packages for mathematical optimization in [MATLAB®](https://mathworks.com/products/matlab.html).
# They are independently developed and are in no way affiliated with JuMP.

# The purpose of this tutorial is to help new users to JuMP who have previously
# used YALMIP or CVX by comparing and contrasting their different features.

# !!! tip
#     If you have not used Julia before, read the [Getting started with Julia](@ref)
#     tutorial.

# ## Namespaces

# Julia has namespaces, which MATLAB lacks. Therefore one needs to either use
# the command:

using JuMP

# in order bring all names exported by JuMP into scope, or:

import JuMP

# in order to merely make the JuMP package available. `import` requires
# prefixing everything you use from JuMP with `JuMP.`. In this tutorial we use
# the former.

# ## Models

# YALMIP and CVX have a single, implicit optimization model that you build by
# defining variables and constraints.

# In JuMP, we create an explicit model first, and then, when you declare
# variables, constraints, or the objective function, you specify to which model
# they are being added.

# Create a new JuMP model with the command:

model = Model()

# ## Variables

# In most cases there is a direct translation between variable declarations.
# The following table shows some common examples:

# | JuMP                                              | YALMIP                                  | CVX                         |
# | :------------------------------------------------ | :-------------------------------------- | :-------------------------- |
# | `@variable(model, x)`                             | `x = sdpvar`                            | `variable x`                |
# | `@variable(model, x, Int)`                        | `x = intvar`                            | `variable x integer`        |
# | `@variable(model, x, Bin)`                        | `x = binvar`                            | `variable x binary`         |
# | `@variable(model, v[1:d])`                        | `v = sdpvar(d, 1)`                      | `variable v(d)`             |
# | `@variable(model, m[1:d, 1:d])`                   | `m = sdpvar(d,d,'full')`                | `variable m(d, d)`          |
# | `@variable(model, m[1:d, 1:d] in ComplexPlane())` | `m = sdpvar(d,d,'full','complex')`      | `variable m(d,d) complex`   |
# | `@variable(model, m[1:d, 1:d], Symmetric)`        | `m = sdpvar(d)`                         | `variable m(d,d) symmetric` |
# | `@variable(model, m[1:d, 1:d], Hermitian)`        | `m = sdpvar(d,d,'hermitian','complex')` | `variable m(d,d) hermitian` |

# Like CVX, but unlike YALMIP, JuMP can also constrain variables upon creation:

# | JuMP                                                  | CVX                                    |
# | :---------------------------------------------------- | :------------------------------------- |
# | `@variable(model, v[1:d] >= 0)`                       | `variable v(d) nonnegative`            |
# | `@variable(model, m[1:d, 1:d], PSD)`                  | `variable m(d,d) semidefinite`         |
# | `@variable(model, m[1:d, 1:d] in PSDCone())`          | `variable m(d,d) semidefinite`         |
# | `@variable(model, m[1:d, 1:d] in HermitianPSDCone())` | `variable m(d,d) complex semidefinite` |

# JuMP can additionally set variable bounds, which may be handled more
# efficiently by a solver than an equivalent linear constraint. For example:

@variable(model, -1 <= x[i in 1:3] <= i)
upper_bound.(x)

# A more interesting case is when you want to declare, for example, `n` real
# symmetric matrices. Both YALMIP and CVX allow you to put the matrices as the
# slices of a 3-dimensional array, via the commands `m = sdpvar(d, d, n)` and
# `variable m(d, d, n) symmetric`, respectively. With JuMP this is not possible.
# Instead, to achieve the same result one needs to declare a vector of `n`
# matrices:

d, n = 3, 2
m = [@variable(model, [1:d, 1:d], Symmetric) for _ in 1:n]

#-

m[1]

#-

m[2]

# The analogous construct in MATLAB would be a cell array containing the
# optimization variables, which every discerning programmer avoids as cell
# arrays are rather slow. This is not a problem in Julia: a vector of matrices
# is almost as fast as a 3-dimensional array.

# ## [Constraints](@id matlab_constraints)

# As in the case of variables, in most cases there is a direct translation
# between the packages:

# | JuMP                                                     | YALMIP               | CVX                                      |
# | :------------------------------------------------------- | :------------------- | :--------------------------------------- |
# | `@constraint(model, v == c)`                             | `v == c`             | `v == c`                                 |
# | `@constraint(model, v >= 0)`                             | `v >= 0`             | `v >= 0`                                 |
# | `@constraint(model, m >= 0, PSDCone())`                  | `m >= 0`             | `m == semidefinite(length(m))`           |
# | `@constraint(model, m >= 0, HermitianPSDCone())`         | `m >= 0`             | `m == hermitian_semidefinite(length(m))` |
# | `@constraint(model, [t; v] in SecondOrderCone())`        | `cone(v, t)`         | `{v, t} == lorentz(length(v))`           |
# | `@constraint(model, [x, y, z] in MOI.ExponentialCone())` | `expcone([x, y, z])` | `{x, y, z} == exponential(1)`            |

# A subtlety appears when declaring equality constraints for matrices. In
# general, JuMP uses `@constraint(model, m .== c)`, with the dot meaning
# broadcasting in Julia, except when `m` is `Symmetric` or `Hermitian`: in this
# case `@constraint(model, m == c)` is allowed, and is much better, as JuMP is
# smart enough to not generate redundant constraints for the lower diagonal and
# the imaginary part of the diagonal (in the complex case). Both YALMIP and CVX
# are also smart enough to do this and the syntax is always just `m == c`.

# Experienced YALMIP users will probably be relieved to see that you must pass
# `PSDCone()` or `HermitianPSDCone()` to make a matrix positive semidefinite, as
# the `>=` ambiguity in YALMIP is common source of bugs.

# ## Setting the objective

# Like CVX, but unlike YALMIP, JuMP has a specific command for setting an
# objective function:

@objective(model, Min, sum(i * x[i] for i in 1:3))

# Here the third argument is any expression you want to optimize, and `Min` is
# an objective sense (the other possibility is `Max`).

# ## Setting solver and options

#  In order to set an optimizer with JuMP, do:

import Clarabel
set_optimizer(model, Clarabel.Optimizer)

# where "Clarabel" is an example solver. See the list of [Supported solvers](@ref)
# for other choices.

# To configure the solver options you use the command:

set_attribute(model, "verbose", true)

# where `verbose` is an option specific to Clarabel.

# A crucial difference is that with JuMP you must explicitly choose a solver
# before optimizing. Both YALMIP and CVX allow you to leave it empty and will
# try to guess an appropriate solver for the problem.

# ## Optimizing

# Like YALMIP, but unlike CVX, with JuMP you need to explicitly start the
# optimization, with the command:

optimize!(model)

# The exclamation mark here is a Julia-ism that means the function is modifying
# its argument, `model`.

# ## Querying solution status

# After the optimization is done, you should check for the solution status to
# see what solution (if any) the solver found.

# Like YALMIP and CVX, JuMP provides a solver-independent way to check it, via
# the command:

is_solved_and_feasible(model)

# If the return value is `false`, you should investigate with [`termination_status`](@ref),
# [`primal_status`](@ref), and [`raw_status`](@ref), See [Solutions](@ref jump_solutions)
# for more details on how to query and interpret solution statuses.

# ## Extracting variables

# Like YALMIP, but unlike CVX, with JuMP you need to explicitly ask for the value
# of your variables after optimization is done, with the function call `value(x)`
# to obtain the value of variable `x`.

value.(m[1][1, 1])

# A subtlety is that, unlike YALMIP, the function `value` is only defined for
# scalars. For vectors and matrices you need to use Julia broadcasting:
# `value.(v)`.

value.(m[1])

# There is also a specialized function for extracting the value of the objective,
# `objective_value(model)`, which is useful if your objective doesn't have a
# convenient expression.

objective_value(model)

# ## Dual variables

# Like YALMIP and CVX, JuMP allows you to recover the dual variables. In order
# to do that, the simplest method is to name the constraint you're interested in,
# for example, `@constraint(model, bob, sum(v) == 1)` and then, after the
# optimzation is done, call `dual(bob)`. See [Duality](@ref) for more details.

# ## Reformulating problems

# Perhaps the biggest difference between JuMP and YALMIP and CVX is how far the
# package is willing to go in reformulating the problems you give to it.

# CVX is happy to reformulate anything it can, even using approximations if your
# solver cannot handle the problem.

# YALMIP will only do exact reformulations, but is still fairly adventurous,
# for example, being willing to reformulate a nonlinear objective in terms of
# conic constraints.

# JuMP does no such thing: it only reformulates objectives into objectives, and
# constraints into constraints, and is fairly conservative at that. As a result,
# you might need to do some reformulations manually, for which a good guide is
# the [Tips and tricks](@ref conic_tips_and_tricks) tutorial.

# ## Vectorization

# In MATLAB, it is absolutely essential to "vectorize" your code to obtain
# acceptable performance. This is because MATLAB is a slow interpreted
# language, which sends your commands to fast libraries. When you "vectorize"
# your code you are minimizing the MATLAB part of the work and sending it to the
# fast libraries instead.

# There's no such duality with Julia.

# Everything you write and most libraries you use will compile down to LLVM, so
# "vectorization" has no effect.

# For example, if you are writing a linear program in MATLAB and instead of the
# usual `constraints = [v >= 0]` you write:
# ```matlab
# for i = 1:n
#    constraints = [constraints, v(i) >= 0];
# end
# ```
# performance will be poor.

# With Julia, on the other hand, there is hardly any difference between
# ```julia
# @constraint(model, v >= 0)
# ```
#  and
# ```julia
# for i in 1:n
#     @constraint(model, v[i] >= 0)
# end
# ```

# ## Symmetric and Hermitian matrices

# Julia has specialized support for symmetric and Hermitian matrices in the
# `LinearAlgebra` package:

import LinearAlgebra

# If you have a matrix that is numerically symmetric:

x = [1 2; 2 3]

#-

LinearAlgebra.issymmetric(x)

# then you can wrap it in a `LinearAlgebra.Symmetric` matrix to tell Julia's
# type system that the matrix is symmetric.

LinearAlgebra.Symmetric(x)

# Using a `Symmetric` matrix lets Julia and JuMP use more efficient algorithms
# when they are working with symmetric matrices.

# If you have a matrix that is nearly but not exactly symmetric:

x = [1.0 2.0; 2.001 3.0]
LinearAlgebra.issymmetric(x)

# then you could, as you might do in MATLAB, make it numerically symmetric as
# follows:

x_sym = 0.5 * (x + x')

# In Julia, you can explicitly choose whether to use the lower or upper triangle
# of the matrix:

x_sym = LinearAlgebra.Symmetric(x, :L)

#-

x_sym = LinearAlgebra.Symmetric(x, :U)

# The same applies for Hermitian matrices, using `LinearAlgebra.Hermitian` and
# `LinearAlgebra.ishermitian`.

# ## Primal versus dual form

# When you translate some optimization problems from YALMIP or CVX to JuMP, you
# might be surprised to see it get much faster or much slower, even if you're
# using exactly the same solver. The most likely reason is that YALMIP will
# always interpret the problem as the dual form, whereas CVX and JuMP will try to
# interpret the problem in the form most appropriate to the solver. If the
# problem is more naturally formulated in the primal form it is likely that
# YALMIP's performance will suffer, or if JuMP gets it wrong, its performance will
# suffer. It might be worth trying both primal and dual forms if you're having
# trouble, which can be done automatically with the package [Dualization.jl](@ref).
#
# For an in-depth explanation of this issue, see the [Dualization](@ref) tutorial.
# ## Rosetta stone

# In this section, we show a complete example of the same optimization problem
# being solved with JuMP, YALMIP, and CVX. It is a semidefinite program that
# computes a lower bound on the random robustness of entanglement using the
# partial transposition criterion.

# The code is complete, apart from the function that does partial transposition.
# With both YALMIP and CVX we use the function `PartialTranspose` from
# [QETLAB](https://github.com/nathanieljohnston/QETLAB). With JuMP, we could use
# the function `Convex.partialtranspose` from [Convex.jl](https://jump.dev/Convex.jl/stable/),
# but we reproduce it here for simplicity:

function partial_transpose(x::AbstractMatrix, sys::Int, dims::Vector)
    @assert size(x, 1) == size(x, 2) == prod(dims)
    @assert 1 <= sys <= length(dims)
    n = length(dims)
    s = n - sys + 1
    p = collect(1:2n)
    p[s], p[n+s] = n + s, s
    r = reshape(x, (reverse(dims)..., reverse(dims)...))
    return reshape(permutedims(r, p), size(x))
end

# ### JuMP

# The JuMP code to solve this problem is:

using JuMP
import Clarabel
import LinearAlgebra

function random_state_pure(d)
    x = randn(Complex{Float64}, d)
    y = x * x'
    return LinearAlgebra.Hermitian(y / LinearAlgebra.tr(y))
end

function robustness_jump(d)
    rho = random_state_pure(d^2)
    id = LinearAlgebra.Hermitian(LinearAlgebra.I(d^2))
    rhoT = LinearAlgebra.Hermitian(partial_transpose(rho, 1, [d, d]))
    model = Model()
    @variable(model, λ)
    @constraint(model, PPT, rhoT + λ * id in HermitianPSDCone())
    @objective(model, Min, λ)
    set_optimizer(model, Clarabel.Optimizer)
    set_attribute(model, "verbose", true)
    optimize!(model)
    if is_solved_and_feasible(model)
        WT = dual(PPT)
        return value(λ), real(LinearAlgebra.dot(WT, rhoT))
    else
        return "Something went wrong: $(raw_status(model))"
    end
end

robustness_jump(3)

# ### YALMIP

# The corresponding YALMIP code is:

# ```matlab
# function robustness_yalmip(d)
#     rho = random_state_pure(d^2);
#     # PartialTranspose from https://github.com/nathanieljohnston/QETLAB
#     rhoT = PartialTranspose(rho, 1, [d d]);
#     lambda = sdpvar;
#     constraints = [(rhoT + lambda*eye(d^2) >= 0):'PPT'];
#     ops = sdpsettings(sdpsettings, 'verbose', 1, 'solver', 'sedumi');
#     sol = optimize(constraints, lambda, ops);
#     if sol.problem == 0
#         WT = dual(constraints('PPT'));
#         value(lambda)
#         real(WT(:).' * rhoT(:))
#     else
#         display(['Something went wrong: ', sol.info])
#     end
# end
#
# function rho = random_state_pure(d)
#     x = randn(d, 1) + 1i * randn(d, 1);
#     y = x * x';
#     rho = y / trace(y);
# end
# ```

# ### CVX

# The corresponding CVX code is:

# ```matlab
# function robustness_cvx(d)
#     rho = random_state_pure(d^2);
#     # PartialTranspose from https://github.com/nathanieljohnston/QETLAB
#     rhoT = PartialTranspose(rho, 1, [d d]);
#     cvx_begin
#         variable lambda
#         dual variable WT
#         WT : rhoT + lambda * eye(d^2) == hermitian_semidefinite(d^2)
#         minimise lambda
#     cvx_end
#     if strcmp(cvx_status, 'Solved')
#         lambda
#         real(WT(:)' * rhoT(:))
#     else
#         display('Something went wrong.')
#     end
# end
#
# function rho = random_state_pure(d)
#     x = randn(d, 1) + 1i * randn(d, 1);
#     y = x * x';
#     rho = y / trace(y);
# end
# ```