# tests/testthat/test-g02-examples.R In anqif/cvxr: Disciplined Convex Optimization

```context("test-g02-examples")
TOL <- 1e-6

ECOS.dims_to_solver_dict <- CVXR:::ECOS.dims_to_solver_dict

test_that("Find the largest Euclidean ball in the polyhedron", {
skip_on_cran()
## Create the input data
a1 <- matrix(c(2,1), nrow = 2, ncol = 1)
a2 <- matrix(c(2,-1), nrow = 2, ncol = 1)
a3 <- matrix(c(-1,2), nrow = 2, ncol = 1)
a4 <- matrix(c(-1,-2), nrow = 2, ncol = 1)
b <- rep(1,4)

## Create and solve the model
r <- Variable(name = "r")
x_c <- Variable(2, name = "x_c")
obj <- Maximize(r)
constraints <- list(
t(a1) %*% x_c + norm(a1,"F")*r <= b[1],
t(a2) %*% x_c + norm(a2,"F")*r <= b[2],
t(a3) %*% x_c + norm(a3,"F")*r <= b[3],
t(a4) %*% x_c + norm(a4,"F")*r <= b[4]
)

p <- Problem(obj, constraints)
result <- solve(p)
expect_equal(result\$value, 0.447214, tolerance = TOL)
expect_equal(result\$getValue(r), result\$value, tolerance = TOL)
expect_equal(result\$getValue(x_c), matrix(c(0,0)), tolerance = TOL)
})

test_that("Test issue with scalars", {
skip_on_cran()
n <- 6
eps <- 1e-6
set.seed(10)
P0 <- matrix(rnorm(n*n), nrow = n, ncol = n)
eye <- diag(n)
P0 <- t(P0) %*% P0 + eps*eye

print(P0)

P1 <- matrix(rnorm(n*n), nrow = n, ncol = n)
P1 <- t(P1) %*% P1
P2 <- matrix(rnorm(n*n), nrow = n, ncol = n)
P2 <- t(P2) %*% P2
P3 <- matrix(rnorm(n*n), nrow = n, ncol = n)
P3 <- t(P3) %*% P3

q0 <- matrix(rnorm(n), nrow = n)
q1 <- matrix(rnorm(n), nrow = n)
q2 <- matrix(rnorm(n), nrow = n)
q3 <- matrix(rnorm(n), nrow = n)

r0 <- matrix(rnorm(1))
r1 <- matrix(rnorm(1))
r2 <- matrix(rnorm(1))
r3 <- matrix(rnorm(1))

slack <- Variable()
# Form the problem
x <- Variable(n)
objective <- Minimize(0.5*quad_form(x, P0) + t(q0) %*% x + r0 + slack)
constraints <- list(0.5*quad_form(x, P1) + t(q1) %*% x + r1 <= slack,
0.5*quad_form(x, P2) + t(q2) %*% x + r2 <= slack,
0.5*quad_form(x, P3) + t(q3) %*% x + r3 <= slack)

# We now find the primal result and compare it to the dual result to check if strong duality holds, i.e., the duality gap is effectively zero
p <- Problem(objective, constraints)
result <- solve(p)

# Note that since our data is random, we may need to run this program multiple times to get a feasible primal.
# When feasible, we can print out the following values.
print(result\$getValue(x))   # Solution
lam1 <- result\$getDualValue(constraints[[1]])
lam2 <- result\$getDualValue(constraints[[2]])
lam3 <- result\$getDualValue(constraints[[3]])
print(class(lam1))

P_lam <- P0 + lam1*P1 + lam2*P2 + lam3*P3
q_lam <- q0 + lam1*q1 + lam2*q2 + lam3*q3
r_lam <- r0 + lam1*r1 + lam2*r2 + lam3*r3
dual_result <- -0.5*t(q_lam) %*% P_lam %*% q_lam + r_lam
print(dim(dual_result))
expect_equal(dim(dual_result), c(1,1))
})

## Solving a problem with different solvers
skip_on_cran()
x <- Variable(2)
obj <- Minimize(x[1] + norm1(x))
constraints <- list(x >= 2)
prob <- Problem(obj, constraints)

## Solve with ECOS
result <- solve(prob, solver = "ECOS")
print(paste("optimal value with ECOS:", result\$value))
expect_equal(result\$value, 6, tolerance = TOL)

## Solve with ECOS_BB
result <- solve(prob, solver = "ECOS_BB")
print(paste("optimal value with ECOS_BB:", result\$value))
expect_equal(result\$value, 6, tolerance = TOL)

## Solve with CVXOPT
if("CVXOPT" %in% installed_solvers()) {
result <- solve(prob, solver = "CVXOPT")
print(paste("optimal value with CVXOPT:", result\$value))
expect_equal(result\$value, 6, tolerance = TOL)
}

## Solve with SCS
result <- solve(prob, solver = "SCS")
print(paste("optimal value with SCS:", result\$value))
expect_equal(result\$value, 6, tolerance = 1e-2)

## Solve with CPLEX
if("CPLEX" %in% installed_solvers()) {
result <- solve(prob, solver = "CPLEX")
print(paste("optimal value with CPLEX:", result\$value))
expect_equal(result\$value, 6, tolerance = TOL)
}

## Solve with GLPK
if("GLPK" %in% installed_solvers()) {
result <- solve(prob, solver = "GLPK")
print(paste("optimal value with GLPK:", result\$value))
expect_equal(result\$value, 6, tolerance = TOL)

result <- solve(prob, solver = "GLPK_MI")
print(paste("optimal value with GLPK_MI:", result\$value))
expect_equal(result\$value, 6, tolerance = TOL)
}

## Solve with GUROBI
if("GUROBI" %in% installed_solvers()) {
result <- solve(prob, solver = "GUROBI")
print(paste("optimal value with GUROBI:", result\$value))
expect_equal(result\$value, 6, tolerance = TOL)
}

print(installed_solvers())
})

test_that("Test log determinant", {
skip_on_cran()
x <- t(data.frame(c(0.55, 0.25, -0.2, -0.25, -0.0, 0.4),
c(0.0, 0.35, 0.2, -0.1, -0.3, -0.2)))
n <- nrow(x)
m <- ncol(x)

## Create and solve the model
A <- Variable(n, n)
b <- Variable(n)
obj <- Maximize(log_det(A))
constraints <- lapply(1:m, function(i) { p_norm(A %*% as.matrix(x[,i]) + b) <= 1 })
p <- Problem(obj, constraints)
result <- solve(p)
expect_equal(result\$value, 1.9746, tolerance = 1e-2)
})

test_that("Test portfolio problem", {
skip_on_cran()
set.seed(5)
n <- 100
m <- 10
pbar <- matrix(1, nrow = n, ncol = 1) * 0.03 +
c(matrix(stats::rnorm(n-1), nrow = n-1, ncol = 1), 0) * 0.12

Fmat <- Matrix::rsparsematrix(m, n, density = 0.01)
Fmat@x <- rep(1, length(Fmat@x))
D <- Matrix::sparseMatrix(1:n, 1:n, x = rep(1,n))
D@x <- stats::rnorm(length(D@x))^2
Z <- matrix(stats::rnorm(m), nrow = m, ncol = 1)
Z <- Z %*% t(Z)

x <- Variable(n)
y <- Fmat %*% x
mu <- 1
ret <- t(pbar) %*% x
## DCP attr causes error because not all the curvature
## matrices are reduced to constants when an atom is scalar.
risk <- p_norm(D %*% x)^2 + (Z %*% y)^2
})

test_that("Test examples from CVXR introduction", {
skip_on_cran()
m <- 30
n <- 20
set.seed(1)
A <- matrix(stats::rnorm(m*n), nrow = m, ncol = n)
b <- matrix(stats::rnorm(m), nrow = m, ncol = 1)

## Construct the problem
x <- Variable(n)
objective <- Minimize(sum_squares(A %*% x - b))
constraints <- list(0 <= x, x <= 1)
prob <- Problem(objective, constraints)

## The optimal objective is returned by solve(prob)
result <- solve(prob)
## The optimal value for x
print(result\$getValue(x))
## The optimal Lagrange multiplier for a constraint is stored in result\$getDualValue(constraints[[1]])
print(result\$getDualValue(constraints[[1]]))

###########################################
## Create two scalar variables
x <- Variable()
y <- Variable()

## Create two constraints
constraints <- list(x + y == 1, x - y >= 1)

## Form objective
obj <- Minimize((x - y)^2)

## Form and solve problem
prob <- Problem(obj, constraints)
result <- solve(prob)
cat("status:", result\$status)
cat("\noptimal value:", result\$value)
cat("\noptimal var:", result\$getValue(x), ", ", result\$getValue(y))

expect_equal(tolower(result\$status), "optimal")
expect_equal(result\$value, 1.0, tolerance = TOL)
expect_equal(result\$getValue(x), 1.0, tolerance = TOL)
expect_equal(result\$getValue(y), 0, tolerance = TOL)

###########################################
## Replace the objective
prob@objective <- Maximize(x + y)
result <- solve(prob)
cat("optimal value:", result\$value)

expect_equal(result\$value, 1.0, tolerance = TOL)

## Replace the constraints (x + y == 1)
prob@constraints[[1]] <- (x + y <= 3)
result <- solve(prob)
cat("optimal value:", result\$value)

expect_equal(result\$value, 3.0, tolerance = TOL)

###########################################
x <- Variable()

## An infeasible problem
prob <- Problem(Minimize(x), list(x >= 1, x <= 0))
# TODO_NARAS_11: result <- solve(prob)
result <- solve(prob, solver = "ECOS")
cat("status:", result\$status)
cat("optimal value:", result\$value)

expect_equal(result\$status, "infeasible")
expect_equal(result\$value, Inf)

## An unbounded problem
prob <- Problem(Minimize(x))
result <- solve(prob)
cat("status:", result\$status)
cat("optimal value:", result\$value)

expect_equal(tolower(result\$status), "unbounded")
expect_equal(result\$value, -Inf)

###########################################
## A scalar variable
a <- Variable()

## Column vector variable of length 5
x <- Variable(5)

## Matrix variable with 4 rows and 7 columns
A <- Variable(4, 7)

###########################################
m <- 10
n <- 5
set.seed(1)
A <- matrix(stats::rnorm(m*n), nrow = m, ncol = n)
b <- matrix(stats::rnorm(m), nrow = m, ncol = 1)

## Construct the problem
x <- Variable(n)
objective <- Minimize(sum((A %*% x - b)^2))
constraints <- list(0 <= x, x <= 1)
prob <- Problem(objective, constraints)

result <- solve(prob)
cat("Optimal value:", result\$value)
cat("Optimal var:", result\$getValue(x))

expect_equal(result\$value, 7.244277, tolerance = TOL)

###########################################
## Positive scalar parameter
m <- Parameter(nonneg = TRUE)

## Column vector parameter with unknown sign (by default)
c <- Parameter(5)

## Matrix parameter with negative entries
G <- Parameter(4, 7, nonpos = TRUE)

## Assigns a constant value to G
value(G) <- -matrix(1, nrow = 4, ncol = 7)

###########################################
## Create parameter, then assign value
rho <- Parameter(nonneg = TRUE)
value(rho) <- 2

## Initialize parameter with a value
rho <- Parameter(nonneg = TRUE, value = 2)

###########################################
## Problem data
n <- 15
m <- 10
set.seed(1)
A <- matrix(stats::rnorm(n*m), nrow = n, ncol = m)
b <- matrix(stats::rnorm(n), nrow = n, ncol = 1)
## gamma must be positive due to DCP rules
gamma <- Parameter(nonneg = TRUE)

## Construct the problem
x <- Variable(m)
error <- sum_squares(A %*% x - b)

## Construct a trade-off curve of ||Ax-b||^2 vs. ||x||_1
sq_penalty <- c()
l1_penalty <- c()
x_values <- list()
gammas <- 10^seq(-4, 6, length.out = 50)
for(gamma in gammas) {
obj <- Minimize(error + gamma*p_norm(x,1))
prob <- Problem(obj)

result <- solve(prob)
sq_penalty <- c(sq_penalty, result\$getValue(error))
l1_penalty <- c(l1_penalty, result\$getValue(norm1(x)))
x_values <- c(x_values, result\$getValue(x))
}

###########################################
X <- Variable(5, 4)
A <- matrix(1, nrow = 3, ncol = 5)

## Use size(expr) to get the dimensions
cat("dimensions of X:", size(X))
cat("\ndimensions of sum_entries(X):", size(sum_entries(X)))
cat("\ndimensions of A %*% X:", size(A %*% X))

## ValueError raised for invalid dimensions
expect_error(A + X)
})

skip_on_cran()
x <- Variable()
prob <- Problem(Minimize(x^2), list(x == 2))

# Get ECOS arguments.
tmp <- get_problem_data(prob, "ECOS")

# Get ECOS_BB arguments.
tmp <- get_problem_data(prob, "ECOS_BB")

# Get CVXOPT arguments.
if("CVXOPT" %in% installed_solvers())
tmp <- get_problem_data(prob, "CVXOPT")

# Get SCS arguments.
tmp <- get_problem_data(prob, "SCS")

# Get ECOS arguments.
library(ECOSolveR)
tmp <- get_problem_data(prob, "ECOS")
data <- tmp[[1]]
chain <- tmp[[2]]
inverse <- tmp[[3]]

# Call ECOS solver.
solution <- ECOSolveR::ECOS_csolve(c = data\$c, G = data\$G, h = data\$h, dims = ECOS.dims_to_solver_dict(data\$dims), A = data\$A, b = data\$b)

# Unpack raw solver output.
unpack_results(prob, solution, chain, inverse)
})
```
anqif/cvxr documentation built on April 1, 2020, 10:30 a.m.