# tests/testthat/test-g01-dgp2dcp.R In anqif/cvxr: Disciplined Convex Optimization

```context("test-g01-dgp2dcp")
TOL <- 1e-6

Solution <- CVXR:::Solution
perform <- CVXR:::perform
reduce <- CVXR:::reduce
invert <- CVXR:::invert
retrieve <- CVXR:::retrieve
unpack_problem <- CVXR:::unpack_problem
cvxr_expr <- CVXR:::expr

Dgp2Dcp <- function(problem = NULL) { new("Dgp2Dcp", problem = problem) }
Dgp2Dcp.mulexpression_canon <- CVXR:::Dgp2Dcp.mulexpression_canon
Dgp2Dcp.prod_canon <- CVXR:::Dgp2Dcp.prod_canon
Dgp2Dcp.trace_canon <- CVXR:::Dgp2Dcp.trace_canon

form_solution <- function(prob, result) {
primal <- list()
for(var in prob@variables)
primal[[as.character(var@id)]] <- result\$getValue(var)

dual <- list()
for(constr in prob@constraints)
dual[[as.character(constr@id)]] <- result\$getDualValue(constr)

attr <- list(solve_time = result\$solve_time, setup_time = result\$setup_time, num_iters = result\$num_iters)
Solution(result\$status, result\$value, primal, dual, attr)
}

test_that("test unconstrained monomial", {
skip_on_cran()
x <- Variable(pos = TRUE)
y <- Variable(pos = TRUE)
prod <- x*y
dgp <- Problem(Minimize(prod))
dgp2dcp <- Dgp2Dcp(dgp)

tmp <- reduce(dgp2dcp)
dgp2dcp <- tmp[[1]]
dcp <- tmp[[2]]

expect_equal(length(cvxr_expr(dcp@objective)@args), 2)
expect_equal(class(cvxr_expr(dcp@objective)@args[[1]])[1], "Variable")
expect_equal(class(cvxr_expr(dcp@objective)@args[[2]])[1], "Variable")
opt <- solve(dcp)

# dcp is solved in log-space, so it is unbounded below
# (since the OPT for dgp is 0 + epsilon).
expect_equal(opt\$value, -Inf)
expect_equal(opt\$status, "unbounded")

solution <- form_solution(dcp, opt)
dgp_unpack <- unpack_problem(dgp, retrieve(dgp2dcp, solution))
expect_equal(dgp_unpack\$value, 0.0)
expect_equal(dgp_unpack\$status, "unbounded")
opt <- solve(dgp, gp = TRUE)
expect_equal(opt\$value, 0.0)
expect_equal(opt\$status, "unbounded")

dgp <- Problem(Maximize(prod))
dgp2dcp <- Dgp2Dcp(dgp)
tmp <- reduce(dgp2dcp)
dgp2dcp <- tmp[[1]]
dcp <- tmp[[2]]
opt <- solve(dcp)
expect_equal(opt\$value, Inf)
expect_equal(opt\$status, "unbounded")

solution <- form_solution(dcp, opt)
dgp_unpack <- unpack_problem(dgp, retrieve(dgp2dcp, solution))
expect_equal(dgp_unpack\$value, Inf)
expect_equal(dgp_unpack\$status, "unbounded")
opt <- solve(dgp, gp = TRUE)
expect_equal(opt\$value, Inf)
expect_equal(opt\$status, "unbounded")
})

test_that("test basic equality constraint", {
skip_on_cran()
x <- Variable(pos = TRUE)
dgp <- Problem(Minimize(x), list(x == 1.0))
dgp2dcp <- Dgp2Dcp(dgp)

tmp <- reduce(dgp2dcp)
dgp2dcp <- tmp[[1]]
dcp <- tmp[[2]]
expect_equal(class(dcp@objective@args[[1]])[1], "Variable")

opt <- solve(dcp)
expect_equal(opt\$value, 0.0, tolerance = TOL)
expect_equal(opt\$getValue(variables(dcp)[[1]]), 0.0, tolerance = TOL)

solution <- form_solution(dcp, opt)
dgp_unpack <- unpack_problem(dgp, retrieve(dgp2dcp, solution))
expect_equal(dgp_unpack\$value, 1.0, tolerance = TOL)
expect_equal(dgp_unpack\$getValue(x), 1.0, tolerance = TOL)
result <- solve(dgp, gp = TRUE)
expect_equal(result\$value, 1.0, tolerance = TOL)
expect_equal(result\$getValue(x), 1.0, tolerance = TOL)
})

test_that("test basic GP", {
skip_on_cran()
x <- Variable(pos = TRUE)
y <- Variable(pos = TRUE)
z <- Variable(pos = TRUE)
constraints <- list(2*x*y + 2*x*z + 2*y*z <= 1.0, x >= 2*y)
problem <- Problem(Minimize(1/(x*y*z)), constraints)
result <- solve(problem, gp = TRUE)
expect_equal(15.59, result\$value, tolerance = 1e-2)
})

test_that("test max_elemwise", {
skip_on_cran()
x <- Variable(pos = TRUE)
y <- Variable(pos = TRUE)
z <- Variable(pos = TRUE)

prod1 <- x*y^(0.5)
prod2 <- 3.0*x*y^(0.5)
obj <- Minimize(max_elemwise(prod1, prod2))
constr <- list(x == 1.0, y == 4.0)

dgp <- Problem(obj, constr)
dgp2dcp <- Dgp2Dcp(dgp)
tmp <- reduce(dgp2dcp)
dgp2dcp <- tmp[[1]]
dcp <- tmp[[2]]
opt <- solve(dcp)
solution <- form_solution(dcp, opt)
dgp_unpack <- unpack_problem(dgp, retrieve(dgp2dcp, solution))
expect_equal(dgp_unpack\$value, 6.0, tolerance = TOL)
expect_equal(dgp_unpack\$getValue(x), 1.0, tolerance = TOL)
expect_equal(dgp_unpack\$getValue(y), 4.0, tolerance = TOL)
result <- solve(dgp, gp = TRUE)
expect_equal(result\$value, 6.0, tolerance = TOL)
expect_equal(result\$getValue(x), 1.0, tolerance = TOL)
})

# TODO_NARAS_5: Test keepdims. Need to edit test to match CVXPY.
test_that("test prod_entries", {
skip_on_cran()
X <- matrix(0:11, nrow = 4, ncol = 3)
expect_equal(prod(X), value(prod_entries(X)))
expect_equal(apply(X, 1, prod), value(prod_entries(X, axis = 1)))
expect_equal(apply(X, 2, prod), value(prod_entries(X, axis = 2)))

prod <- prod_entries(X)
X_canon <- Dgp2Dcp.prod_canon(prod, prod@args)[[1]]
expect_equal(sum(X), value(X_canon))

prod <- prod_entries(X, axis = 1)
X_canon <- Dgp2Dcp.prod_canon(prod, prod@args)[[1]]
expect_equal(apply(X, 1, sum), value(X_canon))

prod <- prod_entries(X, axis = 2)
X_canon <- Dgp2Dcp.prod_canon(prod, prod@args)[[1]]
expect_equal(apply(X, 2, sum), value(X_canon))

X <- matrix(0:11, nrow = 12, ncol = 1)
expect_equal(prod(X), value(prod_entries(X)))

prod <- prod_entries(X)
X_canon <- Dgp2Dcp.prod_canon(prod, prod@args)[[1]]
expect_equal(sum(X), value(X_canon))

x <- Variable(pos = TRUE)
y <- Variable(pos = TRUE)
posy1 <- x*y^(0.5) + 3.0*x*y^(0.5)
posy2 <- x*y^(0.5) + 3.0*x^2*y^(0.5)
expect_true(is_log_log_convex(prod_entries(posy1, posy2)))
expect_false(is_log_log_concave(prod_entries(posy1, posy2)))
expect_false(is_dgp(prod_entries(posy1, 1/posy1)))

m <- x*y^(0.5)
expect_true(is_log_log_affine(prod_entries(m, m)))
expect_true(is_log_log_concave(prod_entries(m, 1/posy1)))
expect_false(is_log_log_convex(prod_entries(m, 1/posy1)))
})

test_that("test max_entries", {
skip_on_cran()
x <- Variable(pos = TRUE)
y <- Variable(pos = TRUE)
z <- Variable(pos = TRUE)

prod1 <- x*y^(0.5)
prod2 <- 3.0*x*y^(0.5)
obj <- Minimize(max_entries(hstack(prod1, prod2)))
constr <- list(x == 1.0, y == 4.0)

dgp <- Problem(obj, constr)
dgp2dcp <- Dgp2Dcp(dgp)
tmp <- reduce(dgp2dcp)
dgp2dcp <- tmp[[1]]
dcp <- tmp[[2]]
opt <- solve(dcp)
solution <- form_solution(dcp, opt)
dgp_unpack <- unpack_problem(dgp, retrieve(dgp2dcp, solution))
expect_equal(dgp_unpack\$value, 6.0)
expect_equal(dgp_unpack\$getValue(x), 1.0)
expect_equal(dgp_unpack\$getValue(y), 4.0)
result <- solve(dgp, gp = TRUE)
expect_equal(result\$value, 6.0, tolerance = TOL)
expect_equal(result\$getValue(x), 1.0, tolerance = TOL)
})

test_that("test min_elemwise", {
skip_on_cran()
x <- Variable(pos = TRUE)
y <- Variable(pos = TRUE)
z <- Variable(pos = TRUE)

prod1 <- x*y^(0.5)
prod2 <- 3.0*x*y^(0.5)
posy <- prod1 + prod2
obj <- Maximize(min_elemwise(prod1, prod2, 1/posy))
constr <- list(x == 1.0, y == 4.0)

dgp <- Problem(obj, constr)
result <- solve(dgp, gp = TRUE)
expect_equal(result\$value, 1.0/(2.0 + 6.0), tolerance = TOL)
expect_equal(result\$getValue(x), 1.0, tolerance = TOL)
expect_equal(result\$getValue(y), 4.0, tolerance = TOL)
})

test_that("test min_entries", {
skip_on_cran()
x <- Variable(pos = TRUE)
y <- Variable(pos = TRUE)
z <- Variable(pos = TRUE)

prod1 <- x*y^(0.5)
prod2 <- 3.0*x*y^(0.5)
posy <- prod1 + prod2
obj <- Maximize(min_entries(hstack(prod1, prod2, 1/posy)))
constr <- list(x == 1.0, y == 4.0)

dgp <- Problem(obj, constr)
result <- solve(dgp, gp = TRUE)
expect_equal(result\$value, 1.0/(2.0 + 6.0), tolerance = TOL)
expect_equal(result\$getValue(x), 1.0, tolerance = TOL)
expect_equal(result\$getValue(y), 4.0, tolerance = TOL)
})

# TODO: Enable test once sum_largest is implemented.
# test_that("test sum_largest", {
#   x <- Variable(4, pos = TRUE)
#   obj <- Minimize(sum_largest(x, 3))
#   constr <- list(x[1]*x[2]*x[3]*x[4] >= 16)
#   dgp <- Problem(obj, constr)
#   dgp2dcp <- Dgp2Dcp(dgp)
#   tmp <- reduce(dgp2dcp)
#   dgp2dcp <- tmp[[1]]
#   dcp <- tmp[[2]]
#   result <- solve(dcp)
#   dgp <- unpack_problem(dgp, retrieve(dgp2dcp, result\$solution))
#   opt <- 6.0
#   expect_equal([email protected], opt, tolerance = TOL)
#   expect_equal(value(x[1]*x[2]*x[3]*x[4]), 16, tolerance = 1e-2)
#   dgp <- clear_solution(dgp)
#   result <- solve(dgp, gp = TRUE)
#   expect_equal(result\$value, opt, tolerance = TOL)
#   expect_equal(value(x[1]*x[2]*x[3]*x[4]), 16, tolerance = 1e-2)
#
#   # An unbounded problem.
#   x <- Variable(4, pos = TRUE)
#   y <- Variable(pos = TRUE)
#   obj <- Minimize(sum_largest(x, 3)*y)
#   constr <- list(x[1]*x[2]*x[3]*x[4] >= 16)
#   dgp <- Problem(obj, constr)
#   dgp2dcp <- Dgp2Dcp(dgp)
#   tmp <- reduce(dgp2dcp)
#   dgp2dcp <- tmp[[1]]
#   dcp <- tmp[[2]]
#   opt <- solve(dcp)
#   expect_equal([email protected], -Inf)
#   dgp <- unpack_problem(dgp, retrieve(dgp2dcp, opt\$solution))
#   expect_equal([email protected], 0.0, tolerance = TOL)
#   expect_equal([email protected], "unbounded")
#   dgp <- clear_solution(dgp)
#   result <- solve(dgp, gp = TRUE)
#   expect_equal(result\$value, 0.0, tolerance = TOL)
#   expect_equal(result\$status, "unbounded")
#
#   # Another unbounded problem.
#   x <- Variable(2, pos = TRUE)
#   obj <- Minimize(sum_largest(x, 1))
#   dgp <- Problem(obj)
#   dgp2dcp <- Dgp2Dcp(dgp)
#   tmp <- reduce(dgp2dcp)
#   dgp2dcp <- tmp[[1]]
#   dcp <- tmp[[2]]
#   opt <- solve(dcp)
#   expect_equal(result\$value, -Inf)
#   dgp <- unpack_problem(dgp, retrieve(dgp2dcp, opt\$solution))
#   expect_equal([email protected], 0.0, tolerance = TOL)
#   expect_equal([email protected], "unbounded")
#   dgp <- clear_solution(dgp)
#   result <- solve(dgp, gp = TRUE)
#   expect_equal([email protected], 0.0, tolerance = TOL)
#   expect_equal([email protected], "unbounded")
#
#   # Composition with posynomials.
#   x <- Variable(4, pos = TRUE)
#   obj <- Minimize(sum_largest(hstack(list(3*x[1]^(0.5) * x[2]^(0.5), x[1]*x[2] + 0.5*x[2]*x[4]^3, x[3])), 2))
#   constr <- list(x[1]*x[2] >= 16)
#   dgp <- Problem(obj, constr)
#   dgp2dcp <- Dgp2Dcp(dgp)
#   dcp <- reduce(dgp2dcp)
#   result <- solve(dcp)
#   dgp <- unpack_problem(dgp, retrieve(dgp2dcp, result\$solution))
#
#   # opt = 3 * sqrt(4) * sqrt(4) + (4 * 4 + 0.5 * 4 * epsilon) = 28
#   opt <- 28.0
#   expect_equal([email protected], opt, tolerance = 1e-2)
#   expect_equal(value(x[1]*x[2]), 16.0, tolerance = 1e-2)
#   expect_equal(value(x[4]), 0.0, tolerance = 1e-2)
#   dgp <- clear_solution(dgp)
#   result <- solve(dgp, gp = TRUE)
#   expect_equal(result\$value, opt, tolerance = 1e-2)
#   expect_equal(value(x[1]*x[2]), 16.0, tolerance = 1e-2)
#   expect_equal(value(x[4]), 0.0, tolerance = 1e-2)
# })

test_that("test div", {
skip_on_cran()
x <- Variable(pos = TRUE)
y <- Variable(pos = TRUE)
p <- Problem(Minimize(x*y), list(y/3 <= x, y >= 1))
result <- solve(p, gp = TRUE)
expect_equal(result\$value, 1.0/3.0, tolerance = TOL)
expect_equal(result\$getValue(y), 1.0, tolerance = TOL)
expect_equal(result\$getValue(x), 1.0/3.0, tolerance = TOL)
})

test_that("test geo_mean", {
skip_on_cran()
x <- Variable(3, pos = TRUE)
p <- c(1, 2, 0.5)
expr <- geo_mean(x, p)
dgp <- Problem(Minimize(expr))
dgp2dcp <- Dgp2Dcp(dgp)
tmp <- reduce(dgp2dcp)
dgp2dcp <- tmp[[1]]
dcp <- tmp[[2]]
result <- solve(dcp)
expect_equal(result\$value, -Inf)

solution <- form_solution(dcp, result)
dgp_unpack <- unpack_problem(dgp, retrieve(dgp2dcp, solution))
expect_equal(dgp_unpack\$value, 0.0)
expect_equal(dgp_unpack\$status, "unbounded")
result <- solve(dgp, gp = TRUE)
expect_equal(result\$value, 0.0)
expect_equal(result\$status, "unbounded")
})

test_that("test solving non-dgp problem raises error", {
skip_on_cran()
problem <- Problem(Minimize(-1.0*Variable()))
expect_error(solve(problem, gp = TRUE))
result <- solve(problem)
expect_equal(result\$status, "unbounded")
expect_equal(result\$value, -Inf)
})

test_that("test solving non-dcp problem raises error", {
skip_on_cran()
problem <- Problem(Minimize(Variable(pos = TRUE) * Variable(pos = TRUE)))
expect_error(solve(problem))
result <- solve(problem, gp = TRUE)
expect_equal(result\$status, "unbounded", tolerance = TOL)
expect_equal(result\$value, 0.0, tolerance = TOL)
})

skip_on_cran()
X <- Constant(rbind(1:3, 4:6))
Y <- Constant(rbind(2:4, 5:7))
Z <- X + Y
canon_matrix <- canon[[1]]
constraints <- canon[[2]]
expect_equal(length(constraints), 0)
expect_equal(dim(canon_matrix), dim(Z))
expected <- log(exp(value(X)) + exp(value(Y)))
expect_equal(expected, value(canon_matrix), tolerance = TOL)

# Test promotion.
X <- Constant(rbind(1:3, 4:6))
y <- Constant(2.0)
Z <- X + y
canon_matrix <- canon[[1]]
constraints <- canon[[2]]
expect_equal(length(constraints), 0)
expect_equal(dim(canon_matrix), dim(Z))
expected <- log(exp(value(X)) + as.numeric(exp(value(y))))
expect_equal(expected, value(canon_matrix), tolerance = TOL)
})

test_that("test matmul_canon", {
skip_on_cran()
X <- Constant(rbind(1:3, 4:6))
Y <- Constant(matrix(1:3, nrow = 3))
Z <- X %*% Y
canon <- Dgp2Dcp.mulexpression_canon(Z, Z@args)
canon_matrix <- canon[[1]]
constraints <- canon[[2]]
expect_equal(length(constraints), 0)
expect_equal(dim(canon_matrix), c(2,1))
first_entry <- log(exp(2.0) + exp(4.0) + exp(6.0))
second_entry <- log(exp(5.0) + exp(7.0) + exp(9.0))
expect_equal(first_entry, as.numeric(value(canon_matrix[1,1])), tolerance = TOL)
expect_equal(second_entry, as.numeric(value(canon_matrix[2,1])), tolerance = TOL)
})

test_that("test trace_canon", {
skip_on_cran()
X <- Constant(rbind(c(1.0, 5.0), c(9.0, 14.0)))
Y <- matrix_trace(X)
canon <- Dgp2Dcp.trace_canon(Y, Y@args)
canon_matrix <- canon[[1]]
constraints <- canon[[2]]
expect_equal(length(constraints), 0)
expect_true(is_scalar(canon_matrix))
expected <- log(exp(1.0) + exp(14.0))
expect_equal(expected, value(canon_matrix), tolerance = TOL)
})

test_that("test one_minus_pos", {
skip_on_cran()
x <- Variable(pos = TRUE)
obj <- Maximize(x)
constr <- list(one_minus_pos(x) >= 0.4)
problem <- Problem(obj, constr)
result <- solve(problem, gp = TRUE)
expect_equal(result\$value, 0.6, tolerance = TOL)
expect_equal(result\$getValue(x), 0.6, tolerance = TOL)
})

test_that("test qp solver not allowed", {
skip_on_cran()
x <- Variable(pos = TRUE)
problem <- Problem(Minimize(x))
expect_error(solve(problem, solver = "OSQP", gp = TRUE))
})

# TODO: Enable test once sum_largest is implemented.
# test_that("test paper example sum_largest", {
#   x <- Variable(4, pos = TRUE)
#   x1 <- x[1]
#   x2 <- x[2]
#   x3 <- x[3]
#   x4 <- x[4]
#   obj <- Minimize(sum_largest(hstack(3*x1^(0.5) * x2^(0.5), x1*x2 + 0.5*x2*x4^3, x3), 2))
#   constr <- list(x1*x2*x3 >= 16)
#   p <- Problem(obj, constr)
#   expect_silent(solve(p, gp = TRUE))   # Smoke test.
# })

test_that("test paper example one_minus_pos", {
skip_on_cran()
x <- Variable(pos = TRUE)
y <- Variable(pos = TRUE)
obj <- Minimize(x*y)
constr <- list((y * one_minus_pos(x/y))^2 >= 1, x >= y/3)
problem <- Problem(obj, constr)
expect_silent(solve(problem, gp = TRUE))   # Smoke test.
})

test_that("test paper example eye_minus_inv", {
skip_on_cran()
X <- Variable(2,2, pos = TRUE)
obj <- Minimize(matrix_trace(eye_minus_inv(X)))
constr <- list(geo_mean(diag(X)) == 0.1)
problem <- Problem(obj, constr)
expect_silent(solve(problem, gp = TRUE, solver = "SCS"))   # Smoke test.
})

test_that("test paper example exp_log", {
skip_on_cran()
x <- Variable(pos = TRUE)
y <- Variable(pos = TRUE)
obj <- Minimize(x*y)
constr <- list(exp(y/x) <= log(y))
problem <- Problem(obj, constr)
expect_silent(solve(problem, gp = TRUE))   # Smoke test.
})

test_that("test pf matrix completion", {
skip_on_cran()
X <- Variable(3,3, pos = TRUE)
obj <- Minimize(pf_eigenvalue(X))
known_indices <- cbind(c(1,1,2,3,3), c(1,3,2,1,2))
constr <- list(X[known_indices] == c(1.0, 1.9, 0.8, 3.2, 5.9),
X[1,2] * X[2,1] * X[2,3] * X[3,3] == 1.0)
problem <- Problem(obj, constr)
expect_silent(solve(problem, gp = TRUE))   # Smoke test.
})

test_that("test rank one nmf", {
skip_on_cran()
X <- Variable(3,3, pos = TRUE)
x <- Variable(3, pos = TRUE)
y <- Variable(3, pos = TRUE)
xy <- hstack(x[1]*y, x[2]*y, x[3]*y)
R <- max_elemwise(multiply(X, (xy)^(-1.0)), multiply(X^(-1.0), xy))
objective <- sum(R)
constraints <- list(X[1,1] == 1.0, X[1,3] == 1.9, X[2,2] == 0.8, X[3,1] == 3.2, X[3,2] == 5.9,
x[1] * x[2] * x[3] == 1.0)
# Smoke test.
prob <- Problem(Minimize(objective), constraints)
expect_silent(solve(prob, gp = TRUE))
})

test_that("test documentation problem", {
skip_on_cran()
x <- Variable(pos = TRUE)
y <- Variable(pos = TRUE)
z <- Variable(pos = TRUE)

objective_fn <- x*y*z
constraints <- list(4*x*y*z + 2*x*z <= 10, x <= 2*y, y <= 2*x, z >= 1)
problem <- Problem(Maximize(objective_fn), constraints)
expect_silent(solve(problem, gp = TRUE))   # Smoke test.
})

test_that("test solver error", {
skip_on_cran()
x <- Variable(pos = TRUE)
y <- Variable(pos = TRUE)
prod <- x*y
dgp <- Problem(Minimize(prod))
dgp2dcp <- Dgp2Dcp()
tmp <- perform(dgp2dcp, dgp)
dgp2dcp <- tmp[[1]]
data <- tmp[[2]]
inverse_data <- tmp[[3]]

soln <- Solution("solver_error", NA_real_, list(), list(), list())
dgp_soln <- invert(dgp2dcp, soln, inverse_data)
expect_equal(dgp_soln@status, "solver_error")
})
```
anqif/cvxr documentation built on April 8, 2020, 5:08 p.m.