# Class V - Fun visualizations in JuMP

This notebook contains a number of fun JuMP examples.


```julia
import Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()
```

## Example 1: minimum bounding ellipse

Given a set of ellipses centered on the origin

$E(A) = \{ u\;|\;u^\top A^{-1} u <= 1 \}$

find a "minimal" ellipse that contains the provided ellipses

We can formulate this as an SDP:

$\begin{align}
     minimize \quad& trace(WX)\\
   subject to \quad& X \ge A_i, \quad i = 1,...,m \\
              &  X \succeq 0
\end{align}$

where $W$ is a positive-definite matrix of weights to choose between different solutions.


```julia
using JuMP, SCS, Plots, LinearAlgebra, Interact

function draw_ellipse(A::Matrix, args...; kwargs...)
    x_values = Float64[]
    y_values = Float64[]
    for angle in 0:0.001π:2π
        u = [cos(angle), sin(angle)]
        z = A * u
        push!(x_values, z[1])
        push!(y_values, z[2])
    end
    plot!(x_values, y_values, args...; kwargs...)
end

function solve_minimum_ellipse_problem(W, A_matrices)
    model = Model(with_optimizer(SCS.Optimizer, eps = 1e-6, verbose = false))
    @variable(model, X[1:2, 1:2], PSD)
    @objective(model, Min, tr(W * X))
    for A in A_matrices
        @SDconstraint(model, X >= A)
    end
    optimize!(model)
    return termination_status(model), value.(X)
end
```

Investigate the model. Here are some things to try:

- What happens if you comment out the first A matrix?

- What happens if you comment out the second A matrix?

- What happens if you comment out the third A matrix?

You can comment lines in Julia using the `#` symbol. As a shortcut, use `[CTRL] + [/]`.


```julia
@manipulate for weight in 1:20
    A_matrices = [
        [2.0 0.0; 0.0 1.0],
        [1.0 0.0; 0.0 3.0],
        [2.3896 1.5433; 1.5433 1.35584]
    ]
    W = [1.0 0.0; 0.0 weight]
    status, X_value = solve_minimum_ellipse_problem(W, A_matrices)
    if status == MOI.OPTIMAL
        plot(legend = false)
        draw_ellipse.(A_matrices, color = "gray")
        draw_ellipse(X_value, color="purple", linewidth=2)
    else
        println("Could not solve. Termination status = $(status)")
    end
end
```

# Example 2: clustering points

See Section 2.3 of:
http://www.optimization-online.org/DB_FILE/2005/04/1114.pdf

Given 2-D data pairs $d_i$, $i=1,\ldots,N$, these points can be partitioned into $k$ clusters by solving the following SDP.

$\begin{align}
     minimize \quad& trace(W * (\mathbf{I} - X))\\
   subject\ to \quad& \Sigma_j X_{ij} = 1,\quad i=1, \ldots,N \\
              & trace(X) = k \\
              &  X \succeq 0
\end{align}$,

where $W_{ij} = e^{\frac{-||d_i - d_j||}{\sigma}}$.


```julia
""""
    calculate_weight_matrix(data::Matrix{Float64})

Calculates the distance between a list of 2-D data points given as
rows in `data` matrix.
"""
function calculate_weight_matrix(data::Matrix{Float64}, σ = 1.0)
    num_points = size(data, 1)
    W = zeros(num_points, num_points)
    for i in 1:num_points
        for j in (i+1):num_points
            dist = exp(-norm(data[i, :] - data[j, :]) / σ)
            W[i, j] = dist
            W[j, i] = dist
        end
    end
    return W
end

function eye(x, y)
    z = zeros(x, y)
    for i in 1:min(x, y)
        z[i, i] = 1.0
    end
    return z
end

function solve_cluster_problem(data::Matrix{Float64}, num_clusters::Int)
    W = calculate_weight_matrix(data)
    num_points = size(data, 1)
    
    model = Model(with_optimizer(SCS.Optimizer, verbose = false))

    @variable(model, X[1:num_points, 1:num_points], PSD)
    @objective(model, Min, tr(W * (eye(num_points, num_points) - X)))
    @constraints(model, begin
        X .>= 0
        [i in 1:num_points], sum(X[i, :]) == 1
        tr(X) == num_clusters
    end)

    optimize!(model)

    X_value = value.(X)

    cluster = zeros(Int, num_points)
    cluster_counter = 0
    for i in 1:num_points
        if cluster[i] == 0
            cluster_counter += 1
            cluster[i] = cluster_counter
            for j in (i+1):num_points
                if norm(X_value[i, j] - X_value[i, i]) <= 1e-6
                    cluster[j] = cluster[i]
                end
            end
        end
    end
    return cluster
end
```

Investigate the model. What goes wrong when?


```julia
data = vcat(
    rand(Float64, (10, 2)) .+ [2 3],
    rand(Float64, (10, 2)) .+ [4 6],
    rand(Float64, (10, 2)) .+ [4 3]
)

@manipulate for num_clusters = 1:4
    which_clusters = solve_cluster_problem(data, num_clusters)
    plot(
        xlabel = "x", ylabel = "y", 
        xlims=(0, 8), ylims = (0, 8),
        legend = :bottomright
    )
    for k in 1:maximum(which_clusters)
        points = which_clusters .== k
        scatter!(data[points, 1], data[points, 2], 
            label = "Cluster $(k)", markersize=10)
    end
    plot!()
end
```

## Example 3 - Ridge regression

In this example, we consider fitting a quadratic function to a noisy signal using [ridge regression](https://en.wikipedia.org/wiki/Tikhonov_regularization).

$\begin{align}
      minimize \quad& ||z - y||_2 + \lambda||\theta||_2\\
   subject\ to \quad& z_i = \theta_1 x_i^2 + \theta_2 x_i + \theta_3
\end{align}$,


```julia
using Ipopt

function solve_ridge_regression(x, y, λ)            
    N = length(x)
    model = Model(with_optimizer(Ipopt.Optimizer, print_level = 0))
    @variable(model, θ[i = 1:4], start = 1.0)
    @variable(model, z[i = 1:N], start = x[i]^2)
    @constraint(model, [i = 1:N], z[i] == θ[1] + θ[2] * x[i] + θ[3] * x[i]^2+ θ[4] * x[i]^3)
    @NLobjective(model, Min, 
        sum((z[i] - y[i])^2 for i in 1:N) + λ * sum(θ[i]^2 for i in 1:4)
    )
    optimize!(model)
    return value.(z)
end
```

Play round with the weight in front of the regularizer term. What changes? What happens if you increase the number of data points?


```julia
N = 10
x = rand(N)
sort!(x)
y = 5 * x.^2 .+ x
yϵ = y .+ (1.5 * rand(N) .- 0.5)

@manipulate for λ in 0:0.01:0.25
    z = solve_ridge_regression(x, yϵ, λ)
    plot(x, y, linestyle=:dash, color="purple")
    scatter!(x, yϵ, 
        legend=false, color="gray", markerstrokewidth=0
    )
    plot!(x, z, linewidth=2, color="purple")
end
```

# Example 4 - feasible regions

This example explores a problem with a nonlinear feasible region.


```julia
@manipulate for λ in 0.05:0.05:0.95
    model = Model(with_optimizer(Ipopt.Optimizer, print_level = 0))
    @variable(model, 0 <= x <= 20)
    @variable(model, 0 <= y <= 20)
    @NLconstraint(model, y >= 0.05x^2)
    @NLconstraint(model, y >= exp(-0.25 * (x - 10)))
    @objective(model, Min, λ * x + (1 - λ) * y)
    optimize!(model)

    b = 1
    a = b * (λ - 1) / λ
    ab = a^2 + b^2
    b /= ab
    a /= ab

    x_value = value(x)
    y_value = value(y)

    z = -5:0.1:25.0
    plot(xlims=(-2.5, 10), ylims=(0, 15), legend=false)
    vline!([0.0], color="gray", linewidth=2)
    plot!(z, 0.05z.^2, color="gray", linewidth=2)
    plot!(z, exp.(-0.25(z.-10)), color="gray", linewidth=2)
    scatter!([x_value], [y_value], color="purple", markersize=10)
    plot!(z, y_value .+ (b / a) * (z .- x_value), linestyle=:dash, color="purple")
end
```

6c4fbe07fc75dfdcc730564fcbd5e72df225a84b
ipynb
Jupyter Notebook
Class II.B - Fun visualizations in JuMP.ipynb
lanl-ansi/tutorial-grid-science-2019
75ea7e77e3819ecf5e8cc1a0ec24465d458e7cfe
MIT
2019-01-07T08:31:59.000Z
2020-03-04T16:42:04.000Z
Class II.B - Fun visualizations in JuMP.ipynb
lanl-ansi/tutorial-grid-science
75ea7e77e3819ecf5e8cc1a0ec24465d458e7cfe
MIT
2019-01-06T15:37:33.000Z
2019-06-22T15:33:18.000Z
Class II.B - Fun visualizations in JuMP.ipynb
lanl-ansi/tutorial-grid-science-2019
75ea7e77e3819ecf5e8cc1a0ec24465d458e7cfe
MIT
2019-01-08T23:31:36.000Z
2020-05-01T22:51:59.000Z
Qwen/Qwen-72B
1. YES
2. YES