text:
@graphs begin
	move_wk_cond = """
	The Witch King is not in an *mobile* army but is able to join or create one.
	"""

	wk_prio = """
	Move the Witch King, place it in a valid region with an army.

	Priority:
	1. Army is *mobile*
	2. Army's *target* is in nation at war
	3. Army becomes *mobile* if the Witch King is added
	4. Free Peoples' army at *target* or on the route to *target* does not contain Gandalf the White
	5. Free Peoples' army at *target* or on the route to *target* does not contain a hobbit
	6. Army is adjacent to a *threat*
	7. Army that is conduction a siege
	8. Army is adjacent to its *target*
	9. Highest *value* Shadow army
	10. Random
	"""

	move_nazgul_cond = """
	No Nazgûls are in the Fellowship's region but they are able to move there.
	Or, a Nazgûl is in a non-*mobile* army but is able to join or create one.
	"""

	nazgul_prio = """
	Gather all Nazgûl and place them one at the time.

	Priority:
	1. One, and only one, in the Fellowship's region
	2. Army with leadership *value* less than the number of army units and less than 5
	3. Army which contain the Which King
	4. Shadow stronghold under siege
	5. *Mobile* army
	6. Army adjacent to *threat*
	7. Army whose *target* nation is active
	8. Army that becomes *mobile* if Nazgûl is added
	9. Army not sieging
	10. Army that is adjacent to its *target*
	11. Highest *value* shadow army
	12. Random
	"""

	move_mos_cond = """
	Mouth of Souron is not in a *mobile* army.
	"""

	mos_prio = """
	Move Mouth of Sauron.

	Priority:
	1. Towards army with leadership *value* less than the number of army units and 5
	2. Towards *mobile* army
	3. Towards army adjacent to its *target*
	4. Army which can be reach with this die
	5. Towards closest army
	6. Random
	"""

	################################################################################
	@node character_army = Start() -> lc_1

	@node lc_1 = BinaryCondition("""
								 An *aggressive* army with the Witch King or maximum leadership is adjacent to its *target*.
								 """) -> [n_true = lc_1_yes, n_false = lc_2]
	@node lc_1_yes = UseActiveDie() -> lc_1_action
	@node lc_1_action = PerformAction("""
									  Attack according to the latest statement. Select army at random if several can perform such an attack.
									  """) -> lc_1_end
	@node lc_1_end = End() -> []

	@node lc_2 = BinaryCondition("""
								 A *mobile* army with leadership and a valid move/attack towards *target* exists.
								 """) -> [n_true = lc_2_yes, n_false = lc_3]
	@node lc_2_yes = JumpToGraph("movement_attack_basic") -> lc_3

	################################################################################
	@node character_move = Start() -> lc_3
	@node lc_3 = BinaryCondition("""
								 A Nazgûl or the Witch King is in play.
								 """) -> [n_true = lc_wk, n_false = lc_3_no]
	@node lc_3_no = JumpToGraph("event_cards_preferred") -> lc_3_return
	@node lc_3_return = ReturnFromGraph() -> []

	########################################
	@node character_which_king = Start() -> lc_wk_yes

	@node lc_wk = BinaryCondition(move_wk_cond) -> [n_true = lc_wk_yes, n_false = lc_naz_1]
	@node lc_wk_yes = UseActiveDie() -> lc_wk_action
	@node lc_wk_action = PerformAction(wk_prio) -> lc_naz_2

	########################################
	@node lc_naz_1 = BinaryCondition(move_nazgul_cond) -> [n_true = lc_naz_1_yes, n_false = lc_mos_1]
	@node lc_naz_1_yes = UseActiveDie() -> lc_naz_1_action
	@node lc_naz_1_action = PerformAction(nazgul_prio) -> lc_mos_2

	@node lc_naz_2 = BinaryCondition(move_nazgul_cond) -> [n_true = lc_naz_2_yes, n_false = lc_mos_3]
	@node lc_naz_2_yes = UseActiveDie() -> lc_naz_2_action
	@node lc_naz_2_action = PerformAction(nazgul_prio) -> lc_mos_4

	########################################
	@node lc_mos_1 = BinaryCondition(move_mos_cond) -> [n_true = lc_mos_1_yes, n_false = lc_play_card]
	@node lc_mos_1_yes = UseActiveDie() -> lc_mos_1_action
	@node lc_mos_1_action = PerformAction(mos_prio) -> lc_mos_end_1

	@node lc_mos_2 = BinaryCondition(move_mos_cond) -> [n_true = lc_mos_2_yes, n_false = lc_mos_end_2]
	@node lc_mos_2_yes = UseActiveDie() -> lc_mos_2_action
	@node lc_mos_2_action = PerformAction(mos_prio) -> lc_mos_end_2

	@node lc_mos_3 = BinaryCondition(move_mos_cond) -> [n_true = lc_mos_3_yes, n_false = lc_mos_end_3]
	@node lc_mos_3_yes = UseActiveDie() -> lc_mos_3_action
	@node lc_mos_3_action = PerformAction(mos_prio) -> lc_mos_end_3

	@node lc_mos_4 = BinaryCondition(move_mos_cond) -> [n_true = lc_mos_4_yes, n_false = lc_mos_end_4]
	@node lc_mos_4_yes = UseActiveDie() -> lc_mos_4_action
	@node lc_mos_4_action = PerformAction(mos_prio) -> lc_mos_end_4

	@node lc_mos_end_1 = End() -> []
	@node lc_mos_end_2 = End() -> []
	@node lc_mos_end_3 = End() -> []
	@node lc_mos_end_4 = End() -> []

	########################################
	@node lc_play_card = JumpToGraph("event_cards_preferred") -> lc_play_card_return
	@node lc_play_card_return = ReturnFromGraph() -> []

	################################################################################
	@node character_wk_prio = Start() -> lc_wk_prio
	@node lc_wk_prio = PerformAction(wk_prio) -> lc_wk_prio_end
	@node lc_wk_prio_end = End() -> []

	################################################################################
	@node character_nazgul_prio = Start() -> lc_nazgul_prio
	@node lc_nazgul_prio = PerformAction(nazgul_prio) -> lc_nazgul_prio_end
	@node lc_nazgul_prio_end = End() -> []

	################################################################################
	@node character_mos_prio = Start() -> lc_mos_prio
	@node lc_mos_prio = PerformAction(mos_prio) -> lc_mos_prio_end
	@node lc_mos_prio_end = End() -> []
end

max_stars_repo_name: mvmorin/queller-bot

max_issues_repo_name: mvmorin/queller-bot

max_forks_repo_name: mvmorin/queller-bot


text:
# Autogenerated wrapper script for Ncurses_jll for armv7l-linux-gnueabihf
export libform, libmenu, libncurses, libpanel

JLLWrappers.@generate_wrapper_header("Ncurses")
JLLWrappers.@declare_library_product(libform, "libformw.so.6")
JLLWrappers.@declare_library_product(libmenu, "libmenuw.so.6")
JLLWrappers.@declare_library_product(libncurses, "libncursesw.so.6")
JLLWrappers.@declare_library_product(libpanel, "libpanelw.so.6")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libform,
        "lib/libform.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmenu,
        "lib/libmenu.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libncurses,
        "lib/libncurses.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libpanel,
        "lib/libpanel.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()

ext: jl


text:
# by Lazaro Alonso
# port from the original code used in:
# L Alonso, et. al. https://doi.org/10.1093/comnet/cnx053

# using CairoMakie # HIDE
using LinearAlgebra, Random, CairoMakie
CairoMakie.activate!()
function RRGAdjacencyM(;lengthBase = 1, radius = 0.1, nodes = 500, rseed = 123)
    Random.seed!(rseed)
    xy = rand(nodes, 2)
    x = lengthBase .* xy[:,1]
    y = xy[:,2] ./ lengthBase
    matrixAdjDiag = Diagonal(√2 * randn(nodes))
    matrixAdj = zeros(nodes, nodes)
    for point in 1:nodes-1
        distance = sqrt.((x[point + 1:end] .- x[point]).^2 .+ (y[point + 1:end] .- y[point]).^2)
        dindx = findall(distance .<= radius) .+ point
        if length(dindx) > 0
            rnd = randn(length(dindx))
            matrixAdj[point, dindx] = rnd
            matrixAdj[dindx, point] = rnd
        end
    end
    return (matrixAdj .+ matrixAdjDiag, x, y)
end
adjacencyM, x, y = RRGAdjacencyM()

function getGraphEdges(adjMatrix, x, y)
    xyos = []
    weights = []
    for i in 1:length(x), j in i+1:length(x)
        if adjMatrix[i,j] != 0.0
            push!(xyos, [x[i], y[i]])
            push!(xyos, [x[j], y[j]])
            push!(weights, adjMatrix[i,j])
            push!(weights, adjMatrix[i,j])
        end
    end
    return (Point2f0.(xyos), Float32.(weights))
end

function plotGraph(adjacencyM, x, y)
    cmap = (:Spectral_11, 0.75)
    adjmin = minimum(adjacencyM)
    adjmax = maximum(adjacencyM)
    diagValues = diag(adjacencyM)
    segm, weights = getGraphEdges(adjacencyM, x, y)

    fig, ax, pltobj = linesegments(segm, color = weights, colormap = cmap,
        linewidth = abs.(weights)/2, colorrange = (adjmin, adjmax),
        figure = (;resolution = (500,400)), axis = (;aspect = DataAspect()))
    scatter!(ax, x, y, color = diagValues, markersize = 3*abs.(diagValues),
        colorrange = (adjmin, adjmax), colormap = cmap)
    limits!(ax, -0.02,1.02,-0.02,1.02)
    Colorbar(fig[1,2], pltobj, label = "weights")
    fig
end
#plotGraph(adjacencyM, x, y)
rrgraph = with_theme(theme_black()) do
    plotGraph(adjacencyM, x, y)
end
save(joinpath(@__DIR__, "output", "RRGraph.png"), rrgraph, px_per_unit = 2.0) # HIDE
using Pkg # HIDE
Pkg.status(["CairoMakie", "LinearAlgebra", "Random"]) # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
# Copyright 2018 IBM Corp.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

module OntologyRDF
export ontology_to_rdf, presentation_to_rdf, annotation_to_rdf,
  wiring_diagram_to_rdf

using Serd
using Catlab

using ..OntologyDBs

# Submodules
############

include("RDFUtils.jl")
include("WiringRDF.jl")
include("ConceptRDF.jl")
include("ConceptPROV.jl")
include("AnnotationRDF.jl")

using .RDFUtils
using .WiringRDF
using .ConceptRDF
using .ConceptPROV
using .AnnotationRDF

# Ontology RDF
##############

""" Convert ontology's concepts and annotations into RDF/OWL ontology.
"""
function ontology_to_rdf(db::OntologyDB, prefix::RDF.Prefix;
                         include_provenance::Bool=true,
                         include_wiring_diagrams::Bool=true)::Vector{<:RDF.Statement}
  # Create RDF triples for concepts.
  function concept_labels(expr, node::RDF.Node)::Vector{<:RDF.Statement}
    # Add RDFS labels for concept.
    doc = concept_document(db, first(expr))
    rdfs_labels(doc, node)
  end
  stmts = presentation_to_rdf(concepts(db), prefix; extra_rdf=concept_labels)

  # Create RDF triples for concept hierarchy based on PROV-O.
  if include_provenance
    append!(stmts, presentation_to_prov(concepts(db), prefix))
  end

  # Create RDF triples for annotations.
  for note in annotations(db)
    append!(stmts, annotation_to_rdf(note, prefix;
      include_provenance=include_provenance,
      include_wiring_diagrams=include_wiring_diagrams))

    # Add RDFS labels for annotation.
    doc = annotation_document(db, note.name)
    node = annotation_rdf_node(note, prefix)
    append!(stmts, rdfs_labels(doc, node))
  end

  return stmts
end

end

max_stars_repo_name: cclauss/semanticflowgraph

max_issues_repo_name: cclauss/semanticflowgraph

max_forks_repo_name: cclauss/semanticflowgraph


text:
using MultivariatePolynomials
using DynamicPolynomials
#using TypedPolynomials
using BenchmarkTools

#   let
#       using MultiPoly
#       x, y = generators(MPoly{Float64}, :x, :y)
#       p = x + y + 2x^2 + y^3
#
#       println(@benchmark(evaluate($p, 1.0, 2.0)))
#   end

let
    @polyvar x y
    p = x + y + 2x^2 + y^3

    varst = (x, y)
    varsv = [x, y]
    valst = (1.0, 2.0)
    valsv = [1.0, 2.0]

    display(@benchmark(($p)(variables($p) => $valst)))
    display(@benchmark(($p)(variables($p) => $valsv)))
    display(@benchmark(($p)($varst => $valst)))
    display(@benchmark(($p)($varsv => $valsv)))

#    vars_out_of_order = [vars[2], vars[1]]
#    println(@benchmark(($p)($vals, $vars_out_of_order)))
end


text:
using SMLMData
using Test

@testset "SMLMData.jl" begin
    # Write your tests here.
end


text:
using Flux, DiffEqFlux, DifferentialEquations, Plots, StatsBase, RecursiveArrayTools

const u0 = Float32[2.; 0.]
const datasize = 1000
const batchsize = 30
const batchtime = 10
const tspan = (0.0f0,25f0)

function trueODEfunc(du,u,p,t)
    true_A = [-0.1 2.0; -2.0 -0.1]
    du .= ((u.^3)'true_A)'
end
t = range(tspan[1],tspan[2],length=datasize)
prob = ODEProblem(trueODEfunc,u0,tspan)
true_sol = solve(prob,Tsit5(),saveat=t)
ode_data = Array(true_sol)

function get_batch()
    batch_0s = sample(1:datasize-batchtime, batchsize, replace=false,ordered=false)
    batch_idxs = range.(batch_0s,batch_0s.+(batchtime-1))
    batch_ts = [t[i] for i in batch_idxs]
    batch_u0 = [true_sol[u0] for u0 in batch_0s]
    batch_u = [true_sol[:,idxs] for idxs in batch_idxs]
    return batch_u0, batch_u, batch_ts
end
batch_u0, batch_u, batch_ts = get_batch();

dudt = Chain(x -> x.^3,
             Dense(2,50,tanh),
             Dense(50,2))
ps = Flux.params(dudt)
function n_ode(batch_u0, batch_t)
    neural_ode(dudt,batch_u0,tspan,Tsit5(),
               saveat=batch_t, reltol=1e-7,abstol=1e-9)
end
#= n_ode = x->neural_ode(dudt,x,tspan,Tsit5(), =#
                      #= saveat=t,reltol=1e-7,abstol=1e-9) =#

n_ode(batch_u0, batch_ts)

pred = n_ode(u0,t) # Get the prediction using the correct initial condition
scatter(t,ode_data[1,:],label="data")
scatter!(t,Flux.data(pred[1,:]),label="prediction")

function predict_n_ode(batch_u0,batch_t)
  n_ode(batch_u0, batch_t)
end
loss_n_ode(batch_u0, batch_u,batch_t) = sum(abs2,batch_u .- predict_n_ode(batch_u0,batch_t))

data = Iterators.repeated(get_batch(), 1000)
opt = ADAM(0.1)
cb = function () #callback function to observe training
    bu0,bu,bt = get_batch()
    display(loss_n_ode(bu0,bu,bt))
  # plot current prediction against data
  #= cur_pred = Flux.data(predict_n_ode(bu0,bt)) =#
  #= pl = scatter(bt,ode_data[1,:],label="data") =#
  #= scatter!(pl,bt,cur_pred[1,:],label="prediction") =#
  #= display(plot(pl)) =#
end

# Display the ODE with the initial parameter values.
cb()

Flux.train!(loss_n_ode, ps, data, opt, cb = cb)

max_stars_repo_name: jessebett/model-zoo

max_issues_repo_name: jessebett/model-zoo

max_forks_repo_name: jessebett/model-zoo


text:
# 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]; numeric_type = T)

    # @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; numeric_type = T)

    # @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; numeric_type = T)

    # # @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; numeric_type = T)

    # @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; numeric_type = T)

    # @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; numeric_type = T)

    # @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}
    x = Variable(Positive())
    y = Variable((3, 3))
    p = minimize(x + y[1, 1], isposdef(y), x >= 1, y[2, 1] == 1; numeric_type = T)

    # @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; numeric_type = T)

    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_nuclear_norm_atom_complex(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    A = [1 2im 3 4; 4im 3im 2 1; 4 5 6 7]
    y = ComplexVariable(3, 4)
    p = minimize(nuclearnorm(y), y == A; numeric_type = T)

    if test
        @test vexity(p) == ConvexVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ sum(svdvals(A)) atol=atol rtol=rtol
        @test evaluate(nuclearnorm(y)) ≈ sum(svdvals(A)) 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; numeric_type = T)

    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_operator_norm_atom_complex(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    A = [1 2im 3 4; 4im 3im 2 1; 4 5 6 7]
    y = ComplexVariable(3, 4)
    p = minimize(opnorm(y), y == A; numeric_type = T)

    if test
        @test vexity(p) == ConvexVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ maximum(svdvals(A)) atol=atol rtol=rtol
        @test evaluate(opnorm(y)) ≈ maximum(svdvals(A)) 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; numeric_type = T)

    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_dual_lambda_max_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    y = Semidefinite(3)
    p = minimize(eigmax(y), y[1,1]>=4; numeric_type = T)

    if test
        @test vexity(p) == ConvexVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 4 atol=atol rtol=rtol
        @test evaluate(eigmax(y)) ≈ 4 atol=atol rtol=rtol
    end

    # https://github.com/jump-dev/Convex.jl/issues/337
    x = ComplexVariable(2, 2)
    p = minimize( eigmax(x), [ x[1,2] == im, x[2,2] == 1.0, x ⪰ - eye(2) ]; numeric_type = T)
    handle_problem!(p)
    if test
        @test p.optval ≈ 1.5 atol=atol rtol=rtol
        @test p.constraints[1].dual ≈ im atol=atol rtol=rtol
        @test p.constraints[2].dual ≈ 0.75 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(eigmin(y), tr(y)<=6; numeric_type = T)

    if test
        @test vexity(p) == ConvexVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 2 atol=atol rtol=rtol
        @test evaluate(eigmin(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); numeric_type = T)

    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), eigmax(P) <= 2, x[1] >= 1; numeric_type = T)

    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; numeric_type = T)

    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) + sumlargesteigs(x,0), [x[i,:] >= i for i=1:3]...; numeric_type = T)

    handle_problem!(p)

    if test
        @test p.optval ≈ 8.4853 atol=atol rtol=rtol
    end

    A = [1   -2im  3   4
         2im  7    3im 5
         3   -3im  2   9
         4    5    9   4]

    x = ComplexVariable(4, 4)
    p = minimize(sumlargesteigs(x, 3), x == A; numeric_type = T)

    handle_problem!(p)

    if test
        @test p.optval ≈ sum(eigvals(A)[2:end]) atol=atol rtol=rtol
    end

    x1 = Semidefinite(3)
    p1 = minimize(eigmax(x1), x1[1,1]>=4; numeric_type = T)

    handle_problem!(p1)

    x2 = Semidefinite(3)
    p2 = minimize(sumlargesteigs(x2, 1), x2[1,1]>=4; numeric_type = T)

    handle_problem!(p2)

    if test
        @test p1.optval ≈ p2.optval atol=atol rtol=rtol
    end

    x1 = Semidefinite(3)
    p1 = minimize(eigmax(x1), [x1[i,:] >= i for i=1:3]...; numeric_type = T)

    handle_problem!(p1)

    x2 = Semidefinite(3)
    p2 = minimize(sumlargesteigs(x2, 1), [x2[i,:] >= i for i=1:3]...; numeric_type = T)

    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; numeric_type = T)

    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; numeric_type = T)

    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; numeric_type = T)

    handle_problem!(p1)

    if test
        x1 = copy(evaluate(x))
    end

    p2 = minimize(sum(x), real(A)*x == real(b), imag(A)*x==imag(b), x>=0; numeric_type = T)

    handle_problem!(p2)
    if test
        x2 = evaluate(x)
        @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; numeric_type = T)

    handle_problem!(p1)

    if test
        x1 = evaluate(x)
    end

    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; numeric_type = T)

    handle_problem!(p2)

    if test
        real_diff = real(x1) - evaluate(xr)
        @test real_diff ≈ zeros(10, 1) atol=atol rtol=rtol

        imag_diff = imag(x1) - evaluate(xi)
        @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; numeric_type = T)

    handle_problem!(p)
    if test
        @test p.status == MOI.OPTIMAL
        @test evaluate(ρ) ≈ [ 1. 0.; 0.  1.] atol=atol rtol=rtol
        @test p.optval ≈ 0 atol=atol rtol=rtol
    end
end

@add_problem sdp function sdp_socp_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; numeric_type = T)

    handle_problem!(p)
    if test
        @test p.optval ≈ 0 atol=atol rtol=rtol
        @test evaluate(objective) ≈ 0 atol=atol rtol=rtol

        real_diff = real(evaluate(x)) - real(a)
        imag_diff = imag(evaluate(x)) - 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_socp_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; numeric_type = T)

    handle_problem!(p)
    if test
        @test p.optval ≈ 0 atol=atol rtol=rtol
        @test evaluate(objective) ≈ 0.0 atol=atol rtol=rtol

        real_diff = real.(evaluate(x)) - real.(a)
        imag_diff = imag.(evaluate(x)) - 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_socp_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; numeric_type = T)

    handle_problem!(p)
    if test
        @test p.optval ≈ 0 atol=atol rtol=rtol
        @test evaluate(objective) ≈ 0.0 atol=atol rtol=rtol

        real_diff = real(evaluate(x)) .- real(a)
        imag_diff = imag(evaluate(x)) .- imag(a)
        @test real_diff ≈ 0.0 atol=atol rtol=rtol
        @test imag_diff ≈ 0.0 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; numeric_type = T)

    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.(evaluate(x)) - real.(posA)
        imag_diff = imag.(evaluate(x)) - imag.(posA)
        @test real_diff ≈ zeros(n, n) atol=atol rtol=rtol
        @test imag_diff ≈ zeros(n, n) atol=atol rtol=rtol
    end
end

@add_problem sdp function sdp_geom_mean_hypocone_argcheck(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    if test
        @test_throws DimensionMismatch GeomMeanHypoCone(   zeros(2,3),    zeros(2,3), 1//2)
        @test_throws DimensionMismatch GeomMeanHypoCone(   zeros(2,3), Variable(2,3), 1//2)
        @test_throws DimensionMismatch GeomMeanHypoCone(Variable(2,3),    zeros(2,3), 1//2)
        @test_throws DimensionMismatch GeomMeanHypoCone(Variable(2,3), Variable(2,3), 1//2)

        @test_throws DimensionMismatch GeomMeanHypoCone(   zeros(2,2),    zeros(3,3), 1//2)
        @test_throws DimensionMismatch GeomMeanHypoCone(   zeros(2,2), Variable(3,3), 1//2)
        @test_throws DimensionMismatch GeomMeanHypoCone(Variable(2,2),    zeros(3,3), 1//2)
        @test_throws DimensionMismatch GeomMeanHypoCone(Variable(2,2), Variable(2,3), 1//2)

        @test_throws DimensionMismatch Variable(2,2) in GeomMeanHypoCone(Variable(3,3), Variable(3,3), 1//2)

        @test_throws DomainError GeomMeanHypoCone(Variable(3,3), Variable(3,3), -1//2)
        @test_throws DomainError GeomMeanHypoCone(Variable(3,3), Variable(3,3), 3//2)
    end
end

@add_problem sdp function sdp_geom_mean_epicone_argcheck(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    if test
        @test_throws DimensionMismatch GeomMeanEpiCone(   zeros(2,3),    zeros(2,3), -1//2)
        @test_throws DimensionMismatch GeomMeanEpiCone(   zeros(2,3), Variable(2,3), -1//2)
        @test_throws DimensionMismatch GeomMeanEpiCone(Variable(2,3),    zeros(2,3), -1//2)
        @test_throws DimensionMismatch GeomMeanEpiCone(Variable(2,3), Variable(2,3), -1//2)

        @test_throws DimensionMismatch GeomMeanEpiCone(   zeros(2,2),    zeros(3,3), -1//2)
        @test_throws DimensionMismatch GeomMeanEpiCone(   zeros(2,2), Variable(3,3), -1//2)
        @test_throws DimensionMismatch GeomMeanEpiCone(Variable(2,2),    zeros(3,3), -1//2)
        @test_throws DimensionMismatch GeomMeanEpiCone(Variable(2,2), Variable(2,3), -1//2)

        @test_throws DimensionMismatch Variable(2,2) in GeomMeanEpiCone(Variable(3,3), Variable(3,3), -1//2)
        @test_throws DomainError GeomMeanEpiCone(Variable(3,3), Variable(3,3), -3//2)
        @test_throws DomainError GeomMeanEpiCone(Variable(3,3), Variable(3,3), 1//2)
        @test_throws DomainError GeomMeanEpiCone(Variable(3,3), Variable(3,3), 5//2)
    end
end

@add_problem sdp function sdp_relative_entropy_argcheck(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    if test
        @test_throws DimensionMismatch RelativeEntropyEpiCone(   zeros(2,3),    zeros(2,3), 3, 3, eye(2))
        @test_throws DimensionMismatch RelativeEntropyEpiCone(   zeros(2,3), Variable(2,3), 3, 3, eye(2))
        @test_throws DimensionMismatch RelativeEntropyEpiCone(Variable(2,3),    zeros(2,3), 3, 3, eye(2))
        @test_throws DimensionMismatch RelativeEntropyEpiCone(Variable(2,3), Variable(2,3), 3, 3, eye(2))

        @test_throws DimensionMismatch RelativeEntropyEpiCone(   zeros(2,2),    zeros(3,3), 3, 3)
        @test_throws DimensionMismatch RelativeEntropyEpiCone(   zeros(2,2), Variable(3,3), 3, 3)
        @test_throws DimensionMismatch RelativeEntropyEpiCone(Variable(2,2),    zeros(3,3), 3, 3)
        @test_throws DimensionMismatch RelativeEntropyEpiCone(Variable(2,2), Variable(2,3), 3, 3)

        @test_throws DimensionMismatch RelativeEntropyEpiCone(   zeros(3,3),    zeros(3,3), 3, 3, zeros(2, 2))
        @test_throws DimensionMismatch RelativeEntropyEpiCone(   zeros(3,3), Variable(3,3), 3, 3, zeros(2, 2))
        @test_throws DimensionMismatch RelativeEntropyEpiCone(Variable(3,3),    zeros(3,3), 3, 3, zeros(2, 2))
        @test_throws DimensionMismatch RelativeEntropyEpiCone(Variable(3,3), Variable(2,3), 3, 3, zeros(2, 2))

        @test_throws DimensionMismatch RelativeEntropyEpiCone(   zeros(3,3),    zeros(3,3), 3, 3, zeros(2, 3))
        @test_throws DimensionMismatch RelativeEntropyEpiCone(   zeros(3,3), Variable(3,3), 3, 3, zeros(2, 3))
        @test_throws DimensionMismatch RelativeEntropyEpiCone(Variable(3,3),    zeros(3,3), 3, 3, zeros(2, 3))
        @test_throws DimensionMismatch RelativeEntropyEpiCone(Variable(3,3), Variable(2,3), 3, 3, zeros(2, 3))

        @test_throws DimensionMismatch RelativeEntropyEpiCone(   zeros(3,3),    zeros(3,3), 3, 3, zeros(2))
        @test_throws DimensionMismatch RelativeEntropyEpiCone(   zeros(3,3), Variable(3,3), 3, 3, zeros(2))
        @test_throws DimensionMismatch RelativeEntropyEpiCone(Variable(3,3),    zeros(3,3), 3, 3, zeros(2))
        @test_throws DimensionMismatch RelativeEntropyEpiCone(Variable(3,3), Variable(2,3), 3, 3, zeros(2))

        @test_throws DimensionMismatch Variable(2,2) in RelativeEntropyEpiCone(Variable(3,3), Variable(3,3), 3, 3)
    end
end

@add_problem sdp function sdp_relative_entropy(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    for cplx in [false, true]
        n = 4
        X = eye(n)
        if cplx
            Y = ComplexVariable(n,n)
        else
            Y = Variable(n,n)
        end

        c1 = eye(n) in RelativeEntropyEpiCone(X, Y)
        objective = real(tr(Y))
        p = minimize(objective, c1; numeric_type = T)

        handle_problem!(p)

        if test
            @test Y.value ≈ eye(n)*exp(-1) atol=atol rtol=rtol
        end
    end
end

@add_problem sdp function sdp_geom_mean_hypocone_real_0(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    n = 4
    A = Variable(n,n)
    B = randn(n,n)
    B = B * B' # now A is positive semidefinite
    B += 0.2 * I # prevent numerical instability

    c1 = eye(n) in GeomMeanHypoCone(A, B, 0)
    objective = tr(A)
    p = minimize(objective, c1; numeric_type = T)

    handle_problem!(p)

    if test
        @test A.value ≈ eye(n) atol=atol rtol=rtol
    end
end

@add_problem sdp function sdp_geom_mean_hypocone_real_1(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    n = 4
    A = Variable(n,n)
    B = randn(n,n)
    B = B * B' # now A is positive semidefinite
    B += 0.2 * I # prevent numerical instability

    c1 = eye(n) in GeomMeanHypoCone(B, A, 1)
    objective = tr(A)
    p = minimize(objective, c1; numeric_type = T)

    handle_problem!(p)

    if test
        @test A.value ≈ eye(n) atol=atol rtol=rtol
    end
end

@add_problem sdp function sdp_geom_mean_hypocone_real_1_2(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    n = 4
    A = randn(n,n)
    A = A * A' # now A is positive semidefinite
    A += 0.2 * I # prevent numerical instability
    B = Variable(n,n)

    c1 = eye(n) in GeomMeanHypoCone(A, B, 1//2)
    objective = tr(B)
    p = minimize(objective, c1; numeric_type = T)

    handle_problem!(p)

    if test
        @test B.value ≈ A^-1 atol=atol rtol=rtol
    end
end

@add_problem sdp function sdp_geom_mean_hypocone_real_3_8(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    n = 4
    A = randn(n,n)
    A = A * A' # now A is positive semidefinite
    A += 0.2 * I # prevent numerical instability
    B = Variable(n,n)

    c1 = eye(n) in GeomMeanHypoCone(A, B, 3//8)
    objective = tr(B)
    p = minimize(objective, c1; numeric_type = T)

    handle_problem!(p)

    # A #_t B = I  =>  B = A^(1-1/t)
    if test
        @test B.value ≈ A^(-5//3) atol=atol rtol=rtol
    end
end

@add_problem sdp function sdp_geom_mean_hypocone_real_3_5(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    n = 4
    A = randn(n,n)
    A = A * A' # now A is positive semidefinite
    A += 0.2 * I # prevent numerical instability
    B = Variable(n,n)

    c1 = eye(n) in GeomMeanHypoCone(A, B, 3//5)
    objective = tr(B)
    p = minimize(objective, c1; numeric_type = T)

    handle_problem!(p)

    # A #_t B = I  =>  B = A^(1-1/t)
    if test
        @test B.value ≈ A^(-2//3) atol=atol rtol=rtol
    end
end

@add_problem sdp function sdp_geom_mean_hypocone_cplx_1_2(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    n = 4
    A = randn(ComplexF64, n, n)
    A = A * A' # now A is positive semidefinite
    A += 0.2 * I # prevent numerical instability
    B = ComplexVariable(n,n)

    c1 = eye(n) in GeomMeanHypoCone(A, B, 1//2)
    objective = real(tr(B))
    p = minimize(objective, c1; numeric_type = T)

    handle_problem!(p)

    if test
        @test real.(B.value) ≈ real.(A^-1) atol=atol rtol=rtol
        @test imag.(B.value) ≈ imag.(A^-1) atol=atol rtol=rtol
    end
end

@add_problem sdp function sdp_geom_mean_hypocone_cplx_3_8(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    n = 4
    A = randn(ComplexF64, n, n)
    A = A * A' # now A is positive semidefinite
    A += 0.2 * I # prevent numerical instability
    B = ComplexVariable(n,n)

    c1 = eye(n) in GeomMeanHypoCone(A, B, 3//8)
    objective = real(tr(B))
    p = minimize(objective, c1; numeric_type = T)

    handle_problem!(p)

    # A #_t B = I  =>  B = A^(1-1/t)
    if test
        @test real.(B.value) ≈ real.(A^(-5//3)) atol=atol rtol=rtol
        @test imag.(B.value) ≈ imag.(A^(-5//3)) atol=atol rtol=rtol
    end
end

@add_problem sdp function sdp_geom_mean_hypocone_cplx_3_5(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    n = 4
    A = randn(ComplexF64, n, n)
    A = A * A' # now A is positive semidefinite
    A += 0.2 * I # prevent numerical instability
    B = ComplexVariable(n,n)

    c1 = eye(n) in GeomMeanHypoCone(A, B, 3//5)
    objective = real(tr(B))
    p = minimize(objective, c1; numeric_type = T)

    handle_problem!(p)

    # A #_t B = I  =>  B = A^(1-1/t)
    if test
        @test real.(B.value) ≈ real.(A^(-2//3)) atol=atol rtol=rtol
        @test imag.(B.value) ≈ imag.(A^(-2//3)) atol=atol rtol=rtol
    end
end

@add_problem sdp function sdp_geom_mean_hypocone_fullhyp(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    # These matrices satisfy A^(1/2) ⪰ B and so (A,eye(2),B) \in hyp_{1/2}
    # A and B however do not satisfy A ⪰ B^2 and so [A B; B eye(2)] not
    # psd, and using fullhyp=false will give an infeasible SDP.
    A = [6.25 0; 0 16]
    B = [2 1; 1 2]
    S = Semidefinite(2)

    p = minimize(0, [B in GeomMeanHypoCone(A, eye(2), 1//2, false)]; numeric_type = T)
    handle_problem!(p)
    if test
        @test p.status == MOI.INFEASIBLE
    end

    p = minimize(0, [B in GeomMeanHypoCone(A, eye(2), 1//2)]; numeric_type = T)
    handle_problem!(p)
    if test
        @test p.status == MOI.OPTIMAL
    end
end

@add_problem sdp function sdp_geom_mean_epicone_real_neg1_optA(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    n = 3
    B = randn(n,n)
    B = B * B' # now B is positive semidefinite
    B += 0.2 * I # prevent numerical instability
    A = Variable(n,n)

    c1 = eye(n) in GeomMeanEpiCone(A, B, -1)
    objective = tr(A)
    p = maximize(objective, c1; numeric_type = T)

    handle_problem!(p)

    # A #_t B = I  =>  B = A^(1-1/t)
    if test
        @test B ≈ A.value^2 atol=atol rtol=rtol
    end
end

@add_problem sdp function sdp_geom_mean_epicone_real_neg1_optB(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    n = 3
    A = randn(n,n)
    A = A * A' # now A is positive semidefinite
    A += 0.2 * I # prevent numerical instability
    A /= tr(A) # solver has problems if B is large
    B = Variable(n,n)

    c1 = eye(n) in GeomMeanEpiCone(A, B, -1)
    objective = tr(B)
    p = minimize(objective, c1; numeric_type = T)

    handle_problem!(p)

    # A #_t B = I  =>  B = A^(1-1/t)
    if test
        @test B.value ≈ A^2 atol=atol rtol=rtol
    end
end

@add_problem sdp function sdp_geom_mean_epicone_real_neg3_5_optA(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    n = 3
    B = randn(n,n)
    B = B * B' # now B is positive semidefinite
    B += 0.2 * I # prevent numerical instability
    A = Variable(n,n)

    c1 = eye(n) in GeomMeanEpiCone(A, B, -3//5)
    objective = tr(A)
    p = maximize(objective, c1; numeric_type = T)

    handle_problem!(p)

    # A #_t B = I  =>  B = A^(1-1/t)
    if test
        @test B ≈ A.value^(8//3) atol=atol rtol=rtol
    end
end

@add_problem sdp function sdp_geom_mean_epicone_real_neg3_5_optB(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    n = 3
    A = randn(n,n)
    A = A * A' # now A is positive semidefinite
    A += 0.2 * I # prevent numerical instability
    A /= tr(A) # solver has problems if B is large
    B = Variable(n,n)

    c1 = eye(n) in GeomMeanEpiCone(A, B, -3//5)
    objective = tr(B)
    p = minimize(objective, c1; numeric_type = T)

    handle_problem!(p)

    # A #_t B = I  =>  B = A^(1-1/t)
    if test
        @test B.value ≈ A^(8//3) atol=atol rtol=rtol
    end
end

@add_problem sdp function sdp_geom_mean_epicone_real_8_5_optA(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    n = 3
    B = randn(n,n)
    B = B * B' # now B is positive semidefinite
    B += 0.2 * I # prevent numerical instability
    B /= tr(B) # solver has problems if B is large
    A = Variable(n,n)

    c1 = eye(n) in GeomMeanEpiCone(A, B, 8//5)
    objective = tr(A)
    p = minimize(objective, c1; numeric_type = T)

    handle_problem!(p)

    # A #_t B = I  =>  B = A^(1-1/t)
    if test
        @test B ≈ A.value^(3//8) atol=atol rtol=rtol
    end
end

@add_problem sdp function sdp_geom_mean_epicone_real_8_5_optB(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    n = 3
    A = randn(n,n)
    A = A * A' # now A is positive semidefinite
    A += 0.2 * I # prevent numerical instability
    B = Variable(n,n)

    c1 = eye(n) in GeomMeanEpiCone(A, B, 8//5)
    objective = tr(B)
    p = maximize(objective, c1; numeric_type = T)

    handle_problem!(p)

    # A #_t B = I  =>  B = A^(1-1/t)
    if test
        @test B.value ≈ A^(3//8) atol=atol rtol=rtol
    end
end

@add_problem sdp function sdp_quantum_relative_entropy_argcheck(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    if test
        @test_throws DimensionMismatch quantum_relative_entropy(   zeros(2,3),    zeros(2,3))
        @test_throws DimensionMismatch quantum_relative_entropy(   zeros(2,3), Variable(2,3))
        @test_throws DimensionMismatch quantum_relative_entropy(Variable(2,3),    zeros(2,3))
        @test_throws DimensionMismatch quantum_relative_entropy(Variable(2,3), Variable(2,3))

        @test_throws DimensionMismatch quantum_relative_entropy(   zeros(2,2),    zeros(3,3))
        @test_throws DimensionMismatch quantum_relative_entropy(   zeros(2,2), Variable(3,3))
        @test_throws DimensionMismatch quantum_relative_entropy(Variable(2,2),    zeros(3,3))
        @test_throws DimensionMismatch quantum_relative_entropy(Variable(2,2), Variable(2,3))

        z = zeros(2,2)
        v = Variable(2,2)
        nh = [ 1 1 ; 0 1 ] # not hermitian
        np = [ 1 2 ; 2 1 ] # not positive semidefinite
        @test_throws DomainError quantum_relative_entropy(z, nh)
        @test_throws DomainError quantum_relative_entropy(z, np)
        @test_throws DomainError quantum_relative_entropy(nh, z)
        @test_throws DomainError quantum_relative_entropy(np, z)
        @test_throws DomainError quantum_relative_entropy(v, nh)
        @test_throws DomainError quantum_relative_entropy(v, np)
        @test_throws DomainError quantum_relative_entropy(nh, v)
        @test_throws DomainError quantum_relative_entropy(np, v)
    end
end

function sdp_quantum_relative_entropy_impl(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}, lowrank::Bool, mode::Integer) where {T, test}
    # Too slow for n>3.  Anything that can be done to speed it up?
    n = 3
    X = randn(ComplexF64, n, n)
    if lowrank
        X[:,1] .= 0
    end
    X = X * X' # now X is positive semidefinite
    X /= tr(X) # make it a quantum state

    A = ComplexVariable(n,n)
    B = ComplexVariable(n,n)
    constraints = Array{Constraint, 1}([
        A == X,
    ])
    if mode != 5
        push!(constraints, tr(B) == 1)
    end

    if mode == 1
        objective = quantum_relative_entropy(A, B)
    elseif mode == 2
        objective = quantum_relative_entropy(X, B)
    elseif mode == 3
        objective = quantum_relative_entropy(B, A)
    elseif mode == 4
        objective = quantum_relative_entropy(B, X)
    elseif mode == 5
        objective = quantum_relative_entropy(X, X)
    end

    p = minimize(objective, constraints; numeric_type = T)
    handle_problem!(p)

    if test
        @test real.(evaluate(A)) ≈ real.(X) atol=atol rtol=rtol
        @test imag.(evaluate(A)) ≈ imag.(X) atol=atol rtol=rtol
        if mode != 5
            @test real.(evaluate(B)) ≈ real.(X) atol=atol rtol=rtol
            @test imag.(evaluate(B)) ≈ imag.(X) atol=atol rtol=rtol
        end
        @test p.optval ≈ 0 atol=atol rtol=rtol
        if mode == 1
            @test p.optval ≈ evaluate(quantum_relative_entropy(A, B)) atol=atol rtol=rtol
        elseif mode == 2
            @test p.optval ≈ evaluate(quantum_relative_entropy(X, B)) atol=atol rtol=rtol
        elseif mode == 3
            @test p.optval ≈ evaluate(quantum_relative_entropy(B, A)) atol=atol rtol=rtol
        elseif mode == 4
            @test p.optval ≈ evaluate(quantum_relative_entropy(B, X)) atol=atol rtol=rtol
        elseif mode == 5
            @test p.optval ≈ evaluate(quantum_relative_entropy(X, X)) atol=atol rtol=rtol
        end
    end
end

@add_problem sdp function sdp_quantum_relative_entropy1_fullrank(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    sdp_quantum_relative_entropy_impl(handle_problem!, Val(test), atol, rtol, T, false, 1)
end

@add_problem sdp function sdp_quantum_relative_entropy2_fullrank(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    sdp_quantum_relative_entropy_impl(handle_problem!, Val(test), atol, rtol, T, false, 2)
end

@add_problem sdp function sdp_quantum_relative_entropy3_fullrank(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    sdp_quantum_relative_entropy_impl(handle_problem!, Val(test), atol, rtol, T, false, 3)
end

@add_problem sdp function sdp_quantum_relative_entropy4_fullrank(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    sdp_quantum_relative_entropy_impl(handle_problem!, Val(test), atol, rtol, T, false, 4)
end

@add_problem sdp function sdp_quantum_relative_entropy5_fullrank(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    sdp_quantum_relative_entropy_impl(handle_problem!, Val(test), atol, rtol, T, false, 5)
end

@add_problem sdp function sdp_quantum_relative_entropy1_lowrank(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    sdp_quantum_relative_entropy_impl(handle_problem!, Val(test), atol, rtol, T, true, 1)
end

@add_problem sdp function sdp_quantum_relative_entropy2_lowrank(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    sdp_quantum_relative_entropy_impl(handle_problem!, Val(test), atol, rtol, T, true, 2)
end

@add_problem sdp function sdp_quantum_relative_entropy3_lowrank(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    sdp_quantum_relative_entropy_impl(handle_problem!, Val(test), atol, rtol, T, true, 3)
end

@add_problem sdp function sdp_quantum_relative_entropy4_lowrank(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    sdp_quantum_relative_entropy_impl(handle_problem!, Val(test), atol, rtol, T, true, 4)
end

@add_problem sdp function sdp_quantum_relative_entropy5_lowrank(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    sdp_quantum_relative_entropy_impl(handle_problem!, Val(test), atol, rtol, T, true, 5)
end

@add_problem sdp function sdp_quantum_relative_entropy_const(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    if test
        @test quantum_relative_entropy(diagm(0=>[0.5, 0.5]), diagm(0=>[0.5, 0.5])) ≈ 0 atol=atol rtol=rtol
        @test quantum_relative_entropy(diagm(0=>[1.0, 0.0]), diagm(0=>[1.0, 0.0])) ≈ 0 atol=atol rtol=rtol
        @test quantum_relative_entropy(diagm(0=>[1.0, 0.0]), diagm(0=>[0.5, 0.5])) ≈ log(2) atol=atol rtol=rtol
        @test isinf(quantum_relative_entropy(diagm(0=>[0.5, 0.5]), diagm(0=>[1.0, 0.0])))
    end
end

@add_problem sdp function sdp_quantum_channel_capacity(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    # Capacity of classical-quantum channel.  See chapter 3 of arxiv:1705.06671,
    # Fawzi, Fawzi "Efficient optimization of the quantum relative entropy".
    # For channels with two symbols, both pure states, we have an exact formula to compare against.

    n = 3

    for cplx in [false, true]
        if cplx
            ψ1 = normalize(randn(ComplexF64, n))
            ψ2 = normalize(randn(ComplexF64, n))
        else
            ψ1 = normalize(randn(Float64, n))
            ψ2 = normalize(randn(Float64, n))
        end
        ρ1 = ψ1 * ψ1'
        ρ2 = ψ2 * ψ2'

        p1 = Variable()
        p2 = Variable()

        objective = (quantum_entropy(p1*ρ1 + p2*ρ2) - p1*quantum_entropy(ρ1) - p2*quantum_entropy(ρ2)) / log(2)
        constraints = [ p1 >= 0, p2 >= 0, p1+p2 == 1 ]
        p = maximize(objective, constraints; numeric_type = T)
        handle_problem!(p)

        ϵ = abs(ψ1' * ψ2)
        q = (1 + ϵ)/2
        r = [q, 1-q]
        v = -r' * log2.(r)

        if test
            @test p.optval ≈ v atol=atol rtol=rtol
            @test evaluate(objective) ≈ v atol=atol rtol=rtol

            @test_throws DimensionMismatch quantum_entropy(Variable(2, 3))
        end
    end
end

@add_problem sdp function sdp_trace_mpower_argcheck(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    if test
        @test_throws DimensionMismatch trace_mpower(   zeros(2,3), 1//2, zeros(2,3))
        @test_throws DimensionMismatch trace_mpower(Variable(2,3), 1//2, zeros(2,3))

        @test_throws DimensionMismatch trace_mpower(   zeros(2,2), 1//2, zeros(3,3))
        @test_throws DimensionMismatch trace_mpower(Variable(2,2), 1//2, zeros(3,3))

        @test_throws DomainError trace_mpower(Variable(3,3),  5//2, zeros(3,3))

        nh = [ 1 1 ; 0 1 ] # not hermitian
        np = [ 1 2 ; 2 1 ] # not positive semidefinite
        @test_throws DomainError trace_mpower(Variable(2,2), 1//2, nh)
        @test_throws DomainError trace_mpower(Variable(2,2), 1//2, np)
    end
end

@add_problem sdp function sdp_trace_mpower_real_2_3(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    # Too slow for n>3.  Anything that can be done to speed it up?
    n = 3
    A = randn(n,n)
    A = A * A' # now A is positive semidefinite
    C = randn(n,n)
    C = C * C' # now C is positive semidefinite
    t = 2//3

    B = Semidefinite(n,n)

    constraints = [
        B ⪯ A
    ]
    objective = trace_mpower(B, t, C)
    p = maximize(objective, constraints; numeric_type = T)

    handle_problem!(p)

    if test
        @test real.(B.value) ≈ real.(A) atol=atol rtol=rtol
        @test imag.(B.value) ≈ imag.(A) atol=atol rtol=rtol
        @test p.optval ≈ tr(C*A^t) atol=atol rtol=rtol
        @test p.optval ≈ trace_mpower(A, t, C) atol=atol rtol=rtol
        @test p.optval ≈ evaluate(trace_mpower(B, t, C)) atol=atol rtol=rtol
    end
end

@add_problem sdp function sdp_trace_mpower_cplx_2_3(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    # Too slow for n>3.  Anything that can be done to speed it up?
    n = 3
    A = randn(ComplexF64, n, n)
    A = A * A' # now A is positive semidefinite
    C = randn(ComplexF64, n, n)
    C = C * C' # now C is positive semidefinite
    t = 2//3

    B = HermitianSemidefinite(n,n)

    constraints = [
        B ⪯ A
    ]
    objective = trace_mpower(B, t, C)
    p = maximize(objective, constraints; numeric_type = T)

    handle_problem!(p)

    if test
        @test real.(B.value) ≈ real.(A) atol=atol rtol=rtol
        @test imag.(B.value) ≈ imag.(A) atol=atol rtol=rtol
        @test p.optval ≈ tr(C*A^t) atol=atol rtol=rtol
        @test p.optval ≈ trace_mpower(A, t, C) atol=atol rtol=rtol
        @test p.optval ≈ evaluate(trace_mpower(B, t, C)) atol=atol rtol=rtol
    end
end

@add_problem sdp function sdp_trace_mpower_real_5_4(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    # Too slow for n>3.  Anything that can be done to speed it up?
    n = 3
    A = randn(n,n)
    A = A * A' # now A is positive semidefinite
    C = randn(n,n)
    C = C * C' # now C is positive semidefinite
    t = 5//4

    B = Semidefinite(n,n)

    constraints = [
        # note: B ⪰ A doesn't work for this unit test because trace_mpower is not operator monotone
        # unless t ∈ [0,1].
        B == A
    ]
    objective = trace_mpower(B, t, C)
    p = minimize(objective, constraints; numeric_type = T)

    handle_problem!(p)

    if test
        @test real.(B.value) ≈ real.(A) atol=atol rtol=rtol
        @test imag.(B.value) ≈ imag.(A) atol=atol rtol=rtol
        @test p.optval ≈ tr(C*A^t) atol=atol rtol=rtol
        @test p.optval ≈ trace_mpower(A, t, C) atol=atol rtol=rtol
        @test p.optval ≈ evaluate(trace_mpower(B, t, C)) atol=atol rtol=rtol
    end
end

@add_problem sdp function sdp_trace_mpower_cplx_5_4(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    # Too slow for n>3.  Anything that can be done to speed it up?
    n = 3
    A = randn(ComplexF64, n, n)
    A = A * A' # now A is positive semidefinite
    C = randn(ComplexF64, n, n)
    C = C * C' # now C is positive semidefinite
    t = 5//4

    B = HermitianSemidefinite(n,n)

    constraints = [
        # note: B ⪰ A doesn't work for this unit test because trace_mpower is not operator monotone
        # unless t ∈ [0,1].
        B == A
    ]
    objective = trace_mpower(B, t, C)
    p = minimize(objective, constraints; numeric_type = T)

    handle_problem!(p)

    if test
        @test real.(B.value) ≈ real.(A) atol=atol rtol=rtol
        @test imag.(B.value) ≈ imag.(A) atol=atol rtol=rtol
        @test p.optval ≈ tr(C*A^t) atol=atol rtol=rtol
        @test p.optval ≈ trace_mpower(A, t, C) atol=atol rtol=rtol
        @test p.optval ≈ evaluate(trace_mpower(B, t, C)) atol=atol rtol=rtol
    end
end

@add_problem sdp function sdp_trace_mpower_real_neg1_4(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    # Too slow for n>3.  Anything that can be done to speed it up?
    n = 3
    A = randn(n,n)
    A = A * A' # now A is positive semidefinite
    A += 0.2 * I # prevent numerical instability
    C = randn(n,n)
    C = C * C' # now C is positive semidefinite
    t = -1//4

    B = Semidefinite(n,n)

    constraints = [
        # -trace_mpower is operator monotone for t ∈ [-1, 0].
        B ⪯ A
    ]
    objective = trace_mpower(B, t, C)
    p = minimize(objective, constraints; numeric_type = T)

    handle_problem!(p)

    if test
        @test real.(B.value) ≈ real.(A) atol=atol rtol=rtol
        @test imag.(B.value) ≈ imag.(A) atol=atol rtol=rtol
        @test p.optval ≈ tr(C*A^t) atol=atol rtol=rtol
        @test p.optval ≈ trace_mpower(A, t, C) atol=atol rtol=rtol
        @test p.optval ≈ evaluate(trace_mpower(B, t, C)) atol=atol rtol=rtol
    end
end

@add_problem sdp function sdp_trace_mpower_cplx_neg1_4(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    # Too slow for n>3.  Anything that can be done to speed it up?
    n = 3
    A = randn(ComplexF64, n, n)
    A = A * A' # now A is positive semidefinite
    A += 0.2 * I # prevent numerical instability
    C = randn(ComplexF64, n, n)
    C = C * C' # now C is positive semidefinite
    t = -1//4

    B = HermitianSemidefinite(n,n)

    constraints = [
        # -trace_mpower is operator monotone for t ∈ [-1, 0].
        B ⪯ A
    ]
    objective = trace_mpower(B, t, C)
    p = minimize(objective, constraints; numeric_type = T)

    handle_problem!(p)

    if test
        @test real.(B.value) ≈ real.(A) atol=atol rtol=rtol
        @test imag.(B.value) ≈ imag.(A) atol=atol rtol=rtol
        @test p.optval ≈ tr(C*A^t) atol=atol rtol=rtol
        @test p.optval ≈ trace_mpower(A, t, C) atol=atol rtol=rtol
        @test p.optval ≈ evaluate(trace_mpower(B, t, C)) atol=atol rtol=rtol
    end
end

@add_problem sdp function sdp_trace_logm_argcheck(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    if test
        @test_throws DimensionMismatch trace_logm(   zeros(2,3), zeros(2,3))
        @test_throws DimensionMismatch trace_logm(Variable(2,3), zeros(2,3))

        @test_throws DimensionMismatch trace_logm(   zeros(2,2), zeros(3,3))
        @test_throws DimensionMismatch trace_logm(Variable(2,2), zeros(3,3))

        nh = [ 1 1 ; 0 1 ] # not hermitian
        np = [ 1 2 ; 2 1 ] # not positive semidefinite
        @test_throws DomainError trace_logm(Variable(2,2), nh)
        @test_throws DomainError trace_logm(Variable(2,2), np)
    end
end

@add_problem sdp function sdp_trace_logm_real(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    n = 3
    C = randn(Float64, n, n)
    C = C * C' # now C is positive semidefinite

    X = Variable(n,n)

    constraints = [
        X ⪯ eye(n)
    ]
    objective = trace_logm(X, C)
    p = maximize(objective, constraints; numeric_type = T)

    handle_problem!(p)

    if test
        @test real.(X.value) ≈ real.(eye(n)) atol=atol rtol=rtol
        @test imag.(X.value) ≈ imag.(eye(n)) atol=atol rtol=rtol
        @test p.optval ≈ tr(C*log(evaluate(X))) atol=atol rtol=rtol
        @test p.optval ≈ trace_logm(evaluate(X), C) atol=atol rtol=rtol
        @test p.optval ≈ evaluate(trace_logm(X, C)) atol=atol rtol=rtol
    end
end

@add_problem sdp function sdp_trace_logm_cplx(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    n = 3
    C = randn(ComplexF64, n, n)
    C = C * C' # now C is positive semidefinite

    X = ComplexVariable(n,n)

    constraints = [
        X ⪯ eye(n)
    ]
    objective = trace_logm(X, C)
    p = maximize(objective, constraints; numeric_type = T)

    handle_problem!(p)

    if test
        @test real.(X.value) ≈ real.(eye(n)) atol=atol rtol=rtol
        @test imag.(X.value) ≈ imag.(eye(n)) atol=atol rtol=rtol
        @test p.optval ≈ tr(C*log(evaluate(X))) atol=atol rtol=rtol
        @test p.optval ≈ trace_logm(evaluate(X), C) atol=atol rtol=rtol
        @test p.optval ≈ evaluate(trace_logm(X, C)) atol=atol rtol=rtol
    end
end

@add_problem sdp function sdp_lieb_ando(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    if test
        @test_throws DomainError lieb_ando(   zeros(3,3),    zeros(3,3), eye(3), -3//2)
        @test_throws DomainError lieb_ando(   zeros(3,3),    zeros(3,3), eye(3), 5//2)
        @test_throws DomainError lieb_ando(   zeros(3,3), Variable(3,3), eye(3), -3//2)
        @test_throws DomainError lieb_ando(   zeros(3,3), Variable(3,3), eye(3), 5//2)
        @test_throws DomainError lieb_ando(Variable(3,3),    zeros(3,3), eye(3), -3//2)
        @test_throws DomainError lieb_ando(Variable(3,3),    zeros(3,3), eye(3), 5//2)
        @test_throws DomainError lieb_ando(Variable(3,3), Variable(3,3), eye(3), -3//2)
        @test_throws DomainError lieb_ando(Variable(3,3), Variable(3,3), eye(3), 5//2)
    end

    for n in [2, 3]
        for t in [1//2, 1//4, 3//4, 1//8, 3//2, 5//4]
            for cplx in [false, true]
                #@show n,t,cplx

                if cplx
                    A = randn(ComplexF64, n, n)
                    B = randn(ComplexF64, n, n)
                    X = ComplexVariable(n, n)
                    Y = ComplexVariable(n, n)
                else
                    A = randn(n, n)
                    B = randn(n, n)
                    X = Variable(n, n)
                    Y = Variable(n, n)
                end
                A = A*A'
                B = B*B'
                A = A/tr(A)
                B = B/tr(B)

                QtAB = lieb_ando(A,B,eye(n),t)

                objective = lieb_ando(X,B,eye(n),t)
                if t >= 0 && t <= 1
                    p = maximize(objective, [X == A]; numeric_type = T)
                else
                    p = minimize(objective, [X == A]; numeric_type = T)
                end
                handle_problem!(p)
                if test
                    @test p.optval ≈ QtAB atol=atol*5 rtol=rtol
                end

                objective = lieb_ando(A,Y,eye(n),t)
                if t >= 0 && t <= 1
                    p = maximize(objective, [Y == B]; numeric_type = T)
                else
                    p = minimize(objective, [Y == B]; numeric_type = T)
                end
                handle_problem!(p)
                if test
                    @test p.optval ≈ QtAB atol=atol*5 rtol=rtol
                end

                objective = lieb_ando(X,Y,eye(n),t)
                if t >= 0 && t <= 1
                    p = maximize(objective, [X == A, Y == B]; numeric_type = T)
                else
                    p = minimize(objective, [X == A, Y == B]; numeric_type = T)
                end
                handle_problem!(p)
                if test
                    @test p.optval ≈ QtAB atol=atol*5 rtol=rtol
                end
            end
        end
    end
end


text:
"""
 Image filling Algorithm

 ```
 filled_img = imfill(img, interval)
 filled_img = imfill(img, interval, value)
 filled_img = imfill(img, interval, value, connectivity)
 ```

 Connected components of an image is found using flood-fill algorithm and returns a copy of
 the original image after filling objects that falls in the range of interval.
 For filling objects, represent the holes(part to be filled) with `true` in your array.

 Parameters:

  -  img            = Input image (Boolean array type)
  -  interval       = objects of size in this range will be filled with `false`
  -  connectivity   = connectivity takes the same values as in label_components (Default value is 1:ndims(img))

 """

function imfill(img::AbstractArray{Bool}, interval::Tuple{Real,Real}, connectivity::Union{Dims, AbstractVector{Int}, BitArray}=1:ndims(img))

    if interval[1] > interval[2] || interval[1] < 0 || interval[2] < 0
        throw(DomainError(interval,"Interval must be non-negative and in format (min_range,max_range)"))
    end

    labels = label_components(img,connectivity)

    count_labels = Dict([(i,count(x->x==i,labels)) for i in unique(labels)])

    new_img = similar(img)
    for ind in eachindex(img)
        if img[ind] == true && interval[1] <= count_labels[labels[ind]] <= interval[2]
            new_img[ind] = false
        else
            new_img[ind] = img[ind]
        end
    end

    return new_img
 end


text:
"""
# module Insort

# Examples

```jldoctest
julia>
```
"""
module Insort

using BisectPy: bisect_left, bisect_right

export insort_left, insort_right

"""
    insort_left(a, x, lo=1, hi=nothing)

Insert item `x` in array `a`, and keep it sorted assuming `a` is sorted.

If `x` is already in `a`, insert it to the left of the leftmost `x`.
Optional args `lo` (default `1`) and `hi` (default `length(a)`) bound the
slice of `a` to be searched.
"""
function insort_left(a, x, lo=1, hi=nothing)
    lo = bisect_left(a, x, lo, hi)
    insert!(a, lo, x)
end

"""
    insort_right(a, x, lo=1, hi=nothing)

Insert item `x` in array `a`, and keep it sorted assuming `a` is sorted.

If `x` is already in `a`, insert it to the right of the rightmost `x`.
Optional args `lo` (default `1`) and `hi` (default `length(a)`) bound the
slice of `a` to be searched.
"""
function insort_right(a, x, lo=1, hi=nothing)
    lo = bisect_right(a, x, lo, hi)
    insert!(a, lo, x)
end

end


text:
module Day12
export day12

parseAngle(arg) = div(arg, 90)

# char => (:name, factor, parsefunction)
const COMMANDS = Dict('R' => (:rotate, -1, parseAngle),
                      'L' => (:rotate,  1, parseAngle),
                      'F' => (:forward, 1, identity),
                      'N' => (:go, im, identity),
                      'S' => (:go, -im, identity),
                      'W' => (:go, -1, identity),
                      'E' => (:go, 1, identity))

function parseLine(line)
    command = line[1]
    arg = parse(Int, line[2:end])
    name, factor, func = COMMANDS[command]
    name, factor*func(arg)
end

function readInput(inputfile="inputs/input12.txt")
    [parseLine(line) for line in eachline(inputfile)]
end

function day12(inputfile="inputs/input12.txt")
    instructions = readInput(inputfile)
    (day12Part1(instructions), day12Part2(instructions))
end

function day12Part1(instructions)
    at = complex(0)
    facing = complex(1)
    for (command, value) in instructions
        if command === :go
            at += value
        elseif command === :rotate
            facing *= (sign(value)*im)^abs(value)
        elseif command === :forward
            at += facing*value
        else
            error("Unknown instruction: $command")
        end
    end
    abs(real(at)) + abs(imag(at))
end

function day12Part2(instructions)
    waypoint = 10+1im
    at = complex(0)
    for (command, value) in instructions
        if command === :go
            waypoint += value
        elseif command === :rotate
            waypoint *= (sign(value)*im)^abs(value)
        elseif command === :forward
            at += waypoint*value
        else
            error("Unknown instruction: $command")
        end
    end
    abs(real(at)) + abs(imag(at))
end

end #module


text:
const psxy  = plot
const psxy! = plot!
const psxyz  = plot3d
const psxyz! = plot3d!

# ---------------------------------------------------------------------------------------------------
function common_plot_xyz(cmd0, arg1, caller, first, is3D, kwargs...)
	global multi_col
	N_args = (arg1 === nothing) ? 0 : 1

	is_ternary = (caller == "ternary") ? true : false
	if (is3D)	        gmt_proggy = "psxyz "
	elseif (is_ternary) gmt_proggy = "psternary "
	else		        gmt_proggy = "psxy "
	end

	(occursin(" -", cmd0)) && return monolitic(gmt_proggy, cmd0, arg1)

	cmd = "";	sub_module = ""			# Will change to "scatter", etc... if called by sub-modules
	if (caller != "")
		if (occursin(" -", caller))		# some sub-modues use this piggy-backed call
			if ((ind = findfirst("|", caller)) !== nothing)	# A mixed case with "caler|partiall_command"
				sub_module = caller[1:ind[1]-1]
				cmd = caller[ind[1]+1:end]
				caller = sub_module		# Because of parse_BJR()
			else
				cmd = caller
				caller = "others"		# It was piggy-backed
			end
		else
			sub_module = caller
		end
	end

	d = KW(kwargs)
    K, O = set_KO(first)		# Set the K O dance

	cmd, opt_B, opt_J, opt_R = parse_BJR(d, cmd, caller, O)
	if (is3D)	cmd,opt_JZ = parse_JZ(cmd, d)	end
	cmd = parse_common_opts(d, cmd, [:a :e :f :g :p :t :yx :params], first)
	cmd = parse_these_opts(cmd, d, [[:D :shift :offset], [:I :intens], [:N :no_clip :noclip]])
	if (is_ternary)  cmd = add_opt(cmd, 'M', d, [:M :no_plot])  end
	opt_UVXY = parse_UVXY("", d)	# Need it separate to not risk to double include it.
	cmd, opt_c = parse_c(cmd, d)	# Need opt_c because we may need to remove it from double calls

	# If a file name sent in, read it and compute a tight -R if this was not provided
	if (opt_R == "" && sub_module == "bar")  opt_R = "///0"  end	# Make sure y_min = 0
	cmd, arg1, opt_R, opt_i = read_data(d, cmd0, cmd, arg1, opt_R, is3D)

	if (is3D && isempty(opt_JZ) && length(collect(eachmatch(r"/", opt_R))) == 5)
		cmd *= " -JZ6c"		# Default -JZ
	end

	cmd = add_opt(cmd, 'A', d, [:A :steps :straight_lines], (x="x", y="y", meridian="m", parallel="p"))
	opt_F = add_opt("", "", d, [:F :conn :connection],
	                (continuous=("c", nothing, 1), net=("n", nothing, 1), network=("n", nothing, 1), refpoint=("r", nothing, 1),  ignore_hdr="_a", single_group="_f", segments="_s", segments_reset="_r", anchor=("", arg2str)))
	if (length(opt_F) > 1 && !occursin("/", opt_F))  opt_F = opt_F[1]  end	# Allow con=:net or con=(1,2)
	if (opt_F != "")  cmd *= " -F" * opt_F  end

	# Error Bars?
	got_Ebars = false
	val, symb = find_in_dict(d, [:E :error :error_bars])
	if (val !== nothing)
		cmd, arg1 = add_opt(add_opt, (cmd, 'E', d, [symb]),
					        (x="|x",y="|y",xy="|xy",X="|X",Y="|Y", asym="_+a", colored="_+c", cline="_+cl", csymbol="_+cf", wiskers="|+n",cap="+w",pen=("+p",add_opt_pen)), false, arg1)
		got_Ebars = true
	end

	# Look for color request. Do it after error bars because they may add a column
	len = length(cmd);	n_prev = N_args;
	cmd, arg1, arg2, N_args = add_opt_cpt(d, cmd, [:C :color :cmap], 'C', N_args, arg1)
	#cmd, arg1, arg2, N_args = add_opt_cpt(d, cmd, [:C :color :cmap], 'C', N_args, arg1, nothing, true, true, "", true)

	# See if we got a CPT. If yes there may be some work to do if no color column provided in input data.
	cmd, arg1, arg2, N_args, mcc = make_color_column(d, cmd, opt_i, len, N_args, n_prev, is3D, got_Ebars, arg1, arg2)

	cmd = add_opt_fill(cmd, d, [:G :fill], 'G')
	opt_Gsymb = add_opt_fill("", d, [:G :markerfacecolor :MarkerFaceColor :mc], 'G')	# Filling of symbols

	# To track a still existing bug in sessions management at GMT lib level
	got_pattern = false
	if (occursin("-Gp", cmd) || occursin("-GP", cmd) || occursin("-Gp", opt_Gsymb) || occursin("-GP", opt_Gsymb))
		got_pattern = true
	end

	if (is_ternary)			# Means we are in the psternary mode
		cmd = add_opt(cmd, 'L', d, [:L :labels])
	else
		cmd = add_opt(cmd, 'L', d, [:L :close],
			(left="_+xl", right="_+xr", x0="+x", bot="_+yb", top="_+yt", y0="+y", sym="_+d", asym="_+D", envelope="_+b", pen=("+p",add_opt_pen)))
		if (occursin("-L", cmd) && !occursin("-G", cmd) && !occursin("+p", cmd))  cmd *= "+p0.5p"  end
	end

	opt_Wmarker = ""
	if ((val = find_in_dict(d, [:markeredgecolor :MarkerEdgeColor])[1]) !== nothing)
		opt_Wmarker = "0.5p," * arg2str(val)		# 0.25p is too thin?
	end

	opt_W = add_opt_pen(d, [:W :pen], "W")
	if (occursin("+c", opt_W) || occursin("+c", cmd) && !occursin("-C", cmd))
		@warn("Color lines (or fill) from a color scale was selected but no color scale provided. Expect ...")
	end

	opt_S = add_opt("", 'S', d, [:S :symbol], (symb="1", size="", unit="1"))
	if (opt_S == "")			# OK, no symbol given via the -S option. So fish in aliases
		marca, arg1, more_cols = get_marker_name(d, [:marker :Marker :shape], is3D, false, arg1)
		if ((val = find_in_dict(d, [:markersize :MarkerSize :ms :size])[1]) !== nothing)
			if (marca == "")  marca = "c"  end		# If a marker name was not selected, defaults to circle
			if (isa(val, AbstractArray))
				if (length(val) == size(arg1,1))
					arg1 = hcat(arg1, val[:])
				else
					error("The size array must have the same number of elements rows in the data")
				end
			else
				marca *= arg2str(val);
			end
			opt_S = " -S" * marca
		elseif (marca != "")		# User only selected a marker name but no size.
			opt_S = " -S" * marca
			# If data comes from a file, then no automatic symbol size is added
			op = lowercase(marca[1])
			if (!more_cols && arg1 !== nothing && !isa(arg1, GMTcpt) && !occursin(op, "bekmrvw"))  opt_S *= "7p"  end
		end
	end

	opt_ML = ""
	if (opt_S != "")
		if ((val = find_in_dict(d, [:markerline :MarkerLine :ml])[1]) !== nothing)
			if (isa(val, Tuple))  opt_ML = " -W" * parse_pen(val) # This can hold the pen, not extended atts
			elseif (isa(val, NamedTuple))  opt_ML = add_opt_pen(nt2dict(val), [:pen], "W")
			else                  opt_ML = " -W" * arg2str(val)
			end
			if (opt_Wmarker != "")
				opt_Wmarker = ""
				@warn("markerline overrides markeredgecolor")
			end
		end
		if (opt_W != "" && opt_ML != "")
			@warn("You cannot use both markerline and W or pen keys.")
		end
	end

	# See if any of the scatter, bar, lines, etc... was the caller and if yes, set sensible defaults.
	cmd = check_caller(d, cmd, opt_S, opt_W, sub_module, O)

	if (opt_W != "" && opt_S == "") 						# We have a line/polygon request
		#cmd = finish_PS(d, cmd * opt_W * opt_UVXY, output, K, O)
		cmd *= opt_W * opt_UVXY

	elseif (opt_W == "" && opt_S != "")						# We have a symbol request
		if (opt_Wmarker != "" && opt_W == "") opt_Gsymb *= " -W" * opt_Wmarker  end		# reuse var name
		if (opt_ML != "")  cmd *= opt_ML  end				# If we have a symbol outline pen
		#cmd = finish_PS(d, cmd * opt_S * opt_Gsymb * opt_UVXY, output, K, O)
		cmd *= opt_S * opt_Gsymb * opt_UVXY

	elseif (opt_W != "" && opt_S != "")						# We have both line/polygon and a symbol
		if (occursin(opt_Gsymb, cmd))  opt_Gsymb = ""  end
		if (opt_S[4] == 'v' || opt_S[4] == 'V' || opt_S[4] == '=')
			#cmd = finish_PS(d, cmd * opt_W * opt_S * opt_Gsymb * opt_UVXY, output, K, O)
			cmd *= opt_W * opt_S * opt_Gsymb * opt_UVXY
		else
			if (opt_Wmarker != "")  opt_Wmarker = " -W" * opt_Wmarker  end		# Set Symbol edge color
			cmd1 = cmd * opt_W * opt_UVXY
			cmd2 = replace(cmd, opt_B => "") * opt_S * opt_Gsymb * opt_Wmarker	# Don't repeat option -B
			if (opt_c != "")  cmd2 = replace(cmd2, opt_c => "")  end			# Not in scond call (subplots)
			if (opt_ML != "")  cmd1 = cmd1 * opt_ML  end	# If we have a symbol outline pen
			#cmd = [finish_PS(d, cmd1, "", true, O); finish_PS(d, cmd2, "", K, true)]
			cmd = [cmd1; cmd2]
		end

	else
		#cmd = finish_PS(d, cmd * opt_UVXY, output, K, O)
		cmd *= opt_UVXY
	end

	# Let matrices with more data columns, and for which Color info was NOT set, plot multiple lines at once
	if (!mcc && opt_S == "" && (caller == "lines" || caller == "plot") && isa(arg1, Array{<:Number,2}) && size(arg1,2) > 2+is3D && size(arg1,1) > 1 && (multi_col || haskey(d, :multicol)) )
		multi_col = false							# Reset because this is a use-only-once option
		penC = "";		penS = "";	cycle=:cycle
		# But if we have a color in opt_W (idiotic) let it overrule the automatic color cycle in mat2ds()
		if (opt_W != "")  penT, penC, penS = break_pen(scan_opt(opt_W, "-W"))  end
		if (penC  != "")  cycle = [penC]  end
		arg1 = mat2ds(arg1, color=cycle, ls=penS)	# Convert to multi-segment GMTdataset
		D = gmt("gmtinfo -C", arg1)					# But now also need to update the -R string
		if (isa(cmd, Array))						# Replace old -R by the new one
			cmd[1] = replace(cmd[1], opt_R => " -R" * arg2str(round_wesn(D[1].data)))
		else
			cmd = replace(cmd, opt_R => " -R" * arg2str(round_wesn(D[1].data)))
		end
	end

	put_in_legend_bag(d, cmd, arg1)

	r = finish_PS_module(d, gmt_proggy .* cmd, "", K, O, true, arg1, arg2)
	if (got_pattern || occursin("-Sk", opt_S))  gmt("destroy")  end 	# Apparently patterns are screweing the session
	return r
end

# ---------------------------------------------------------------------------------------------------
function make_color_column(d, cmd, opt_i, len, N_args, n_prev, is3D, got_Ebars, arg1, arg2)
	# See if we got a CPT. If yes, there is quite some work to do if no color column provided in input data.
	# N_ARGS will be == n_prev+1 when a -Ccpt was used. Otherwise they are equal.

	if (arg1 === nothing || isa(arg1, GMT.GMTcpt))  return cmd, arg1, arg2, N_args, false  end		# Play safe

	mz, the_kw = find_in_dict(d, [:zcolor :markerz :mz])
	if (!(N_args > n_prev || len < length(cmd)) && mz === nothing)	# No color request, so return right away
		return cmd, arg1, arg2, N_args, false
	end

	if (isa(arg1, Array{<:Number}))  n_rows, n_col = size(arg1)
	elseif (isa(arg1,GMTdataset))    n_rows, n_col = size(arg1.data)
	else                             n_rows, n_col = size(arg1[1].data)
	end

	warn1 = string("Probably color column in ", the_kw, " has incorrect dims. Ignoring it.")
	warn2 = "Plotting with color table requires adding one more column to the dataset but your -i
	option didn't do it, so you won't get what you expect. Try -i0-1,1 for 2D or -i0-2,2 for 3D plots"

	if (n_col <= 2+is3D)
		if (mz !== nothing)
			if (length(mz) != n_rows)  @warn(warn1); @goto noway  end
			if (isa(arg1, Array))          arg1 = hcat(arg1, mz[:])
			elseif (isa(arg1,GMTdataset))  arg1.data = hcat(arg1.data, mz[:])
			else                           arg1[1].data = hcat(arg1[1].data, mz[:])
			end
		else
			if (opt_i != "")  @warn(warn2);		@goto noway		end
			cmd = @sprintf("%s -i0-%d,%d", cmd, 1+is3D, 1+is3D)
			if ((val = find_in_dict(d, [:markersize :ms :size])[1]) !== nothing)
				cmd = @sprintf("%s-%d", cmd, 2+is3D)	# Because we know that an extra col will be added later
			end
		end
	else
		if (mz !== nothing)			# Here we must insert the color col right after the coords
			if (length(mz) != n_rows)  @warn(warn1); @goto noway  end
			if     (isa(arg1, Array))      arg1 = hcat(arg1[:,1:2+is3D], mz[:], arg1[:,3+is3D:end])
			elseif (isa(arg1,GMTdataset))  arg1.data = hcat(arg1.data[:,1:2+is3D], mz[:], arg1.data[:,3+is3D:end])
			else                           arg1[1].data = hcat(arg1[1].data[:,1:2+is3D], mz[:], arg1[1].data[:,3+is3D:end])
			end
		elseif (!got_Ebars)		# The Error bars case is very multi. Don't try to guess then.
			if (opt_i == "")  cmd = @sprintf("%s -i0-%d,%d,%d-%d", cmd, 1+is3D, 1+is3D, 2+is3D, n_col-1)
			else              @warn(warn2);		@goto noway
			end
		end
	end

	if (N_args == n_prev)		# No cpt transmitted, so need to compute one
		if (GMTver >= 7)		# 7 because this solution is currently still bugged
			#=
			if (mz !== nothing)
				arg2 = gmt("makecpt -E " * cmd[len+2:end], mz[:])
			else
				if (isa(arg1, Array))  arg2 = gmt("makecpt -E " * cmd[len+2:end], arg1)
				else                   arg2 = gmt("makecpt -E " * cmd[len+2:end], arg1.data)
				end
			end
			=#
		else
			if (mz !== nothing)               mi, ma = extrema(mz)
			else
				if     (isa(arg1, Array))     mi, ma = extrema(view(arg1, :, 2+is3D))
				elseif (isa(arg1,GMTdataset)) mi, ma = extrema(view(arg1.data, :, 2+is3D))
				else                          mi, ma = extrema(view(arg1[1].data, :, 2+is3D))
				end
			end
			just_C = cmd[len+2:end];	reset_i = ""
			if ((ind = findfirst(" -i", just_C)) !== nothing)
				reset_i = just_C[ind[1]:end]
				just_C  = just_C[1:ind[1]-1]
			end
			arg2 = gmt(string("makecpt -T", mi-0.001*abs(mi), '/', ma+0.001*abs(ma), " ", just_C))
			if (occursin(" -C", cmd))  cmd = cmd[1:len+3]  end		# Strip the cpt name
			if (reset_i != "")  cmd *= reset_i  end		# Reset -i, in case it existed
		end
		if (!occursin(" -C", cmd))  cmd *= " -C"  end	# Need to inform that there is a cpt to use
		N_args = 2
	end

	@label noway

	return cmd, arg1, arg2, N_args, true
end

# ---------------------------------------------------------------------------------------------------
function get_marker_name(d::Dict, symbs, is3D, del=false, arg1=nothing)
	marca = "";		N = 0
	for symb in symbs
		if (haskey(d, symb))
			t = d[symb]
			if (isa(t, Tuple))				# e.g. marker=(:r, [2 3])
				msg = "";	cst = false
				o = string(t[1])
				if     (startswith(o, "E"))  opt = "E";  N = 3; cst = true
				elseif (startswith(o, "e"))  opt = "e";  N = 3
				elseif (o == "J" || startswith(o, "Rot"))  opt = "J";  N = 3; cst = true
				elseif (o == "j" || startswith(o, "rot"))  opt = "j";  N = 3
				elseif (o == "M" || startswith(o, "Mat"))  opt = "M";  N = 3
				elseif (o == "m" || startswith(o, "mat"))  opt = "m";  N = 3
				elseif (o == "R" || startswith(o, "Rec"))  opt = "R";  N = 3
				elseif (o == "r" || startswith(o, "rec"))  opt = "r";  N = 2
				elseif (o == "V" || startswith(o, "Vec"))  opt = "V";  N = 2
				elseif (o == "v" || startswith(o, "vec"))  opt = "v";  N = 2
				elseif (o == "w" || o == "pie" || o == "web" || o == "wedge")  opt = "w";  N = 2
				elseif (o == "W" || o == "Pie" || o == "Web" || o == "Wedge")  opt = "W";  N = 2
				end
				if (N > 0)  marca, arg1, msg = helper_markers(opt, t[2], arg1, N, cst)  end
				if (msg != "")  error(msg)  end
				if (length(t) == 3 && isa(t[3], NamedTuple))
					if (marca == "w" || marca == "W")	# Ex (spiderweb): marker=(:pie, [...], (inner=1,))
						marca *= add_opt(t[3], (inner="/", arc="+a", radial="+r", size=("", arg2str, 1), pen=("+p", add_opt_pen)) )
					elseif (marca == "m" || marca == "M")
						marca *= vector_attrib(t[3])
					end
				end
			elseif (isa(t, NamedTuple))		# e.g. marker=(pie=true, inner=1, ...)
				key = keys(t)[1];	opt = ""
				if     (key == :w || key == :pie || key == :web || key == :wedge)  opt = "w"
				elseif (key == :W || key == :Pie || key == :Web || key == :Wedge)  opt = "W"
				elseif (key == :b || key == :bar)     opt = "b"
				elseif (key == :B || key == :HBar)    opt = "B"
				elseif (key == :l || key == :letter)  opt = "l"
				elseif (key == :K || key == :Custom)  opt = "K"
				elseif (key == :k || key == :custom)  opt = "k"
				elseif (key == :M || key == :Matang)  opt = "M"
				elseif (key == :m || key == :matang)  opt = "m"
				end
				if (opt == "w" || opt == "W")
					marca = opt * add_opt(t, (size=("", arg2str, 1), inner="/", arc="+a", radial="+r", pen=("+p", add_opt_pen)))
				elseif (opt == "b" || opt == "B")
					marca = opt * add_opt(t, (size=("", arg2str, 1), base="+b", Base="+B"))
				elseif (opt == "l")
					marca = opt * add_opt(t, (size=("", arg2str, 1), letter="+t", justify="+j", font=("+f", font)))
				elseif (opt == "m" || opt == "M")
					marca = opt * add_opt(t, (size=("", arg2str, 1), arrow=("", vector_attrib)))
				elseif (opt == "k" || opt == "K")
					marca = opt * add_opt(t, (custom="", size="/"))
				end
			else
				if (isa(t, Symbol))	t = string(t)	end
				if     (t == "-" || t == "x-dash")    marca = "-"
				elseif (t == "+" || t == "plus")      marca = "+"
				elseif (t == "a" || t == "*" || t == "star")  marca = "a"
				elseif (t == "k" || t == "custom")    marca = "k"
				elseif (t == "x" || t == "cross")     marca = "x"
				elseif (t[1] == 'c' || t[1] == 'C')   marca = "c"
				elseif (t == "d" || t == "diamond")   marca = "d"
				elseif (t == "g" || t == "octagon")   marca = "g"
				elseif (t == "h" || t == "hexagon")   marca = "h"
				elseif (t == "i" || t == "v" || t == "inverted_tri")  marca = "i"
				elseif (t == "l" || t == "letter")    marca = "l"
				elseif (t == "n" || t == "pentagon")  marca = "n"
				elseif (t == "p" || t == "." || t == "point")  marca = "p"
				elseif (t == "s" || t == "square")    marca = "s"
				elseif (t == "t" || t == "^" || t == "triangle")  marca = "t"
				elseif (t == "T" || t == "Triangle")  marca = "T"
				elseif (is3D && (t == "u" || t == "cube"))  marca = "u"
				elseif (t == "y" || t == "y-dash")    marca = "y"
				end
				# Still need to check the simpler forms of these
				if (marca == "")  marca = helper2_markers(t, ["e" "ellipse"])   end
				if (marca == "")  marca = helper2_markers(t, ["E" "Ellipse"])   end
				if (marca == "")  marca = helper2_markers(t, ["j" "rotrect"])   end
				if (marca == "")  marca = helper2_markers(t, ["J" "RotRect"])   end
				if (marca == "")  marca = helper2_markers(t, ["m" "matangle"])  end
				if (marca == "")  marca = helper2_markers(t, ["M" "Matangle"])  end
				if (marca == "")  marca = helper2_markers(t, ["r" "rectangle"])   end
				if (marca == "")  marca = helper2_markers(t, ["R" "RRectangle"])  end
				if (marca == "")  marca = helper2_markers(t, ["v" "vector"])  end
				if (marca == "")  marca = helper2_markers(t, ["V" "Vector"])  end
				if (marca == "")  marca = helper2_markers(t, ["w" "pie" "web"])  end
				if (marca == "")  marca = helper2_markers(t, ["W" "Pie" "Web"])  end
			end
			if (del)  delete!(d, symb)  end
			break
		end
	end
	return marca, arg1, N > 0
end

function helper_markers(opt, ext, arg1, N, cst)
	# Helper function to deal with the cases where one sends marker's extra columns via command
	# Example that will land and be processed here:  marker=(:Ellipse, [30 10 15])
	# N is the number of extra columns
	marca = "";		msg = ""
	if (size(ext,2) == N && arg1 !== nothing)
		S = Symbol(opt)
		marca, arg1 = add_opt(add_opt, (opt, "", Dict(S => (par=ext,)), [S]), (par="|",), true, arg1)
	elseif (cst && length(ext) == 1)
		marca = opt * "-" * string(ext)
	else
		msg = string("Wrong number of extra columns for marker (", opt, "). Got ", size(ext,2), " but expected ", N)
	end
	return marca, arg1, msg
end

function helper2_markers(opt, alias)
	marca = ""
	if (opt == alias[1])			# User used only the one letter syntax
		marca = alias[1]
	else
		for k = 2:length(alias)		# Loop because of cases like ["w" "pie" "web"]
			o2 = alias[k][1:min(2,length(alias[k]))]	# check the first 2 chars and Ro, Rotrect or RotRec are all good
			if (startswith(opt, o2))  marca = alias[1]; break  end		# Good when, for example, marker=:Pie
		end
	end

	# If we still have found nothing, assume that OPT is a full GMT opt string (e.g. W/5+a30+r45+p2,red)
	if (marca == "" && opt[1] == alias[1][1])  marca = opt  end
	return marca
end

# ---------------------------------------------------------------------------------------------------
function check_caller(d::Dict, cmd::String, opt_S::String, opt_W::String, caller::String, O::Bool)
	# Set sensible defaults for the sub-modules "scatter" & "bar"
	if (caller == "scatter")
		if (opt_S == "")  cmd *= " -Sc7p"  end
	elseif (caller == "scatter3")
		if (opt_S == "")  cmd *= " -Su7p"  end
	elseif (caller == "lines")
		if (!occursin("+p", cmd) && opt_W == "") cmd *= " -W0.25p"  end # Do not leave without a pen specification
	elseif (caller == "bar")
		if (opt_S == "")
			if (haskey(d, :bar))
				cmd = GMT.parse_bar_cmd(d, :bar, cmd, "Sb")
			elseif (haskey(d, :hbar))
				cmd = GMT.parse_bar_cmd(d, :hbar, cmd, "SB")
			else
				opt = add_opt("", "",  d, [:width])		# No need to purge because width is not a psxy option
				if (opt == "")	opt = "0.8"	end			# The default
				cmd *= " -Sb" * opt * "u"

				optB = add_opt("", "",  d, [:base])
				if (optB == "")	optB = "0"	end
				cmd *= "+b" * optB
			end
		end
		if (!occursin(" -G", cmd) && !occursin(" -C", cmd))  cmd *= " -G0/115/190"	end		# Default color
	elseif (caller == "bar3")
		if (haskey(d, :noshade) && occursin("-So", cmd))
			cmd = replace(cmd, "-So" => "-SO", count=1)
		end
		if (!occursin(" -G", cmd) && !occursin(" -C", cmd))  cmd *= " -G0/115/190"	end
		if (!occursin(" -J", cmd))  cmd *= " -JX12c/0"  end
	end

	if (occursin('3', caller))
		if (!occursin(" -B", cmd) && !O)  cmd *= def_fig_axes3  end	# For overlays default is no axes
		if (!occursin(" -p", cmd))  cmd *= " -p200/30"  end
	#else
		#if (!occursin(" -B", cmd) && !O)  cmd *= def_fig_axes   end
	end

	return cmd
end

# ---------------------------------------------------------------------------------------------------
function parse_bar_cmd(d::Dict, key::Symbol, cmd::String, optS::String, no_u=false)
	# Deal with parsing the 'bar' & 'hbar' keywors of psxy. Also called by plot/bar3. For this
	# later module if input is not a string or NamedTuple the scatter options must be processed in bar3().
	# KEY is either :bar or :hbar
	# OPTS is either "Sb", "SB" or "So"
	# NO_U if true means to NO automatic adding of flag 'u'
	opt =""
	if (haskey(d, key))
		if (isa(d[key], String))
			cmd *= " -" * optS * d[key]
		elseif (isa(d[key], NamedTuple))
			opt = add_opt("", optS, d, [key], (width="",unit="1",base="+b",height="+B",nbands="+z",Nbands="+Z"))
		else
			error("Argument of the *bar* keyword can be only a string or a NamedTuple.")
		end
	end

	if (opt != "")				# Still need to finish parsing this
		flag_u = no_u ? "" : 'u'
		if ((ind = findfirst("+", opt)) !== nothing)	# See if need to insert a 'u'
			if (!isletter(opt[ind[1]-1]))  opt = opt[1:ind[1]-1] * flag_u * opt[ind[1]:end]  end
		else
			pb = (optS != "So") ? "+b0" : ""		# The default for bar3 (So) is set in the bar3() fun
			if (!isletter(opt[end]))  opt *= flag_u	  end	# No base set so default to ...
			opt *= pb
		end
		cmd *= opt
	end
	return cmd
end


text:
using MathProgBase

include("mpb_model.jl")

export MathProgNLSModel

mutable struct MathProgNLSModel <: AbstractNLSModel
  meta :: NLPModelMeta
  nls_meta :: NLSMeta
  Fmodel :: MathProgModel
  cmodel :: MathProgModel
  counters :: NLSCounters      # Evaluation counters.

  Fjrows :: Vector{Int}      # Jacobian sparsity pattern.
  Fjcols :: Vector{Int}
  Fjvals :: Vector{Float64}  # Room for the constraints Jacobian.

  Fhrows :: Vector{Int}      # Hessian sparsity pattern.
  Fhcols :: Vector{Int}
  Fhvals :: Vector{Float64}  # Room for the Lagrangian Hessian.

  cjrows :: Vector{Int}      # Jacobian sparsity pattern.
  cjcols :: Vector{Int}
  cjvals :: Vector{Float64}  # Room for the constraints Jacobian.

  chrows :: Vector{Int}      # Hessian sparsity pattern.
  chcols :: Vector{Int}
  chvals :: Vector{Float64}  # Room for the Lagrangian Hessian.
end

"Construct a `MathProgNLSModel` from two `MathProgModel`s."
function MathProgNLSModel(Fmodel :: MathProgModel,
                          cmodel :: MathProgModel;
                          name :: String="Generic")

  nvar = cmodel.numVar
  nequ = Fmodel.numConstr
  lvar = cmodel.lvar
  uvar = cmodel.uvar

  nlin = length(cmodel.eval.m.linconstr)         # Number of linear constraints.
  nquad = length(cmodel.eval.m.quadconstr)       # Number of quadratic constraints.
  nnln = length(cmodel.eval.m.nlpdata.nlconstr)  # Number of nonlinear constraints.
  ncon = cmodel.numConstr                        # Total number of constraints.
  lcon = cmodel.lcon
  ucon = cmodel.ucon

  Fjrows, Fjcols = MathProgBase.jac_structure(Fmodel.eval)
  cjrows, cjcols = MathProgBase.jac_structure(cmodel.eval)
  Fhrows, Fhcols = MathProgBase.hesslag_structure(Fmodel.eval)
  chrows, chcols = MathProgBase.hesslag_structure(cmodel.eval)

  meta = NLPModelMeta(nvar,
                      x0=cmodel.x,
                      lvar=lvar,
                      uvar=uvar,
                      ncon=ncon,
                      y0=zeros(ncon),
                      lcon=lcon,
                      ucon=ucon,
                      nnzj=length(cjrows),
                      nnzh=length(chrows),
                      lin=collect(1:nlin),  # linear constraints appear first in MPB
                      nln=collect(nlin+1:ncon),
                      name=name,
                      )

  return MathProgNLSModel(meta,
                          NLSMeta(nequ, nvar, nnzj=length(Fjrows), nnzh=length(Fhrows)),
                          Fmodel,
                          cmodel,
                          NLSCounters(),
                          Fjrows,
                          Fjcols,
                          zeros(length(Fjrows)),  # Fjvals
                          Fhrows,
                          Fhcols,
                          zeros(length(Fhrows)),  # Fhvals
                          cjrows,
                          cjcols,
                          zeros(length(cjrows)),  # cjvals
                          chrows,
                          chcols,
                          zeros(length(chrows)),  # chvals
                         )
end

function NLPModels.residual!(nls :: MathProgNLSModel, x :: AbstractVector, Fx :: AbstractVector)
  increment!(nls, :neval_residual)
  MathProgBase.eval_g(nls.Fmodel.eval, Fx, x)
  return Fx
end

function NLPModels.jac_residual(nls :: MathProgNLSModel, x :: AbstractVector)
  increment!(nls, :neval_jac_residual)
  m, n = nls.nls_meta.nequ, nls.meta.nvar
  MathProgBase.eval_jac_g(nls.Fmodel.eval, nls.Fjvals, x)
  return sparse(nls.Fjrows, nls.Fjcols, nls.Fjvals, m, n)
end

function NLPModels.jac_coord_residual(nls :: MathProgNLSModel, x :: AbstractVector)
  increment!(nls, :neval_jac_residual)
  m, n = nls.nls_meta.nequ, nls.meta.nvar
  MathProgBase.eval_jac_g(nls.Fmodel.eval, nls.Fjvals, x)
  return nls.Fjrows, nls.Fjcols, nls.Fjvals
end

function NLPModels.jprod_residual!(nls :: MathProgNLSModel, x :: AbstractVector, v :: AbstractVector, Jv :: AbstractVector)
  nls.counters.neval_jac_residual -= 1
  increment!(nls, :neval_jprod_residual)
  Jv[:] = jac_residual(nls, x) * v
  return Jv
end

function NLPModels.jtprod_residual!(nls :: MathProgNLSModel, x :: AbstractVector, v :: AbstractVector, Jtv :: AbstractVector)
  nls.counters.neval_jac_residual -= 1
  increment!(nls, :neval_jtprod_residual)
  Jtv[:] = jac_residual(nls, x)' * v
  return Jtv
end

function NLPModels.hess_residual(nls :: MathProgNLSModel, x :: AbstractVector, v :: AbstractVector)
  increment!(nls, :neval_hess_residual)
  n = nls.meta.nvar
  MathProgBase.eval_hesslag(nls.Fmodel.eval, nls.Fhvals, x, 0.0, v)
  return sparse(nls.Fhrows, nls.Fhcols, nls.Fhvals, n, n)
end

function NLPModels.hess_coord_residual(nls :: MathProgNLSModel, x :: AbstractVector, v :: AbstractVector)
  increment!(nls, :neval_hess_residual)
  n = nls.meta.nvar
  MathProgBase.eval_hesslag(nls.Fmodel.eval, nls.Fhvals, x, 0.0, v)
  return nls.Fhrows, nls.Fhcols, nls.Fhvals
end

function NLPModels.jth_hess_residual(nls :: MathProgNLSModel, x :: AbstractVector, i :: Int)
  increment!(nls, :neval_jhess_residual)
  y = [j == i ? 1.0 : 0.0 for j = 1:nls.nls_meta.nequ]
  n = nls.meta.nvar
  MathProgBase.eval_hesslag(nls.Fmodel.eval, nls.Fhvals, x, 0.0, y)
  return sparse(nls.Fhrows, nls.Fhcols, nls.Fhvals, n, n)
end

function NLPModels.hprod_residual!(nls :: MathProgNLSModel, x :: AbstractVector, i :: Int, v :: AbstractVector, Hiv :: AbstractVector)
  increment!(nls, :neval_hprod_residual)
  y = [j == i ? 1.0 : 0.0 for j = 1:nls.nls_meta.nequ]
  MathProgBase.eval_hesslag_prod(nls.Fmodel.eval, Hiv, x, v, 0.0, y)
  return Hiv
end

function NLPModels.obj(nls :: MathProgNLSModel, x :: AbstractVector)
  increment!(nls, :neval_obj)
  return MathProgBase.eval_f(nls.cmodel.eval, x)
end

function NLPModels.grad(nls :: MathProgNLSModel, x :: AbstractVector)
  g = zeros(nls.meta.nvar)
  return grad!(nls, x, g)
end

function NLPModels.grad!(nls :: MathProgNLSModel, x :: AbstractVector, g :: AbstractVector)
  increment!(nls, :neval_grad)
  MathProgBase.eval_grad_f(nls.cmodel.eval, g, x)
  return g
end

function NLPModels.cons(nls :: MathProgNLSModel, x :: AbstractVector)
  c = zeros(nls.meta.ncon)
  return cons!(nls, x, c)
end

function NLPModels.cons!(nls :: MathProgNLSModel, x :: AbstractVector, c :: AbstractVector)
  increment!(nls, :neval_cons)
  MathProgBase.eval_g(nls.cmodel.eval, c, x)
  return c
end

function NLPModels.jac_coord(nls :: MathProgNLSModel, x :: AbstractVector)
  increment!(nls, :neval_jac)
  MathProgBase.eval_jac_g(nls.cmodel.eval, nls.cjvals, x)
  return (nls.cjrows, nls.cjcols, nls.cjvals)
end

function NLPModels.jac(nls :: MathProgNLSModel, x :: AbstractVector)
  return sparse(jac_coord(nls, x)..., nls.meta.ncon, nls.meta.nvar)
end

function NLPModels.jprod(nls :: MathProgNLSModel, x :: AbstractVector, v :: AbstractVector)
  Jv = zeros(nls.meta.ncon)
  return jprod!(nls, x, v, Jv)
end

function NLPModels.jprod!(nls :: MathProgNLSModel,
                x :: AbstractVector,
                v :: AbstractVector,
                Jv :: AbstractVector)
  nls.counters.counters.neval_jac -= 1
  increment!(nls, :neval_jprod)
  Jv[:] = jac(nls, x) * v
  return Jv
end

function NLPModels.jtprod(nls :: MathProgNLSModel, x :: AbstractVector, v :: AbstractVector)
  Jtv = zeros(nls.meta.nvar)
  return jtprod!(nls, x, v, Jtv)
end

function NLPModels.jtprod!(nls :: MathProgNLSModel,
                x :: AbstractVector,
                v :: AbstractVector,
                Jtv :: AbstractVector)
  nls.counters.counters.neval_jac -= 1
  increment!(nls, :neval_jtprod)
  Jtv[1:nls.meta.nvar] = jac(nls, x)' * v
  return Jtv
end

function NLPModels.hess_coord(nls :: MathProgNLSModel, x :: AbstractVector;
    obj_weight :: Float64=1.0, y :: AbstractVector=zeros(nls.meta.ncon))
  increment!(nls, :neval_hess)
  MathProgBase.eval_hesslag(nls.cmodel.eval, nls.chvals, x, obj_weight, y)
  return (nls.chrows, nls.chcols, nls.chvals)
end

function NLPModels.hess(nls :: MathProgNLSModel, x :: AbstractVector;
    obj_weight :: Float64=1.0, y :: AbstractVector=zeros(nls.meta.ncon))
  return sparse(hess_coord(nls, x, y=y, obj_weight=obj_weight)...,
                nls.meta.nvar, nls.meta.nvar)
end

function NLPModels.hprod(nls :: MathProgNLSModel, x :: AbstractVector, v :: AbstractVector;
    obj_weight :: Float64=1.0, y :: AbstractVector=zeros(nls.meta.ncon))
  hv = zeros(nls.meta.nvar)
  return hprod!(nls, x, v, hv, obj_weight=obj_weight, y=y)
end

#=
# Removed due to bug https://github.com/JuliaOpt/JuMP.jl/issues/1204
function NLPModels.hprod!(nls :: MathProgNLSModel, x :: AbstractVector, v :: AbstractVector,
    hv :: AbstractVector;
    obj_weight :: Float64=1.0, y :: AbstractVector=zeros(nls.meta.ncon))
  MathProgBase.eval_hesslag_prod(nls.cmodel.eval, hv, x, v, obj_weight, y)
  return hv
end
=#

function NLPModels.hprod!(nls :: MathProgNLSModel, x :: AbstractVector, v :: AbstractVector,
    hv :: AbstractVector;
    obj_weight :: Float64=1.0, y :: AbstractVector=zeros(nls.meta.ncon))
  increment!(nls, :neval_hprod)
  # See bug https://github.com/JuliaOpt/JuMP.jl/issues/1204
  MathProgBase.eval_hesslag_prod(nls.cmodel.eval, hv, x, v, 0.0, y)
  n = nls.meta.nvar
  if obj_weight != 0.0
    Fx = residual(nls, x)
    Jv = jprod_residual(nls, x, v)
    w = jtprod_residual(nls, x, Jv)
    hv[1:n] .+= w
    m = length(Fx)
    for i = 1:m
      hprod_residual!(nls, x, i, v, w)
      @views hv[1:n] .= hv[1:n] .+ Fx[i] * w
    end
    hv[1:n] .*= obj_weight
  end
  return hv
end


text: #include("./sampleTest.jl")