diff --git a/LICENSE.md b/LICENSE.md index b8af352a1..626a48bb5 100644 --- a/LICENSE.md +++ b/LICENSE.md @@ -24,9 +24,9 @@ The Convex.jl package is licensed under the Simplified "2-clause" BSD License: > (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE > OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -The file benchmark/benchmarks/benchmark.jl contains some utilities copied -from the MathOptInterface.jl package (https://github.com/JuliaOpt/MathOptInterface.jl) -which is licensed under the MIT License: +The file src/problem_depot/problem_depot.jl contains some utilities copied from +the MathOptInterface.jl package (https://github.com/JuliaOpt/MathOptInterface.jl) +which is licensed under the following MIT License: >Copyright (c) 2017: Miles Lubin and contributors Copyright (c) 2017: Google Inc. > diff --git a/Project.toml b/Project.toml index 536586376..27ca76071 100644 --- a/Project.toml +++ b/Project.toml @@ -4,10 +4,12 @@ version = "0.12.5" [deps] AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c" +BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MathProgBase = "fdba3010-5040-5b88-9595-932c9decdf73" OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] AbstractTrees = "^0.2.1" diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl index 063b7148f..c56b4dc86 100644 --- a/benchmark/benchmarks.jl +++ b/benchmark/benchmarks.jl @@ -2,20 +2,32 @@ using Pkg tempdir = mktempdir() Pkg.activate(tempdir) Pkg.develop(PackageSpec(path=joinpath(@__DIR__, ".."))) -Pkg.add(["BenchmarkTools", "SCS", "ECOS", "PkgBenchmark"]) +Pkg.add(["BenchmarkTools", "PkgBenchmark"]) Pkg.resolve() -using Convex: Convex -using SCS: SCSSolver -using ECOS: ECOSSolver +using Convex: Convex, ProblemDepot using BenchmarkTools -include("benchmarks/benchmarks.jl") # defines module Benchmarks - const SUITE = BenchmarkGroup() -SUITE["SCS"] = Benchmarks.suite(p -> Convex.solve!(p, SCSSolver(verbose=0))) -SUITE["ECOS"] = Benchmarks.suite(p -> Convex.solve!(p, ECOSSolver(verbose=0)); - exclude = [r"sdp", r"SDP"]) -SUITE["formulation"] = Benchmarks.suite(Convex.conic_problem) +problems = [ + "constant_fix!_with_complex_numbers", + "affine_dot_multiply_atom", + "affine_hcat_atom", + "affine_trace_atom", + "exp_entropy_atom", + "exp_log_perspective_atom", + "socp_norm_2_atom", + "socp_quad_form_atom", + "socp_sum_squares_atom", + "lp_norm_inf_atom", + "lp_maximum_atom", + "sdp_and_exp_log_det_atom", + "sdp_norm2_atom", + "sdp_lambda_min_atom", + "sdp_sum_largest_eigs", + "mip_integer_variables", + ] + +SUITE["formulation"] = ProblemDepot.benchmark_suite(Convex.conic_problem, problems) diff --git a/benchmark/benchmarks/affine_benchmarks.jl b/benchmark/benchmarks/affine_benchmarks.jl deleted file mode 100644 index a64112191..000000000 --- a/benchmark/benchmarks/affine_benchmarks.jl +++ /dev/null @@ -1,164 +0,0 @@ -@add_benchmark function affine_negate_atom(handle_problem) - x = Variable() - p = minimize(-x, [x <= 0]) - handle_problem(p) -end - -@add_benchmark function affine_kron_atom(handle_problem) - x = ComplexVariable(3, 3) - y = [1.0 2.0; 3.0 4.0] - p = satisfy(kron(x,y) == kron(eye(3), y)) - handle_problem(p) -end - -@add_benchmark function affine_mult_atom(handle_problem) - probs=[ - let - x = Variable(1) - p = minimize(2.0 * x, [x >= 2, x <= 4]) - end, - let - x = Variable(2) - A = 1.5 * eye(2) - p = minimize([2 2] * x, [A * x >= [1.1; 1.1]]) - end, - let - y = Variable(1) - x = Variable(3) - z = [1.0, 2.0, 3.0] * y - k = -y * [1.0, 2.0, 3.0] - c = [y <= 3.0, y >= 0.0, x >= ones(3), k <= x, x <= z] - o = 3 * y - p = Problem(:minimize, o, c) - end, - let - x = ComplexVariable(2,2) - p = minimize( real( [1.0im, 0.0]' * x * [1.0im, 0.0] ), - [ x == [1.0 0.0; 0.0 1.0] ] ) - end - ] - handle_problem.(probs) -end - -@add_benchmark function affine_dot_atom(handle_problem) - probs = [ - let - x = Variable(2) - p = minimize(dot([2.0; 2.0], x), x >= [1.1; 1.1]) - end, - let - x = Variable(2,2) - p = minimize(dot(fill(2.0, (2,2)), x), x >= 1.1) - end - ] - - handle_problem.(probs) -end - - -@add_benchmark function affine_add_atom(handle_problem) - - probs = [ - let - x = Variable(1) - y = Variable(1) - p = minimize(x + y, [x >= 3, y >= 2]) - end, - let - x = Variable(1) - p = minimize(x, [eye(2) + x >= eye(2)]) - end, - let - y = Variable() - p = minimize(y - 5, y >= -1) - end ] - - handle_problem.(probs) -end - - - -@add_benchmark function affine_transpose_atom(handle_problem) - probs = [ - let - x = Variable(2) - c = ones(2, 1) - p = minimize(x' * c, x >= 1) - end, - let - X = Variable(2, 2) - c = ones(2, 1) - p = minimize(c' * X' * c, [X >= ones(2, 2)]) - end, - let - rows = 2 - cols = 3 - r = rand(rows, cols) - r_2 = rand(cols, rows) - x = Variable(rows, cols) - c = ones(1, cols) - d = ones(rows, 1) - p = minimize(c * x' * d + d' * x * c' + (c * x''''' * d)', - [x' >= r_2, x >= r, x''' >= r_2, x'' >= r]) - end ] - - handle_problem.(probs) -end - -@add_benchmark function affine_index_atom(handle_problem) - probs = [ - let - x = Variable(2) - p = minimize(x[1] + x[2], [x >= 1]) - end, - let - x = Variable(3) - I = [true true false] - p = minimize(sum(x[I]), [x >= 1]) - end, - let - rows = 6 - cols = 8 - n = 2 - X = Variable(rows, cols) - A = randn(rows, cols) - c = rand(1, n) - p = minimize(c * X[1:n, 5:5+n-1]' * c', X >= A) - end] - handle_problem.(probs) -end - -@add_benchmark function affine_sum_atom(handle_problem) - probs = [ - let - x = Variable(2,2) - p = minimize(sum(x) - 2*x[1,1], x>=1, x[1,1]<=2) - end, - let - x = Variable(2,2) - p = minimize(sum(x) - 2*x[1,1], x>=1, x[1,1]<=2) - end - ] - return handle_problem.(probs) -end - -@add_benchmark function affine_diag_atom(handle_problem) - probs = [ - let - x = Variable(2,2) - p = minimize(sum(diag(x,1)), x >= 1) - end, - let - x = Variable(4, 4) - p = minimize(sum(diag(x)), x >= 2) - end, - ] - return handle_problem.(probs) - -end - -@add_benchmark function affine_tr_atom(handle_problem) - x = Variable(2,2) - p = minimize(tr(x), x >= 1) - handle_problem(p) -end diff --git a/benchmark/benchmarks/benchmarks.jl b/benchmark/benchmarks/benchmarks.jl deleted file mode 100644 index ec22a274c..000000000 --- a/benchmark/benchmarks/benchmarks.jl +++ /dev/null @@ -1,56 +0,0 @@ -# `BENCHMARKS`, `suite`, and `add_benchmark` were taken from MathOptInterface -# which is available under an MIT license (see LICENSE). -module Benchmarks - -using BenchmarkTools -using Convex, LinearAlgebra - -const BENCHMARKS = Dict{String, Function}() - -""" - suite( - handle_problem::Function; - exclude::Vector{Regex} = Regex[] - ) - -Create a suite of benchmarks. `handle_problem` should be a function that takes one -argument, a Convex.jl `Problem` and processes it (e.g. `solve!` the problem with -a specific solver). - -Use `exclude` to exclude a subset of benchmarks. - -### Examples - -```julia -suite() do p - solve!(p, GLPK.Optimizer()) -end -``` -""" -function suite(handle_problem::Function; exclude::Vector{Regex} = Regex[]) - group = BenchmarkGroup() - for (name, func) in BENCHMARKS - any(occursin.(exclude, Ref(name))) && continue - group[name] = @benchmarkable $func($handle_problem) setup=Convex.clearmemory() - end - return group -end - -### -### Benchmarks -### - -macro add_benchmark(f) - name = f.args[1].args[1] - return quote - $(esc(f)) - BENCHMARKS[String($(Base.Meta.quot(name)))] = $(esc(name)) - end -end - -eye(n) = Matrix(1.0I, n, n) -include("constraint_benchmarks.jl") -include("affine_benchmarks.jl") -include("sdp_benchmarks.jl") - -end diff --git a/benchmark/benchmarks/sdp_benchmarks.jl b/benchmark/benchmarks/sdp_benchmarks.jl deleted file mode 100644 index 670124401..000000000 --- a/benchmark/benchmarks/sdp_benchmarks.jl +++ /dev/null @@ -1,58 +0,0 @@ -@add_benchmark function sdp_nuclear_norm_atom(handle_problem) - y = Semidefinite(3) - p = minimize(nuclearnorm(y), y[2,1]<=4, y[2,2]>=3, y[3,3]<=2) - handle_problem(p) -end - -@add_benchmark function sdp_operator_norm_atom(handle_problem) - y = Variable((3,3)) - p = minimize(opnorm(y), y[2,1]<=4, y[2,2]>=3, sum(y)>=12) - handle_problem(p) -end - -@add_benchmark function sdp_sigma_max_atom(handle_problem) - y = Variable((3,3)) - p = minimize(sigmamax(y), y[2,1]<=4, y[2,2]>=3, sum(y)>=12) - handle_problem(p) -end - -@add_benchmark function sdp_lambda_max_atom(handle_problem) - y = Semidefinite(3) - p = minimize(lambdamax(y), y[1,1]>=4) - handle_problem(p) -end - -@add_benchmark function sdp_lambda_min_atom(handle_problem) - y = Semidefinite(3) - p = maximize(lambdamin(y), tr(y)<=6) - handle_problem(p) -end - - -@add_benchmark function sdp_matrix_frac_atom(handle_problem) - probs = [ - let - x = [1, 2, 3] - P = Variable(3, 3) - p = minimize(matrixfrac(x, P), P <= 2*eye(3), P >= 0.5 * eye(3)) - end, - let - x = Variable(3) - P = Variable(3, 3) - p = minimize(matrixfrac(x, P), lambdamax(P) <= 2, x[1] >= 1) - end - ] - - handle_problem.(probs) -end - -@add_benchmark function sdp_sum_squares_atom(handle_problem) - n = 10 - A = rand(n,n) + im*rand(n,n) - A = A + A' # now A is hermitian - x = ComplexVariable(n,n) - objective = sumsquares(A - x) - c1 = x in :SDP - p = minimize(objective, c1) - handle_problem(p) -end diff --git a/docs/make.jl b/docs/make.jl index 742b796d9..358fdb010 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -15,6 +15,7 @@ makedocs(; "FAQ" => "faq.md", "Optimizing in a Loop" => "loop.md", "Advanced" => "advanced.md", + "Problem Depot" => "problem_depot.md", "Contributing" => "contributing.md", "Credits" => "credits.md" ], diff --git a/docs/src/contributing.md b/docs/src/contributing.md index f9cf4a058..73e46b168 100644 --- a/docs/src/contributing.md +++ b/docs/src/contributing.md @@ -5,11 +5,11 @@ We'd welcome contributions to the Convex.jl package. Here are some short instructions on how to get started. If you don't know what you'd like to contribute, you could -> - take a look at the current -> [issues](https://github.com/JuliaOpt/Convex.jl/issues) and pick -> one. (Feature requests are probably the easiest to tackle.) -> - add a [usage -> example](https://github.com/JuliaOpt/Convex.jl/tree/master/examples). + - take a look at the current + [issues](https://github.com/JuliaOpt/Convex.jl/issues) and pick + one. (Feature requests are probably the easiest to tackle.) + - add a [usage + example](https://github.com/JuliaOpt/Convex.jl/tree/master/examples). Then submit a pull request (PR). (Let us know if it's a work in progress by putting \[WIP\] in the name of the PR.) @@ -17,13 +17,13 @@ progress by putting \[WIP\] in the name of the PR.) Adding examples --------------- -> - Take a look at our exising [usage -> examples](https://github.com/JuliaOpt/Convex.jl/tree/master/examples) -> and add another in similar style. -> - Submit a PR. (Let us know if it's a work in progress by putting -> \[WIP\] in the name of the PR.) -> - We'll look it over, fix up anything that doesn't work, and merge -> it! + - Take a look at our exising [usage + examples](https://github.com/JuliaOpt/Convex.jl/tree/master/examples) + and add another in similar style. + - Submit a PR. (Let us know if it's a work in progress by putting + \[WIP\] in the name of the PR.) + - We'll look it over, fix up anything that doesn't work, and merge + it! Adding atoms ------------ @@ -31,31 +31,32 @@ Adding atoms Here are the steps to add a new function or operation (atom) to Convex.jl. Let's say you're adding the new function $f$. -> - Take a look at the [nuclear norm -> atom](https://github.com/JuliaOpt/Convex.jl/blob/master/src/atoms/sdp_cone/nuclearnorm.jl) -> for an example of how to construct atoms, and see the [norm -> atom](https://github.com/JuliaOpt/Convex.jl/blob/master/src/atoms/second_order_cone/norm.jl) -> for an example of an atom that depends on a parameter. -> - Copy paste (eg) the nuclear norm file, replace anything saying -> nuclear norm with the name of the atom $f$, fill in monotonicity, -> curvature, etc. Save it in the appropriate subfolder of -> `src/atoms/`. -> - Add as a comment a description of what the atom does and its -> parameters. -> - The most mathematically interesting part is the `conic_form!` -> function. Following the example in the nuclear norm atom, you'll -> see that you can just construct the problem whose optimal value is -> $f(x)$, introducing any auxiliary variables you need, exactly as -> you would normally in Convex.jl, and then call `cache_conic_form!` -> on that problem. -> - Add a test for the atom so we can verify it works in -> `test/test_`, where `` matches the subfolder of -> `src/atoms`. -> - Submit a PR, including a description of what the atom does and its -> parameters. (Let us know if it's a work in progress by putting -> \[WIP\] in the name of the PR.) -> - We'll look it over, fix up anything that doesn't work, and merge -> it! + - Take a look at the [nuclear norm + atom](https://github.com/JuliaOpt/Convex.jl/blob/master/src/atoms/sdp_cone/nuclearnorm.jl) + for an example of how to construct atoms, and see the [norm + atom](https://github.com/JuliaOpt/Convex.jl/blob/master/src/atoms/second_order_cone/norm.jl) + for an example of an atom that depends on a parameter. + - Copy paste (eg) the nuclear norm file, replace anything saying + nuclear norm with the name of the atom $f$, fill in monotonicity, + curvature, etc. Save it in the appropriate subfolder of + `src/atoms/`. + - Add as a comment a description of what the atom does and its + parameters. + - The most mathematically interesting part is the `conic_form!` + function. Following the example in the nuclear norm atom, you'll + see that you can just construct the problem whose optimal value is + $f(x)$, introducing any auxiliary variables you need, exactly as + you would normally in Convex.jl, and then call `cache_conic_form!` + on that problem. + - Add a test for the atom so we can verify it works in + `src/problem_depot/problem/`, where `` matches the subfolder of + `src/atoms`. See [How to write a ProblemDepot problem](@ref) for details + on how to write the tests. + - Submit a PR, including a description of what the atom does and its + parameters. (Let us know if it's a work in progress by putting + \[WIP\] in the name of the PR.) + - We'll look it over, fix up anything that doesn't work, and merge + it! Fixing the guts --------------- @@ -69,18 +70,18 @@ helpful as well; you can find the ipython notebook presented in the talk Then read the conic form code: -> - We define data structures for conic objectives and conic -> constraints, and simple ways of combining them, in -> [conic\_form.jl](https://github.com/JuliaOpt/Convex.jl/blob/master/src/conic_form.jl) -> - We convert the internal conic form representation into the -> [standard form for conic -> solvers](http://mathprogbasejl.readthedocs.io/en/latest/conic.html) -> in the function -> [conic\_problem](https://github.com/JuliaOpt/Convex.jl/blob/master/src/problems.jl#L97). -> - We solve problems (that is, pass the standard form of the problem -> to a solver, and put the solution back into the values of the -> appropriate variables) in -> [solution.jl](https://github.com/JuliaOpt/Convex.jl/blob/master/src/solution.jl#L8). + - We define data structures for conic objectives and conic + constraints, and simple ways of combining them, in + [conic\_form.jl](https://github.com/JuliaOpt/Convex.jl/blob/master/src/conic_form.jl) + - We convert the internal conic form representation into the + [standard form for conic + solvers](http://mathprogbasejl.readthedocs.io/en/latest/conic.html) + in the function + [conic\_problem](https://github.com/JuliaOpt/Convex.jl/blob/master/src/problems.jl#L97). + - We solve problems (that is, pass the standard form of the problem + to a solver, and put the solution back into the values of the + appropriate variables) in + [solution.jl](https://github.com/JuliaOpt/Convex.jl/blob/master/src/solution.jl#L8). You're now armed and dangerous. Go ahead and open an issue (or comment on a previous one) if you can't figure something out, or submit a PR if diff --git a/docs/src/problem_depot.md b/docs/src/problem_depot.md new file mode 100644 index 000000000..2de7151f7 --- /dev/null +++ b/docs/src/problem_depot.md @@ -0,0 +1,97 @@ +Problem Depot +============= + +Convex.jl has a submodule, `ProblemDepot` which holds a collection of convex optimization problems. The problems are used by Convex itself to test and benchmark its code, but can also be used by solvers to test and benchmark their code. + +ProblemDepot has two main methods for accessing these problems: `Convex.ProblemDepot.run_tests` and `Convex.ProblemDepot.benchmark_suite`. + +For example, to test the solver SCS on all the problems of the depot except the mixed-integer problems (which it cannot handle), run + +```julia +using Convex, SCS, Test +@testset "SCS" begin + Convex.ProblemDepot.run_tests(; exclude=[r"mip"]) do p + solve!(p, SCSSolver(verbose=0, eps=1e-6)) + end +end +``` + +How to write a ProblemDepot problem +----------------------------------- + +The problems are organized into folders in `src/problem_depot/problems`. Each is written as a function, annotated by `@add_problem`, and a name, which is used to group the problems. For example, here is a simple problem: + +```julia +@add_problem affine function affine_negate_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + x = Variable() + p = minimize(-x, [x <= 0]) + if test + @test vexity(p) == AffineVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 0 atol=atol rtol=rtol + @test evaluate(-x) ≈ 0 atol=atol rtol=rtol + end +end +``` + +The `@add_problem` call adds the problem to the registry of problems in [`Convex.ProblemDepot.PROBLEMS`](@ref), which in turn is used by [`run_tests`](@ref) and [`benchmark_suite`](@ref). Next, `affine` is the grouping of the problem; this problem came from one of the affine tests, and in particular is testing the negation atom. Next is the function signature: + +```julia +function affine_negate_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} +``` + +this should be the same for every problem, except for the name, which is a description of the problem. It should include what kind of atoms it uses (`affine` in this case), so that certain kinds of atoms can be ruled out by the `exclude` keyword to [`run_tests`](@ref) and [`benchmark_suite`](@ref); for example, many solvers cannot solve mixed-integer problems, so `mip` is included in the name of such problems. + +Then begins the body of the problem. It is setup like any other Convex.jl problem, only `handle_problem!` is called instead of `solve!`. This allows particular solvers to be used (via e.g. choosing `handle_problem! = p -> solve!(p, solver)`), or for any other function of the problem (e.g. `handle_problem! = p -> Convex.conic_problem(p)` which is used for benchmarking problem formulation speed.) Tests should be included and gated behind `if test` blocks, so that tests can be skipped for benchmarking, or in the case that the problem is not in fact solved during `handle_problem!`. + +The fact that the problems may not be solved during `handle_problem!` brings with it a small complication: any command that assumes the problem has been solved should be behind an `if test` check. For example, in some of the problems, `real(x.value)` is used, for a variable `x`; perhaps as + +```julia +x_re = real(x.value) +if test + @test x_re = ... +end +``` + +However, if the problem `x` is used in has not been solved, then `x.value === nothing`, and `real(nothing)` throws an error. So instead, this should be rewritten as + +```julia +if test + x_re = real(x.value) + @test x_re = ... +end +``` + +Benchmark-only problems +----------------------- + +To add problems for benchmarking without tests, place problems in `src/problem_depot/problems/benchmark`, and include `benchmark` in the name. These problems will be automatically skipped during `run_tests` calls. For example, to benchmark the time it takes to add an SDP constraint, we have the problem + +```julia +@add_problem constraints_benchmark function sdp_constraint(handle_problem!, args...) + p = satisfy() + x = Variable(44, 44) # 990 vectorized entries + push!(p.constraints, x ⪰ 0) + handle_problem!(p) + nothing +end +``` + +However, this "problem" has no tests or interesting content for testing, so we skip it during testing. +Note, we use `args...` in the function signature so that it may be called with the standard function signature + +```julia +f(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} +``` + +Reference +--------- + +```@docs +Convex.ProblemDepot.run_tests +Convex.ProblemDepot.benchmark_suite +Convex.ProblemDepot.foreach_problem +Convex.ProblemDepot.PROBLEMS +``` diff --git a/src/Convex.jl b/src/Convex.jl index f18ec509f..ffd0d899b 100644 --- a/src/Convex.jl +++ b/src/Convex.jl @@ -87,6 +87,7 @@ include("utilities/tree_interface.jl") include("utilities/show.jl") include("utilities/iteration.jl") include("utilities/broadcast.jl") +include("problem_depot/problem_depot.jl") # Deprecated workaround for memory leak (https://github.com/JuliaOpt/Convex.jl/issues/83) function clearmemory() diff --git a/src/problem_depot/problem_depot.jl b/src/problem_depot/problem_depot.jl new file mode 100644 index 000000000..6bf0f89a5 --- /dev/null +++ b/src/problem_depot/problem_depot.jl @@ -0,0 +1,201 @@ +# Some code in `src/problem_depot` was modified from MathOptInterface +# which is available under an MIT license (see LICENSE). + +module ProblemDepot +using BenchmarkTools, Test + +using Convex +using LinearAlgebra +using LinearAlgebra: eigen, I, opnorm + +randperm(d) = sortperm(rand(d)) +shuffle(x) = x[randperm(length(x))] +mean(x) = sum(x) / length(x) +eye(n, T) = Matrix{T}(I, n, n) +eye(n) = Matrix{Float64}(I, n, n) + +""" + const PROBLEMS = Dict{String, Dict{String, Function}}() + +A "depot" of Convex.jl problems, subdivided into categories. +Each problem is stored as a function with the signature + + f(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + +where `handle_problem!` specifies what to do with the `Problem` instance +(e.g., `solve!` it with a chosen solver), an option `test` to choose +whether or not to test the values (assuming it has been solved), +tolerances for the tests, and a numeric type in which the problem +should be specified (currently, this is not respected and all +problems are specified in `Float64` precision). + +See also [`run_tests`](@ref) and [`benchmark_suite`](@ref) for helpers +to use these problems in testing or benchmarking. + +### Examples + +```julia +julia> PROBLEMS["affine"]["affine_diag_atom"] +affine_diag_atom (generic function with 1 method) +``` +""" +const PROBLEMS = Dict{String, Dict{String, Function}}() + +""" + foreach_problem(apply::Function, [class::String], + problems::Union{Nothing, Vector{String}, Vector{Regex}} = nothing; + exclude::Vector{Regex} = Regex[]) + +Provides a convience method for iterating over problems in [`PROBLEMS`](@ref). +For each problem in [`PROBLEMS`](@ref), apply the function `apply`, which +takes two arguments: the name of the function associated to the problem, +and the function associated to the problem itself. + +Optionally, pass a second argument `class` to only iterate over a class of +problems (`class` should satsify `class ∈ keys(PROBLEMS)`), and pass third +argument `problems` to only allow certain problems (specified by exact names or +regex). Use the `exclude` keyword argument to exclude problems by regex. +""" +function foreach_problem( apply::Function, + problems::Union{Nothing, Vector{String}, Vector{Regex}} = nothing; + exclude::Vector{Regex} = Regex[]) + for class in keys(PROBLEMS) + any(occursin.(exclude, Ref(class))) && continue + foreach_problem(apply, class, problems; exclude=exclude) + end +end + +function foreach_problem( apply::Function, + class::String, + problems::Union{Nothing, Vector{String}, Vector{Regex}} = nothing; + exclude::Vector{Regex} = Regex[]) + for (name, func) in PROBLEMS[class] + any(occursin.(exclude, Ref(name))) && continue + if problems !== nothing + problems isa Vector{String} && !(name ∈ problems) && continue + problems isa Vector{Regex} && !any(occursin.(problems, Ref(name))) && continue + end + apply(name, func) + end +end + + +""" + run_tests( + handle_problem!::Function; + problems::Union{Nothing, Vector{String}, Vector{Regex}} = nothing; + exclude::Vector{Regex} = Regex[], + T=Float64, atol=1e-3, rtol=0.0, + ) + +Run a set of tests. `handle_problem!` should be a function that takes one +argument, a Convex.jl `Problem` and processes it (e.g. `solve!` the problem with +a specific solver). + +Use `exclude` to exclude a subset of sets; automatically excludes +`r"benchmark"`. Optionally, pass a second argument `problems` to only allow certain problems +(specified by exact names or regex). The test tolerances specified by `atol` and +`rtol`. Set `T` to choose a numeric type for the problem. Currently +this is only used for choosing the type parameter of the underlying +MathOptInterface model, but not for the actual problem data. + +### Examples + +```julia +run_tests(exclude=[r"mip"]) do p + solve!(p, SCSSolver(verbose=0)) +end +``` +""" +function run_tests( handle_problem!::Function, + problems::Union{Nothing, Vector{String}, Vector{Regex}} = nothing; + exclude::Vector{Regex} = Regex[], T=Float64, atol=1e-3, rtol=0.0) + push!(exclude, r"benchmark") + for class in keys(PROBLEMS) + any(occursin.(exclude, Ref(class))) && continue + @testset "$class" begin + foreach_problem(class, problems; exclude=exclude) do name, problem_func + @testset "$name" begin + problem_func(handle_problem!, Val(true), atol, rtol, T) + end + end + end + end +end + + +""" + benchmark_suite( + handle_problem!::Function, + problems::Union{Nothing, Vector{String}, Vector{Regex}} = nothing; + exclude::Vector{Regex} = Regex[], + test = Val(false), + T=Float64, atol=1e-3, rtol=0.0, + ) + +Create a benchmark_suite of benchmarks. `handle_problem!` should be a function +that takes one argument, a Convex.jl `Problem` and processes it (e.g. `solve!` +the problem with a specific solver). Pass a second argument `problems` to specify +run benchmarks only with certain problems (specified by exact names or regex). + +Use `exclude` to exclude a subset of benchmarks. Optionally, pass a second +argument `problems` to only allow certain problems (specified by exact names or +regex). Set `test=true` to also check the answers, with tolerances specified by +`atol` and `rtol`. Set `T` to choose a numeric type for the problem. Currently +this is only used for choosing the type parameter of the underlying +MathOptInterface model, but not for the actual problem data. + +### Examples + +```julia +benchmark_suite(exclude=[r"mip"]) do p + solve!(p, SCSSolver(verbose=0)) +end +``` +""" +function benchmark_suite(handle_problem!::Function, + problems::Union{Nothing, Vector{String}, Vector{Regex}} = nothing; + exclude::Vector{Regex} = Regex[], + T=Float64, atol=1e-3, rtol=0.0, test = Val(false)) + group = BenchmarkGroup() + for class in keys(PROBLEMS) + any(occursin.(exclude, Ref(class))) && continue + group[class] = BenchmarkGroup() + foreach_problem(class, problems; exclude=exclude) do name, problem_func + group[class][name] = @benchmarkable $problem_func($handle_problem!, $test, $atol, $rtol, $T) + end + end + return group +end + + +macro add_problem(prefix, q) + @assert prefix isa Symbol + if q.head == :block + f = q.args[2] + elseif q.head == :function + f = q + else + error("head $(q.head) unexpected") + end + name = f.args[1].args[1] + if name isa Expr + name = name.args[1] + end + return quote + $(esc(f)) + dict = get!(PROBLEMS, String($(Base.Meta.quot(prefix))), Dict{String,Function}()) + dict[String($(Base.Meta.quot(name)))] = $(esc(name)) + end +end + +include("problems/affine.jl") +include("problems/constant.jl") +include("problems/exp.jl") +include("problems/lp.jl") +include("problems/mip.jl") +include("problems/sdp_and_exp.jl") +include("problems/sdp.jl") +include("problems/socp.jl") + +end diff --git a/src/problem_depot/problems/affine.jl b/src/problem_depot/problems/affine.jl new file mode 100644 index 000000000..4d197e026 --- /dev/null +++ b/src/problem_depot/problems/affine.jl @@ -0,0 +1,644 @@ +@add_problem affine function affine_negate_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + x = Variable() + p = minimize(-x, [x <= 0]) + if test + @test vexity(p) == AffineVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 0 atol=atol rtol=rtol + @test evaluate(-x) ≈ 0 atol=atol rtol=rtol + end +end + +@add_problem affine function affine_kron_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + x = ComplexVariable(3, 3) + y = [1.0 2.0; 3.0 4.0] + if test + @test size(kron(x, y)) == (6, 6) + @test size(kron(y, x)) == (6, 6) + end +end + +@add_problem affine function affine_multiply_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + x = Variable(1) + p = minimize(2.0 * x, [x >= 2, x <= 4]) + if test + @test vexity(p) == AffineVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 4 atol=atol rtol=rtol + @test (evaluate(2.0x))[1] ≈ 4 atol=atol rtol=rtol + end + + x = Variable(2) + A = 1.5 * eye(2) + p = minimize([2 2] * x, [A * x >= [1.1; 1.1]]) + if test + @test vexity(p) == AffineVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 2.93333 atol=atol rtol=rtol + @test (evaluate([2 2] * x))[1] ≈ 2.93333 atol=atol rtol=rtol + @test vec(evaluate(A * x)) ≈ [1.1; 1.1] atol=atol rtol=rtol + end + + y = Variable(1) + x = Variable(3) + z = [1.0, 2.0, 3.0] * y + k = -y * [1.0, 2.0, 3.0] + c = [y <= 3.0, y >= 0.0, x >= ones(3), k <= x, x <= z] + o = 3 * y + p = Problem(:minimize, o, c) + if test + @test vexity(p) == AffineVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 3 atol=atol rtol=rtol + end + + p = Problem(:minimize, o, c...) + if test + @test vexity(p) == AffineVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 3 atol=atol rtol=rtol + end + + # Check #274 + x = ComplexVariable(2,2) + p = minimize( real( [1.0im, 0.0]' * x * [1.0im, 0.0] ), [ x == [1.0 0.0; 0.0 1.0] ]) + handle_problem!(p) + if test + @test p.optval ≈ 1.0 atol=atol rtol=rtol + end +end + +@add_problem affine function affine_dot_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + x = Variable(2) + p = minimize(dot([2.0; 2.0], x), x >= [1.1; 1.1]) + if test + @test vexity(p) == AffineVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 4.4 atol=atol rtol=rtol + @test (evaluate(dot([2.0; 2.0], x)))[1] ≈ 4.4 atol=atol rtol=rtol + end +end + +@add_problem affine function affine_dot_atom_for_matrix_variables(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + x = Variable(2,2) + p = minimize(dot(fill(2.0, (2,2)), x), x >= 1.1) + if test + @test vexity(p) == AffineVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 8.8 atol=atol rtol=rtol + @test (evaluate(dot(fill(2.0, (2, 2)), x)))[1] ≈ 8.8 atol=atol rtol=rtol + end +end + +@add_problem affine function affine_add_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + x = Variable(1) + y = Variable(1) + p = minimize(x + y, [x >= 3, y >= 2]) + if test + @test vexity(p) == AffineVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 5 atol=atol rtol=rtol + @test evaluate(x + y) ≈ 5 atol=atol rtol=rtol + end + + x = Variable(1) + p = minimize(x, [eye(2) + x >= eye(2)]) + if test + @test vexity(p) == AffineVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 0 atol=atol rtol=rtol + @test evaluate(eye(2) + x) ≈ eye(2) atol=atol rtol=rtol + end + + y = Variable() + p = minimize(y - 5, y >= -1) + if test + @test vexity(p) == AffineVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ -6 atol=atol rtol=rtol + @test evaluate(y - 5) ≈ -6 atol=atol rtol=rtol + end +end + +@add_problem affine function affine_transpose_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + x = Variable(2) + c = ones(2, 1) + p = minimize(x' * c, x >= 1) + if test + @test vexity(p) == AffineVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 2 atol=atol rtol=rtol + @test (evaluate(x' * c))[1] ≈ 2 atol=atol rtol=rtol + end + + X = Variable(2, 2) + c = ones(2, 1) + p = minimize(c' * X' * c, [X >= ones(2, 2)]) + if test + @test vexity(p) == AffineVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 4 atol=atol rtol=rtol + @test (evaluate(c' * X' * c))[1] ≈ 4 atol=atol rtol=rtol + end + + rows = 2 + cols = 3 + r = rand(rows, cols) + r_2 = rand(cols, rows) + x = Variable(rows, cols) + c = ones(1, cols) + d = ones(rows, 1) + p = minimize(c * x' * d + d' * x * c' + (c * x''''' * d)', + [x' >= r_2, x >= r, x''' >= r_2, x'' >= r]) + if test + @test vexity(p) == AffineVexity() + end + handle_problem!(p) + s = sum(max.(r, r_2')) * 3 + if test + @test p.optval ≈ s atol=atol rtol=rtol + @test (evaluate(c * x' * d + d' * x * c' + (c * ((((x')')')')' * d)'))[1] ≈ s atol=atol rtol=rtol + end +end + +@add_problem affine function affine_index_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + x = Variable(2) + p = minimize(x[1] + x[2], [x >= 1]) + if test + @test vexity(p) == AffineVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 2 atol=atol rtol=rtol + @test (evaluate(x[1] + x[2]))[1] ≈ 2 atol=atol rtol=rtol + end + + x = Variable(3) + I = [true true false] + p = minimize(sum(x[I]), [x >= 1]) + if test + @test vexity(p) == AffineVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 2 atol=atol rtol=rtol + @test (evaluate(sum(x[I])))[1] ≈ 2 atol=atol rtol=rtol + end + + rows = 6 + cols = 8 + n = 2 + X = Variable(rows, cols) + A = randn(rows, cols) + c = rand(1, n) + p = minimize(c * X[1:n, 5:5+n-1]' * c', X >= A) + if test + @test vexity(p) == AffineVexity() + end + handle_problem!(p) + s = c * A[1:n, 5:5+n-1]' * c' + if test + @test p.optval ≈ s[1] atol=atol rtol=rtol + @test evaluate(c * (X[1:n, 5:(5 + n) - 1])' * c') ≈ s atol=atol rtol=rtol + end +end + +@add_problem affine function affine_sum_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + x = Variable(2,2) + p = minimize(sum(x), x>=1) + if test + @test vexity(p) == AffineVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 4 atol=atol rtol=rtol + @test evaluate(sum(x)) ≈ 4 atol=atol rtol=rtol + end + + x = Variable(2,2) + p = minimize(sum(x) - 2*x[1,1], x>=1, x[1,1]<=2) + if test + @test vexity(p) == AffineVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 1 atol=atol rtol=rtol + @test (evaluate(sum(x) - 2 * x[1, 1]))[1] ≈ 1 atol=atol rtol=rtol + end + + x = Variable(10) + a = rand(10, 1) + p = maximize(sum(x[2:6]), x <= a) + if test + @test vexity(p) == AffineVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ sum(a[2:6]) atol=atol rtol=rtol + @test evaluate(sum(x[2:6])) ≈ sum(a[2:6]) atol=atol rtol=rtol + end +end + +@add_problem affine function affine_diag_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + x = Variable(2,2) + p = minimize(sum(diag(x,1)), x >= 1) + if test + @test vexity(p) == AffineVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 1 atol=atol rtol=rtol + @test evaluate(sum(diag(x, 1))) ≈ 1 atol=atol rtol=rtol + end + + x = Variable(4, 4) + p = minimize(sum(diag(x)), x >= 2) + if test + @test vexity(p) == AffineVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 8 atol=atol rtol=rtol + @test evaluate(sum(diag(x))) ≈ 8 atol=atol rtol=rtol + end +end + +@add_problem affine function affine_trace_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + x = Variable(2,2) + p = minimize(tr(x), x >= 1) + if test + @test vexity(p) == AffineVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 2 atol=atol rtol=rtol + @test evaluate(tr(x)) ≈ 2 atol=atol rtol=rtol + end +end + +@add_problem affine function affine_dot_multiply_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + x = Variable(3) + p = maximize(sum(dot(*)(x,[1,2,3])), x<=1) + if test + @test vexity(p) == AffineVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 6 atol=atol rtol=rtol + @test evaluate(sum((dot(*))(x, [1, 2, 3]))) ≈ 6 atol=atol rtol=rtol + end + + x = Variable(3, 3) + p = maximize(sum(dot(*)(x,eye(3))), x<=1) + if test + @test vexity(p) == AffineVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 3 atol=atol rtol=rtol + @test evaluate(sum((dot(*))(x, eye(3)))) ≈ 3 atol=atol rtol=rtol + end + + x = Variable(5, 5) + p = minimize(x[1, 1], dot(*)(3,x) >= 3) + if test + @test vexity(p) == AffineVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 1 atol=atol rtol=rtol + @test (evaluate(x[1, 1]))[1] ≈ 1 atol=atol rtol=rtol + end + + x = Variable(3,1) + p = minimize(sum(dot(*)(ones(3,3), x)), x>=1) + if test + @test vexity(p) == AffineVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 9 atol=atol rtol=rtol + @test (evaluate(x[1, 1]))[1] ≈ 1 atol=atol rtol=rtol + end + + x = Variable(1,3) + p = minimize(sum(dot(*)(ones(3,3), x)), x>=1) + if test + @test vexity(p) == AffineVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 9 atol=atol rtol=rtol + @test (evaluate(x[1, 1]))[1] ≈ 1 atol=atol rtol=rtol + end + + x = Variable(1, 3, Positive()) + p = maximize(sum(dot(/)(x,[1 2 3])), x<=1) + if test + @test vexity(p) == AffineVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 11 / 6 atol=atol rtol=rtol + @test evaluate(sum((dot(/))(x, [1 2 3]))) ≈ 11 / 6 atol=atol rtol=rtol + end + + # Broadcast fusion works + x = Variable(5, 5) + a = 2.0 .* x .* ones(Int, 5) + if test + @test a isa Convex.DotMultiplyAtom + end +end + +@add_problem affine function affine_reshape_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + A = rand(2, 3) + X = Variable(3, 2) + c = rand() + p = minimize(sum(reshape(X, 2, 3) + A), X >= c) + if test + @test vexity(p) == AffineVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ sum(A .+ c) atol=atol rtol=rtol + @test evaluate(sum(reshape(X, 2, 3) + A)) ≈ sum(A .+ c) atol=atol rtol=rtol + end + + b = rand(6) + p = minimize(sum(vec(X) + b), X >= c) + if test + @test vexity(p) == AffineVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ sum(b .+ c) atol=atol rtol=rtol + @test evaluate(sum(vec(X) + b)) ≈ sum(b .+ c) atol=atol rtol=rtol + end + + x = Variable(4, 4) + c = ones(16, 1) + reshaped = reshape(x, 16, 1) + a = collect(1:16) + p = minimize(c' * reshaped, reshaped >= a) + if test + @test vexity(p) == AffineVexity() + end + handle_problem!(p) + # TODO: why is accuracy lower here? + if test + @test p.optval ≈ 136 atol=10atol atol=atol rtol=rtol + @test (evaluate(c' * reshaped))[1] ≈ 136 atol=10atol atol=atol rtol=rtol + end +end + +@add_problem affine function affine_hcat_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + x = Variable(4, 4) + y = Variable(4, 6) + p = maximize(sum(x) + sum([y fill(4.0, 4)]), [x y fill(2.0, (4, 2))] <= 2) + if test + @test vexity(p) == AffineVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 96 atol=atol rtol=rtol + @test evaluate(sum(x) + sum([y fill(4.0, 4)])) ≈ 96 atol=atol rtol=rtol + @test evaluate([x y fill(2.0, (4, 2))]) ≈ fill(2.0, (4, 12)) atol=atol rtol=rtol + end +end + +@add_problem affine function affine_vcat_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + x = Variable(4, 4) + y = Variable(4, 6) + + # TODO: fix dimension mismatch [y 4*eye(4); x -ones(4, 6)] + p = maximize(sum(x) + sum([y 4*eye(4); x -ones(4, 6)]), [x;y'] <= 2) + if test + @test vexity(p) == AffineVexity() + end + handle_problem!(p) + # TODO: why is accuracy lower here? + if test + @test p.optval ≈ 104 atol=10atol atol=atol rtol=rtol + @test evaluate(sum(x) + sum([y 4 * eye(4); x -(ones(4, 6))])) ≈ 104 atol=10atol atol=atol rtol=rtol + @test evaluate([x; y']) ≈ 2 * ones(10, 4) atol=atol rtol=rtol + end +end + +@add_problem affine function affine_Diagonal_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + x = Variable(2, 2) + if test + @test_throws ArgumentError Diagonal(x) + end + + x = Variable(4) + p = minimize(sum(Diagonal(x)), x == [1, 2, 3, 4]) + if test + @test vexity(p) == AffineVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 10 atol=atol rtol=rtol + @test all(abs.(evaluate(Diagonal(x)) - Diagonal([1, 2, 3, 4])) .<= atol) + end + + x = Variable(3) + c = [1; 2; 3] + p = minimize(c' * Diagonal(x) * c, x >= 1, sum(x) == 10) + if test + @test vexity(p) == AffineVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 21 atol=atol rtol=rtol + end + + x = Variable(3) + p = minimize(sum(x), x >= 1, Diagonal(x)[1, 2] == 1) + output = handle_problem!(p) + if test + @test output === nothing + @test p.status != :Optimal + end +end + +@add_problem affine function affine_conv_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + x = Variable(3) + h = [1, -1] + p = minimize(sum(conv(h, x)) + sum(x), x >= 1, x <= 2) + if test + @test vexity(p) == AffineVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 3 atol=atol rtol=rtol + @test evaluate(sum(conv(h, x))) ≈ 0 atol=atol rtol=rtol + end + + x = Variable(3) + h = [1, -1] + p = minimize(sum(conv(x, h)) + sum(x), x >= 1, x <= 2) + if test + @test vexity(p) == AffineVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 3 atol=atol rtol=rtol + @test evaluate(sum(conv(h, x))) ≈ 0 atol=atol rtol=rtol + end + +end + +@add_problem affine function affine_satisfy_problems(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + x = Variable() + p = satisfy(x >= 0) + add_constraints!(p, x >= 1) + add_constraints!(p, [x >= -1, x <= 4]) + handle_problem!(p) + if test + @test p.status == :Optimal + end + + p = satisfy([x >= 0, x >= 1, x <= 2]) + handle_problem!(p) + if test + @test p.status == :Optimal + end + + p = maximize(1, [x >= 1, x <= 2]) + handle_problem!(p) + if test + @test p.status == :Optimal + end + + constr = x >= 0 + constr += x >= 1 + constr += x <= 10 + constr2 = x >= 0 + constr2 += [x >= 2, x <= 3] + constr + p = satisfy(constr) + handle_problem!(p) + if test + @test p.status == :Optimal + end +end + +@add_problem affine function affine_dual(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + x = Variable() + p = minimize(x, x >= 0) + handle_problem!(p) + if test + if p.solution.has_dual + println("Solution object has dual value, checking for dual correctness.") + @test p.constraints[1].dual ≈ 1 atol=atol rtol=rtol + end + end + + x = Variable() + p = maximize(x, x <= 0) + handle_problem!(p) + if test + if p.solution.has_dual + println("Solution object has dual value, checking for dual correctness.") + @test p.constraints[1].dual ≈ 1 atol=atol rtol=rtol + end + end + + x = Variable() + p = minimize(x, x >= 0, x == 2) + handle_problem!(p) + if test + if p.solution.has_dual + println("Solution object has dual value, checking for dual correctness.") + @test p.constraints[1].dual ≈ 0 atol=atol rtol=rtol + @test abs.(p.constraints[2].dual) ≈ 1 atol=atol rtol=rtol + end + end + + x = Variable(2) + A = 1.5 * eye(2) + p = minimize(dot([2.0; 2.0], x), [A * x >= [1.1; 1.1]]) + handle_problem!(p) + if test + if p.solution.has_dual + println("Solution object has dual value, checking for dual correctness.") + dual = [4/3; 4/3] + @test all(abs.(p.constraints[1].dual - dual) .<= atol) + end + end +end + +@add_problem affine function affine_Partial_transpose(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + dims = [2,3,4] + d = prod(dims) + A = rand(ComplexF64,2,2) + B = rand(ComplexF64,3,3) + C = rand(ComplexF64,4,4) + M = kron(A,B,C) + Mt1 = kron(transpose(A),B,C) + Mt2 = kron(A,transpose(B),C) + Mt3 = kron(A,B,transpose(C)) + + Rt1 = ComplexVariable(d,d) + Rt2 = ComplexVariable(d,d) + Rt3 = ComplexVariable(d,d) + S = rand(ComplexF64,d,d) + handle_problem!(satisfy(partialtranspose(Rt1, 1, dims) == S )) + handle_problem!(satisfy(partialtranspose(Rt2, 2, dims) == S )) + handle_problem!(satisfy(partialtranspose(Rt3, 3, dims) == S )) + + + if test + @test partialtranspose(M,1,dims) ≈ Mt1 atol=atol rtol=rtol + @test partialtranspose(M,2,dims) ≈ Mt2 atol=atol rtol=rtol + @test partialtranspose(M,3,dims) ≈ Mt3 atol=atol rtol=rtol + @test partialtranspose(S,1,dims) ≈ evaluate(Rt1) atol=atol rtol=rtol + @test partialtranspose(S,2,dims) ≈ evaluate(Rt2) atol=atol rtol=rtol + @test partialtranspose(S,3,dims) ≈ evaluate(Rt3) atol=atol rtol=rtol + end + + if test + @test_throws ArgumentError partialtrace(rand(6, 6), 3, [2, 3]) + @test_throws ArgumentError partialtrace(rand(6, 6), 1, [2, 4]) + @test_throws ArgumentError partialtrace(rand(3, 4), 1, [2, 3]) + end +end + +@add_problem affine function affine_permuteddims_matrix(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} +#this function is used in the partial transpose + for n in (2, 3, 4, 5) + dims = ntuple( i -> rand(2:5), n) + d = prod(dims) + v = rand(d) + p = randperm(n) + out1 = vec(permutedims(reshape(v, dims), p)) + out2 = Convex.permutedims_matrix(dims, p) * v + if test + @test out1 ≈ out2 atol=atol rtol=rtol + end + end +end diff --git a/benchmark/benchmarks/constraint_benchmarks.jl b/src/problem_depot/problems/benchmark.jl similarity index 50% rename from benchmark/benchmarks/constraint_benchmarks.jl rename to src/problem_depot/problems/benchmark.jl index e00bb377d..bcf2a295e 100644 --- a/benchmark/benchmarks/constraint_benchmarks.jl +++ b/src/problem_depot/problems/benchmark.jl @@ -1,64 +1,72 @@ -@add_benchmark function LT_constraints(handle_problem) +@add_problem constraints_benchmark function LT_constraints(handle_problem!, args...) p = satisfy() x = [Variable() for _ = 1:1000] for (i, xi) in enumerate(x) push!(p.constraints, xi <= 1.0 * i) end - return handle_problem(p) + handle_problem!(p) + nothing end -@add_benchmark function LT_constraint(handle_problem) +@add_problem constraints_benchmark function LT_constraint(handle_problem!, args...) p = satisfy() x = Variable(1000) push!(p.constraints, x <= collect(1.0:1000.0)) - return handle_problem(p) + handle_problem!(p) + nothing end -@add_benchmark function GT_constraints(handle_problem) +@add_problem constraints_benchmark function GT_constraints(handle_problem!, args...) p = satisfy() x = [Variable() for _ = 1:1000] for (i, xi) in enumerate(x) push!(p.constraints, xi >= 1.0 * i) end - return handle_problem(p) + handle_problem!(p) + nothing end -@add_benchmark function GT_constraint(handle_problem) +@add_problem constraints_benchmark function GT_constraint(handle_problem!, args...) p = satisfy() x = Variable(1000) push!(p.constraints, x >= collect(1.0:1000.0)) - return handle_problem(p) + handle_problem!(p) + nothing end -@add_benchmark function equality_constraints(handle_problem) +@add_problem constraints_benchmark function equality_constraints(handle_problem!, args...) p = satisfy() x = [Variable() for _ = 1:1000] for (i, xi) in enumerate(x) push!(p.constraints, xi == 1.0 * i) end - return handle_problem(p) + handle_problem!(p) + nothing end -@add_benchmark function equality_constraint(handle_problem) +@add_problem constraints_benchmark function equality_constraint(handle_problem!, args...) p = satisfy() x = Variable(1000) push!(p.constraints, x == collect(1.0:1000.0)) - return handle_problem(p) + handle_problem!(p) + nothing end -@add_benchmark function SDP_constraint(handle_problem) +@add_problem constraints_benchmark function sdp_constraint(handle_problem!, args...) p = satisfy() x = Variable(44, 44) # 990 vectorized entries push!(p.constraints, x ⪰ 0) - return handle_problem(p) + handle_problem!(p) + nothing end -@add_benchmark function SDP_constraints(handle_problem) +@add_problem constraints_benchmark function sdp_constraints(handle_problem!, args...) p = satisfy() x = [Variable(4, 4) for _ = 1:100] # 1000 total vectorized entries for v in x push!(p.constraints, v ⪰ 0) end - return handle_problem(p) + handle_problem!(p) + nothing end diff --git a/src/problem_depot/problems/constant.jl b/src/problem_depot/problems/constant.jl new file mode 100644 index 000000000..2f7bca287 --- /dev/null +++ b/src/problem_depot/problems/constant.jl @@ -0,0 +1,138 @@ +@add_problem constant function constant_Issue_166(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + # Issue #166 + α = Variable(5) + fix!(α, ones(5,1)) + + # has const vexity, but not at the head + c = (rand(5,5) * α) * ones(1,5) + + β = Variable(5) + β.value = ones(5) + + problem = minimize(sum(c * β), [β >= 0]) + handle_problem!(problem) + if test + @test problem.optval ≈ evaluate(sum(c * β)) atol=atol rtol=rtol + @test problem.optval ≈ 0.0 atol=atol rtol=rtol + @test β.value ≈ zeros(5) atol=atol rtol=rtol + end +end + +@add_problem constant function constant_Issue_228(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + x = Variable(2) + y = Variable(2) + fix!(x, [1 1]') + prob = minimize(y'*(x+[2 2]'), [y>=0]) + handle_problem!(prob) + if test + @test prob.optval ≈ 0.0 atol=atol rtol=rtol + end + + prob = minimize(x'*y, [y>=0]) + handle_problem!(prob) + if test + @test prob.optval ≈ 0.0 atol=atol rtol=rtol + end +end + +@add_problem constant function constant_Test_double_fix!(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + x = Variable() + y = Variable() + fix!(x, 1.0) + prob = minimize(y*x, [y >= x, x >= 0.5]) + handle_problem!(prob) + if test + @test prob.optval ≈ 1.0 atol=atol rtol=rtol + end + + fix!(x, 2.0) + handle_problem!(prob) + if test + @test prob.optval ≈ 4.0 atol=atol rtol=rtol + end + + free!(x) + fix!(y, 1.0) + handle_problem!(prob) + if test + @test prob.optval ≈ 0.5 atol=atol rtol=rtol + end +end + +@add_problem constant function constant_fix!_and_multiply(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + p = Variable() + fix!(p, 1.0) + x = Variable(2,2) + prob = minimize( tr(p*x), [ x >= 1 ]) + handle_problem!(prob) + if test + @test prob.optval ≈ 2.0 atol=atol rtol=rtol + @test evaluate( tr(p*x) ) ≈ 2.0 atol=atol rtol=rtol + end +end + +@add_problem constant function constant_fix!_with_complex_numbers(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + x = ComplexVariable() + fix!(x, 1.0 + im*1.0) + y = Variable() + prob = minimize( real(x*y), [ y >= .5, real(x) >= .5, imag(x) >= 0]) + handle_problem!(prob) + if test + @test prob.optval ≈ .5 atol=atol rtol=rtol + @test evaluate(real(x*y)) ≈ .5 atol=atol rtol=rtol + @test evaluate(y) ≈ 0.5 atol=atol rtol=rtol + end + + free!(x) + if test + fix!(y) + else # if we haven't solved the problem, we need to fix! it to a number + fix!(y, 0.5) + end + handle_problem!(prob) + if test + @test prob.optval ≈ 0.25 atol=atol rtol=rtol + @test evaluate(real(x*y)) ≈ 0.25 atol=atol rtol=rtol + @test real(evaluate(x)) ≈ 0.5 atol=atol rtol=rtol + @test evaluate(y) ≈ 0.5 atol=atol rtol=rtol + end + + if test + @test_throws DimensionMismatch fix!(x, rand(2)) + @test_throws DimensionMismatch fix!(x, rand(2,2)) + end +end + +@add_problem constant function constant_fix!_with_vectors(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + x = ComplexVariable(5) + fix!(x, ones(5) + im*ones(5)) + y = Variable() + prob = minimize( real(y*sum(x)), [ y >= .5, real(x) >= .5, imag(x) >= 0]) + handle_problem!(prob) + if test + @test prob.optval ≈ 2.5 atol=atol rtol=rtol + @test evaluate(real(y*sum(x))) ≈ 2.5 atol=atol rtol=rtol + @test evaluate(y) ≈ 0.5 atol=atol rtol=rtol + end + + free!(x) + if test + fix!(y) + else # if we haven't solved the problem, we need to fix! it to a number + fix!(y, 0.5) + end + + handle_problem!(prob) + if test + @test prob.optval ≈ 1.25 atol=atol rtol=rtol + @test evaluate(real(y*sum(x))) ≈ 1.25 atol=atol rtol=rtol + @test real(evaluate(x)) ≈ 0.5*ones(5) atol=atol rtol=rtol + @test evaluate(y) ≈ 0.5 atol=atol rtol=rtol + end + + if test + @test_throws DimensionMismatch fix!(x, rand(5,5)) + @test_throws DimensionMismatch fix!(x, rand(4)) + end + +end diff --git a/src/problem_depot/problems/exp.jl b/src/problem_depot/problems/exp.jl new file mode 100644 index 000000000..b74f85896 --- /dev/null +++ b/src/problem_depot/problems/exp.jl @@ -0,0 +1,140 @@ +@add_problem exp function exp_exp_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + y = Variable() + p = minimize(exp(y), y>=0) + if test + @test vexity(p) == ConvexVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 1 atol=atol rtol=rtol + @test evaluate(exp(y)) ≈ 1 atol=atol rtol=rtol + end + + y = Variable() + p = minimize(exp(y), y>=1) + if test + @test vexity(p) == ConvexVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ exp(1) atol=atol rtol=rtol + @test evaluate(exp(y)) ≈ exp(1) atol=atol rtol=rtol + end + + y = Variable(5) + p = minimize(sum(exp(y)), y>=0) + if test + @test vexity(p) == ConvexVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 5 atol=atol rtol=rtol + @test evaluate(sum(exp(y))) ≈ 5 atol=atol rtol=rtol + end + + y = Variable(5) + p = minimize(sum(exp(y)), y>=0) + if test + @test vexity(p) == ConvexVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 5 atol=atol rtol=rtol + end +end + +@add_problem exp function exp_log_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + y = Variable() + p = maximize(log(y), y<=1) + if test + @test vexity(p) == ConvexVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 0 atol=atol rtol=rtol + end + + y = Variable() + p = maximize(log(y), y<=2) + if test + @test vexity(p) == ConvexVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ log(2) atol=atol rtol=rtol + end + + y = Variable() + p = maximize(log(y), [y<=2, exp(y)<=10]) + if test + @test vexity(p) == ConvexVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ log(2) atol=atol rtol=rtol + end +end + +@add_problem exp function exp_log_sum_exp_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + y = Variable(5) + p = minimize(logsumexp(y), y>=1) + if test + @test vexity(p) == ConvexVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ log(exp(1) * 5) atol=atol rtol=rtol + end +end + +@add_problem exp function exp_logistic_loss_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + y = Variable(5) + p = minimize(logisticloss(y), y>=1) + if test + @test vexity(p) == ConvexVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ log(exp(1) + 1) * 5 atol=atol rtol=rtol + end +end + +@add_problem exp function exp_entropy_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + y = Variable(5, Positive()) + p = maximize(entropy(y), sum(y)<=1) + if test + @test vexity(p) == ConvexVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ -(log(1 / 5)) atol=atol rtol=rtol + end +end + +@add_problem exp function exp_relative_entropy_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + x = Variable(1) + y = Variable(1) + # x log (x/y) + p = minimize(relative_entropy(x,y), y==1, x >= 2) + if test + @test vexity(p) == ConvexVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 2 * log(2) atol=atol rtol=rtol + end +end + +@add_problem exp function exp_log_perspective_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + x = Variable(1) + y = Variable(1) + # y log (x/y) + p = maximize(log_perspective(x,y), y==5, x <= 10) + if test + @test vexity(p) == ConvexVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 5 * log(2) atol=atol rtol=rtol + end +end diff --git a/src/problem_depot/problems/lp.jl b/src/problem_depot/problems/lp.jl new file mode 100644 index 000000000..367039e98 --- /dev/null +++ b/src/problem_depot/problems/lp.jl @@ -0,0 +1,243 @@ +@add_problem lp function lp_dual_abs_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + x = Variable() + p = minimize(abs(x), x<=-1) + if test + @test vexity(p) == ConvexVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 1 atol=atol rtol=rtol + @test evaluate(abs(x)) ≈ 1 atol=atol rtol=rtol + end + + x = Variable(2,2) + p = minimize(sum(abs(x)), x[2,2]>=1, x[1,1]>=1, x>=0) + if test + @test vexity(p) == ConvexVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 2 atol=atol rtol=rtol + @test evaluate(sum(abs(x))) ≈ 2 atol=atol rtol=rtol + @test p.constraints[1].dual ≈ 1 atol=atol rtol=rtol + @test p.constraints[2].dual ≈ 1 atol=atol rtol=rtol + @test p.constraints[3].dual[1,1] ≈ 0 atol=atol rtol=rtol + @test p.constraints[3].dual[2,2] ≈ 0 atol=atol rtol=rtol + @test p.constraints[3].dual[1,2] ≈ p.constraints[3].dual[2,1] atol=atol rtol=rtol + end +end + +@add_problem lp function lp_maximum_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + x = Variable(10) + a = shuffle(collect(0.1:0.1:1.0)) + p = minimize(maximum(x), x >= a) + if test + @test vexity(p) == ConvexVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ maximum(a) atol=atol rtol=rtol + @test evaluate(maximum(x)) ≈ maximum(a) atol=atol rtol=rtol + end +end + +@add_problem lp function lp_minimum_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + x = Variable(1) + a = reshape(shuffle(collect(0.01:0.01:1.0)), (10, 10)) + p = maximize(minimum(x), x <= a) + if test + @test vexity(p) == ConvexVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ minimum(a) atol=atol rtol=rtol + @test evaluate(minimum(x)) ≈ minimum(a) atol=atol rtol=rtol + end + + x = Variable(4, 4) + y = Variable(4, 6) + z = Variable(1) + c = ones(4, 1) + d = fill(2.0, (6, 1)) + constraints = [[x y] <= 2, z <= 0, z <= x, 2z >= -1] + objective = sum(x + z) + minimum(y) + c' * y * d + p = maximize(objective, constraints) + if test + @test vexity(p) == ConvexVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 130 atol=atol rtol=rtol + @test (evaluate(objective))[1] ≈ 130 atol=atol rtol=rtol + end +end + +@add_problem lp function lp_max_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + x = Variable(10, 10) + y = Variable(10, 10) + a = reshape(shuffle(collect(0.01:0.01:1.0)), (10, 10)) + b = reshape(shuffle(collect(0.01:0.01:1.0)), (10, 10)) + p = minimize(maximum(max(x, y)), [x >= a, y >= b]) + if test + @test vexity(p) == ConvexVexity() + end + handle_problem!(p) + max_a = maximum(a) + max_b = maximum(b) + if test + @test p.optval ≈ max(max_a, max_b) atol=10atol atol=atol rtol=rtol + @test evaluate(maximum(max(x, y))) ≈ max(max_a, max_b) atol=10atol atol=atol rtol=rtol + end +end + +@add_problem lp function lp_min_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + x = Variable(10, 10) + y = Variable(10, 10) + a = reshape(shuffle(collect(0.01:0.01:1.0)), (10, 10)) + b = reshape(shuffle(collect(0.01:0.01:1.0)), (10, 10)) + p = maximize(minimum(min(x, y)), [x <= a, y <= b]) + if test + @test vexity(p) == ConvexVexity() + end + handle_problem!(p) + min_a = minimum(a) + min_b = minimum(b) + if test + @test p.optval ≈ min(min_a, min_b) atol=10atol atol=atol rtol=rtol + @test evaluate(minimum(min(x, y))) ≈ min(min_a, min_b) atol=10atol atol=atol rtol=rtol + end +end + +@add_problem lp function lp_pos_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + x = Variable(3) + a = [-2; 1; 2] + p = minimize(sum(pos(x)), [x >= a, x <= 2]) + if test + @test vexity(p) == ConvexVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 3 atol=atol rtol=rtol + @test evaluate(sum(pos(x))) ≈ 3 atol=atol rtol=rtol + end +end + +@add_problem lp function lp_neg_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + x = Variable(3) + p = minimize(1, [x >= -2, x <= -2, neg(x) <= 3]) + if test + @test vexity(p) == ConvexVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 1 atol=atol rtol=rtol + @test evaluate(sum(neg(x))) ≈ 6 atol=atol rtol=rtol + end +end + +@add_problem lp function lp_sumlargest_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + x = Variable(2) + p = minimize(sumlargest(x, 2), x >= [1; 1]) + if test + @test vexity(p) == ConvexVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 2 atol=atol rtol=rtol + @test evaluate(sumlargest(x, 2)) ≈ 2 atol=atol rtol=rtol + end + + x = Variable(4, 4) + p = minimize(sumlargest(x, 3), x >= eye(4), x[1, 1] >= 1.5, x[2, 3] >= 2.1) + if test + @test vexity(p) == ConvexVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 4.6 atol=atol rtol=rtol + @test evaluate(sumlargest(x, 2)) ≈ 3.6 atol=atol rtol=rtol + end +end + +@add_problem lp function lp_sumsmallest_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + x = Variable(4, 4) + p = minimize(sumlargest(x, 2), sumsmallest(x, 4) >= 1) + if test + @test vexity(p) == ConvexVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 0.5 atol=atol rtol=rtol + @test evaluate(sumsmallest(x, 4)) ≈ 1 atol=atol rtol=rtol + end + + x = Variable(3, 2) + p = maximize(sumsmallest(x, 3), x >= 2, x <= 5, sumlargest(x, 3) <= 12) + if test + @test vexity(p) == ConvexVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 12 atol=atol rtol=rtol + @test evaluate(sumsmallest(x, 3)) ≈ 12 atol=atol rtol=rtol + end +end + +@add_problem lp function lp_dotsort_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + x = Variable(4, 1) + p = minimize(dotsort(x, [1, 2, 3, 4]), sum(x) >= 7, x >= 0, x <= 2, x[4] <= 1) + if test + @test vexity(p) == ConvexVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 19 atol=atol rtol=rtol + @test vec(x.value) ≈ [2; 2; 2; 1] atol=atol rtol=rtol + @test evaluate(dotsort(x, [1, 2, 3, 4])) ≈ 19 atol=atol rtol=rtol + end + + x = Variable(2, 2) + p = minimize(dotsort(x, [1 2; 3 4]), sum(x) >= 7, x >= 0, x <= 2, x[2, 2] <= 1) + if test + @test vexity(p) == ConvexVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 19 atol=atol rtol=rtol + @test evaluate(dotsort(x, [1, 2, 3, 4])) ≈ 19 atol=atol rtol=rtol + end +end + +@add_problem lp function lp_hinge_loss_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + # TODO: @davidlizeng. We should finish this someday. +end + +@add_problem lp function lp_dual_norm_inf_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + x = Variable(3) + p = minimize(norm_inf(x), [-2 <= x, x <= 1]) + if test + @test vexity(p) == ConvexVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 0 atol=atol rtol=rtol + @test evaluate(norm_inf(x)) ≈ 0 atol=atol rtol=rtol + @test norm(p.constraints[1].dual) ≈ 0 atol=atol rtol=rtol + @test norm(p.constraints[2].dual) ≈ 0 atol=atol rtol=rtol + end +end + +@add_problem lp function lp_dual_norm_1_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + x = Variable(3) + p = minimize(norm_1(x), [-2 <= x, x <= 1]) + if test + @test vexity(p) == ConvexVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 0 atol=atol rtol=rtol + @test evaluate(norm_1(x)) ≈ 0 atol=atol rtol=rtol + @test norm(p.constraints[1].dual) ≈ 0 atol=atol rtol=rtol + @test norm(p.constraints[2].dual) ≈ 0 atol=atol rtol=rtol + end +end diff --git a/src/problem_depot/problems/mip.jl b/src/problem_depot/problems/mip.jl new file mode 100644 index 000000000..2d454aaeb --- /dev/null +++ b/src/problem_depot/problems/mip.jl @@ -0,0 +1,106 @@ +@add_problem mip function mip_lp_fallback_interface(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + x = Variable() + p = minimize(x, x>=4.3) + if test + @test vexity(p) == AffineVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 4.3 atol=atol rtol=rtol + end + + x = Variable(2) + p = minimize(norm(x,1), x[1]>=4.3) + if test + @test vexity(p) == ConvexVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 4.3 atol=atol rtol=rtol + end +end + +@add_problem mip function mip_integer_variables(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + x = Variable(:Int) + p = minimize(x, x>=4.3) + if test + @test vexity(p) == AffineVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 5 atol=atol rtol=rtol + end + + x = Variable(2, :Int) + p = minimize(sum(x), x>=4.3) + if test + @test vexity(p) == AffineVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 10 atol=atol rtol=rtol + end + + x = Variable(:Int) + y = Variable() + p = minimize(sum(x + y), x>=4.3, y>=7) + if test + @test vexity(p) == AffineVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 12 atol=atol rtol=rtol + end + + x = Variable(2, :Int) + p = minimize(norm(x, 1), x[1]>=4.3) + if test + @test vexity(p) == ConvexVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 5 atol=atol rtol=rtol + end + + x = Variable(2, :Int) + p = minimize(sum(x), x[1]>=4.3, x>=0) + if test + @test vexity(p) == AffineVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 5 atol=atol rtol=rtol + end + + x = Variable(2, :Int) + p = minimize(sum(x), x>=.5) + if test + @test vexity(p) == AffineVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 2 atol=atol rtol=rtol + end +end + +@add_problem mip function mip_binary_variables(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + x = Variable(2, :Bin) + p = minimize(sum(x), x>=.5) + if test + @test vexity(p) == AffineVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 2 atol=atol rtol=rtol + end + + x = Variable(2, :Bin) + p = minimize(sum(x), x[1]>=.5, x>=0) + if test + @test vexity(p) == AffineVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 1 atol=atol rtol=rtol + end +end diff --git a/src/problem_depot/problems/sdp.jl b/src/problem_depot/problems/sdp.jl new file mode 100644 index 000000000..0c84be4da --- /dev/null +++ b/src/problem_depot/problems/sdp.jl @@ -0,0 +1,394 @@ +# TODO: uncomment vexity checks once SDP on vars/constraints changes vexity of problem +@add_problem sdp function sdp_sdp_variables(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + y = Variable((2,2), :Semidefinite) + p = minimize(y[1,1]) + # @fact vexity(p) --> ConvexVexity() + handle_problem!(p) + if test + @test p.optval ≈ 0 atol=atol rtol=rtol + end + + y = Variable((3,3), :Semidefinite) + p = minimize(y[1,1], y[2,2]==1) + # @fact vexity(p) --> ConvexVexity() + handle_problem!(p) + if test + @test p.optval ≈ 0 atol=atol rtol=rtol + end + + # Solution is obtained as y[2,2] -> infinity + # This test fails on Mosek. See + # https://github.com/JuliaOpt/Mosek.jl/issues/29 + # y = Variable((2, 2), :Semidefinite) + # p = minimize(y[1, 1], y[1, 2] == 1) + # # @fact vexity(p) --> ConvexVexity() + # handle_problem!(p) + # @fact p.optval --> roughly(0, atol) + + y = Semidefinite(3) + p = minimize(sum(diag(y)), y[1, 1] == 1) + # @fact vexity(p) --> ConvexVexity() + handle_problem!(p) + if test + @test p.optval ≈ 1 atol=atol rtol=rtol + end + + y = Variable((3, 3), :Semidefinite) + p = minimize(tr(y), y[2,1]<=4, y[2,2]>=3) + # @fact vexity(p) --> ConvexVexity() + handle_problem!(p) + if test + @test p.optval ≈ 3 atol=atol rtol=rtol + end + + x = Variable(Positive()) + y = Semidefinite(3) + p = minimize(y[1, 2], y[2, 1] == 1) + # @fact vexity(p) --> ConvexVexity() + handle_problem!(p) + if test + @test p.optval ≈ 1 atol=atol rtol=rtol + end +end + +@add_problem sdp function sdp_sdp_constraints(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + # This test fails on Mosek + x = Variable(Positive()) + y = Variable((3, 3)) + p = minimize(x + y[1, 1], isposdef(y), x >= 1, y[2, 1] == 1) + # @fact vexity(p) --> ConvexVexity() + handle_problem!(p) + if test + @test p.optval ≈ 1 atol=atol rtol=rtol + end +end + +@add_problem sdp function sdp_nuclear_norm_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + y = Semidefinite(3) + p = minimize(nuclearnorm(y), y[2,1]<=4, y[2,2]>=3, y[3,3]<=2) + if test + @test vexity(p) == ConvexVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 3 atol=atol rtol=rtol + @test evaluate(nuclearnorm(y)) ≈ 3 atol=atol rtol=rtol + end +end + +@add_problem sdp function sdp_operator_norm_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + y = Variable((3,3)) + p = minimize(opnorm(y), y[2,1]<=4, y[2,2]>=3, sum(y)>=12) + if test + @test vexity(p) == ConvexVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 4 atol=atol rtol=rtol + @test evaluate(opnorm(y)) ≈ 4 atol=atol rtol=rtol + end +end + +@add_problem sdp function sdp_sigma_max_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + y = Variable((3,3)) + p = minimize(sigmamax(y), y[2,1]<=4, y[2,2]>=3, sum(y)>=12) + if test + @test vexity(p) == ConvexVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 4 atol=atol rtol=rtol + @test evaluate(sigmamax(y)) ≈ 4 atol=atol rtol=rtol + end +end + +@add_problem sdp function sdp_lambda_max_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + y = Semidefinite(3) + p = minimize(lambdamax(y), y[1,1]>=4) + if test + @test vexity(p) == ConvexVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 4 atol=atol rtol=rtol + @test evaluate(lambdamax(y)) ≈ 4 atol=atol rtol=rtol + end +end + +@add_problem sdp function sdp_lambda_min_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + y = Semidefinite(3) + p = maximize(lambdamin(y), tr(y)<=6) + if test + @test vexity(p) == ConvexVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 2 atol=atol rtol=rtol + @test evaluate(lambdamin(y)) ≈ 2 atol=atol rtol=rtol + end +end + +@add_problem sdp function sdp_matrix_frac_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + x = [1, 2, 3] + P = Variable(3, 3) + p = minimize(matrixfrac(x, P), P <= 2*eye(3), P >= 0.5 * eye(3)) + if test + @test vexity(p) == ConvexVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 7 atol=atol rtol=rtol + @test (evaluate(matrixfrac(x, P)))[1] ≈ 7 atol=atol rtol=rtol + end +end + +@add_problem sdp function sdp_matrix_frac_atom_both_arguments_variable(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + x = Variable(3) + P = Variable(3, 3) + p = minimize(matrixfrac(x, P), lambdamax(P) <= 2, x[1] >= 1) + if test + @test vexity(p) == ConvexVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 0.5 atol=atol rtol=rtol + @test (evaluate(matrixfrac(x, P)))[1] ≈ 0.5 atol=atol rtol=rtol + end +end + +@add_problem sdp function sdp_sum_largest_eigs(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + x = Semidefinite(3) + p = minimize(sumlargesteigs(x, 2), x >= 1) + handle_problem!(p) + + if test + @test p.optval ≈ 3 atol=atol rtol=rtol + @test evaluate(x) ≈ ones(3, 3) atol=atol rtol=rtol + end + + x = Semidefinite(3) + p = minimize(sumlargesteigs(x, 2), [x[i,:] >= i for i=1:3]...) + handle_problem!(p) + + if test + @test p.optval ≈ 8.4853 atol=atol rtol=rtol + end + + x1 = Semidefinite(3) + p1 = minimize(lambdamax(x1), x1[1,1]>=4) + handle_problem!(p1) + + x2 = Semidefinite(3) + p2 = minimize(sumlargesteigs(x2, 1), x2[1,1]>=4) + handle_problem!(p2) + + if test + @test p1.optval ≈ p2.optval atol=atol rtol=rtol + end + + x1 = Semidefinite(3) + p1 = minimize(lambdamax(x1), [x1[i,:] >= i for i=1:3]...) + handle_problem!(p1) + + x2 = Semidefinite(3) + p2 = minimize(sumlargesteigs(x2, 1), [x2[i,:] >= i for i=1:3]...) + handle_problem!(p2) + + if test + @test p1.optval ≈ p2.optval atol=atol rtol=rtol + end + +end + +@add_problem sdp function sdp_kron_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + id = eye(4) + X = Semidefinite(4) + W = kron(id, X) + p = maximize(tr(W), tr(X) ≤ 1) + if test + @test vexity(p) == AffineVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 4 atol=atol rtol=rtol + end +end + +@add_problem sdp function sdp_Partial_trace(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + A = Semidefinite(2) + B = [1 0; 0 0] + ρ = kron(B, A) + constraints = [partialtrace(ρ, 1, [2; 2]) == [0.09942819 0.29923607; 0.29923607 0.90057181], ρ in :SDP] + p = satisfy(constraints) + handle_problem!(p) + if test + @test evaluate(ρ) ≈ [0.09942819 0.29923607 0 0; 0.299237 0.900572 0 0; 0 0 0 0; 0 0 0 0] atol=atol rtol=rtol + @test evaluate(partialtrace(ρ, 1, [2; 2])) ≈ [0.09942819 0.29923607; 0.29923607 0.90057181] atol=atol rtol=rtol + end + + function rand_normalized(n) + A = 5*randn(n, n) + im*5*randn(n, n) + A / tr(A) + end + + As = [ rand_normalized(3) for _ = 1:5] + Bs = [ rand_normalized(2) for _ = 1:5] + p = rand(5) + + AB = sum(i -> p[i]*kron(As[i],Bs[i]), 1:5) + if test + @test partialtrace(AB, 2, [3, 2]) ≈ sum( p .* As ) atol=atol rtol=rtol + @test partialtrace(AB, 1, [3, 2]) ≈ sum( p .* Bs ) atol=atol rtol=rtol + end + + A, B, C = rand(5,5), rand(4,4), rand(3,3) + ABC = kron(kron(A, B), C) + if test + @test kron(A,B)*tr(C) ≈ partialtrace(ABC, 3, [5, 4, 3]) atol=atol rtol=rtol + end + + # Test 281 + A = rand(6,6) + expr = partialtrace(Constant(A), 1, [2, 3]) + if test + @test size(expr) == size(evaluate(expr)) + + @test_throws ArgumentError partialtrace(rand(6, 6), 3, [2, 3]) + @test_throws ArgumentError partialtrace(rand(6, 6), 1, [2, 4]) + @test_throws ArgumentError partialtrace(rand(3, 4), 1, [2, 3]) + end +end + + +@add_problem sdp function sdp_Real_Variables_with_complex_equality_constraints(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + n = 10 # variable dimension (parameter) + m = 5 # number of constraints (parameter) + xo = rand(n) + A = randn(m,n) + im*randn(m,n) + b = A * xo + x = Variable(n) + p1 = minimize(sum(x), A*x == b, x>=0) + handle_problem!(p1) + x1 = x.value + + p2 = minimize(sum(x), real(A)*x == real(b), imag(A)*x==imag(b), x>=0) + handle_problem!(p2) + x2 = x.value + if test + @test x1 == x2 + end +end + +@add_problem sdp function sdp_Complex_Variable_with_complex_equality_constraints(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + n = 10 # variable dimension (parameter) + m = 5 # number of constraints (parameter) + xo = rand(n)+im*rand(n) + A = randn(m,n) + im*randn(m,n) + b = A * xo + x = ComplexVariable(n) + p1 = minimize(real(sum(x)), A*x == b, real(x)>=0, imag(x)>=0) + handle_problem!(p1) + x1 = x.value + + xr = Variable(n) + xi = Variable(n) + p2 = minimize(sum(xr), real(A)*xr-imag(A)*xi == real(b), imag(A)*xr+real(A)*xi == imag(b), xr>=0, xi>=0) + handle_problem!(p2) + #x2 = xr.value + im*xi.value + + if test + real_diff = real(x1) - xr.value + @test real_diff ≈ zeros(10, 1) atol=atol rtol=rtol + + imag_diff = imag(x1) - xi.value + @test imag_diff ≈ zeros(10, 1) atol=atol rtol=rtol + end +end + +@add_problem sdp function sdp_Issue_198(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + ρ = HermitianSemidefinite(2) + constraints = [ρ == [ 1. 0.; 0. 1.]] + p = satisfy(constraints) + handle_problem!(p) + if test + @test p.status == :Optimal + @test p.solution.primal ≈ [0.; 1.; 0.; 0.; 1.; zeros(4)] atol=atol rtol=rtol + @test p.optval ≈ 0 atol=atol rtol=rtol + end +end + +@add_problem sdp function sdp_norm2_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + a = 2+4im + x = ComplexVariable() + objective = norm2(a-x) + c1 = real(x)>=0 + p = minimize(objective,c1) + handle_problem!(p) + if test + @test p.optval ≈ 0 atol=atol rtol=rtol + @test evaluate(objective) ≈ 0 atol=atol rtol=rtol + + real_diff = real(x.value) - real(a) + imag_diff = imag(x.value) - imag(a) + @test real_diff ≈ 0 atol=atol rtol=rtol + @test imag_diff ≈ 0 atol=atol rtol=rtol + end +end + +@add_problem sdp function sdp_sumsquares_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + a = [2+4im;4+6im] + x = ComplexVariable(2) + objective = sumsquares(a-x) + c1 = real(x)>=0 + p = minimize(objective,c1) + handle_problem!(p) + if test + @test p.optval ≈ 0 atol=atol rtol=rtol + @test evaluate(objective) ≈ zeros(1, 1) atol=atol rtol=rtol + + real_diff = real.(x.value) - real.(a) + imag_diff = imag.(x.value) - imag.(a) + @test real_diff ≈ zeros(2, 1) atol=atol rtol=rtol + @test imag_diff ≈ zeros(2, 1) atol=atol rtol=rtol + end +end + +@add_problem sdp function sdp_abs_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + a = [5-4im] + x = ComplexVariable() + objective = abs(a-x) + c1 = real(x)>=0 + p = minimize(objective,c1) + handle_problem!(p) + if test + @test p.optval ≈ 0 atol=atol rtol=rtol + @test evaluate(objective) ≈ zeros(1) atol=atol rtol=rtol + + real_diff = real(x.value) .- real(a) + imag_diff = imag(x.value) .- imag(a) + @test real_diff ≈ zeros(1) atol=atol rtol=rtol + @test imag_diff ≈ zeros(1) atol=atol rtol=rtol + end +end + +@add_problem sdp function sdp_Complex_Semidefinite_constraint(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + n = 10 + A = rand(n,n) + im*rand(n,n) + A = A + A' # now A is hermitian + x = ComplexVariable(n,n) + objective = sumsquares(A - x) + c1 = x in :SDP + p = minimize(objective, c1) + handle_problem!(p) + # test that X is approximately equal to posA: + l,v = eigen(A) + posA = v*Diagonal(max.(l,0))*v' + + + if test + real_diff = real.(x.value) - real.(posA) + imag_diff = imag.(x.value) - imag.(posA) + @test real_diff ≈ zeros(n, n) atol=atol rtol=rtol + @test imag_diff ≈ zeros(n, n) atol=atol rtol=rtol + end +end diff --git a/src/problem_depot/problems/sdp_and_exp.jl b/src/problem_depot/problems/sdp_and_exp.jl new file mode 100644 index 000000000..ef11c2a8f --- /dev/null +++ b/src/problem_depot/problems/sdp_and_exp.jl @@ -0,0 +1,12 @@ +@add_problem sdp_and_exp function sdp_and_exp_log_det_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + x = Variable(2, 2) + p = maximize(logdet(x), [x[1, 1] == 1, x[2, 2] == 1]) + if test + @test vexity(p) == ConvexVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 0 atol=10atol rtol=rtol + @test evaluate(logdet(x)) ≈ 0 atol=10atol rtol=rtol + end +end diff --git a/src/problem_depot/problems/socp.jl b/src/problem_depot/problems/socp.jl new file mode 100644 index 000000000..5269a6d6e --- /dev/null +++ b/src/problem_depot/problems/socp.jl @@ -0,0 +1,424 @@ +@add_problem socp function socp_norm_2_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + x = Variable(2, 1) + A = [1 2; 2 1; 3 4] + b = [2; 3; 4] + p = minimize(norm2(A * x + b)) + if test + @test vexity(p) == ConvexVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 0.64888 atol=atol rtol=rtol + @test evaluate(norm2(A * x + b)) ≈ 0.64888 atol=atol rtol=rtol + end + + x = Variable(2, 1) + A = [1 2; 2 1; 3 4] + b = [2; 3; 4] + lambda = 1 + p = minimize(norm2(A * x + b) + lambda * norm2(x), x >= 1) + if test + @test vexity(p) == ConvexVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 14.9049 atol=atol rtol=rtol + @test evaluate(norm2(A * x + b) + lambda * norm2(x)) ≈ 14.9049 atol=atol rtol=rtol + end + + x = Variable(2) + + p = minimize(norm2([x[1] + 2x[2] + 2; 2x[1] + x[2] + 3; 3x[1]+4x[2] + 4]) + lambda * norm2(x), x >= 1) + if test + @test vexity(p) == ConvexVexity() + end + + handle_problem!(p) + if test + @test p.optval ≈ 14.9049 atol=atol rtol=rtol + @test evaluate(norm2(A * x + b) + lambda * norm2(x)) ≈ 14.9049 atol=atol rtol=rtol + end + + x = Variable(2, 1) + A = [1 2; 2 1; 3 4] + b = [2; 3; 4] + lambda = 1 + p = minimize(norm2(A * x + b) + lambda * norm_1(x), x >= 1) + if test + @test vexity(p) == ConvexVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 15.4907 atol=atol rtol=rtol + @test evaluate(norm2(A * x + b) + lambda * norm_1(x)) ≈ 15.4907 atol=atol rtol=rtol + end +end + +@add_problem socp function socp_frobenius_norm_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + m = Variable(4, 5) + c = [m[3, 3] == 4, m >= 1] + p = minimize(norm(vec(m), 2), c) + if test + @test vexity(p) == ConvexVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ sqrt(35) atol=atol rtol=rtol + @test evaluate(norm(vec(m), 2)) ≈ sqrt(35) atol=atol rtol=rtol + end +end + +@add_problem socp function socp_quad_over_lin_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + x = Variable(3, 1) + A = [2 -3 5; -2 9 -3; 5 -8 3] + b = [-3; 9; 5] + c = [3 2 4] + d = -3 + p = minimize(quadoverlin(A*x + b, c*x + d)) + if test + @test vexity(p) == ConvexVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 17.7831 atol=atol rtol=rtol + @test (evaluate(quadoverlin(A * x + b, c * x + d)))[1] ≈ 17.7831 atol=atol rtol=rtol + end +end + +@add_problem socp function socp_sum_squares_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + x = Variable(2, 1) + A = [1 2; 2 1; 3 4] + b = [2; 3; 4] + p = minimize(sumsquares(A*x + b)) + if test + @test vexity(p) == ConvexVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 0.42105 atol=atol rtol=rtol + @test (evaluate(sumsquares(A * x + b)))[1] ≈ 0.42105 atol=atol rtol=rtol + end +end + +@add_problem socp function socp_square_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + x = Variable(2, 1) + A = [1 2; 2 1; 3 4] + b = [2; 3; 4] + p = minimize(sum(square(A*x + b))) + if test + @test vexity(p) == ConvexVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 0.42105 atol=atol rtol=rtol + @test evaluate(sum(square(A * x + b))) ≈ 0.42105 atol=atol rtol=rtol + end + + x = Variable(2, 1) + A = [1 2; 2 1; 3 4] + b = [2; 3; 4] + expr = A * x + b + p = minimize(sum(dot(^)(expr,2))) # elementwise ^ + if test + @test vexity(p) == ConvexVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 0.42105 atol=atol rtol=rtol + @test evaluate(sum(broadcast(^, expr, 2))) ≈ 0.42105 atol=atol rtol=rtol + end + + p = minimize(sum(dot(*)(expr, expr))) # elementwise * + if test + @test vexity(p) == ConvexVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 0.42105 atol=atol rtol=rtol + @test evaluate(sum((dot(*))(expr, expr))) ≈ 0.42105 atol=atol rtol=rtol + end +end + +@add_problem socp function socp_inv_pos_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + x = Variable(4) + p = minimize(sum(invpos(x)), invpos(x) < 2, x > 1, x == 2, 2 == x) + if test + @test vexity(p) == ConvexVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 2 atol=atol rtol=rtol + @test evaluate(sum(invpos(x))) ≈ 2 atol=atol rtol=rtol + end + + x = Variable(3) + p = minimize(sum(dot(/)([3,6,9], x)), x<=3) + handle_problem!(p) + if test + @test x.value ≈ fill(3.0, (3, 1)) atol=atol rtol=rtol + @test p.optval ≈ 6 atol=atol rtol=rtol + @test evaluate(sum((dot(/))([3, 6, 9], x))) ≈ 6 atol=atol rtol=rtol + end + + x = Variable() + p = minimize(sum([3,6,9]/x), x<=3) + handle_problem!(p) + if test + @test x.value ≈ 3 atol=atol rtol=rtol + @test p.optval ≈ 6 atol=atol rtol=rtol + @test evaluate(sum([3, 6, 9] / x)) ≈ 6 atol=atol rtol=rtol + end +end + +@add_problem socp function socp_geo_mean_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + x = Variable(2) + y = Variable(2) + p = minimize(geomean(x, y), x >= 1, y >= 2) + # not DCP compliant + if test + @test vexity(p) == ConcaveVexity() + end + p = maximize(geomean(x, y), 1 < x, x < 2, y < 2) + # Just gave it a vector as an objective, not okay + if test + @test_throws Exception handle_problem!(p) + end + + p = maximize(sum(geomean(x, y)), 1 < x, x < 2, y < 2) + handle_problem!(p) + if test + @test p.optval ≈ 4 atol=atol rtol=rtol + @test evaluate(sum(geomean(x, y))) ≈ 4 atol=atol rtol=rtol + end +end + +@add_problem socp function socp_sqrt_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + x = Variable() + p = maximize(sqrt(x), 1 >= x) +end + +@add_problem socp function socp_quad_form_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + x = Variable(3, 1) + A = [0.8608 0.3131 0.5458; 0.3131 0.8584 0.5836; 0.5458 0.5836 1.5422] + p = minimize(quadform(x, A), [x >= 1]) + if test + @test vexity(p) == ConvexVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 6.1464 atol=atol rtol=rtol + @test (evaluate(quadform(x, A)))[1] ≈ 6.1464 atol=atol rtol=rtol + end + + x = Variable(3, 1) + A = -1.0*[0.8608 0.3131 0.5458; 0.3131 0.8584 0.5836; 0.5458 0.5836 1.5422] + c = [3 2 4] + p = maximize(c*x , [quadform(x, A) >= -1]) + if test + @test vexity(p) == ConvexVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 3.7713 atol=atol rtol=rtol + @test (evaluate(quadform(x, A)))[1] ≈ -1 atol=atol rtol=rtol + end +end + +@add_problem socp function socp_huber_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + x = Variable(3) + p = minimize(sum(huber(x, 1)), x >= 2) + if test + @test vexity(p) == ConvexVexity() + end + handle_problem!(p) + if test + @test p.optval ≈ 9 atol=atol rtol=rtol + @test evaluate(sum(huber(x, 1))) ≈ 9 atol=atol rtol=rtol + end +end + +@add_problem socp function socp_rational_norm_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + A = [1 2 3; -1 2 3] + b = A * ones(3) + x = Variable(3) + p = minimize(norm(x, 4.5), [A * x == b]) + if test + @test vexity(p) == ConvexVexity() + end + # Solution is approximately x = [1, .93138, 1.04575] + handle_problem!(p) + if test + @test p.optval ≈ 1.2717 atol=atol rtol=rtol + @test evaluate(norm(x, 4.5)) ≈ 1.2717 atol=atol rtol=rtol + end +end + +@add_problem socp function socp_rational_norm_dual_norm(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + v = [0.463339, 0.0216084, -2.07914, 0.99581, 0.889391] + x = Variable(5) + q = 1.379; # q norm constraint that generates many inequalities + qs = q / (q - 1); # Conjugate to q + p = minimize(x' * v) + p.constraints += (norm(x, q) <= 1) + if test + @test vexity(p) == ConvexVexity() + end + handle_problem!(p) # Solution is -norm(v, q / (q - 1)) + if test + @test p.optval ≈ -2.144087 atol=atol rtol=rtol + @test sum(evaluate(x' * v)) ≈ -2.144087 atol=atol rtol=rtol + @test evaluate(norm(x, q)) ≈ 1 atol=atol rtol=rtol + @test sum(evaluate(x' * v)) ≈ -(sum(abs.(v) .^ qs) ^ (1 / qs)) atol=atol rtol=rtol + end +end + +@add_problem socp function socp_rational_norm_atom_sum(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + A = [-0.719255 -0.229089 + -1.33632 -1.37121 + 0.703447 -1.4482] + b = [-1.82041, -1.67516, -0.866884] + q = 1.5 + xvar = Variable(2) + p = minimize(.5 * sumsquares(xvar) + norm(A * xvar - b, q)) + if test + @test vexity(p) == ConvexVexity() + end + handle_problem!(p) + + if test + # Compute gradient, check it is zero(ish) + x_opt = xvar.value + margins = A * x_opt - b + qs = q / (q - 1); # Conjugate + denom = sum(abs.(margins).^q)^(1/qs) + g = x_opt + A' * (abs.(margins).^(q-1) .* sign.(margins)) / denom + @test p.optval ≈ 1.7227 atol=atol rtol=rtol + @test norm(g, 2) ^ 2 ≈ 0 atol=atol rtol=rtol + end +end + +@add_problem socp function socp_norm_consistent_with_Base_for_matrix_variables(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + A = randn(4, 4) + x = Variable(4, 4) + x.value = A + # Matrix norm + if test + @test evaluate(opnorm(x)) ≈ opnorm(A) atol=atol rtol=rtol + @test evaluate(opnorm(x, 1)) ≈ opnorm(A, 1) atol=atol rtol=rtol + @test evaluate(opnorm(x, 2)) ≈ opnorm(A, 2) atol=atol rtol=rtol + @test evaluate(opnorm(x, Inf)) ≈ opnorm(A, Inf) atol=atol rtol=rtol + end + # Vector norm + # TODO: Once the deprecation for norm on matrices is removed, remove the `vec` calls + if test + @test evaluate(norm(vec(x), 1)) ≈ norm(vec(A), 1) atol=atol rtol=rtol + @test evaluate(norm(vec(x), 2)) ≈ norm(vec(A), 2) atol=atol rtol=rtol + @test evaluate(norm(vec(x), 7)) ≈ norm(vec(A), 7) atol=atol rtol=rtol + @test evaluate(norm(vec(x), Inf)) ≈ norm(vec(A), Inf) atol=atol rtol=rtol + end +end + + +@add_problem socp function socp_fix_and_free_addition(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + x = Variable() + y = Variable() + + p = minimize(x+y, x>=0, y>=0) + handle_problem!(p) + if test + @test p.optval ≈ 0 atol=atol rtol=rtol + end + + y.value = 4 + fix!(y) + handle_problem!(p) + if test + @test p.optval ≈ 4 atol=atol rtol=rtol + end + + free!(y) + handle_problem!(p) + if test + @test p.optval ≈ 0 atol=atol rtol=rtol + end +end + +@add_problem socp function socp_fix_multiplication(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + a = [1,2,3,2,1] + x = Variable(length(a)) + gamma = Variable(Positive()) + fix!(gamma, 0.7) + + p = minimize(norm(x-a) + gamma*norm(x[1:end-1] - x[2:end])) + handle_problem!(p) + if test + o1 = p.optval + # x should be very close to a + @test o1 ≈ 0.7 * norm(a[1:end - 1] - a[2:end]) atol=atol rtol=rtol + end + # increase regularization + fix!(gamma, 1.0) + handle_problem!(p) + + if test + o2 = p.optval + # x should be very close to mean(a) + @test o2 ≈ norm(a .- mean(a)) atol=atol rtol=rtol + end + + if test + @test o1 <= o2 + end +end + + +@add_problem socp function socp_dual_minimal_norm_solutions(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test} + x = Variable(2) + A = [1 2; 2 4]; + b = [3, 6]; + p = minimize(norm(x, 1), A*x==b) + + if test + @test vexity(p) == ConvexVexity() + end + + handle_problem!(p) + if test + @test p.optval ≈ 1.5 atol=atol rtol=rtol + @test evaluate(x) ≈ [0, 1.5] atol=atol rtol=rtol + @test evaluate(norm(x, 1)) ≈ p.optval atol=atol rtol=rtol + @test dot(b, p.constraints[1].dual) ≈ p.optval atol=atol rtol=rtol + end + + x = Variable(2) + A = [1 2; 2 4]; + b = [3, 6]; + p = minimize(norm(x, 2), A*x==b) + + test && @test vexity(p) == ConvexVexity() + + handle_problem!(p) + + if test + @test p.optval ≈ 3/sqrt(5) atol=atol rtol=rtol + @test evaluate(x) ≈ [3/5, 6/5] atol=atol rtol=rtol + @test evaluate(norm(x, 2)) ≈ p.optval atol=atol rtol=rtol + @test dot(b, p.constraints[1].dual) ≈ p.optval atol=atol rtol=rtol + end + + x = Variable(2) + A = [1 2; 2 4]; + b = [3, 6]; + p = minimize(norm(x, Inf), A*x==b) + + test && @test vexity(p) == ConvexVexity() + + handle_problem!(p) + + if test + @test p.optval ≈ 1.0 atol=atol rtol=rtol + @test evaluate(x) ≈ [1, 1] atol=atol rtol=rtol + @test evaluate(norm(x, Inf)) ≈ p.optval atol=atol rtol=rtol + @test dot(b, p.constraints[1].dual) ≈ p.optval atol=atol rtol=rtol + end +end diff --git a/test/runtests.jl b/test/runtests.jl index 8f15f9461..0d9e542ad 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,45 +1,43 @@ using Convex -using Convex: DotMultiplyAtom +using Convex.ProblemDepot: run_tests using Test -using ECOS -using SCS -using GLPKMathProgInterface -using Random - -import LinearAlgebra.eigen -import LinearAlgebra.I -import LinearAlgebra.opnorm -import Random.shuffle -import Statistics.mean +using SCS, ECOS, GLPKMathProgInterface -TOL = 1e-3 -eye(n) = Matrix(1.0I, n, n) # Seed random number stream to improve test reliability +using Random Random.seed!(2) -solvers = Any[] - -push!(solvers, ECOSSolver(verbose=0)) -push!(solvers, GLPKSolverMIP()) -push!(solvers, SCSSolver(verbose=0, eps=1e-6)) - -# If Gurobi is installed, uncomment to test with it: -#using Gurobi -#push!(solvers, GurobiSolver(OutputFlag=0)) - -# If Mosek is installed, uncomment to test with it: -#using Mosek -#push!(solvers, MosekSolver(LOG=0)) +@testset "ProblemDepot" begin + @testset "Problems can run without `solve!`ing if `test==false`" begin + Convex.ProblemDepot.foreach_problem() do name, func + @testset "$name" begin + # We want to check to make sure this does not throw + func(Convex.conic_problem, Val(false), 0.0, 0.0, Float64) + @test true + end + end + end +end @testset "Convex" begin include("test_utilities.jl") - include("test_const.jl") - include("test_affine.jl") - include("test_lp.jl") - include("test_socp.jl") - include("test_sdp.jl") - include("test_exp.jl") - include("test_sdp_and_exp.jl") - include("test_mip.jl") + + @testset "SCS" begin + run_tests(; exclude=[r"mip"]) do p + solve!(p, SCSSolver(verbose=0, eps=1e-6)) + end + end + + @testset "ECOS" begin + run_tests(; exclude=[r"mip", r"sdp"]) do p + solve!(p, ECOSSolver(verbose=0)) + end + end + + @testset "GLPK MIP" begin + run_tests(; exclude=[r"socp", r"sdp", r"exp", r"dual"]) do p + solve!(p, GLPKSolverMIP()) + end + end end diff --git a/test/test_affine.jl b/test/test_affine.jl deleted file mode 100644 index 70c68793b..000000000 --- a/test/test_affine.jl +++ /dev/null @@ -1,465 +0,0 @@ -@testset "Affine Atoms: $solver" for solver in solvers - @testset "negate atom" begin - x = Variable() - p = minimize(-x, [x <= 0]) - @test vexity(p) == AffineVexity() - solve!(p, solver) - @test p.optval ≈ 0 atol=TOL - @test evaluate(-x) ≈ 0 atol=TOL - end - - @testset "kron atom" begin - x = ComplexVariable(3, 3) - y = [1.0 2.0; 3.0 4.0] - @test size(kron(x, y)) == (6, 6) - @test size(kron(y, x)) == (6, 6) - end - - @testset "multiply atom" begin - x = Variable(1) - p = minimize(2.0 * x, [x >= 2, x <= 4]) - @test vexity(p) == AffineVexity() - solve!(p, solver) - @test p.optval ≈ 4 atol=TOL - @test (evaluate(2.0x))[1] ≈ 4 atol=TOL - - x = Variable(2) - A = 1.5 * eye(2) - p = minimize([2 2] * x, [A * x >= [1.1; 1.1]]) - @test vexity(p) == AffineVexity() - solve!(p, solver) - @test p.optval ≈ 2.93333 atol=TOL - @test (evaluate([2 2] * x))[1] ≈ 2.93333 atol=TOL - @test vec(evaluate(A * x)) ≈ [1.1; 1.1] atol=TOL - - y = Variable(1) - x = Variable(3) - z = [1.0, 2.0, 3.0] * y - k = -y * [1.0, 2.0, 3.0] - c = [y <= 3.0, y >= 0.0, x >= ones(3), k <= x, x <= z] - o = 3 * y - p = Problem(:minimize, o, c) - @test vexity(p) == AffineVexity() - solve!(p, solver) - @test p.optval ≈ 3 atol=TOL - - p = Problem(:minimize, o, c...) - @test vexity(p) == AffineVexity() - solve!(p, solver) - @test p.optval ≈ 3 atol=TOL - - # Check #274 - x = ComplexVariable(2,2) - p = minimize( real( [1.0im, 0.0]' * x * [1.0im, 0.0] ), [ x == [1.0 0.0; 0.0 1.0] ]) - solve!(p, solver) - @test p.optval ≈ 1.0 - end - - @testset "dot atom" begin - x = Variable(2) - p = minimize(dot([2.0; 2.0], x), x >= [1.1; 1.1]) - @test vexity(p) == AffineVexity() - solve!(p, solver) - @test p.optval ≈ 4.4 atol=TOL - @test (evaluate(dot([2.0; 2.0], x)))[1] ≈ 4.4 atol=TOL - end - - @testset "dot atom for matrix variables" begin - x = Variable(2,2) - p = minimize(dot(fill(2.0, (2,2)), x), x >= 1.1) - @test vexity(p) == AffineVexity() - solve!(p, solver) - @test p.optval ≈ 8.8 atol=TOL - @test (evaluate(dot(fill(2.0, (2, 2)), x)))[1] ≈ 8.8 atol=TOL - end - - @testset "add atom" begin - x = Variable(1) - y = Variable(1) - p = minimize(x + y, [x >= 3, y >= 2]) - @test vexity(p) == AffineVexity() - solve!(p, solver) - @test p.optval ≈ 5 atol=TOL - @test evaluate(x + y) ≈ 5 atol=TOL - - x = Variable(1) - p = minimize(x, [eye(2) + x >= eye(2)]) - @test vexity(p) == AffineVexity() - solve!(p, solver) - @test p.optval ≈ 0 atol=TOL - @test evaluate(eye(2) + x) ≈ eye(2) atol=TOL - - y = Variable() - p = minimize(y - 5, y >= -1) - @test vexity(p) == AffineVexity() - solve!(p, solver) - @test p.optval ≈ -6 atol=TOL - @test evaluate(y - 5) ≈ -6 atol=TOL - end - - @testset "transpose atom" begin - x = Variable(2) - c = ones(2, 1) - p = minimize(x' * c, x >= 1) - @test vexity(p) == AffineVexity() - solve!(p, solver) - @test p.optval ≈ 2 atol=TOL - @test (evaluate(x' * c))[1] ≈ 2 atol=TOL - - X = Variable(2, 2) - c = ones(2, 1) - p = minimize(c' * X' * c, [X >= ones(2, 2)]) - @test vexity(p) == AffineVexity() - solve!(p, solver) - @test p.optval ≈ 4 atol=TOL - @test (evaluate(c' * X' * c))[1] ≈ 4 atol=TOL - - rows = 2 - cols = 3 - r = rand(rows, cols) - r_2 = rand(cols, rows) - x = Variable(rows, cols) - c = ones(1, cols) - d = ones(rows, 1) - p = minimize(c * x' * d + d' * x * c' + (c * x''''' * d)', - [x' >= r_2, x >= r, x''' >= r_2, x'' >= r]) - @test vexity(p) == AffineVexity() - solve!(p, solver) - s = sum(max.(r, r_2')) * 3 - @test p.optval ≈ s atol=TOL - @test (evaluate(c * x' * d + d' * x * c' + (c * ((((x')')')')' * d)'))[1] ≈ s atol=TOL - end - - @testset "index atom" begin - x = Variable(2) - p = minimize(x[1] + x[2], [x >= 1]) - @test vexity(p) == AffineVexity() - solve!(p, solver) - @test p.optval ≈ 2 atol=TOL - @test (evaluate(x[1] + x[2]))[1] ≈ 2 atol=TOL - - x = Variable(3) - I = [true true false] - p = minimize(sum(x[I]), [x >= 1]) - @test vexity(p) == AffineVexity() - solve!(p, solver) - @test p.optval ≈ 2 atol=TOL - @test (evaluate(sum(x[I])))[1] ≈ 2 atol=TOL - - rows = 6 - cols = 8 - n = 2 - X = Variable(rows, cols) - A = randn(rows, cols) - c = rand(1, n) - p = minimize(c * X[1:n, 5:5+n-1]' * c', X >= A) - @test vexity(p) == AffineVexity() - solve!(p, solver) - s = c * A[1:n, 5:5+n-1]' * c' - @test p.optval ≈ s[1] atol=TOL - @test evaluate(c * (X[1:n, 5:(5 + n) - 1])' * c') ≈ s atol=TOL - end - - @testset "sum atom" begin - x = Variable(2,2) - p = minimize(sum(x), x>=1) - @test vexity(p) == AffineVexity() - solve!(p, solver) - @test p.optval ≈ 4 atol=TOL - @test evaluate(sum(x)) ≈ 4 atol=TOL - - x = Variable(2,2) - p = minimize(sum(x) - 2*x[1,1], x>=1, x[1,1]<=2) - @test vexity(p) == AffineVexity() - solve!(p, solver) - @test p.optval ≈ 1 atol=TOL - @test (evaluate(sum(x) - 2 * x[1, 1]))[1] ≈ 1 atol=TOL - - x = Variable(10) - a = rand(10, 1) - p = maximize(sum(x[2:6]), x <= a) - @test vexity(p) == AffineVexity() - solve!(p, solver) - @test p.optval ≈ sum(a[2:6]) atol=TOL - @test evaluate(sum(x[2:6])) ≈ sum(a[2:6]) atol=TOL - end - - @testset "diag atom" begin - x = Variable(2,2) - p = minimize(sum(diag(x,1)), x >= 1) - @test vexity(p) == AffineVexity() - solve!(p, solver) - @test p.optval ≈ 1 atol=TOL - @test evaluate(sum(diag(x, 1))) ≈ 1 atol=TOL - - x = Variable(4, 4) - p = minimize(sum(diag(x)), x >= 2) - @test vexity(p) == AffineVexity() - solve!(p, solver) - @test p.optval ≈ 8 atol=TOL - @test evaluate(sum(diag(x))) ≈ 8 atol=TOL - end - - @testset "trace atom" begin - x = Variable(2,2) - p = minimize(tr(x), x >= 1) - @test vexity(p) == AffineVexity() - solve!(p, solver) - @test p.optval ≈ 2 atol=TOL - @test evaluate(tr(x)) ≈ 2 atol=TOL - end - - @testset "dot multiply atom" begin - x = Variable(3) - p = maximize(sum(dot(*)(x,[1,2,3])), x<=1) - @test vexity(p) == AffineVexity() - solve!(p, solver) - @test p.optval ≈ 6 atol=TOL - @test evaluate(sum((dot(*))(x, [1, 2, 3]))) ≈ 6 atol=TOL - - x = Variable(3, 3) - p = maximize(sum(dot(*)(x,eye(3))), x<=1) - @test vexity(p) == AffineVexity() - solve!(p, solver) - @test p.optval ≈ 3 atol=TOL - @test evaluate(sum((dot(*))(x, eye(3)))) ≈ 3 atol=TOL - - x = Variable(5, 5) - p = minimize(x[1, 1], dot(*)(3,x) >= 3) - @test vexity(p) == AffineVexity() - solve!(p, solver) - @test p.optval ≈ 1 atol=TOL - @test (evaluate(x[1, 1]))[1] ≈ 1 atol=TOL - - x = Variable(3,1) - p = minimize(sum(dot(*)(ones(3,3), x)), x>=1) - @test vexity(p) == AffineVexity() - solve!(p, solver) - @test p.optval ≈ 9 atol=TOL - @test (evaluate(x[1, 1]))[1] ≈ 1 atol=TOL - - x = Variable(1,3) - p = minimize(sum(dot(*)(ones(3,3), x)), x>=1) - @test vexity(p) == AffineVexity() - solve!(p, solver) - @test p.optval ≈ 9 atol=TOL - @test (evaluate(x[1, 1]))[1] ≈ 1 atol=TOL - - x = Variable(1, 3, Positive()) - p = maximize(sum(dot(/)(x,[1 2 3])), x<=1) - @test vexity(p) == AffineVexity() - solve!(p, solver) - @test p.optval ≈ 11 / 6 atol=TOL - @test evaluate(sum((dot(/))(x, [1 2 3]))) ≈ 11 / 6 atol=TOL - - # Broadcast fusion works - x = Variable(5, 5) - a = 2.0 .* x .* ones(Int, 5) - @test a isa DotMultiplyAtom - end - - @testset "reshape atom" begin - A = rand(2, 3) - X = Variable(3, 2) - c = rand() - p = minimize(sum(reshape(X, 2, 3) + A), X >= c) - @test vexity(p) == AffineVexity() - solve!(p, solver) - @test p.optval ≈ sum(A .+ c) atol=TOL - @test evaluate(sum(reshape(X, 2, 3) + A)) ≈ sum(A .+ c) atol=TOL - - b = rand(6) - p = minimize(sum(vec(X) + b), X >= c) - @test vexity(p) == AffineVexity() - solve!(p, solver) - @test p.optval ≈ sum(b .+ c) atol=TOL - @test evaluate(sum(vec(X) + b)) ≈ sum(b .+ c) atol=TOL - - x = Variable(4, 4) - c = ones(16, 1) - reshaped = reshape(x, 16, 1) - a = collect(1:16) - p = minimize(c' * reshaped, reshaped >= a) - @test vexity(p) == AffineVexity() - solve!(p, solver) - # TODO: why is accuracy lower here? - @test p.optval ≈ 136 atol=10TOL - @test (evaluate(c' * reshaped))[1] ≈ 136 atol=10TOL - end - - @testset "hcat atom" begin - x = Variable(4, 4) - y = Variable(4, 6) - p = maximize(sum(x) + sum([y fill(4.0, 4)]), [x y fill(2.0, (4, 2))] <= 2) - @test vexity(p) == AffineVexity() - solve!(p, solver) - @test p.optval ≈ 96 atol=TOL - @test evaluate(sum(x) + sum([y fill(4.0, 4)])) ≈ 96 atol=TOL - @test evaluate([x y fill(2.0, (4, 2))]) ≈ fill(2.0, (4, 12)) atol=TOL - end - - @testset "vcat atom" begin - x = Variable(4, 4) - y = Variable(4, 6) - - # TODO: fix dimension mismatch [y 4*eye(4); x -ones(4, 6)] - p = maximize(sum(x) + sum([y 4*eye(4); x -ones(4, 6)]), [x;y'] <= 2) - @test vexity(p) == AffineVexity() - solve!(p, solver) - # TODO: why is accuracy lower here? - @test p.optval ≈ 104 atol=10TOL - @test evaluate(sum(x) + sum([y 4 * eye(4); x -(ones(4, 6))])) ≈ 104 atol=10TOL - @test evaluate([x; y']) ≈ 2 * ones(10, 4) atol=TOL - end - - @testset "Diagonal atom" begin - x = Variable(2, 2) - @test_throws ArgumentError Diagonal(x) - - x = Variable(4) - p = minimize(sum(Diagonal(x)), x == [1, 2, 3, 4]) - @test vexity(p) == AffineVexity() - solve!(p, solver) - @test p.optval ≈ 10 atol=TOL - @test all(abs.(evaluate(Diagonal(x)) - Diagonal([1, 2, 3, 4])) .<= TOL) - - x = Variable(3) - c = [1; 2; 3] - p = minimize(c' * Diagonal(x) * c, x >= 1, sum(x) == 10) - @test vexity(p) == AffineVexity() - solve!(p, solver) - @test p.optval ≈ 21 atol=TOL - - x = Variable(3) - p = minimize(sum(x), x >= 1, Diagonal(x)[1, 2] == 1) - @test solve!(p, solver) === nothing - @test p.status != :Optimal - end - - @testset "conv atom" begin - x = Variable(3) - h = [1, -1] - p = minimize(sum(conv(h, x)) + sum(x), x >= 1, x <= 2) - @test vexity(p) == AffineVexity() - solve!(p, solver) - @test p.optval ≈ 3 atol=TOL - @test evaluate(sum(conv(h, x))) ≈ 0 atol=TOL - - x = Variable(3) - h = [1, -1] - p = minimize(sum(conv(x, h)) + sum(x), x >= 1, x <= 2) - @test vexity(p) == AffineVexity() - solve!(p, solver) - @test p.optval ≈ 3 atol=TOL - @test evaluate(sum(conv(h, x))) ≈ 0 atol=TOL - - end - - @testset "satisfy problems" begin - x = Variable() - p = satisfy(x >= 0) - add_constraints!(p, x >= 1) - add_constraints!(p, [x >= -1, x <= 4]) - solve!(p, solver) - @test p.status == :Optimal - - p = satisfy([x >= 0, x >= 1, x <= 2]) - solve!(p, solver) - @test p.status == :Optimal - - p = maximize(1, [x >= 1, x <= 2]) - solve!(p, solver) - @test p.status == :Optimal - - constr = x >= 0 - constr += x >= 1 - constr += x <= 10 - constr2 = x >= 0 - constr2 += [x >= 2, x <= 3] + constr - p = satisfy(constr) - solve!(p, solver) - @test p.status == :Optimal - end - - @testset "dual" begin - x = Variable() - p = minimize(x, x >= 0) - solve!(p, solver) - if p.solution.has_dual - println("Solution object has dual value, checking for dual correctness.") - @test p.constraints[1].dual ≈ 1 atol=TOL - end - - x = Variable() - p = maximize(x, x <= 0) - solve!(p, solver) - if p.solution.has_dual - println("Solution object has dual value, checking for dual correctness.") - @test p.constraints[1].dual ≈ 1 atol=TOL - end - - x = Variable() - p = minimize(x, x >= 0, x == 2) - solve!(p, solver) - if p.solution.has_dual - println("Solution object has dual value, checking for dual correctness.") - @test p.constraints[1].dual ≈ 0 atol=TOL - @test abs.(p.constraints[2].dual) ≈ 1 atol=TOL - end - - x = Variable(2) - A = 1.5 * eye(2) - p = minimize(dot([2.0; 2.0], x), [A * x >= [1.1; 1.1]]) - solve!(p, solver) - if p.solution.has_dual - println("Solution object has dual value, checking for dual correctness.") - dual = [4/3; 4/3] - @test all(abs.(p.constraints[1].dual - dual) .<= TOL) - end - end - - @testset "Partial transpose" begin - dims = [2,3,4] - d = prod(dims) - A = rand(ComplexF64,2,2) - B = rand(ComplexF64,3,3) - C = rand(ComplexF64,4,4) - M = kron(A,B,C) - Mt1 = kron(transpose(A),B,C) - Mt2 = kron(A,transpose(B),C) - Mt3 = kron(A,B,transpose(C)) - - Rt1 = ComplexVariable(d,d) - Rt2 = ComplexVariable(d,d) - Rt3 = ComplexVariable(d,d) - S = rand(ComplexF64,d,d) - solve!(satisfy(partialtranspose(Rt1, 1, dims) == S ),solver) - solve!(satisfy(partialtranspose(Rt2, 2, dims) == S ),solver) - solve!(satisfy(partialtranspose(Rt3, 3, dims) == S ),solver) - - - @test partialtranspose(M,1,dims) ≈ Mt1 atol = TOL - @test partialtranspose(M,2,dims) ≈ Mt2 atol = TOL - @test partialtranspose(M,3,dims) ≈ Mt3 atol = TOL - @test partialtranspose(S,1,dims) ≈ evaluate(Rt1) atol = TOL - @test partialtranspose(S,2,dims) ≈ evaluate(Rt2) atol = TOL - @test partialtranspose(S,3,dims) ≈ evaluate(Rt3) atol = TOL - - @test_throws ArgumentError partialtrace(rand(6, 6), 3, [2, 3]) - @test_throws ArgumentError partialtrace(rand(6, 6), 1, [2, 4]) - @test_throws ArgumentError partialtrace(rand(3, 4), 1, [2, 3]) - end - - @testset "permuteddims_matrix" begin - #this function is used in the partial transpose - for n in (2, 3, 4, 5) - dims = ntuple( i -> rand(2:5), n) - d = prod(dims) - v = rand(d) - p = randperm(n) - out1 = vec(permutedims(reshape(v, dims), p)) - out2 = Convex.permutedims_matrix(dims, p) * v - @test out1 ≈ out2 - end - end -end diff --git a/test/test_const.jl b/test/test_const.jl deleted file mode 100644 index 696c68203..000000000 --- a/test/test_const.jl +++ /dev/null @@ -1,107 +0,0 @@ -@testset "Constant variables: $solver" for solver in solvers - - @testset "Issue #166" begin - # Issue #166 - α = Variable(5) - fix!(α, ones(5,1)) - - # has const vexity, but not at the head - c = (rand(5,5) * α) * ones(1,5) - - β = Variable(5) - β.value = ones(5) - - problem = minimize(sum(c * β), [β >= 0]) - solve!(problem, solver) - @test problem.optval ≈ evaluate(sum(c * β)) atol=TOL - @test problem.optval ≈ 0.0 atol=TOL - @test β.value ≈ zeros(5) atol=TOL - end - - @testset "Issue #228" begin - x = Variable(2) - y = Variable(2) - fix!(x, [1 1]') - prob = minimize(y'*(x+[2 2]'), [y>=0]) - solve!(prob, solver) - @test prob.optval ≈ 0.0 atol = TOL - - prob = minimize(x'*y, [y>=0]) - solve!(prob, solver) - @test prob.optval ≈ 0.0 atol = TOL - end - - @testset "Test double `fix!`" begin - x = Variable() - y = Variable() - fix!(x, 1.0) - prob = minimize(y*x, [y >= x, x >= 0.5]) - solve!(prob, solver) - @test prob.optval ≈ 1.0 atol = TOL - - fix!(x, 2.0) - solve!(prob, solver) - @test prob.optval ≈ 4.0 atol = TOL - - free!(x) - fix!(y, 1.0) - solve!(prob, solver) - @test prob.optval ≈ 0.5 atol = TOL - end - - @testset "fix! and multiply" begin - p = Variable() - fix!(p, 1.0) - x = Variable(2,2) - prob = minimize( tr(p*x), [ x >= 1 ]) - solve!(prob, solver) - @test prob.optval ≈ 2.0 atol = TOL - @test evaluate( tr(p*x) ) ≈ 2.0 atol = TOL - end - - @testset "fix! with complex numbers" begin - x = ComplexVariable() - fix!(x, 1.0 + im*1.0) - y = Variable() - prob = minimize( real(x*y), [ y >= .5, real(x) >= .5, imag(x) >= 0]) - solve!(prob, solver) - @test prob.optval ≈ .5 atol=TOL - @test evaluate(real(x*y)) ≈ .5 atol=TOL - @test evaluate(y) ≈ 0.5 atol=TOL - - free!(x) - fix!(y) - solve!(prob, solver) - @test prob.optval ≈ 0.25 atol=TOL - @test evaluate(real(x*y)) ≈ 0.25 atol=TOL - @test real(evaluate(x)) ≈ 0.5 atol=TOL - @test evaluate(y) ≈ 0.5 atol=TOL - - @test_throws DimensionMismatch fix!(x, rand(2)) - @test_throws DimensionMismatch fix!(x, rand(2,2)) - end - - @testset "fix! with vectors" begin - x = ComplexVariable(5) - fix!(x, ones(5) + im*ones(5)) - y = Variable() - prob = minimize( real(y*sum(x)), [ y >= .5, real(x) >= .5, imag(x) >= 0]) - solve!(prob, solver) - @test prob.optval ≈ 2.5 atol=TOL - @test evaluate(real(y*sum(x))) ≈ 2.5 atol=TOL - @test evaluate(y) ≈ 0.5 atol=TOL - - free!(x) - fix!(y) - solve!(prob, solver) - @test prob.optval ≈ 1.25 atol=TOL - @test evaluate(real(y*sum(x))) ≈ 1.25 atol=TOL - @test real(evaluate(x)) ≈ 0.5*ones(5) atol=TOL - @test evaluate(y) ≈ 0.5 atol=TOL - - @test_throws DimensionMismatch fix!(x, rand(5,5)) - @test_throws DimensionMismatch fix!(x, rand(4)) - - end - -end \ No newline at end of file diff --git a/test/test_exp.jl b/test/test_exp.jl deleted file mode 100644 index 4354f1a60..000000000 --- a/test/test_exp.jl +++ /dev/null @@ -1,97 +0,0 @@ - -@testset "Exp Atoms: $solver" for solver in solvers - if can_solve_exp(solver) - @testset "exp atom" begin - y = Variable() - p = minimize(exp(y), y>=0) - @test vexity(p) == ConvexVexity() - solve!(p, solver) - @test p.optval ≈ 1 atol=TOL - @test evaluate(exp(y)) ≈ 1 atol=TOL - - y = Variable() - p = minimize(exp(y), y>=1) - @test vexity(p) == ConvexVexity() - solve!(p, solver) - @test p.optval ≈ exp(1) atol=TOL - @test evaluate(exp(y)) ≈ exp(1) atol=TOL - - y = Variable(5) - p = minimize(sum(exp(y)), y>=0) - @test vexity(p) == ConvexVexity() - solve!(p, solver) - @test p.optval ≈ 5 atol=TOL - @test evaluate(sum(exp(y))) ≈ 5 atol=TOL - - y = Variable(5) - p = minimize(sum(exp(y)), y>=0) - @test vexity(p) == ConvexVexity() - solve!(p, solver) - @test p.optval ≈ 5 atol=TOL - end - - @testset "log atom" begin - y = Variable() - p = maximize(log(y), y<=1) - @test vexity(p) == ConvexVexity() - solve!(p, solver) - @test p.optval ≈ 0 atol=TOL - - y = Variable() - p = maximize(log(y), y<=2) - @test vexity(p) == ConvexVexity() - solve!(p, solver) - @test p.optval ≈ log(2) atol=TOL - - y = Variable() - p = maximize(log(y), [y<=2, exp(y)<=10]) - @test vexity(p) == ConvexVexity() - solve!(p, solver) - @test p.optval ≈ log(2) atol=TOL - end - - @testset "log sum exp atom" begin - y = Variable(5) - p = minimize(logsumexp(y), y>=1) - @test vexity(p) == ConvexVexity() - solve!(p, solver) - @test p.optval ≈ log(exp(1) * 5) atol=TOL - end - - @testset "logistic loss atom" begin - y = Variable(5) - p = minimize(logisticloss(y), y>=1) - @test vexity(p) == ConvexVexity() - solve!(p, solver) - @test p.optval ≈ log(exp(1) + 1) * 5 atol=TOL - end - - @testset "entropy atom" begin - y = Variable(5, Positive()) - p = maximize(entropy(y), sum(y)<=1) - @test vexity(p) == ConvexVexity() - solve!(p, solver) - @test p.optval ≈ -(log(1 / 5)) atol=TOL - end - - @testset "relative entropy atom" begin - x = Variable(1) - y = Variable(1) - # x log (x/y) - p = minimize(relative_entropy(x,y), y==1, x >= 2) - @test vexity(p) == ConvexVexity() - solve!(p, solver) - @test p.optval ≈ 2 * log(2) atol=TOL - end - - @testset "log perspective atom" begin - x = Variable(1) - y = Variable(1) - # y log (x/y) - p = maximize(log_perspective(x,y), y==5, x <= 10) - @test vexity(p) == ConvexVexity() - solve!(p, solver) - @test p.optval ≈ 5 * log(2) atol=TOL - end - end -end diff --git a/test/test_lp.jl b/test/test_lp.jl deleted file mode 100644 index c2a511fca..000000000 --- a/test/test_lp.jl +++ /dev/null @@ -1,190 +0,0 @@ -@testset "LP Atoms: $solver" for solver in solvers - @testset "abs atom" begin - x = Variable() - p = minimize(abs(x), x<=-1) - @test vexity(p) == ConvexVexity() - solve!(p, solver) - @test p.optval ≈ 1 atol=TOL - @test evaluate(abs(x)) ≈ 1 atol=TOL - if p.solution.has_dual - println("Solution object has dual value, checking for dual correctness.") - @test p.constraints[1].dual ≈ 1 atol=TOL - end - - x = Variable(2,2) - p = minimize(sum(abs(x)), x[2,2]>=1, x[1,1]>=1, x>=0) - @test vexity(p) == ConvexVexity() - solve!(p, solver) - @test p.optval ≈ 2 atol=TOL - @test evaluate(sum(abs(x))) ≈ 2 atol=TOL - if p.solution.has_dual - println("Solution object has dual value, checking for dual correctness.") - @test p.constraints[1].dual ≈ 1 atol=TOL - @test p.constraints[2].dual ≈ 1 atol=TOL - @test p.constraints[3].dual[1,1] ≈ 0 atol=TOL - @test p.constraints[3].dual[2,2] ≈ 0 atol=TOL - @test p.constraints[3].dual[1,2] ≈ p.constraints[3].dual[2,1] atol=TOL - end - end - - @testset "maximum atom" begin - x = Variable(10) - a = shuffle(collect(0.1:0.1:1.0)) - p = minimize(maximum(x), x >= a) - @test vexity(p) == ConvexVexity() - solve!(p, solver) - @test p.optval ≈ maximum(a) atol=TOL - @test evaluate(maximum(x)) ≈ maximum(a) atol=TOL - end - - @testset "minimum atom" begin - x = Variable(1) - a = reshape(shuffle(collect(0.01:0.01:1.0)), (10, 10)) - p = maximize(minimum(x), x <= a) - @test vexity(p) == ConvexVexity() - solve!(p, solver) - @test p.optval ≈ minimum(a) atol=TOL - @test evaluate(minimum(x)) ≈ minimum(a) atol=TOL - - x = Variable(4, 4) - y = Variable(4, 6) - z = Variable(1) - c = ones(4, 1) - d = fill(2.0, (6, 1)) - constraints = [[x y] <= 2, z <= 0, z <= x, 2z >= -1] - objective = sum(x + z) + minimum(y) + c' * y * d - p = maximize(objective, constraints) - @test vexity(p) == ConvexVexity() - solve!(p, solver) - @test p.optval ≈ 130 atol=TOL - @test (evaluate(objective))[1] ≈ 130 atol=TOL - end - - @testset "max atom" begin - x = Variable(10, 10) - y = Variable(10, 10) - a = reshape(shuffle(collect(0.01:0.01:1.0)), (10, 10)) - b = reshape(shuffle(collect(0.01:0.01:1.0)), (10, 10)) - p = minimize(maximum(max(x, y)), [x >= a, y >= b]) - @test vexity(p) == ConvexVexity() - solve!(p, solver) - max_a = maximum(a) - max_b = maximum(b) - @test p.optval ≈ max(max_a, max_b) atol=10TOL - @test evaluate(maximum(max(x, y))) ≈ max(max_a, max_b) atol=10TOL - end - - @testset "min atom" begin - x = Variable(10, 10) - y = Variable(10, 10) - a = reshape(shuffle(collect(0.01:0.01:1.0)), (10, 10)) - b = reshape(shuffle(collect(0.01:0.01:1.0)), (10, 10)) - p = maximize(minimum(min(x, y)), [x <= a, y <= b]) - @test vexity(p) == ConvexVexity() - solve!(p, solver) - min_a = minimum(a) - min_b = minimum(b) - @test p.optval ≈ min(min_a, min_b) atol=10TOL - @test evaluate(minimum(min(x, y))) ≈ min(min_a, min_b) atol=10TOL - end - - @testset "pos atom" begin - x = Variable(3) - a = [-2; 1; 2] - p = minimize(sum(pos(x)), [x >= a, x <= 2]) - @test vexity(p) == ConvexVexity() - solve!(p, solver) - @test p.optval ≈ 3 atol=TOL - @test evaluate(sum(pos(x))) ≈ 3 atol=TOL - end - - @testset "neg atom" begin - x = Variable(3) - p = minimize(1, [x >= -2, x <= -2, neg(x) <= 3]) - @test vexity(p) == ConvexVexity() - solve!(p, solver) - @test p.optval ≈ 1 atol=TOL - @test evaluate(sum(neg(x))) ≈ 6 atol=TOL - end - - @testset "sumlargest atom" begin - x = Variable(2) - p = minimize(sumlargest(x, 2), x >= [1; 1]) - @test vexity(p) == ConvexVexity() - solve!(p, solver) - @test p.optval ≈ 2 atol=TOL - @test evaluate(sumlargest(x, 2)) ≈ 2 atol=TOL - - x = Variable(4, 4) - p = minimize(sumlargest(x, 3), x >= eye(4), x[1, 1] >= 1.5, x[2, 3] >= 2.1) - @test vexity(p) == ConvexVexity() - solve!(p, solver) - @test p.optval ≈ 4.6 atol=TOL - @test evaluate(sumlargest(x, 2)) ≈ 3.6 atol=TOL - end - - @testset "sumsmallest atom" begin - x = Variable(4, 4) - p = minimize(sumlargest(x, 2), sumsmallest(x, 4) >= 1) - @test vexity(p) == ConvexVexity() - solve!(p, solver) - @test p.optval ≈ 0.5 atol=TOL - @test evaluate(sumsmallest(x, 4)) ≈ 1 atol=TOL - - x = Variable(3, 2) - p = maximize(sumsmallest(x, 3), x >= 2, x <= 5, sumlargest(x, 3) <= 12) - @test vexity(p) == ConvexVexity() - solve!(p, solver) - @test p.optval ≈ 12 atol=TOL - @test evaluate(sumsmallest(x, 3)) ≈ 12 atol=TOL - end - - @testset "dotsort atom" begin - x = Variable(4, 1) - p = minimize(dotsort(x, [1, 2, 3, 4]), sum(x) >= 7, x >= 0, x <= 2, x[4] <= 1) - @test vexity(p) == ConvexVexity() - solve!(p, solver) - @test p.optval ≈ 19 atol=TOL - @test vec(x.value) ≈ [2; 2; 2; 1] atol=TOL - @test evaluate(dotsort(x, [1, 2, 3, 4])) ≈ 19 atol=TOL - - x = Variable(2, 2) - p = minimize(dotsort(x, [1 2; 3 4]), sum(x) >= 7, x >= 0, x <= 2, x[2, 2] <= 1) - @test vexity(p) == ConvexVexity() - solve!(p, solver) - @test p.optval ≈ 19 atol=TOL - @test evaluate(dotsort(x, [1, 2, 3, 4])) ≈ 19 atol=TOL - end - - @testset "hinge loss atom" begin - # TODO: @davidlizeng. We should finish this someday. - end - - @testset "norm inf atom" begin - x = Variable(3) - p = minimize(norm_inf(x), [-2 <= x, x <= 1]) - @test vexity(p) == ConvexVexity() - solve!(p, solver) - @test p.optval ≈ 0 atol=TOL - @test evaluate(norm_inf(x)) ≈ 0 atol=TOL - if p.solution.has_dual - println("Solution object has dual value, checking for dual correctness.") - @test norm(p.constraints[1].dual) ≈ 0 atol=TOL - @test norm(p.constraints[2].dual) ≈ 0 atol=TOL - end - end - - @testset "norm 1 atom" begin - x = Variable(3) - p = minimize(norm_1(x), [-2 <= x, x <= 1]) - @test vexity(p) == ConvexVexity() - solve!(p, solver) - @test p.optval ≈ 0 atol=TOL - @test evaluate(norm_1(x)) ≈ 0 atol=TOL - if p.solution.has_dual - println("Solution object has dual value, checking for dual correctness.") - @test norm(p.constraints[1].dual) ≈ 0 atol=TOL - @test norm(p.constraints[2].dual) ≈ 0 atol=TOL - end - end -end diff --git a/test/test_mip.jl b/test/test_mip.jl deleted file mode 100644 index 8d99cd41e..000000000 --- a/test/test_mip.jl +++ /dev/null @@ -1,72 +0,0 @@ -@testset "Mixed Integer Programs: $solver" for solver in solvers - if can_solve_mip(solver) - mip_solver = !can_solve_mip(solver) ? GLPKSolverMIP() : solver - - @testset "lp fallback interface" begin - x = Variable() - p = minimize(x, x>=4.3) - @test vexity(p) == AffineVexity() - solve!(p, solver) - @test p.optval ≈ 4.3 atol=TOL - - x = Variable(2) - p = minimize(norm(x,1), x[1]>=4.3) - @test vexity(p) == ConvexVexity() - solve!(p, solver) - @test p.optval ≈ 4.3 atol=TOL - end - - @testset "integer variables" begin - x = Variable(:Int) - p = minimize(x, x>=4.3) - @test vexity(p) == AffineVexity() - solve!(p, mip_solver) - @test p.optval ≈ 5 atol=TOL - - x = Variable(2, :Int) - p = minimize(sum(x), x>=4.3) - @test vexity(p) == AffineVexity() - solve!(p, mip_solver) - @test p.optval ≈ 10 atol=TOL - - x = Variable(:Int) - y = Variable() - p = minimize(sum(x + y), x>=4.3, y>=7) - @test vexity(p) == AffineVexity() - solve!(p, mip_solver) - @test p.optval ≈ 12 atol=TOL - - x = Variable(2, :Int) - p = minimize(norm(x, 1), x[1]>=4.3) - @test vexity(p) == ConvexVexity() - solve!(p, mip_solver) - @test p.optval ≈ 5 atol=TOL - - x = Variable(2, :Int) - p = minimize(sum(x), x[1]>=4.3, x>=0) - @test vexity(p) == AffineVexity() - solve!(p, mip_solver) - @test p.optval ≈ 5 atol=TOL - - x = Variable(2, :Int) - p = minimize(sum(x), x>=.5) - @test vexity(p) == AffineVexity() - solve!(p, mip_solver) - @test p.optval ≈ 2 atol=TOL - end - - @testset "binary variables" begin - x = Variable(2, :Bin) - p = minimize(sum(x), x>=.5) - @test vexity(p) == AffineVexity() - solve!(p, mip_solver) - @test p.optval ≈ 2 atol=TOL - - x = Variable(2, :Bin) - p = minimize(sum(x), x[1]>=.5, x>=0) - @test vexity(p) == AffineVexity() - solve!(p, mip_solver) - @test p.optval ≈ 1 atol=TOL - end - end -end diff --git a/test/test_sdp.jl b/test/test_sdp.jl deleted file mode 100644 index 9f641a026..000000000 --- a/test/test_sdp.jl +++ /dev/null @@ -1,321 +0,0 @@ - -# TODO: uncomment vexity checks once SDP on vars/constraints changes vexity of problem -@testset "SDP Atoms: $solver" for solver in solvers - if can_solve_sdp(solver) - @testset "sdp variables" begin - y = Variable((2,2), :Semidefinite) - p = minimize(y[1,1]) - # @fact vexity(p) --> ConvexVexity() - solve!(p, solver) - @test p.optval ≈ 0 atol=TOL - - y = Variable((3,3), :Semidefinite) - p = minimize(y[1,1], y[2,2]==1) - # @fact vexity(p) --> ConvexVexity() - solve!(p, solver) - @test p.optval ≈ 0 atol=TOL - - # Solution is obtained as y[2,2] -> infinity - # This test fails on Mosek. See - # https://github.com/JuliaOpt/Mosek.jl/issues/29 - # y = Variable((2, 2), :Semidefinite) - # p = minimize(y[1, 1], y[1, 2] == 1) - # # @fact vexity(p) --> ConvexVexity() - # solve!(p, solver) - # @fact p.optval --> roughly(0, TOL) - - y = Semidefinite(3) - p = minimize(sum(diag(y)), y[1, 1] == 1) - # @fact vexity(p) --> ConvexVexity() - solve!(p, solver) - @test p.optval ≈ 1 atol=TOL - - y = Variable((3, 3), :Semidefinite) - p = minimize(tr(y), y[2,1]<=4, y[2,2]>=3) - # @fact vexity(p) --> ConvexVexity() - solve!(p, solver) - @test p.optval ≈ 3 atol=TOL - - x = Variable(Positive()) - y = Semidefinite(3) - p = minimize(y[1, 2], y[2, 1] == 1) - # @fact vexity(p) --> ConvexVexity() - solve!(p, solver) - @test p.optval ≈ 1 atol=TOL - end - - @testset "sdp constraints" begin - # This test fails on Mosek - x = Variable(Positive()) - y = Variable((3, 3)) - p = minimize(x + y[1, 1], isposdef(y), x >= 1, y[2, 1] == 1) - # @fact vexity(p) --> ConvexVexity() - solve!(p, solver) - @test p.optval ≈ 1 atol=TOL - end - - @testset "nuclear norm atom" begin - y = Semidefinite(3) - p = minimize(nuclearnorm(y), y[2,1]<=4, y[2,2]>=3, y[3,3]<=2) - @test vexity(p) == ConvexVexity() - solve!(p, solver) - @test p.optval ≈ 3 atol=TOL - @test evaluate(nuclearnorm(y)) ≈ 3 atol=TOL - end - - @testset "operator norm atom" begin - y = Variable((3,3)) - p = minimize(opnorm(y), y[2,1]<=4, y[2,2]>=3, sum(y)>=12) - @test vexity(p) == ConvexVexity() - solve!(p, solver) - @test p.optval ≈ 4 atol=TOL - @test evaluate(opnorm(y)) ≈ 4 atol=TOL - end - - @testset "sigma max atom" begin - y = Variable((3,3)) - p = minimize(sigmamax(y), y[2,1]<=4, y[2,2]>=3, sum(y)>=12) - @test vexity(p) == ConvexVexity() - solve!(p, solver) - @test p.optval ≈ 4 atol=TOL - @test evaluate(sigmamax(y)) ≈ 4 atol=TOL - end - - @testset "lambda max atom" begin - y = Semidefinite(3) - p = minimize(lambdamax(y), y[1,1]>=4) - @test vexity(p) == ConvexVexity() - solve!(p, solver) - @test p.optval ≈ 4 atol=TOL - @test evaluate(lambdamax(y)) ≈ 4 atol=TOL - end - - @testset "lambda min atom" begin - y = Semidefinite(3) - p = maximize(lambdamin(y), tr(y)<=6) - @test vexity(p) == ConvexVexity() - solve!(p, solver) - @test p.optval ≈ 2 atol=TOL - @test evaluate(lambdamin(y)) ≈ 2 atol=TOL - end - - @testset "matrix frac atom" begin - x = [1, 2, 3] - P = Variable(3, 3) - p = minimize(matrixfrac(x, P), P <= 2*eye(3), P >= 0.5 * eye(3)) - @test vexity(p) == ConvexVexity() - solve!(p, solver) - @test p.optval ≈ 7 atol=TOL - @test (evaluate(matrixfrac(x, P)))[1] ≈ 7 atol=TOL - end - - @testset "matrix frac atom both arguments variable" begin - x = Variable(3) - P = Variable(3, 3) - p = minimize(matrixfrac(x, P), lambdamax(P) <= 2, x[1] >= 1) - @test vexity(p) == ConvexVexity() - solve!(p, solver) - @test p.optval ≈ 0.5 atol=TOL - @test (evaluate(matrixfrac(x, P)))[1] ≈ 0.5 atol=TOL - end - - @testset "sum largest eigs" begin - x = Semidefinite(3) - p = minimize(sumlargesteigs(x, 2), x >= 1) - solve!(p, solver) - @test p.optval ≈ 3 atol=TOL - @test evaluate(x) ≈ ones(3, 3) atol=TOL - - x = Semidefinite(3) - p = minimize(sumlargesteigs(x, 2), [x[i,:] >= i for i=1:3]...) - solve!(p, solver) - @test p.optval ≈ 8.4853 atol=TOL - - x1 = Semidefinite(3) - p1 = minimize(lambdamax(x1), x1[1,1]>=4) - solve!(p1, solver) - - x2 = Semidefinite(3) - p2 = minimize(sumlargesteigs(x2, 1), x2[1,1]>=4) - solve!(p2, solver) - - @test p1.optval ≈ p2.optval atol=TOL - - x1 = Semidefinite(3) - p1 = minimize(lambdamax(x1), [x1[i,:] >= i for i=1:3]...) - solve!(p1, solver) - - x2 = Semidefinite(3) - p2 = minimize(sumlargesteigs(x2, 1), [x2[i,:] >= i for i=1:3]...) - solve!(p2, solver) - - @test p1.optval ≈ p2.optval atol=TOL - - println(p1.optval) - end - - @testset "kron atom" begin - id = eye(4) - X = Semidefinite(4) - W = kron(id, X) - p = maximize(tr(W), tr(X) ≤ 1) - @test vexity(p) == AffineVexity() - solve!(p, solver) - @test p.optval ≈ 4 atol=TOL - end - - @testset "Partial trace" begin - A = Semidefinite(2) - B = [1 0; 0 0] - ρ = kron(B, A) - constraints = [partialtrace(ρ, 1, [2; 2]) == [0.09942819 0.29923607; 0.29923607 0.90057181], ρ in :SDP] - p = satisfy(constraints) - solve!(p, solver) - @test evaluate(ρ) ≈ [0.09942819 0.29923607 0 0; 0.299237 0.900572 0 0; 0 0 0 0; 0 0 0 0] atol=TOL - @test evaluate(partialtrace(ρ, 1, [2; 2])) ≈ [0.09942819 0.29923607; 0.29923607 0.90057181] atol=TOL - - function rand_normalized(n) - A = 5*randn(n, n) + im*5*randn(n, n) - A / tr(A) - end - - As = [ rand_normalized(3) for _ = 1:5] - Bs = [ rand_normalized(2) for _ = 1:5] - p = rand(5) - - AB = sum(i -> p[i]*kron(As[i],Bs[i]), 1:5) - @test partialtrace(AB, 2, [3, 2]) ≈ sum( p .* As ) - @test partialtrace(AB, 1, [3, 2]) ≈ sum( p .* Bs ) - - A, B, C = rand(5,5), rand(4,4), rand(3,3) - ABC = kron(kron(A, B), C) - @test kron(A,B)*tr(C) ≈ partialtrace(ABC, 3, [5, 4, 3]) - - # Test 281 - A = rand(6,6) - expr = partialtrace(Constant(A), 1, [2, 3]) - @test size(expr) == size(evaluate(expr)) - - @test_throws ArgumentError partialtrace(rand(6, 6), 3, [2, 3]) - @test_throws ArgumentError partialtrace(rand(6, 6), 1, [2, 4]) - @test_throws ArgumentError partialtrace(rand(3, 4), 1, [2, 3]) - end - - @testset "Optimization with complex variables" begin - @testset "Real Variables with complex equality constraints" begin - n = 10 # variable dimension (parameter) - m = 5 # number of constraints (parameter) - xo = rand(n) - A = randn(m,n) + im*randn(m,n) - b = A * xo - x = Variable(n) - p1 = minimize(sum(x), A*x == b, x>=0) - solve!(p1, solver) - x1 = x.value - - p2 = minimize(sum(x), real(A)*x == real(b), imag(A)*x==imag(b), x>=0) - solve!(p2, solver) - x2 = x.value - @test x1 == x2 - end - - @testset "Complex Variable with complex equality constraints" begin - n = 10 # variable dimension (parameter) - m = 5 # number of constraints (parameter) - xo = rand(n)+im*rand(n) - A = randn(m,n) + im*randn(m,n) - b = A * xo - x = ComplexVariable(n) - p1 = minimize(real(sum(x)), A*x == b, real(x)>=0, imag(x)>=0) - solve!(p1, solver) - x1 = x.value - - xr = Variable(n) - xi = Variable(n) - p2 = minimize(sum(xr), real(A)*xr-imag(A)*xi == real(b), imag(A)*xr+real(A)*xi == imag(b), xr>=0, xi>=0) - solve!(p2, solver) - #x2 = xr.value + im*xi.value - real_diff = real(x1) - xr.value - - @test real_diff ≈ zeros(10, 1) atol=TOL - imag_diff = imag(x1) - xi.value - @test imag_diff ≈ zeros(10, 1) atol=TOL - #@fact x1==x2 --> true - end - - @testset "Issue #198" begin - ρ = HermitianSemidefinite(2) - constraints = [ρ == [ 1. 0.; 0. 1.]] - p = satisfy(constraints) - solve!(p, SCSSolver()) - @test p.status == :Optimal - @test p.solution.primal ≈ [0.; 1.; 0.; 0.; 1.; zeros(4)] atol=TOL - @test p.optval ≈ 0 atol=TOL - end - - @testset "norm2 atom" begin - a = 2+4im - x = ComplexVariable() - objective = norm2(a-x) - c1 = real(x)>=0 - p = minimize(objective,c1) - solve!(p, solver) - @test p.optval ≈ 0 atol=TOL - @test evaluate(objective) ≈ 0 atol=TOL - real_diff = real(x.value) - real(a) - imag_diff = imag(x.value) - imag(a) - @test real_diff ≈ 0 atol=TOL - @test imag_diff ≈ 0 atol=TOL - end - - @testset "sumsquares atom" begin - a = [2+4im;4+6im] - x = ComplexVariable(2) - objective = sumsquares(a-x) - c1 = real(x)>=0 - p = minimize(objective,c1) - solve!(p, solver) - @test p.optval ≈ 0 atol=TOL - @test evaluate(objective) ≈ zeros(1, 1) atol=TOL - real_diff = real.(x.value) - real.(a) - imag_diff = imag.(x.value) - imag.(a) - @test real_diff ≈ zeros(2, 1) atol=TOL - @test imag_diff ≈ zeros(2, 1) atol=TOL - end - - @testset "abs atom" begin - a = [5-4im] - x = ComplexVariable() - objective = abs(a-x) - c1 = real(x)>=0 - p = minimize(objective,c1) - solve!(p, solver) - @test p.optval ≈ 0 atol=TOL - @test evaluate(objective) ≈ zeros(1) atol=TOL - real_diff = real(x.value) .- real(a) - imag_diff = imag(x.value) .- imag(a) - @test real_diff ≈ zeros(1) atol=TOL - @test imag_diff ≈ zeros(1) atol=TOL - end - - @testset "Complex Semidefinite constraint" begin - n = 10 - A = rand(n,n) + im*rand(n,n) - A = A + A' # now A is hermitian - x = ComplexVariable(n,n) - objective = sumsquares(A - x) - c1 = x in :SDP - p = minimize(objective, c1) - solve!(p, solver) - # test that X is approximately equal to posA: - l,v = eigen(A) - posA = v*Diagonal(max.(l,0))*v' - - real_diff = real.(x.value) - real.(posA) - imag_diff = imag.(x.value) - imag.(posA) - @test real_diff ≈ zeros(n, n) atol=TOL - @test imag_diff ≈ zeros(n, n) atol=TOL - end - end - end -end diff --git a/test/test_sdp_and_exp.jl b/test/test_sdp_and_exp.jl deleted file mode 100644 index b7178ebc5..000000000 --- a/test/test_sdp_and_exp.jl +++ /dev/null @@ -1,13 +0,0 @@ -@testset "SDP and Exp Atoms: $solver" for solver in solvers - if can_solve_sdp(solver) && can_solve_exp(solver) - tol = 1e-2 - @testset "log det atom" begin - x = Variable(2, 2) - p = maximize(logdet(x), [x[1, 1] == 1, x[2, 2] == 1]) - @test vexity(p) == ConvexVexity() - solve!(p, solver) - @test p.optval ≈ 0 atol=tol - @test evaluate(logdet(x)) ≈ 0 atol=tol - end - end -end diff --git a/test/test_socp.jl b/test/test_socp.jl deleted file mode 100644 index 3ce701441..000000000 --- a/test/test_socp.jl +++ /dev/null @@ -1,344 +0,0 @@ -@testset "SOCP Atoms: $solver" for solver in solvers - if can_solve_socp(solver) - @testset "norm 2 atom" begin - x = Variable(2, 1) - A = [1 2; 2 1; 3 4] - b = [2; 3; 4] - p = minimize(norm2(A * x + b)) - @test vexity(p) == ConvexVexity() - solve!(p, solver) - @test p.optval ≈ 0.64888 atol=TOL - @test evaluate(norm2(A * x + b)) ≈ 0.64888 atol=TOL - - x = Variable(2, 1) - A = [1 2; 2 1; 3 4] - b = [2; 3; 4] - lambda = 1 - p = minimize(norm2(A * x + b) + lambda * norm2(x), x >= 1) - @test vexity(p) == ConvexVexity() - solve!(p, solver) - @test p.optval ≈ 14.9049 atol=TOL - @test evaluate(norm2(A * x + b) + lambda * norm2(x)) ≈ 14.9049 atol=TOL - if p.solution.has_dual - println("Solution object has dual value, checking for dual correctness.") - dual = [4.4134, 5.1546] - @test p.constraints[1].dual ≈ dual atol=TOL - end - - x = Variable(2) - p = minimize(norm2([x[1] + 2x[2] + 2; 2x[1] + x[2] + 3; 3x[1]+4x[2] + 4]) + lambda * norm2(x), x >= 1) - @test vexity(p) == ConvexVexity() - solve!(p, solver) - @test p.optval ≈ 14.9049 atol=TOL - @test evaluate(norm2(A * x + b) + lambda * norm2(x)) ≈ 14.9049 atol=TOL - if p.solution.has_dual - println("Solution object has dual value, checking for dual correctness.") - dual = [4.4134, 5.1546] - @test p.constraints[1].dual ≈ dual atol=TOL - end - - x = Variable(2, 1) - A = [1 2; 2 1; 3 4] - b = [2; 3; 4] - lambda = 1 - p = minimize(norm2(A * x + b) + lambda * norm_1(x), x >= 1) - @test vexity(p) == ConvexVexity() - solve!(p, solver) - @test p.optval ≈ 15.4907 atol=TOL - @test evaluate(norm2(A * x + b) + lambda * norm_1(x)) ≈ 15.4907 atol=TOL - if p.solution.has_dual - println("Solution object has dual value, checking for dual correctness.") - dual = [4.7062, 5.4475] - @test p.constraints[1].dual ≈ dual atol=TOL - end - end - - @testset "frobenius norm atom" begin - m = Variable(4, 5) - c = [m[3, 3] == 4, m >= 1] - p = minimize(norm(vec(m), 2), c) - @test vexity(p) == ConvexVexity() - solve!(p, solver) - @test p.optval ≈ sqrt(35) atol=TOL - @test evaluate(norm(vec(m), 2)) ≈ sqrt(35) atol=TOL - if p.solution.has_dual - println("Solution object has dual value, checking for dual correctness.") - @test p.constraints[1].dual ≈ 0.6761 atol=TOL - dual = 0.1690 .* ones(4, 5); dual[3, 3] = 0 - @test p.constraints[2].dual ≈ dual atol=TOL - end - end - - @testset "quad over lin atom" begin - x = Variable(3, 1) - A = [2 -3 5; -2 9 -3; 5 -8 3] - b = [-3; 9; 5] - c = [3 2 4] - d = -3 - p = minimize(quadoverlin(A*x + b, c*x + d)) - @test vexity(p) == ConvexVexity() - solve!(p, solver) - @test p.optval ≈ 17.7831 atol=TOL - @test (evaluate(quadoverlin(A * x + b, c * x + d)))[1] ≈ 17.7831 atol=TOL - end - - @testset "sum squares atom" begin - x = Variable(2, 1) - A = [1 2; 2 1; 3 4] - b = [2; 3; 4] - p = minimize(sumsquares(A*x + b)) - @test vexity(p) == ConvexVexity() - solve!(p, solver) - @test p.optval ≈ 0.42105 atol=TOL - @test (evaluate(sumsquares(A * x + b)))[1] ≈ 0.42105 atol=TOL - end - - @testset "square atom" begin - x = Variable(2, 1) - A = [1 2; 2 1; 3 4] - b = [2; 3; 4] - p = minimize(sum(square(A*x + b))) - @test vexity(p) == ConvexVexity() - solve!(p, solver) - @test p.optval ≈ 0.42105 atol=TOL - @test evaluate(sum(square(A * x + b))) ≈ 0.42105 atol=TOL - - x = Variable(2, 1) - A = [1 2; 2 1; 3 4] - b = [2; 3; 4] - expr = A * x + b - p = minimize(sum(dot(^)(expr,2))) # elementwise ^ - @test vexity(p) == ConvexVexity() - solve!(p, solver) - @test p.optval ≈ 0.42105 atol=TOL - @test evaluate(sum(broadcast(^, expr, 2))) ≈ 0.42105 atol=TOL - - p = minimize(sum(dot(*)(expr, expr))) # elementwise * - @test vexity(p) == ConvexVexity() - solve!(p, solver) - @test p.optval ≈ 0.42105 atol=TOL - @test evaluate(sum((dot(*))(expr, expr))) ≈ 0.42105 atol=TOL - end - - @testset "inv pos atom" begin - x = Variable(4) - p = minimize(sum(invpos(x)), invpos(x) < 2, x > 1, x == 2, 2 == x) - @test vexity(p) == ConvexVexity() - solve!(p, solver) - @test p.optval ≈ 2 atol=TOL - @test evaluate(sum(invpos(x))) ≈ 2 atol=TOL - - x = Variable(3) - p = minimize(sum(dot(/)([3,6,9], x)), x<=3) - solve!(p, solver) - @test x.value ≈ fill(3.0, (3, 1)) atol=TOL - @test p.optval ≈ 6 atol=TOL - @test evaluate(sum((dot(/))([3, 6, 9], x))) ≈ 6 atol=TOL - - x = Variable() - p = minimize(sum([3,6,9]/x), x<=3) - solve!(p, solver) - @test x.value ≈ 3 atol=TOL - @test p.optval ≈ 6 atol=TOL - @test evaluate(sum([3, 6, 9] / x)) ≈ 6 atol=TOL - end - - @testset "geo mean atom" begin - x = Variable(2) - y = Variable(2) - p = minimize(geomean(x, y), x >= 1, y >= 2) - # not DCP compliant - @test vexity(p) == ConcaveVexity() - p = maximize(geomean(x, y), 1 < x, x < 2, y < 2) - # Just gave it a vector as an objective, not okay - @test_throws Exception solve!(p, solver) - - p = maximize(sum(geomean(x, y)), 1 < x, x < 2, y < 2) - solve!(p, solver) - @test p.optval ≈ 4 atol=TOL - @test evaluate(sum(geomean(x, y))) ≈ 4 atol=TOL - end - - @testset "sqrt atom" begin - x = Variable() - p = maximize(sqrt(x), 1 >= x) - end - - @testset "quad form atom" begin - x = Variable(3, 1) - A = [0.8608 0.3131 0.5458; 0.3131 0.8584 0.5836; 0.5458 0.5836 1.5422] - p = minimize(quadform(x, A), [x >= 1]) - @test vexity(p) == ConvexVexity() - solve!(p, solver) - @test p.optval ≈ 6.1464 atol=TOL - @test (evaluate(quadform(x, A)))[1] ≈ 6.1464 atol=TOL - - x = Variable(3, 1) - A = -1.0*[0.8608 0.3131 0.5458; 0.3131 0.8584 0.5836; 0.5458 0.5836 1.5422] - c = [3 2 4] - p = maximize(c*x , [quadform(x, A) >= -1]) - @test vexity(p) == ConvexVexity() - solve!(p, solver) - @test p.optval ≈ 3.7713 atol=TOL - @test (evaluate(quadform(x, A)))[1] ≈ -1 atol=TOL - end - - @testset "huber atom" begin - x = Variable(3) - p = minimize(sum(huber(x, 1)), x >= 2) - @test vexity(p) == ConvexVexity() - solve!(p, solver) - @test p.optval ≈ 9 atol=TOL - @test evaluate(sum(huber(x, 1))) ≈ 9 atol=TOL - end - - @testset "rational norm atom" begin - A = [1 2 3; -1 2 3] - b = A * ones(3) - x = Variable(3) - p = minimize(norm(x, 4.5), [A * x == b]) - @test vexity(p) == ConvexVexity() - # Solution is approximately x = [1, .93138, 1.04575] - solve!(p, solver) - @test p.optval ≈ 1.2717 atol=TOL - @test evaluate(norm(x, 4.5)) ≈ 1.2717 atol=TOL - end - - @testset "rational norm dual norm" begin - v = [0.463339, 0.0216084, -2.07914, 0.99581, 0.889391] - x = Variable(5) - q = 1.379; # q norm constraint that generates many inequalities - qs = q / (q - 1); # Conjugate to q - p = minimize(x' * v) - p.constraints += (norm(x, q) <= 1) - @test vexity(p) == ConvexVexity() - solve!(p, solver) # Solution is -norm(v, q / (q - 1)) - @test p.optval ≈ -2.144087 atol=TOL - @test sum(evaluate(x' * v)) ≈ -2.144087 atol=TOL - @test evaluate(norm(x, q)) ≈ 1 atol=TOL - @test sum(evaluate(x' * v)) ≈ -(sum(abs.(v) .^ qs) ^ (1 / qs)) atol=TOL - end - - @testset "rational norm atom sum" begin - A = [-0.719255 -0.229089 - -1.33632 -1.37121 - 0.703447 -1.4482] - b = [-1.82041, -1.67516, -0.866884] - q = 1.5 - xvar = Variable(2) - p = minimize(.5 * sumsquares(xvar) + norm(A * xvar - b, q)) - @test vexity(p) == ConvexVexity() - solve!(p, solver) - # Compute gradient, check it is zero(ish) - x_opt = xvar.value - margins = A * x_opt - b - qs = q / (q - 1); # Conjugate - denom = sum(abs.(margins).^q)^(1/qs) - g = x_opt + A' * (abs.(margins).^(q-1) .* sign.(margins)) / denom - @test p.optval ≈ 1.7227 atol=TOL - @test norm(g, 2) ^ 2 ≈ 0 atol=TOL - end - - @testset "norm consistent with Base for matrix variables" begin - A = randn(4, 4) - x = Variable(4, 4) - x.value = A - # Matrix norm - @test evaluate(opnorm(x)) ≈ opnorm(A) atol=TOL - @test evaluate(opnorm(x, 1)) ≈ opnorm(A, 1) atol=TOL - @test evaluate(opnorm(x, 2)) ≈ opnorm(A, 2) atol=TOL - @test evaluate(opnorm(x, Inf)) ≈ opnorm(A, Inf) atol=TOL - # Vector norm - # TODO: Once the deprecation for norm on matrices is removed, remove the `vec` calls - @test evaluate(norm(vec(x), 1)) ≈ norm(vec(A), 1) atol=TOL - @test evaluate(norm(vec(x), 2)) ≈ norm(vec(A), 2) atol=TOL - @test evaluate(norm(vec(x), 7)) ≈ norm(vec(A), 7) atol=TOL - @test evaluate(norm(vec(x), Inf)) ≈ norm(vec(A), Inf) atol=TOL - end - - @testset "Fixed and freed variables" begin - @testset "fix and free addition" begin - x = Variable() - y = Variable() - - p = minimize(x+y, x>=0, y>=0) - solve!(p, solver) - @test p.optval ≈ 0 atol=TOL - - y.value = 4 - fix!(y) - solve!(p, solver) - @test p.optval ≈ 4 atol=TOL - - free!(y) - solve!(p, solver) - @test p.optval ≈ 0 atol=TOL - end - - @testset "fix multiplication" begin - a = [1,2,3,2,1] - x = Variable(length(a)) - gamma = Variable(Positive()) - fix!(gamma, 0.7) - - p = minimize(norm(x-a) + gamma*norm(x[1:end-1] - x[2:end])) - solve!(p, solver) - o1 = p.optval - # x should be very close to a - @test o1 ≈ 0.7 * norm(a[1:end - 1] - a[2:end]) atol=TOL - # increase regularization - fix!(gamma, 1.0) - solve!(p, solver) - o2 = p.optval - # x should be very close to mean(a) - @test o2 ≈ norm(a .- mean(a)) atol=TOL - - @test o1 <= o2 - end - end - - @testset "minimal norm solutions" begin - x = Variable(2) - A = [1 2; 2 4]; - b = [3, 6]; - p = minimize(norm(x, 1), A*x==b) - @test vexity(p) == ConvexVexity() - solve!(p, solver) - @test p.optval ≈ 1.5 atol=TOL - @test evaluate(x) ≈ [0, 1.5] atol=TOL - @test evaluate(norm(x, 1)) ≈ p.optval atol=TOL - if p.solution.has_dual - println("Solution object has dual value, checking for dual problem correctness.") - @test dot(b, p.constraints[1].dual) ≈ p.optval atol=TOL - end - - x = Variable(2) - A = [1 2; 2 4]; - b = [3, 6]; - p = minimize(norm(x, 2), A*x==b) - @test vexity(p) == ConvexVexity() - solve!(p, solver) - @test p.optval ≈ 3/sqrt(5) atol=TOL - @test evaluate(x) ≈ [3/5, 6/5] atol=TOL - @test evaluate(norm(x, 2)) ≈ p.optval atol=TOL - if p.solution.has_dual - println("Solution object has dual value, checking for dual problem correctness.") - @test dot(b, p.constraints[1].dual) ≈ p.optval atol=TOL - end - - x = Variable(2) - A = [1 2; 2 4]; - b = [3, 6]; - p = minimize(norm(x, Inf), A*x==b) - @test vexity(p) == ConvexVexity() - solve!(p, solver) - @test p.optval ≈ 1.0 atol=TOL - @test evaluate(x) ≈ [1, 1] atol=TOL - @test evaluate(norm(x, Inf)) ≈ p.optval atol=TOL - if p.solution.has_dual - println("Solution object has dual value, checking for dual problem correctness.") - @test dot(b, p.constraints[1].dual) ≈ p.optval atol=TOL - end - end - end -end