text:
#!/usr/bin/env julia

using ArgParse

s = ArgParseSettings()
@add_arg_table s begin
    "-M"
    help = "Print out the list of files this program will create."
    action = :store_true
    "config"
    help = "a YAML configuration file"
    default = "config.yaml"
end

parsed_args = parse_args(ARGS, s)

import YAML
config = YAML.load(open(parsed_args["config"]))

using DiskJockey.model
using DiskJockey.constants

species = config["species"]
transition = config["transition"]
lam0 = lam0s[species*transition]
model = config["model"]
pars = convert_dict(config["parameters"], config["model"])

if parsed_args["M"]
    args = "velocity.png temperature.png scale_height.png surface_density.png density.png grid_topgrid.png"
    # If we have vertical gradient, add some extra to this
    println(args)
    quit()
end

import PyPlot.plt
using LaTeXStrings

# Plot looking down, in polar coordinates
function plot_topgrid(pars::AbstractParameters, grid::Grid)
    fig = plt[:figure](figsize=(8,8))
    ax = fig[:add_subplot](111, polar=true)

    # Something to span the circle
    phis = linspace(0, 2pi, 100)

    # Plot the grid cell edges in radius and phi
    for R in grid.Rs
        rr = R * ones(phis)
        ax[:plot](phis, rr ./AU, "b", lw=0.1)
    end

    fig[:subplots_adjust](left=0.15, bottom=0.15, right=0.85)
    plt[:savefig]("grid_topgrid.png", dpi=300)
end

function plot_sidegrid(pars::AbstractParameters, grid::Grid)
end

# velocity structure
function plot_vel(pars::AbstractParameters, grid::Grid)
    vels = DiskJockey.model.velocity(grid.rs, pars) .* 1e-5 # convert from cm/s to km/s

    fig = plt[:figure]()
    ax = fig[:add_subplot](111)
    ax[:semilogx](rr, vels)

    # Now, go overlay small grey lines vertically
    for cell_edge in grid.Rs/AU
        ax[:axvline](cell_edge, color="0.5", lw=0.4)
    end

    ax[:set_ylabel](L"$v_\phi$ [km/s]")
    ax[:set_xlabel](L"$r$ [AU]")
    fig[:subplots_adjust](left=0.15, bottom=0.15, right=0.85)

    plt[:savefig]("velocity.png")
end

# Instead of a 1D slice, we'll need to do a 2D field
function plot_vel(pars::ParametersVertical, grid::Grid)
    # Instead of spherical coordinates, do this with cartesian
    nz = 64
    zs = linspace(0.0, 20 * AU, nz)
    zz = zs./AU

    nr = grid.nr
    rs = grid.rs

    xx = Array(Float64, (nz, nr))
    yy = Array(Float64, (nz, nr))
    vels = Array(Float64, (nz, nr))

    for i=1:nz
        xx[i, :] = rr
    end

    for j=1:nr
        yy[:, j] = zz
    end

    for i=1:nz
        for j=1:nr
            vels[i,j] = DiskJockey.model.velocity(grid.rs[j], zs[i], pars)
        end
    end

    levels = Float64[0.0, 1.0, 5.0, 10.0, 20.0, 40.]

    fig = plt[:figure]()
    ax = fig[:add_subplot](111)

    ax[:set_ylabel](L"$z$ [AU]")
    ax[:set_xlabel](L"$r$ [AU]")

    #ticks = np.linspace(0, np.max(cov), num=6)
    #cb.set_ticks(ticks)
    #cb.set_ticks(MaxNLocator(nbins=5))

    #Plot the contoured density in cylindrical coordinates, then plot the spherical grid on top of it?
    # Do this by plotting a bunch of lines
    # First, the radial lines

    # Basically, each of these originates from x = 0, z = 0, and has a slope of theta
    # for theta in (pi/2 - grid.Thetas)
    #     slope = tan(theta)
    #     ax[:plot](rr, slope .* rr, "k", lw=0.1)
    # end
    #
    # # Arcs
    # xs = Array(Float64, (grid.ntheta + 1))
    # ys = Array(Float64, (grid.ntheta + 1))
    #
    # for r in rr
    #     for (i, theta) in enumerate(pi/2 - grid.Thetas)
    #         xs[i] = cos(theta) * r
    #         ys[i] = sin(theta) * r
    #     end
    #     ax[:plot](xs, ys, "k", lw=0.1)
    # end
    #
    # ax[:set_ylim](0, maximum(zz))

    # convert from cm/s to km/s
    img = ax[:contourf](xx, yy, 1e-5 * vels, levels=levels)

    ax[:set_xscale]("log")

    # # Now, go overlay small grey lines vertically
    # for cell_edge in grid.Rs/AU
    #     ax[:axvline](cell_edge, color="0.5", lw=0.4)
    # end

    fig[:subplots_adjust](left=0.15, bottom=0.15, right=0.77)

    cax = fig[:add_axes]([0.82, 0.22, 0.03, 0.65])
    cb = fig[:colorbar](img, cax=cax)

    plt[:savefig]("velocity.png")
end

# temperature structure
function plot_temp(pars::AbstractParameters, grid::Grid)
    temps = DiskJockey.model.temperature(grid.rs, pars)

    fig = plt[:figure]()
    ax = fig[:add_subplot](111)
    ax[:semilogx](rr, temps)

    # Now, go overlay small grey lines vertically
    for cell_edge in grid.Rs/AU
        ax[:axvline](cell_edge, color="0.5", lw=0.4)
    end

    ax[:set_ylabel](L"$T$ [K]")
    ax[:set_xlabel](L"$r$ [AU]")
    fig[:subplots_adjust](left=0.15, bottom=0.15, right=0.85)

    plt[:savefig]("temperature.png")
end

function plot_temp(pars::ParametersVertical, grid::Grid)
    # Instead of spherical coordinates, do this with cartesian
    nz = 64
    zs = linspace(0.0, 20 * AU, nz)
    zz = zs./AU

    nr = grid.nr
    rs = grid.rs

    xx = Array(Float64, (nz, nr))
    yy = Array(Float64, (nz, nr))
    temps = Array(Float64, (nz, nr))

    for i=1:nz
        xx[i, :] = rr
    end

    for j=1:nr
        yy[:, j] = zz
    end

    for i=1:nz
        for j=1:nr
            temps[i,j] = DiskJockey.model.temperature(grid.rs[j], zs[i], pars)
        end
    end

    levels = Float64[0.0, 5.0, 10.0, 20.0, 30., 40., 50.0, 100.0]

    fig = plt[:figure]()
    ax = fig[:add_subplot](111)

    ax[:set_ylabel](L"$z$ [AU]")
    ax[:set_xlabel](L"$r$ [AU]")

    #ticks = np.linspace(0, np.max(cov), num=6)
    #cb.set_ticks(ticks)
    #cb.set_ticks(MaxNLocator(nbins=5))

    #Plot the contoured density in cylindrical coordinates, then plot the spherical grid on top of it?
    # Do this by plotting a bunch of lines
    # First, the radial lines

    # Basically, each of these originates from x = 0, z = 0, and has a slope of theta
    # for theta in (pi/2 - grid.Thetas)
    #     slope = tan(theta)
    #     ax[:plot](rr, slope .* rr, "k", lw=0.1)
    # end
    #
    # # Arcs
    # xs = Array(Float64, (grid.ntheta + 1))
    # ys = Array(Float64, (grid.ntheta + 1))
    #
    # for r in rr
    #     for (i, theta) in enumerate(pi/2 - grid.Thetas)
    #         xs[i] = cos(theta) * r
    #         ys[i] = sin(theta) * r
    #     end
    #     ax[:plot](xs, ys, "k", lw=0.1)
    # end
    #
    # ax[:set_ylim](0, maximum(zz))

    img = ax[:contourf](xx, yy, temps, levels=levels)

    ax[:set_xscale]("log")

    # # Now, go overlay small grey lines vertically
    # for cell_edge in grid.Rs/AU
    #     ax[:axvline](cell_edge, color="0.5", lw=0.4)
    # end

    fig[:subplots_adjust](left=0.15, bottom=0.15, right=0.77)

    cax = fig[:add_axes]([0.82, 0.22, 0.03, 0.65])
    cb = fig[:colorbar](img, cax=cax)

    plt[:savefig]("temperature.png")
end

# scale height
function plot_height(pars::AbstractParameters, grid::Grid)
    heights = DiskJockey.model.Hp(grid.rs, pars) ./ AU

    fig = plt[:figure]()
    ax = fig[:add_subplot](111)
    ax[:semilogx](rr, heights)
    ax[:set_ylabel](L"$H_p$ [AU]")
    ax[:set_xlabel](L"$r$ [AU]")
    fig[:subplots_adjust](left=0.15, bottom=0.15, right=0.85)

    plt[:savefig]("scale_height.png")
end

function plot_surface_density(pars::AbstractParameters, grid::Grid)

    fig = plt[:figure]()
    ax = fig[:add_subplot](111)

    Sigmas = Array(Float64, grid.nr)

    for i=1:grid.nr
        Sigmas[i] = DiskJockey.model.Sigma(grid.rs[i], pars)
    end

    # ax[:loglog](rr, Sigmas)
    ax[:semilogy](rr, Sigmas)

    # Now, go overlay small grey lines vertically for the radial cells
    for cell_edge in grid.Rs/AU
        ax[:axvline](cell_edge, color="0.5", lw=0.4)
    end

    ax[:set_ylabel](L"$\Sigma\, [\mathrm{g/cm}^2]$")
    ax[:set_xlabel](L"$r$ [AU]")
    fig[:subplots_adjust](left=0.15, bottom=0.15, right=0.85)

    plt[:savefig]("surface_density.png")
end

# density structure
function plot_dens(pars::AbstractParameters, grid)

    # Instead of spherical coordinates, do this with cartesian
    nz = 64
    zs = linspace(0.0, 150 * AU, nz)
    zz = zs./AU

    nr = grid.nr
    rs = grid.rs

    xx = Array(Float64, (nz, nr))
    yy = Array(Float64, (nz, nr))
    rhos = Array(Float64, (nz, nr))

    for i=1:nz
        xx[i, :] = rr
    end

    for j=1:nr
        yy[:, j] = zz
    end

    for i=1:nz
        for j=1:nr
            rhos[i,j] = DiskJockey.model.rho_gas(grid.rs[j], zs[i], pars)
        end
    end

    nlog = log10(rhos/(mu_gas * m_H))

    levels = Float64[0.0, 1.0, 2.0, 3.0, 4.0, 5, 6, 7, 8, 9]

    fig = plt[:figure]()
    ax = fig[:add_subplot](111)

    ax[:set_ylabel](L"$z$ [AU]")
    ax[:set_xlabel](L"$r$ [AU]")

    #ticks = np.linspace(0, np.max(cov), num=6)
    #cb.set_ticks(ticks)
    #cb.set_ticks(MaxNLocator(nbins=5))

    #Plot the contoured density in cylindrical coordinates, then plot the spherical grid on top of it?
    # Do this by plotting a bunch of lines
    # First, the radial lines

    # Basically, each of these originates from x = 0, z = 0, and has a slope of theta
    # for theta in (pi/2 - grid.Thetas)
    #     slope = tan(theta)
    #     ax[:plot](rr, slope .* rr, "k", lw=0.1)
    # end
    #
    # # Arcs
    # xs = Array(Float64, (grid.ntheta + 1))
    # ys = Array(Float64, (grid.ntheta + 1))
    #
    # for r in rr
    #     for (i, theta) in enumerate(pi/2 - grid.Thetas)
    #         xs[i] = cos(theta) * r
    #         ys[i] = sin(theta) * r
    #     end
    #     ax[:plot](xs, ys, "k", lw=0.1)
    # end
    #
    # ax[:set_ylim](0, maximum(zz))

    img = ax[:contourf](xx, yy, nlog, levels=levels)

    ax[:set_xscale]("log")

    # Now, go overlay small grey lines vertically
    for cell_edge in grid.Rs/AU
        ax[:axvline](cell_edge, color="0.5", lw=0.4)
    end

    fig[:subplots_adjust](left=0.15, bottom=0.15, right=0.77)

    cax = fig[:add_axes]([0.82, 0.22, 0.03, 0.65])
    cb = fig[:colorbar](img, cax=cax)

    plt[:savefig]("density.png")
end

# Plot a 1D slice of the vertical density structure at a particular radius
function plot_density_column(pars::ParametersVertical, grid::Grid)
    r = 10 * AU

    # Calculate the normalizing constant and height of dissocation
    norm_rho, z_phot = DiskJockey.model.rho_norm_and_phot(r, pars)

    # Calculate a range of zs
    nz = 64
    zs = linspace(0, 4 * AU, nz)

    un_lnrhos = Array(Float64, nz)
    for i=1:nz
        un_lnrhos[i] = DiskJockey.model.un_lnrho(r, zs[i], pars)
    end

    # calculate un_rhos for these zs
    rhos = norm_rho .* exp(un_lnrhos)

    fig = plt[:figure]()
    ax = fig[:add_subplot](111)

    ax[:semilogy](zs./AU, rhos/(mu_gas * m_H))

    # Now, go overlay small grey lines vertically for the radial cells
    # for cell_edge in grid.Rs/AU
    #     ax[:axvline](cell_edge, color="0.5", lw=0.4)
    # end

    ax[:set_ylabel](L"$\rho_\mathrm{gas} \, [\mathrm{n/cm}^3]$")
    ax[:set_xlabel](L"$r$ [AU]")
    fig[:subplots_adjust](left=0.15, bottom=0.15, right=0.85)

    plt[:savefig]("density_column.png")

end

function plot_density_column_CO(pars::ParametersVertical, grid::Grid)
    r = 50 * AU

    # Gas density [g/cm^3]
    zs, rhos = DiskJockey.model.rho_column(r, pars)

    # Gas density of CO [g/cm^3]
    zs, rhos = DiskJockey.model.rho_column_CO(r, zs, rhos, pars)

    fig = plt[:figure]()
    ax = fig[:add_subplot](111)

    ax[:semilogy](zs./AU, DiskJockey.constants.f_12CO * DiskJockey.constants.X_H2 * rhos/(mu_gas * m_H))

    # Now, go overlay small grey lines vertically for the radial cells
    # for cell_edge in grid.Rs/AU
    #     ax[:axvline](cell_edge, color="0.5", lw=0.4)
    # end

    ax[:set_ylabel](L"$\rho_\mathrm{CO} \, [\mathrm{n/cm}^3]$")
    ax[:set_xlabel](L"$z$ [AU]")
    fig[:subplots_adjust](left=0.15, bottom=0.15, right=0.85)

    plt[:savefig]("density_column_CO.png")

end

function test_dens(pars::ParametersVertical, grid::Grid)
    tic()
    rho_gas = DiskJockey.model.rho_gas_interpolator(pars, grid)
    println("Constructing interpolator")
    toc()
    #
    # r_in = grid.rs[1]
    # println("r_in ", r_in)
    # println(rho_gas(r_in, 0.0, pars))

    r2 = 0.12 * AU
    println("r2 ", r2)
    println("un_lnrho ", DiskJockey.model.un_lnrho(r2, 0.0, pars))
    println("norm + zphot ", DiskJockey.model.rho_norm_and_phot(r2, pars))
    tic()
    println(rho_gas(r2, 0.0, pars))
    toc()
end

function plot_norm(pars::ParametersVertical, grid::Grid)

    # DiskJockey.model.rho_norm(10 * AU, pars)
    # DiskJockey.model.z_phot(10 * AU, pars)
    #
    # quit()

    fig = plt[:figure]()
    ax = fig[:add_subplot](111)

    norms = Array(Float64, grid.nr)

    for i=1:grid.nr
        norms[i] = DiskJockey.model.rho_norm(grid.rs[i], pars)
    end

    fig, ax = plt[:subplots](figsize=(6,6))

    # ax[:loglog](rr, Sigmas)
    ax[:loglog](rr, norms)

    ax[:set_ylabel](L"$norm$")

    ax[:set_xlabel](L"$r$ [AU]")
    fig[:subplots_adjust](left=0.15, bottom=0.15, right=0.85)

    plt[:savefig]("norm.png")
end

function plot_ztop(pars::ParametersVertical, grid::Grid)

    fig = plt[:figure]()
    ax = fig[:add_subplot](111)

    ztops = Array(Float64, grid.nr)

    for i=1:grid.nr
        ztops[i] = DiskJockey.model.z_top(grid.rs[i], pars)
    end

    fig, ax = plt[:subplots](figsize=(6,6))

    # ax[:loglog](rr, Sigmas)
    ax[:semilogx](rr, ztops ./ AU)

    ax[:set_ylabel](L"$z_\mathrm{top}$ [AU]")

    ax[:set_xlabel](L"$r$ [AU]")
    fig[:subplots_adjust](left=0.15, bottom=0.15, right=0.85)

    plt[:savefig]("ztop.png")
end

function plot_dens(pars::ParametersVertical, grid::Grid)
    # Set up the interpolator

    DiskJockey.model.rho_gas(10 * AU, 0.1 * AU, pars)

    tic()
    DiskJockey.model.rho_gas(10 * AU, 0.1 * AU, pars)
    println("Typical evaluation time")
    toc()

    # Instead of spherical coordinates, do this with cartesian
    nz = 128
    zs = linspace(0.0, 100 * AU, nz)
    zz = zs./AU

    nr = grid.nr
    rs = grid.rs

    xx = Array(Float64, (nz, nr))
    yy = Array(Float64, (nz, nr))
    rhos = Array(Float64, (nz, nr))

    for i=1:nz
        xx[i, :] = rr
    end

    for j=1:nr
        yy[:, j] = zz
    end

    for i=1:nz
        for j=1:nr
            temp = DiskJockey.model.temperature(grid.rs[j], zs[i], pars)
            rhos[i,j] = DiskJockey.model.X_freeze(temp, pars) * DiskJockey.model.rho_gas(grid.rs[j], zs[i], pars)
        end
    end

    ztops = Array(Float64, grid.nr)
    for i=1:grid.nr
        ztops[i] = DiskJockey.model.z_top(grid.rs[i], pars)
    end

    println("rho extrema ", extrema(rhos))

    nlog = log10(rhos/(mu_gas * m_H))

    levels = Float64[0.0, 1.0, 2.0, 3.0, 4.0, 5, 6, 7, 8, 9, 10, 11, 12]

    fig = plt[:figure]()
    ax = fig[:add_subplot](111)

    ax[:set_ylabel](L"$z$ [AU]")
    ax[:set_xlabel](L"$r$ [AU]")

    img = ax[:contourf](xx, yy, nlog, levels=levels)

    ax[:plot](rr, ztops ./ AU)

    ax[:set_xscale]("log")
    ax[:set_ylim](0, 100)

    # Now, go overlay small grey lines vertically
    # for cell_edge in grid.Rs/AU
    #     ax[:axvline](cell_edge, color="0.5", lw=0.4)
    # end

    fig[:subplots_adjust](left=0.15, bottom=0.15, right=0.77)

    cax = fig[:add_axes]([0.82, 0.22, 0.03, 0.65])
    cb = fig[:colorbar](img, cax=cax)

    plt[:savefig]("density.png")
end

grid = Grid(config["grid"])

# The cell centers for plotting purposes
const global rr = grid.rs ./ AU # convert to AU

if config["model"] == "vertical"
    plot_ztop(pars, grid)
end

plot_topgrid(pars, grid)
plot_vel(pars, grid)
plot_temp(pars, grid)
plot_height(pars, grid)
plot_surface_density(pars, grid)
plot_dens(pars, grid)

max_stars_repo_name: elnjensen/DiskJockey

max_issues_repo_name: elnjensen/DiskJockey

max_forks_repo_name: elnjensen/DiskJockey


text:
@shorthands bivariate
@shorthands bivariatelegend
@shorthands trivariate
@shorthands trivariatelegend

"""
test 1
"""
@recipe function plot(layer::T) where {T<:SimpleSDMLayer}
    seriestype --> :heatmap
    if get(plotattributes, :seriestype, :heatmap) in [:heatmap, :contour]
        eltype(layer) <: AbstractFloat || throw(
            ArgumentError(
                "This plot is only supported for layers with number values ($(eltype(layer)))",
            ),
        )
        aspect_ratio --> 1
        xlims --> extrema(longitudes(layer))
        ylims --> extrema(latitudes(layer))
        lg = copy(layer.grid)
        replace!(lg, nothing => NaN)
        lg = convert(Matrix{Float64}, lg)
        longitudes(layer), latitudes(layer), lg
    elseif get(plotattributes, :seriestype, :histogram) in [:histogram, :density]
        collect(layer)
    end
end

"""
test 2
"""
@recipe function plot(
    l1::FT,
    l2::ST;
    classes::Int=3,
    p0=colorant"#e8e8e8ff",
    p1=colorant"#64acbeff",
    p2=colorant"#c85a5aff",
    quantiles=true,
) where {FT<:SimpleSDMLayer,ST<:SimpleSDMLayer}
    eltype(l1) <: Number || throw(
        ArgumentError(
            "Plotting is only supported for layers with number values ($(eltype(l1)))"
        ),
    )
    eltype(l2) <: Number || throw(
        ArgumentError(
            "Plotting is only supported for layers with number values ($(eltype(l2)))"
        ),
    )
    seriestype --> :scatter
    if get(plotattributes, :seriestype, :scatter) in [:scatter, :histogram2d]
        SimpleSDMLayers._layers_are_compatible(l1, l2)
        valid_i = findall(.!isnothing.(l1.grid) .& .!isnothing.(l2.grid))
        l1.grid[valid_i], l2.grid[valid_i]
    elseif get(plotattributes, :seriestype, :bivariate) in [:bivariate]
        SimpleSDMLayers._layers_are_compatible(l1, l2)
        c1 = LinRange(p0, p1, classes)
        c2 = LinRange(p0, p2, classes)
        breakpoints = LinRange(0.0, 1.0, classes + 1)
        if quantiles
            q1 = rescale(l1, collect(LinRange(0.0, 1.0, 10classes)))
            q2 = rescale(l2, collect(LinRange(0.0, 1.0, 10classes)))
        else
            q1 = rescale(l1, (0.0, 1.0))
            q2 = rescale(l2, (0.0, 1.0))
        end
        classified = similar(l1, Int)
        cols = typeof(p0)[]
        for i in 1:classes
            if isequal(classes)(i)
                fi = (v) -> breakpoints[i] < v <= breakpoints[i + 1]
            else
                fi = (v) -> breakpoints[i] <= v < breakpoints[i + 1]
            end
            m1 = broadcast(fi, q1)
            for j in 1:classes
                if isequal(classes)(j)
                    fj = (v) -> breakpoints[j] < v <= breakpoints[j + 1]
                else
                    fj = (v) -> breakpoints[j] <= v < breakpoints[j + 1]
                end
                m2 = broadcast(fj, q2)
                push!(cols, ColorBlendModes.BlendMultiply(c1[i], c2[j]))
                m = reduce(*, [m1, m2])
                replace!(m, false => nothing)
                if length(m) > 0
                    classified[keys(m)] = fill(length(cols), length(m))
                end
            end
        end
        replace!(classified, 0 => 1)
        @series begin
            seriescolor := vec(cols)
            seriestype := :heatmap
            subplot := 1
            legend --> false
            clims --> (1, classes^2)
            convert(Float16, classified)
        end
    elseif get(plotattributes, :seriestype, :bivariatelegend) in [:bivariatelegend]
        c1 = LinRange(p0, p1, classes)
        c2 = LinRange(p0, p2, classes)
        grid --> :none
        ticks --> :none
        legend --> false
        frametype --> :none
        xlims --> (1 - 0.5, classes + 0.5)
        ylims --> (1 - 0.5, classes + 0.5)
        aspect_ratio --> 1
        cols = Vector{typeof(p0)}(undef, classes^2)
        class = 1
        m = zeros(Float64, classes, classes)
        for i in 1:classes
            for j in 1:classes
                m[j, i] = class
                cols[class] = ColorBlendModes.BlendMultiply(c1[i], c2[j])
                class += 1
            end
        end
        @series begin
            seriescolor := cols
            seriestype := :heatmap
            m
        end
    end
end

"""
test 2
"""
@recipe function plot(
    x::FT, y::ST, z::TT; quantiles=true, simplex=false, red="", green="", blue=""
) where {FT<:SimpleSDMLayer,ST<:SimpleSDMLayer,TT<:SimpleSDMLayer}
    eltype(x) <: Number || throw(
        ArgumentError(
            "Plotting is only supported for layers with number values ($(eltype(x)))"
        ),
    )
    eltype(y) <: Number || throw(
        ArgumentError(
            "Plotting is only supported for layers with number values ($(eltype(y)))"
        ),
    )
    eltype(z) <: Number || throw(
        ArgumentError(
            "Plotting is only supported for layers with number values ($(eltype(z)))"
        ),
    )
    SimpleSDMLayers._layers_are_compatible(x, y)
    SimpleSDMLayers._layers_are_compatible(x, z)
    SimpleSDMLayers._layers_are_compatible(y, z)
    seriestype --> :scatter
    if get(plotattributes, :seriestype, :trivariate) in [:trivariate]
        void = colorant"#ffffff00"
        if quantiles
            X = rescale(x, collect(LinRange(0.0, 1.0, 100)))
            Y = rescale(y, collect(LinRange(0.0, 1.0, 100)))
            Z = rescale(z, collect(LinRange(0.0, 1.0, 100)))
        else
            X = rescale(x, (0.0, 1.0))
            Y = rescale(y, (0.0, 1.0))
            Z = rescale(z, (0.0, 1.0))
        end
        trip = fill(void, size(X))
        for i in CartesianIndices(trip)
            if !isnothing(X[i])
                r = X[i]
                g = Y[i]
                b = Z[i]
                if simplex
                    S = r + g + b
                    if S > zero(typeof(S))
                        r = r / S
                        g = g / S
                        b = b / S
                    end
                end
                trip[i] = RGBA(r, g, b, 1.0)
            end
        end
        @series begin
            seriestype := :heatmap
            subplot := 1
            legend --> false
            yflip --> false
            longitudes(X), reverse(latitudes(X)), trip
        end
    elseif get(plotattributes, :seriestype, :trivariatelegend) in [:trivariatelegend]
        # Legend
        steps = 6
        a = 1 / steps
        h = a * sqrt(3) / 2
        aspect_ratio --> 1
        xlims --> (0, 1)
        ylims --> (-0.1, sqrt(3) / 2 + 0.1)
        legend --> false
        ticks --> :none
        framestyle --> :none
        shapes_x = []
        shapes_y = []
        shapes_colors = []
        for (row_number, triangles_to_do) in enumerate(reverse(1:steps))
            for triangle_number in 1:triangles_to_do
                x0 = (triangle_number - 1) * a + (row_number - 1) * (a / 2)
                y0 = (row_number - 1) * h
                push!(shapes_x, [x0, x0 + a / 2, x0 + a])
                push!(shapes_y, [y0, y0 + h, y0])
                push!(shapes_colors, rgb_from_xy(x0 + a / 2, y0 + h / 2; simplex=simplex))
                if row_number > 1
                    push!(shapes_x, [x0, x0 + a / 2, x0 + a])
                    push!(shapes_y, [y0, y0 - h, y0])
                    push!(shapes_colors, rgb_from_xy(x0 + a / 2, y0 - h / 2; simplex=simplex))
                end
            end
        end
        _fs = get(plotattributes, :annotationfontsize, 14)
        @series begin
            seriestype := :shape
            annotations := [(-0.05, 0.0, (red, _fs, :left, 60.0)), (1.0, -0.05, (green, _fs, :right, 0.0)), (0.5+0.05, sqrt(3)/2, (blue, _fs, :left, -60.0))]
            markersize := 0
            seriescolor := hcat(shapes_colors...)
            shapes_x, shapes_y
        end
    end
end

function rgb_from_xy(x, y; simplex=false)
    b = 2y / sqrt(3)
    g = (2x - b) / 2
    r = 1 - (b + g)
    if !simplex
        m = maximum([r, g, b])
        b = b * (1 / m)
        g = g * (1 / m)
        r = r * (1 / m)
    end
    return RGB(r, g, b)
end

ext: jl


text:
using Flux: Flux
using Metalhead: ResNet, VGG19
using Images

#=
This script is used for submitting issue that shows behavior of ResNet on preference time
Do not expect this will work fine.
=#
#=
functions center_crop and resize_smallest_dimension is
taken from https://github.com/FluxML/Metalhead.jl/blob/master/src/utils.jl
=#

# Take the len-by-len square of pixels at the center of image `im`
function center_crop(im, len)
    l2 = div(len,2)
    adjust = len % 2 == 0 ? 1 : 0
    return im[div(end,2)-l2:div(end,2)+l2-adjust,div(end,2)-l2:div(end,2)+l2-adjust]
end

# Resize an image such that its smallest dimension is the given length
function resize_smallest_dimension(im, len)
    reduction_factor = len/minimum(size(im)[1:2])
    new_size = size(im)
    new_size = (
            round(Int, size(im,1)*reduction_factor),
            round(Int, size(im,2)*reduction_factor),
    )
    if reduction_factor < 1.0
        # Images.jl's imresize() needs to first lowpass the image, it won't do it for us
        im = imfilter(im, KernelFactors.gaussian(0.75/reduction_factor), Inner())
    end
    return imresize(im, new_size)
end

function preprocess(im::AbstractMatrix{<:AbstractRGB})
    # Resize such that smallest edge is 256 pixels long
    im = resize_smallest_dimension(im, 256)

    # Center-crop to 224x224
    im = center_crop(im, 224)

    # Convert to channel view and normalize (these coefficients taken
    # from PyTorch's ImageNet normalization code)
    μ = [0.485, 0.456, 0.406]
    σ = [0.229, 0.224, 0.225]
    im = (channelview(im) .- μ)./σ

    # Convert from CHW (Image.jl's channel ordering) to WHCN (Flux.jl's ordering)
    # and enforce Float64
    return Float64.(permutedims(im, (3, 2, 1))[:,:,:,:].*255)
end

function get_label()
    imagenet_labels = String[]
    label_txt = expanduser("~/.julia/packages/Metalhead/rGGAv/datasets/meta/ILSVRC_synset_mappings.txt")
    for (idx, line) in enumerate(eachline(label_txt))
        synset = line[1:9]
        label = line[11:end]
        push!(imagenet_labels, label)
    end
    return imagenet_labels
end

function main()
    model = ResNet()
    Flux.testmode!(model)
    imagenet_labels=get_label()
    preprocessed = preprocess(load("elephant.jpeg"))
    result = imagenet_labels[Flux.onecold(model(preprocessed))][1]
    @show result
end

main()

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/AbaqusReader.jl/blob/master/LICENSE

import Base.parse
import Base: getindex, length

### Model definitions for ABAQUS data model

abstract type AbstractMaterial end
abstract type AbstractMaterialProperty end
abstract type AbstractProperty end
abstract type AbstractStep end
abstract type AbstractBoundaryCondition end
abstract type AbstractOutputRequest end

mutable struct Mesh
    nodes :: Dict{Int, Vector{Float64}}
    node_sets :: Dict{String, Vector{Integer}}
    elements :: Dict{Int, Vector{Int}}
    element_types :: Dict{Int, Symbol}
    element_sets :: Dict{String, Vector{Integer}}
    surface_sets :: Dict{String, Vector{Tuple{Int, Symbol}}}
    surface_types :: Dict{String, Symbol}
end

function Mesh(d::Dict{String, Dict})
    return Mesh(d["nodes"], d["node_sets"], d["elements"],
                d["element_types"], d["element_sets"],
                d["surface_sets"], d["surface_types"])
end

mutable struct Model
    path :: String
    name :: String
    mesh :: Mesh
    materials :: Dict{Symbol, AbstractMaterial}
    properties :: Vector{AbstractProperty}
    boundary_conditions :: Vector{AbstractBoundaryCondition}
    steps :: Vector{AbstractStep}
end

mutable struct SolidSection <: AbstractProperty
    element_set :: Symbol
    material_name :: Symbol
end

mutable struct Material <: AbstractMaterial
    name :: Symbol
    properties :: Vector{AbstractMaterialProperty}
end

mutable struct Elastic <: AbstractMaterialProperty
    E :: Float64
    nu :: Float64
end

mutable struct Step <: AbstractStep
    kind :: Nullable{Symbol} # STATIC, ...
    boundary_conditions :: Vector{AbstractBoundaryCondition}
    output_requests :: Vector{AbstractOutputRequest}
end

mutable struct BoundaryCondition <: AbstractBoundaryCondition
    kind :: Symbol # BOUNDARY, CLOAD, DLOAD, DSLOAD, ...
    data :: Vector
    options :: Dict
end

mutable struct OutputRequest <: AbstractOutputRequest
    kind :: Symbol # NODE, EL, SECTION, ...
    data :: Vector
    options :: Dict
    target :: Symbol # PRINT, FILE
end

### Utility functions to parse ABAQUS .inp file to data model

mutable struct Keyword
    name :: String
    options :: Vector{Union{String, Pair}}
end

mutable struct AbaqusReaderState
    section :: Nullable{Keyword}
    material :: Nullable{AbstractMaterial}
    property :: Nullable{AbstractProperty}
    step :: Nullable{AbstractStep}
    data :: Vector{String}
end

function get_data(state::AbaqusReaderState)
    data = []
    for row in state.data
        row = strip(row, [' ', ','])
        col = split(row, ',')
        col = map(Meta.parse, col)
        push!(data, col)
    end
    return data
end

function get_options(state::AbaqusReaderState)
    return Dict(get(state.section).options)
end

function get_option(state::AbaqusReaderState, what::String)
    return get_options(state)[what]
end

function length(state::AbaqusReaderState)
    return length(state.data)
end

function is_comment(line)
    return startswith(line, "**")
end

function is_keyword(line)
    return startswith(line, "*") && !is_comment(line)
end

function parse_keyword(line; uppercase_keyword=true)
    args = split(line, ",")
    args = map(String, map(strip, args))
    keyword_name = strip(args[1], '*')
    if uppercase_keyword
        keyword_name = uppercase(keyword_name)
    end
    keyword = Keyword(keyword_name, [])
    for option in args[2:end]
        pair = map(String, split(option, "="))
        if uppercase_keyword
            pair[1] = uppercase(pair[1])
        end
        if length(pair) == 1
            push!(keyword.options, pair[1])
        elseif length(pair) == 2
            push!(keyword.options, pair[1] => pair[2])
        else
            error("Keyword failure: $line, $option, $pair")
        end
    end
    return keyword
end

function is_new_section(line)
    is_keyword(line) || return false
    section = parse_keyword(line)
    is_abaqus_keyword_registered(section.name) || return false
    return true
end

# open_section! is called right after keyword is found
function open_section! end

# close_section! is called at the end or section or before new keyword
function close_section! end

function maybe_close_section!(model, state)
    global close_section!
    isnull(state.section) && return
    section_name = get(state.section).name
    @info("Close section: $section_name")
    args = Tuple{Model, AbaqusReaderState, Type{Val{Symbol(section_name)}}}
    if hasmethod(close_section!, args)
        close_section!(model, state, Val{Symbol(section_name)})
    else
        @warn("no close_section! found for $section_name")
    end
    state.section = nothing
end

function maybe_open_section!(model, state)
    global open_section!
    section_name = get(state.section).name
    section_options = get(state.section).options
    @info("New section: $section_name with options $section_options")
    args = Tuple{Model, AbaqusReaderState, Type{Val{Symbol(section_name)}}}
    if hasmethod(open_section!, args)
        open_section!(model, state, Val{Symbol(section_name)})
    else
        @warn("no open_section! found for $section_name")
    end
end

function new_section!(model, state, line::String)
    maybe_close_section!(model, state)
    state.data = []
    state.section = parse_keyword(line)
    maybe_open_section!(model, state)
end

function process_line!(model, state, line::String)
    if isnull(state.section)
        @info("section = nothing! line = $line")
        return
    end
    if is_keyword(line)
        @warn("missing keyword? line = $line")
        # close section, this is probably keyword and collecting data should stop.
        maybe_close_section!(model, state)
        return
    end
    push!(state.data, line)
end

"""
    abaqus_read_model(filename::String)

Read ABAQUS model from file. Include also boundary conditions, load steps
and so on. If only mesh is needed, it's better to use `abaqus_read_mesh`
insted.
"""
function abaqus_read_model(fn::String)

    model_path = dirname(fn)
    model_name = first(splitext(basename(fn)))
    mesh = Mesh(open(parse_abaqus, fn))
    materials = Dict()
    model = Model(model_path, model_name, mesh, materials, [], [], [])

    state = AbaqusReaderState(nothing, nothing, nothing, nothing, [])

    fid = open(fn)
    for line in eachline(fid)
        line = convert(String, strip(line))
        is_comment(line) && continue
        if is_new_section(line)
            new_section!(model, state, line)
        else
            process_line!(model, state, line)
        end
    end
    close(fid)
    maybe_close_section!(model, state)

    return model
end

### Code to parse ABAQUS .inp to data model

# add here only keywords when planning to define open_section! and/or
# close_section!, i.e. actually parse keyword to model. little bit of
# magic is happening here, but after calling macro there is typealias
# defined i.e. typealias SOLID_SECTION Type{Val{Symbol("SOLID_SECTION")}}
# and also is_keyword_registered("SOLID SECTION") returns true after
# registration, also notice underscoring

SOLID_SECTION = register_abaqus_keyword("SOLID SECTION")

MATERIAL = register_abaqus_keyword("MATERIAL")
ELASTIC = register_abaqus_keyword("ELASTIC")

STEP = register_abaqus_keyword("STEP")
STATIC = register_abaqus_keyword("STATIC")
END_STEP = register_abaqus_keyword("END STEP")

BOUNDARY = register_abaqus_keyword("BOUNDARY")
CLOAD = register_abaqus_keyword("CLOAD")
DLOAD = register_abaqus_keyword("DLOAD")
DSLOAD = register_abaqus_keyword("DSLOAD")
const BOUNDARY_CONDITIONS = Union{BOUNDARY, CLOAD, DLOAD, DSLOAD}

NODE_PRINT = register_abaqus_keyword("NODE PRINT")
EL_PRINT = register_abaqus_keyword("EL PRINT")
SECTION_PRINT = register_abaqus_keyword("SECTION PRINT")
const OUTPUT_REQUESTS = Union{NODE_PRINT, EL_PRINT, SECTION_PRINT}

## Properties

function open_section!(model, state, ::SOLID_SECTION)
    element_set = Symbol(get_option(state, "ELSET"))
    material_name = Symbol(get_option(state, "MATERIAL"))
    property = SolidSection(element_set, material_name)
    state.property = property
    push!(model.properties, property)
end

function close_section!(model, state, ::SOLID_SECTION)
    name = model.name
    state.property = nothing
end

## Materials

function open_section!(model, state, ::MATERIAL)
    material_name = Symbol(get_option(state, "NAME"))
    material = Material(material_name, [])
    state.material = material
    model.materials[material_name] = material
end

function close_section!(model, state, ::ELASTIC)
    name = model.name
    @assert length(state) == 1
    E, nu = first(get_data(state))
    material_property = Elastic(E, nu)
    material = get(state.material)
    push!(material.properties, material_property)
end

## Steps

function open_section!(model, state, ::STEP)
    name = model.name
    step_ = Step(nothing, Vector(), Vector())
    state.step = step_
    push!(model.steps, step_)
end

function open_section!(model, state, ::STATIC)
    name = model.name
    isnull(state.step) && error("*STATIC outside *STEP ?")
    get(state.step).kind = :STATIC
end

function open_section!(model, state, ::END_STEP)
    name = model.name
    state.step = nothing
end

## Steps -- boundary conditions

function close_section!(model, state, ::BOUNDARY_CONDITIONS)
    kind = Symbol(get(state.section).name)
    data = get_data(state)
    options = get_options(state)
    bc = BoundaryCondition(kind, data, options)
    if isnull(state.step)
        push!(model.boundary_conditions, bc)
    else
        step_ = get(state.step)
        push!(step_.boundary_conditions, bc)
    end
end

## Steps -- output requests

function close_section!(model, state, ::OUTPUT_REQUESTS)
    name = model.name
    kind, target = map(Meta.parse, split(get(state.section).name, " "))
    data = get_data(state)
    options = get_options(state)
    request = OutputRequest(Symbol(kind), data, options, Symbol(target))
    step_ = get(state.step)
    push!(step_.output_requests, request)
end


text:
using Calculus
using Gadfly
using LaTeXStrings # for L"" strings
using Compat
using QuadGK
using SpecialFunctions

# WARNING: integrate(f,a,b) is deprecated, use (quadgk(f,a,b))[1] instead
# f(x) = integrate(z -> besselj(1, z), 0.0, x)

f(x) = QuadGK.quadgk(z -> besselj(1, z), 0.0, x)[1]

function my_draw(f)
    my_plot_pgf = plot(f, 0, 100
                , Guide.XLabel("r"
                        , orientation=:horizontal)
                , Guide.YLabel(String(L"\int_{z=0}^{z=\textbf{r}}J_1(z)dz") # seems Gadfly doesn't support LaTeX L"" anymore
                                                                    # workaround: edit manually PGF .tex file
                #, Guide.YLabel("∫₀ʳJ₁(z)dz" # it is very ugly
                        , orientation=:vertical))

    draw(PGF("images/Int_J1_tex2.tex", 5inch, 3inch), my_plot_pgf)
end

# Gadfly used to support LaTeX labels but doesn't support anymore
# [Does Gadfly support LaTeX in plot titles and labels?](http://stackoverflow.com/questions/36880531/does-gadfly-support-latex-in-plot-titles-and-labels)
# see also [PGF LaTeX output seems to be broken #753](https://github.com/dcjones/Gadfly.jl/issues/753)

@time my_draw(f)
@time my_draw(f)

#using FastAnonymous
#y = @anon z -> besselj(1, z)
# FastAnonymous doesn't work for 0.5, 0.6

y = z -> besselj(1, z)
g(x) = QuadGK.quadgk(y, 0.0, x)[1]

@time my_draw(g)
@time my_draw(g)

max_stars_repo_name: alsam/julia-samples

max_issues_repo_name: alsam/julia-samples

max_forks_repo_name: alsam/julia-samples


text:
using CSV
using Test
using Logging
using Dates
using DataFrames
using YAML

using SLiDE

import InfrastructureSystems
const IS = InfrastructureSystems

SRC_DIR = joinpath(SLIDE_DIR, "src");

max_stars_repo_name: matbesancon/SLiDE

max_issues_repo_name: matbesancon/SLiDE

max_forks_repo_name: matbesancon/SLiDE


text:
#Eulers

function dyn_update_Euler(f,Δt;par)
    return(par .+ f(par).*Δt)
end

function dyn_update_RT4(f,Δt;par)
    #todo

end

function dyn_update_advMidP(f,Δt;par)
    #todo
end

max_stars_repo_path: Eulers Method.jl

max_issues_repo_path: Eulers Method.jl

max_forks_repo_path: Eulers Method.jl


text:
abstract type AbstractVariable <: JuMP.AbstractVariable end

abstract type HintPoint end

# It will not really be the center and the center for z = 0 is the the same as for <h, x> = 0
struct CenterPoint{T} <: HintPoint
    h::Vector{T}
end

struct InteriorPoint{T} <: HintPoint
    h::Vector{T}
end
_β(model, h::InteriorPoint) = @variable(model)
_b(model, h::InteriorPoint) = @variable(model, [1:length(h.h)])

function polar_perspective_ellipsoid(ell, point::HintPoint, z::SpaceVariable,
                                     x::Vector{SpaceVariable})
    y = [z; x]
    H = Sets._householder(point.h)
    p = y' * H * Sets._HPH(ell) * H * y
    return Sets.perspective_dual(Sets.Householder(ell, p, point.h, z, x))
end
function polar_perspective_ellipsoid(model, Q::Symmetric{JuMP.VariableRef},
                                     point::CenterPoint, z, x)
    psd_constraint(model, Q)
    ell = Sets.Ellipsoid(Q)
    return polar_perspective_ellipsoid(ell, point, z, x)
end
function polar_perspective_ellipsoid(model, Q::Symmetric{JuMP.VariableRef},
                                     point::InteriorPoint, z, x)
    n = LinearAlgebra.checksquare(Q)
    @assert n == length(point.h)
    β = @variable(model, base_name="β")
    b = @variable(model, [1:length(point.h)], base_name="b")
    psd_constraint(model, Symmetric([β+1 b'; b Q]))
    ell = Sets.ShiftedEllipsoid(Q, b, β)
    return polar_perspective_ellipsoid(ell, point, z, x)
end

function perspective_dual_polyset(set, point::HintPoint, z::SpaceVariable,
                                  x::Vector{SpaceVariable})
    y = [z; x]
    H = Sets._householder(point.h)
    p = Sets.perspective_gauge0(set)(y => H * y)
    return Sets.perspective_dual(Sets.Householder(set, p, point.h, z, x))
end
# For `CenterPoint`, `q` should be non-perspective, need refactoring
function perspective_dual_polyset(degree, q, point::InteriorPoint, z, x)
    set = Sets.ConvexPolynomialSet(degree, q, z, x)
    perspective_dual_polyset(set, point, z, x)
end

### Polytope ###
struct Polytope <: AbstractVariable
    symmetric::Bool
    dimension::Union{Nothing, Int}
    piecewise::Union{Polyhedra.Rep, Nothing}
end
function Polytope(;
    symmetric::Bool=false,
    dimension::Union{Int, Nothing}=nothing,
    piecewise::Union{Polyhedra.Rep, Nothing}=nothing,
)
    function update_dim(object, dim_fun)
        if object !== nothing
            d = dim_fun(object)
            if dimension === nothing
                dimension = d
            elseif dimension != d
                throw(DimensionMismatch())
            end
        end
    end
    update_dim(piecewise, Polyhedra.fulldim)
    return Polytope(symmetric, dimension, piecewise)
end
Sets.space_variables(::Polytope) = nothing

function variable_set(model::JuMP.AbstractModel, ell::Polytope, space::Space,
                      space_dimension, space_polyvars)
    n = space_dimension
    if ell.symmetric
        if ell.piecewise === nothing
            set = Sets.PolarPoint(@variable(model, [1:n], base_name = "a"))
        else
            hashyperplanes(ell.piecewise) && error("hyperplanes not supported for piecewise")
            a = @variable(model, [1:nhalfspaces(ell.piecewise), 1:n], base_name="a")
            sets = [Sets.PolarPoint(a[i, :]) for i in 1:nhalfspaces(ell.piecewise)]
            set = Sets.Piecewise(sets, ell.piecewise)
            for i in eachindex(set.graph)
                for (j, v) in set.graph[i]
                    if i < j # The constraints are the same for (i, j) and (j, i)
                        λ = @variable(model, base_name="λc[$i,$j]")
                        # Ensures continuity
                        @constraint(model, sets[i].a - sets[j].a .== λ * v)
                        # Ensures convexity
                        # We need to ensure that `dot(sets[i].a, x)` is maximal
                        # over when `piece[i]`; see Proposition 5 of [R21].
                        #
                        # [R21] Raković, S. V.
                        # *Control Minkowski–Lyapunov functions*
                        # Automatica, Elsevier BV, 2021, 128, 109598
                        #
                        # In fact, we only need this to hold locally to prove
                        # convexity of each point so we can just check for every
                        # neighbor.
                        # Because <ai, x> = <aj, x> by the previous constraint,
                        # we can check the gradient.
                        @constraint(model, dot(sets[i].a, v) <= dot(sets[j].a, v), base_name="conv[$i,$j]")
                    end
                end
            end
        end
        if space == PrimalSpace
            return set
        else
            @assert space == DualSpace
            return Polyhedra.polar(set)
        end
    else
        error("Non-symmetric polytope not supported yet")
    end
end
function JuMP.value(h::Sets.PolarPoint)
    return Sets.PolarPoint(JuMP.value.(h.a))
end

### Ellipsoid ###
struct Ellipsoid <: AbstractVariable
    point::Union{Nothing, HintPoint}
    symmetric::Bool
    dimension::Union{Nothing, Int}
    guaranteed_psd::Bool # Is it already guaranteed that it is PSD ? e.g. by nth_root
    superset::Union{Sets.Ellipsoid, Nothing}
    piecewise::Union{Polyhedra.Rep, Nothing}
end
function Ellipsoid(; point::Union{Nothing, HintPoint}=nothing,
                   symmetric::Bool=false,
                   dimension::Union{Int, Nothing}=nothing,
                   superset::Union{Sets.Ellipsoid, Nothing}=nothing,
                   piecewise::Union{Polyhedra.Rep, Nothing}=nothing)
    function update_dim(object, dim_fun)
        if object !== nothing
            d = dim_fun(object)
            if dimension === nothing
                dimension = d
            elseif dimension != d
                throw(DimensionMismatch())
            end
        end
    end
    update_dim(point, point -> length(point.h))
    update_dim(superset, Sets.dimension)
    update_dim(piecewise, Polyhedra.fulldim)
    return Ellipsoid(point, symmetric, dimension, false, superset, piecewise)
end
Sets.space_variables(::Ellipsoid) = nothing

function variable_set(model::JuMP.AbstractModel, ell::Ellipsoid, space::Space,
                      space_dimension, space_polyvars)
    n = space_dimension
    function new_Q()
        # TODO, we should use constraiend variable instead in case direct mode is used.
        Q = @variable(model, [1:n, 1:n], Symmetric, base_name="Q")
        if !ell.guaranteed_psd
            psd_constraint(model, Q)
        end
        return Q
    end
    if ell.symmetric
        function new_piece()
            if space == PrimalSpace
                if ell.superset !== nothing
                    Q = Symmetric(new_Q() + ell.superset.Q)
                else
                    Q = new_Q()
                end
                return Sets.Ellipsoid(Q)
            else
                ell.superset === nothing || error("superset not supported in dual space")
                @assert space == DualSpace
                return Sets.Ellipsoid(new_Q())
            end
        end
        if ell.piecewise === nothing
            set = new_piece()
        else
            hashyperplanes(ell.piecewise) && error("hyperplanes not supported for piecewise")
            sets = [new_piece() for i in 1:nhalfspaces(ell.piecewise)]
            set = Sets.Piecewise(sets, ell.piecewise)
            DynamicPolynomials.@polyvar x[1:n]
            q = [quad_form(set.Q, x) for set in sets]
            for i in eachindex(set.graph)
                for (j, v) in set.graph[i]
                    if i < j # The constraints are the same for (i, j) and (j, i)
                        @constraint(model, q[i] == q[j], domain = @set x'v == 0)
                        # v corresponds to `-n_ij` in LCSS paper
                        Δ = sets[i].Q * v - sets[j].Q * v
                        inter = set.pieces[i] ∩ set.pieces[j]
                        h = HalfSpace(Δ, zero(eltype(Δ)))
                        @constraint(model, inter ⊆ h)
                    end
                end
            end
        end
        if space == PrimalSpace
            return set
        else
            return Polyhedra.polar(set)
        end
    else
        ell.superset === nothing || error("superset not supported for non-symmetric Ellipsoid")
        ell.piecewise === nothing || error("piecewise not supported for non-symmetric Ellipsoid")
        if space == PrimalSpace
            error("Non-symmetric ellipsoid non implemented yet, use `Ellipsoid(symmetric=true)`.")
        else
            @assert space == DualSpace
            if ell.point === nothing
                throw(ArgumentError("Specify a point for nonsymmetric ellipsoid, e.g. `Ellipsoid(point=InteriorPoint([1.0, 0.0]))"))
            end
            return polar_perspective_ellipsoid(model, new_Q(), ell.point,
                                               data(model).perspective_polyvar,
                                               space_polyvars)
        end
    end
end
function JuMP.value(ell::Sets.Ellipsoid)
    return Sets.Ellipsoid(Symmetric(JuMP.value.(ell.Q)))
end

### PolySet ###
struct PolySet <: AbstractVariable
    point::Union{Nothing, HintPoint}
    symmetric::Bool
    degree::Int
    dimension::Union{Nothing, Int}
    convex::Bool
    variables::Union{Nothing, Vector{SpaceVariable}}
    superset::Union{Sets.PolySet, Nothing}
    basis::Type
end
function PolySet(; point::Union{Nothing, HintPoint}=nothing,
                 symmetric::Bool=false,
                 degree::Union{Int, Nothing}=nothing,
                 dimension::Union{Int, Nothing}=nothing,
                 convex::Bool=false,
                 variables::Union{Vector{SpaceVariable}, Nothing}=nothing,
                 superset::Union{Sets.PolySet, Nothing}=nothing,
                 basis::Type=MultivariateBases.MonomialBasis)
    if degree === nothing
        error("Degree of PolySet not specified, use PolySet(degree=..., ...)")
    end
    if isodd(degree)
        throw(ArgumentError("Degree of PolySet not even"))
    end
    if dimension === nothing
        if point !== nothing
            dimension = length(point.h)
        end
    end
    if superset !== nothing
        if dimension === nothing
            dimension = Sets.dimension(superset)
        elseif dimension != Sets.dimension(superset)
            throw(DimensionMismatch())
        end
        if variables === nothing
            variables = Sets.space_variables(superset)
        elseif variables != Sets.space_variables(superset)
            error("Space variables set does not correspond to superset space variables.")
        end
    end
    return PolySet(point, symmetric, degree, dimension, convex, variables, superset, basis)
end
Sets.space_variables(p::PolySet) = p.variables

function constrain_convex(model, p, vars)
    hessian = MP.differentiate(p, vars, 2)
    # We do not just do `@constraint(model, p in SOSConvex())` as we would
    # like to have access to the PSD matrix of variables for the det volume heuristic
    y = [MP.similarvariable(eltype(hessian), gensym()) for i in 1:LinearAlgebra.checksquare(hessian)]
    q = dot(y, hessian * y)
    X = SumOfSquares.Certificate.monomials_half_newton_polytope(MP.monomials(q), (y,))
    # If `X` is empty, we will need the following bridge
    JuMP.add_bridge(model, SumOfSquares.Bridges.Constraint.EmptyBridge)
    # If `length(X)` is 2, we will need the following bridge
    JuMP.add_bridge(model, SumOfSquares.Bridges.Constraint.PositiveSemidefinite2x2Bridge)
    set = SumOfSquares.matrix_cone(MOI.PositiveSemidefiniteConeTriangle,
                                   length(X))
    Q = @variable(model, [1:MOI.dimension(set)])
    @constraint(model, Q in set)
    s = SumOfSquares.build_gram_matrix(
        Q, MonomialBasis(X), MOI.PositiveSemidefiniteConeTriangle, Float64)
    @constraint(model, q == s)
    return MultivariateMoments.getmat(s)
end

function variable_set(model::JuMP.AbstractModel, set::PolySet, space::Space,
                      space_dimension, space_polyvars)
    n = space_dimension
    d = data(model)
    # General all monomials of degree `degree`, we don't want monomials of
    # lower degree as the polynomial is homogeneous
    @assert iseven(set.degree)
    if set.symmetric
        monos = MP.monomials(space_polyvars, div(set.degree, 2))
    else
        monos = MP.monomials(lift_space_variables(d, space_polyvars),
                             div(set.degree, 2))
    end
    basis = MultivariateBases.basis_covering_monomials(set.basis, monos)
    # TODO If `set.convex` and `set.symmetric`, no need for the poly to be SOS, see Lemma 6.33 of [BPT12]
    p = @variable(model, variable_type=SOSPoly(basis))
    if set.convex
        set.superset === nothing || error("superset not supported for convex PolySet")
        if set.symmetric
            convexity_proof = constrain_convex(model, p, space_polyvars)
            if space == PrimalSpace
                return Sets.ConvexPolySet(set.degree, p, convexity_proof)
            else
                @assert space == DualSpace
                return Polyhedra.polar(Sets.ConvexPolySet(set.degree, p, convexity_proof))
            end
        else
            constrain_convex(model, MP.subs(p, d.perspective_polyvar => 1),
                             space_polyvars)
            if space == PrimalSpace
                error("Non-symmetric PolySet in PrimalSpace not implemented yet")
            else
                @assert space == DualSpace
                if set.point === nothing
                    throw(ArgumentError("Specify a point for nonsymmetric polyset, e.g. `PolySet(point=InteriorPoint([1.0, 0.0]))"))
                end
                return perspective_dual_polyset(set.degree, p, set.point, d.perspective_polyvar, space_polyvars)
            end
        end
    else
        if set.symmetric
            if space == PrimalSpace
                if set.superset !== nothing
                    p = SetProg.SumOfSquares.gram_operate(+, set.superset.p, p)
                end
                return Sets.PolySet(set.degree, p)
            else
                error("Non-convex PolySet not supported in $space")
            end
        else
            error("Non-convex nonsymmetric PolySet not implemented yet.")
        end
    end
end
_value(convexity_proof::Nothing) = nothing
function _value(convexity_proof::MultivariateMoments.SymMatrix)
    return MultivariateMoments.SymMatrix(JuMP.value.(convexity_proof.Q),
                                         convexity_proof.n)
end
function JuMP.value(set::Sets.PolySet)
    return Sets.PolySet(set.degree, JuMP.value(set.p))
end
function JuMP.value(set::Sets.ConvexPolySet)
    return Sets.ConvexPolySet(set.degree, JuMP.value(set.p), _value(set.convexity_proof))
end
function JuMP.value(set::Sets.Polar)
    return Polyhedra.polar(JuMP.value(Polyhedra.polar(set)))
end
function JuMP.value(set::Sets.PerspectiveDual)
    return Sets.perspective_dual(JuMP.value(Sets.perspective_dual(set)))
end
function JuMP.value(set::Sets.Householder)
    return Sets.Householder(JuMP.value(set.set), JuMP.value(set.p), set.h,
                            set.z, set.x)
end
function JuMP.value(set::Sets.ShiftedEllipsoid)
    return Sets.ShiftedEllipsoid(Symmetric(JuMP.value.(set.Q)),
                                 JuMP.value.(set.b), JuMP.value(set.β))
end
function JuMP.value(set::Sets.ConvexPolynomialSet)
    return Sets.ConvexPolynomialSet(set.degree, JuMP.value(set.q), set.z, set.x)
end
function JuMP.value(set::Sets.Piecewise)
    return Sets.Piecewise(JuMP.value.(set.sets), set.polytope, set.pieces, set.graph)
end

### SetVariableRef ###
mutable struct SetVariableRef{M <: JuMP.AbstractModel,
                              S <: AbstractVariable} <: JuMP.AbstractVariableRef
    model::M
    set::S
    name::String
    # `variable` is typically `Sets.AbstractSet{JuMP.VariableRef}` but it can
    # also be `Sets.AbstractSet{JuMP.AffExpr}` with `PolySet(superset = ...)`.
    variable::Union{Nothing, Sets.AbstractSet}
    space_index::Union{Nothing, SpaceIndex}
end
JuMP.name(vref::SetVariableRef) = vref.name
function JuMP.build_variable(_error::Function, info::JuMP.VariableInfo, set::AbstractVariable)
    @assert !info.has_lb && !info.has_ub && !info.has_fix && !info.binary && !info.integer && !info.has_start
    return set
end
function JuMP.add_variable(model::JuMP.AbstractModel, set::AbstractVariable, name::String)
    vref = SetVariableRef(model, set, name, nothing, nothing)
    d = data(model)
    @assert d.state == Modeling
    push!(d.variables, vref)
    return vref
end
JuMP.value(vref::SetVariableRef) = JuMP.value(vref.variable)

function clear_spaces(vref::SetVariableRef)
    vref.space_index = nothing
end
Sets.space_variables(::SetVariableRef) = nothing
Sets.perspective_variable(::SetVariableRef) = nothing
function create_spaces(vref::SetVariableRef, spaces::Spaces)
    if vref.space_index === nothing
        if Sets.space_variables(vref.set) === nothing
            if vref.set.dimension === nothing
                vref.space_index = new_space(spaces)
            else
                vref.space_index = new_space(spaces, vref.set.dimension)
            end
        else
            vref.space_index = new_space(spaces, Sets.space_variables(vref.set))
        end
    end
    return vref.space_index
end
space_index(vref::SetVariableRef) = vref.space_index

function load(model::JuMP.AbstractModel, vref::SetVariableRef)
    d = data(model)
    vref.variable = variable_set(model, vref.set, d.space,
                                 space_dimension(d.spaces, vref.space_index),
                                 space_polyvars(d.spaces, vref.space_index))
end

variablify(p::Sets.Projection) = Polyhedra.project(variablify(p.set), p.indices)


text:
module Chords # Musical

export Chord, Major, Minor
export notes, root, third, minorthird, fifth # , pitches
export Am
import ..Musical.Notations: Notation, A, B, C, D, E, F, G

abstract type Quality end

# M
struct Major <: Quality
end

# m
struct Minor <: Quality
end

# aug
struct Augmented <: Quality
end

# dim
struct Diminished{N<:Int} <: Quality
end

struct HalfDiminished <: Quality
end

struct Dominant <: Quality
end

struct Chord
    notation::Notation
    q::Type{Q} where {Q<:Quality}
end

function notes(chord::Chord)::Tuple
    notes(chord.notation, chord.q)
end

function notes(notation::Notation, q::Type{Major})::Tuple
    (C, E, G)
end

function notes(notation::Notation, q::Type{Minor})::Tuple
    (A, C, E)
end

root(chord::Chord) = notes(chord)[1]
third(chord::Chord) = notes(chord)[2]
minorthird(chord::Chord) = notes(chord)[2]
fifth(chord::Chord) = notes(chord)[3]
# pitches(chord::Chord) =

const Am = Chord(A, Minor)

end # module Musical.Chords


text:
using JuMP

# Inspired from Joey Huchette's test in ConvexHull.jl
function ex1test(lib::Polyhedra.Library)
    m = Model()
    @variable(m, x)
    @variable(m, y)
    @constraints(m, begin
        12 + 2x - y ≥ 0
        -6 - x + 2y ≥ 0
        -3 + x + y ≥ 0
        1 + x ≥ 0
    end)

    poly = polyhedron(m, lib)

    V = [ 0  3;
         -1  4;
         -1 10]

    R = [1 2;
         2 1]

    generator_fulltest(poly, V, R)
end


text:
"""
    CKKSEncoder

A `CKKSEncoder` provides functionality to convert raw data such as scalars and vectors into
`Plaintext` instances using `encode!`, and to convert `Plaintext` elements back to raw data using
`decode!`.

See also: [`Plaintext`](@ref), [`encode!`](@ref), [`decode!`](@ref)
"""
mutable struct CKKSEncoder <: SEALObject
  handle::Ptr{Cvoid}
  context::SEALContext

  function CKKSEncoder(context)
    handleref = Ref{Ptr{Cvoid}}(C_NULL)
    retval = ccall((:CKKSEncoder_Create, libsealc), Clong,
                   (Ptr{Cvoid}, Ref{Ptr{Cvoid}}),
                   context, handleref)
    @check_return_value retval
    return CKKSEncoder(handleref[], context)
  end

  function CKKSEncoder(handle::Ptr{Cvoid}, context)
    object = new(handle, context)
    finalizer(destroy!, object)
    return object
  end
end

function destroy!(object::CKKSEncoder)
  if isallocated(object)
    @check_return_value ccall((:CKKSEncoder_Destroy, libsealc), Clong, (Ptr{Cvoid},), object)
    sethandle!(object, C_NULL)
  end

  return nothing
end

"""
    slot_count(encoder)

Return the number of available slots for a given encoder, i.e., how many raw data values can be
stored and processed simultaneously with the given encryption setup.
"""
function slot_count(encoder::CKKSEncoder)
  count = Ref{UInt64}(0)
  retval = ccall((:CKKSEncoder_SlotCount, libsealc), Clong,
                 (Ptr{Cvoid}, Ref{UInt64}),
                 encoder, count)
  @check_return_value retval
  return Int(count[])
end

"""
    encode!(destination, data::DenseVector{Float64}, scale, encoder)
    encode!(destination, data::Float64, scale, encoder)

Use `CKKSEncoder` instance `encoder` to encode raw `data`, which can either be a scalar or a dense
vector. The result is stored in the `Plaintext` instance `destination` using encoding precision
`scale`. Note that if `data` is a vector, it must have at least as many elements as there are slots
available.

See also: [`slot_count`](@ref)
"""
function encode! end

function encode!(destination::Plaintext, data::DenseVector{Float64}, scale, encoder::CKKSEncoder)
  value_count = UInt64(length(data))
  parms_id = first_parms_id(encoder.context)
  retval = ccall((:CKKSEncoder_Encode1, libsealc), Clong,
                 (Ptr{Cvoid}, UInt64, Ref{Cdouble}, Ref{UInt64}, Float64, Ptr{Cvoid}, Ptr{Cvoid}),
                 encoder, value_count, data, parms_id, scale, destination, C_NULL)
  @check_return_value retval
  return destination
end

function encode!(destination::Plaintext, data::Float64, scale, encoder::CKKSEncoder)
  parms_id = first_parms_id(encoder.context)
  retval = ccall((:CKKSEncoder_Encode3, libsealc), Clong,
                 (Ptr{Cvoid}, Float64, Ref{UInt64}, Float64, Ptr{Cvoid}, Ptr{Cvoid}),
                 encoder, data, parms_id, scale, destination, C_NULL)
  @check_return_value retval
  return destination
end

function encode!(destination::Plaintext, data::Integer, encoder::CKKSEncoder)
  parms_id = first_parms_id(encoder.context)
  retval = ccall((:CKKSEncoder_Encode5, libsealc), Clong,
                 (Ptr{Cvoid}, Int64, Ref{UInt64}, Ptr{Cvoid}),
                 encoder, data, parms_id, destination)
  @check_return_value retval
  return destination
end

"""
    decode!(destination, plain, encoder)

Use `CKKSEncoder` instance `encoder` to convert the `Plaintext` instance `plain` back to raw data.
The result is stored in the dense vector `destination`, which must have at least as many elements as
there are slots available.

See also: [`slot_count`](@ref)
"""
function decode!(destination::DenseVector{Float64}, plain::Plaintext, encoder::CKKSEncoder)
  value_count = UInt64(length(destination))
  retval = ccall((:CKKSEncoder_Decode1, libsealc), Clong,
                 (Ptr{Cvoid}, Ptr{Cvoid}, Ref{UInt64}, Ref{Cdouble}, Ptr{Cvoid}),
                 encoder, plain, value_count, destination, C_NULL)
  @check_return_value retval
  return destination
end


text:
type ADMetadata
  valueSize::Tuple{Int64, Int64}
  opNode::Bool
  gradSize::Tuple{Int64, Int64}
end
type ADNodesMetadata
  id::String
  name::String
  metadata::ADMetadata
  group::Int64

  ADNodesMetadata(node::AD, id::String) = new(id, node.name, ADMetadata(size(node.value), node.opNode, size(node.grad)), node.opNode?1:0)
  ADNodesMetadata(node::AD, id::String, grp::Int64) = new(id, node.name, ADMetadata(size(node.value), node.opNode, size(node.grad)), grp)
end
type ADLinksMetadata
  source::String
  target::String
  value::Int64
end
type D3Data
  nodes::Array{ADNodesMetadata}
  links::Array{ADLinksMetadata}
end

# traverse the graph top-down
# convert AD to ADMetadata
function getGraphMetadata(graph::AD)
  uniqueNodes = Dict()
  nodes = Array(ADNodesMetadata,0)
  links = Array(ADLinksMetadata,0)

  bfs = [graph]
  while length(bfs) != 0
    current = pop!(bfs)
    currentId = string(pointer(current.value))
    if !haskey(uniqueNodes, currentId)
      push!(nodes, ADNodesMetadata(current, currentId))
      uniqueNodes[currentId] = true
    end
    currentPtr = string(pointer(current.value))
    numParents = length(current.parents)
    for i=1:numParents
      push!(links, ADLinksMetadata(currentPtr, string(pointer(current.parents[i].value)), current.level))
      push!(bfs, current.parents[i])
    end
  end

  return D3Data(nodes, links)
end

function d3ForcedDirectedGraphHTML(data::D3Data, supportsHTML::Bool)
  d3HTML =
    """
      <!DOCTYPE html>
      <meta charset="utf-8">
      <style>

        .links line {
          stroke: #999;
          stroke-opacity: 0.6;
        }

        .nodes text {
          pointer-events: none;
          font: 10px sans-serif;
        }

        text{
		      visibility: visible;
	      }

        .nodes circle {
          stroke: #fff;
          stroke-width: 1.5px;
        }

      </style>
      <svg width="960" height="600"></svg>
      <script type=text/javascript>
    """

  if supportsHTML
    # mbostock comment: https://github.com/mpld3/mpld3/pull/34
    d3HTML = d3HTML *
      """
        require.config({
          paths: {
            d3: '//d3js.org/d3.v4.min',
            mathjax: '//edx-static.s3.amazonaws.com/mathjax-MathJax-727332c/MathJax.js?config=TeX-MML-AM_HTMLorMML-full'
          },
          shim: {
            'd3': {
              exports: 'd3',
              init: function() {
                window.d3 = d3;
              }
            },
          }
        });
        require(["d3", "mathjax"], function(d3, mathjax) {
      """
  end

  d3HTML = d3HTML *
    """
      var svg = d3.select("svg"),
      width = +svg.attr("width"),
      height = +svg.attr("height");

      var color = d3.scaleOrdinal(d3.schemeCategory20);
      var simulation = d3.forceSimulation()
                          .force("link", d3.forceLink().id(function(d) { return d.id; }))
                          .force("charge", d3.forceManyBody())
                          .force("center", d3.forceCenter(width / 2, height / 2));

      var graph = $(json(data));
      var link = svg.append("g")
                    .attr("class", "links")
                    .selectAll("line")
                    .data(graph.links)
                    .enter().append("line")
                    .attr("stroke-width", function(d) { return Math.sqrt(d.value); });
      var node = svg.append("g")
                    .attr("class", "nodes")
                    .selectAll("circle")
                    .data(graph.nodes)
                    .enter()
                    .append("circle")
                    .attr("r", 5)
                    .attr("fill", function(d) { return color(d.group); })
                    .call(d3.drag()
                    .on("start", dragstarted)
                    .on("drag", dragged)
                    .on("end", dragended));

      node.append("title")
        .text(function(d) { return JSON.stringify(d.metadata); });
      node.append("text")
        .attr("dx", 12)
        .attr("dy", ".35em")
        .text(function(d) { return d.name });

      simulation
        .nodes(graph.nodes)
        .on("tick", ticked);

      simulation.force("link")
        .links(graph.links);

      setTimeout(textDisplay, 500);

  	  function textDisplay(){
  		    svg.append("g")
  		      .selectAll("text")
  		      .data(d3.selectAll("circle")._groups[0])
            .enter()
  		      .append("text")
  		      .attr("x", function(d,i){ return d.cx.animVal.value+5;})
  		      .attr("y", function(d,i){ return d.cy.animVal.value-5;})
  		      .text(function(d) { return d.__data__.name })
  		      .attr("style","")
  		      .attr("font-family", "sans-serif")
  		      .attr("font-size", "15px");
  	  }

      function ticked() {
        link
          .attr("x1", function(d) { return d.source.x; })
          .attr("y1", function(d) { return d.source.y; })
          .attr("x2", function(d) { return d.target.x; })
          .attr("y2", function(d) { return d.target.y; });
        node
          .attr("cx", function(d) { return d.x; })
          .attr("cy", function(d) { return d.y; });
      }

      function dragstarted(d) {
        if (!d3.event.active) simulation.alphaTarget(0.3).restart();
        d.fx = d.x;
        d.fy = d.y;
        svg.selectAll("text").remove()
	      textDisplay();
      }

      function dragged(d) {
        d.fx = d3.event.x;
        d.fy = d3.event.y;
        svg.selectAll("text").remove()
	      textDisplay();
      }

      function dragended(d) {
        if (!d3.event.active) simulation.alphaTarget(0);
        d.fx = null;
        d.fy = null;
        svg.selectAll("text").remove()
	      textDisplay();
      }
    """

  if supportsHTML
    d3HTML = d3HTML *
      """
        });
      """
  end

  d3HTML = d3HTML *
    """
      </script>
    """

  return d3HTML
end

# D3 operations on graph metadata
function getD3GraphMetadataHTML(graph::AD, supportsHTML::Bool)
  d3Data = getGraphMetadata(graph)
  d3HTML = d3ForcedDirectedGraphHTML(d3Data, supportsHTML)
  return d3HTML
end

# display in a new window (Electron) in case
# text/html mime type is not supported
function blink_show(graph::AD)
    w = Window()
    try
      d3 = joinpath(dirname(@__FILE__), "..", "deps", "d3.min.js")
      for file in (d3,)
        load!(w, file)
      end
    catch
      error("Could not load the file: d3.min.js from the deps directory.")
    end
    body = getD3GraphMetadataHTML(graph, false)
    body!(w, body)
end

# plot the graph
function plot(graph::AD)
  if displayable("text/html")
    display("text/html", getD3GraphMetadataHTML(graph, true))
  else
    #try
      blink_show(graph)
    #catch
    #  error("MIME type text/html is not available.")
    #end
  end
end


text:
using Cbc, MathProgBase, SimpleGraphs, SimpleGraphAlgorithms, Test

function old_max_matching(G::SimpleGraph)
    m = NE(G)
    if m == 0
      T = vertex_type(G)
      S = Tuple{T,T}
      return Set{S}()
    end
    A = incidence(G,false)
    c = ones(m)

    X = mixintprog(-c,A,'<',1,:Int,0,1,CbcSolver())

    indices = findall(x->x!=0,round.(Int,X.sol))
    EE = elist(G)
    result = Set(EE[indices])
    # cache_save(G,:max_matching,result)
    return result
end

function old_max_indep_set(G::SimpleGraph)
    if NV(G)==0
      T = vertex_type(G)
      return Set{T}()
    end
    # if cache_check(G,:max_indep_set)
    #   return cache_recall(G,:max_indep_set)
    # end
    n = NV(G)
    A = incidence(G,false)'
    c = ones(n)

    X = mixintprog(-c,A,'<',1,:Int,0,1,CbcSolver())

    indices = findall(x->x!=0,round.(Int,X.sol))
    VV = vlist(G)
    result = Set(VV[indices])
    # cache_save(G,:max_indep_set,result)
    return result
end

function old_min_edge_cover(G::SimpleGraph)
    # if cache_check(G,:min_edge_cover)
    #   return cache_recall(G,:min_edge_cover)
    # end

    if NV(G) == 0
      T = vertex_type(G)
      S = Tuple{T,T}
      return Set{S}()
    end

    if minimum(deg(G))==0
        error("Graph has an isolated vertex; no minimum edge cover exists.")
    end
    m = NE(G)
    M = incidence(G,false)
    c = ones(m)

    X = mixintprog(c,M,'>',1,:Int,0,1,CbcSolver())

    indices = findall(x->x!=0,round.(Int,X.sol))
    EE = elist(G)
    result = Set(EE[indices])
    # cache_save(G,:min_edge_cover,result)
    return result
end

function old_min_dom_set(G::SimpleGraph)
    # if cache_check(G,:min_dom_set)
    #   return cache_recall(G,:min_dom_set)
    # end
    n = NV(G)
    if n==0
      T = vertex_type(G)
      return Set{T}()
    end

    # A = adjacency(G)+eye(Int,n)

    A = adjacency(G)
    for j=1:n
        A[j,j] = 1
    end
    c = ones(n)

    X = mixintprog(c,A,'>',1,:Int,0,1,CbcSolver())

    indices = findall(x->x!=0,round.(Int,X.sol))
    VV = vlist(G)
    result = Set(VV[indices])
    # cache_save(G,:min_dom_set,result)
    return result
end

function func_test(G::SimpleGraph, old_f::Function, new_f::Function)::Bool
    A = old_f(G)
    B = new_f(G)
    return length(A)==length(B)
end

function quad_test(G::SimpleGraph)
    func_test(G,max_matching,old_max_matching) &&
    func_test(G,max_indep_set,old_max_indep_set) &&
    func_test(G,old_min_edge_cover,min_edge_cover) &&
    func_test(G,old_min_dom_set,min_dom_set)
end

@test quad_test(BuckyBall())
@test quad_test(Paley(17))
@test quad_test(RandomTree(20))
@test quad_test(RandomRegular(20,3))


text:
include("abstract_trajectory.jl")
include("trajectory.jl")
include("vectorial_trajectory.jl")
include("circular_trajectory.jl")
include("vectorial_compact_SARTSA_buffer.jl")
include("circular_compact_SARTSA_buffer.jl")
include("episodic_compact_SARTSA_buffer.jl")
include("common.jl")
include("circular_compact_PSARTSA_buffer.jl")


text:
function test_query(q, text)
    res = run(q)

    @fact res.query --> q
    @fact isempty(res.scores) --> false
    @fact isempty(res.matches) --> false

    # Get the raw docstring.
    doc = getdoc(res)

    @fact doc --> text
end

getdoc(res) = data(docs(first(first(reverse(sort([(a, b) for (a, b) in res.scores])))[2])))

facts("Query results.") do

    context("Functions.") do

        test_query(@query(f), "f")
        test_query(@query(g), "g")

        test_query(@query(f, 1), "f")
        test_query(@query(g, 1), "g")

        test_query(@query(A.f), "A.f")
        test_query(@query(A.g), "A.g")

        test_query(@query(A.f, 1), "A.f")
        test_query(@query(A.g, 1), "A.g")

        test_query(@query(A.B.f), "A.B.f")
        test_query(@query(A.B.g), "A.B.g")

        test_query(@query(A.B.f, 1), "A.B.f")
        test_query(@query(A.B.g, 1), "A.B.g")

    end

    context("Methods.") do

        test_query(@query(f()), "f/0")
        test_query(@query(g()), "g/0")

        test_query(@query(f(1)), "f/1")
        test_query(@query(g(1)), "g/1")

        test_query(@query(f(), 1), "f/0")
        test_query(@query(g(), 1), "g/0")

        test_query(@query(f(1), 1), "f/1")
        test_query(@query(g(1), 1), "g/1")

        test_query(@query(A.f()), "A.f/0")
        test_query(@query(A.g()), "A.g/0")

        test_query(@query(A.f(1)), "A.f/1")
        test_query(@query(A.g(1)), "A.g/1")

        test_query(@query(A.f(), 1), "A.f/0")
        test_query(@query(A.g(), 1), "A.g/0")

        test_query(@query(A.f(1), 1), "A.f/1")
        test_query(@query(A.g(1), 1), "A.g/1")

        test_query(@query(A.B.f()), "A.B.f/0")
        test_query(@query(A.B.g()), "A.B.g/0")

        test_query(@query(A.B.f(1)), "A.B.f/1")
        test_query(@query(A.B.g(1)), "A.B.g/1")

        test_query(@query(A.B.f(), 1), "A.B.f/0")
        test_query(@query(A.B.g(), 1), "A.B.g/0")

        test_query(@query(A.B.f(1), 1), "A.B.f/1")
        test_query(@query(A.B.g(1), 1), "A.B.g/1")

    end

    context("Macros.") do

        test_query(@query(@m), "@m")
        test_query(@query(@m(), 1), "@m")
        test_query(@query(A.@m), "A.@m")

        test_query(@query(A.@m(), 1), "A.@m")
        test_query(@query(A.B.@m), "A.B.@m")
        test_query(@query(A.B.@m(), 1), "A.B.@m")

    end

    context("Types.") do

        test_query(@query(T), "T")
        test_query(@query(A.T), "A.T")
        test_query(@query(A.B.T), "A.B.T")

        test_query(@query(T, 1), "T")
        test_query(@query(A.T, 1), "A.T")
        test_query(@query(A.B.T, 1), "A.B.T")

        test_query(@query(S), "S")
        test_query(@query(A.S), "A.S")
        test_query(@query(A.B.S), "A.B.S")

        test_query(@query(S, 1), "S")
        test_query(@query(A.S, 1), "A.S")
        test_query(@query(A.B.S, 1), "A.B.S")

    end

    context("Globals.") do

        test_query(@query(A.K), "A.K")
        test_query(@query(A.B.K), "A.B.K")

        test_query(@query(A.K, 1), "A.K")
        test_query(@query(A.B.K, 1), "A.B.K")

    end

    context("Partial signature matching.") do

        @fact getdoc(query(f, (Any,))) --> "f/1"
        @fact getdoc(query(g, (Any,))) --> "g/1"

        @fact getdoc(query(A.f, (Any,))) --> "A.f/1"
        @fact getdoc(query(A.g, (Any,))) --> "A.g/1"

        @fact getdoc(query(A.B.f, (Any,))) --> "A.B.f/1"
        @fact getdoc(query(A.B.g, (Any,))) --> "A.B.g/1"

    end

    context("Full text search.") do

        test_query(@query("A.f/0"), "A.f/0")
        test_query(@query("A.f/1"), "A.f/1")

        test_query(@query("A.g/0"), "A.g/0")
        test_query(@query("A.g/1"), "A.g/1")

        test_query(@query("A.B.f/0"), "A.B.f/0")
        test_query(@query("A.B.f/1"), "A.B.f/1")

        test_query(@query("A.B.g/0"), "A.B.g/0")
        test_query(@query("A.B.g/1"), "A.B.g/1")

    end

    context("REPL Help.") do

        Lexicon.help("")
        Lexicon.help("fft")
        Lexicon.help("1")
        Lexicon.help("foobar")
        Lexicon.help("foobar 1")
        Lexicon.help("Lexicon")
        Lexicon.help("Docile.@document")
        Lexicon.help("Docile.Entry")
        Lexicon.help("Docile.Entry 1")
        Lexicon.help("Docile.Entry(\"...\")")

        if VERSION < v"0.4-"
            Lexicon.help("...")
            Lexicon.help("do")
        end

    end
end


text:
"""
    reflect(v, n)

Reflect direction `v` at plane with normal `n`.
"""
function reflect(v, n)
    @assert(abs(norm(n) - 1) < 1e-11, "surface normal must be normalized")
    return v - 2*dot(v, n)*n
end

"""
    refract(v, n, ni_over_nt)

Compute direction of refracted ray according to Snell's law,
or return `nothing` if no solution exists.
"""
function refract(v, n, ni_over_nt)
    @assert(abs(norm(v) - 1) < 1e-11, "input ray direction must be normalized")
    @assert(abs(norm(n) - 1) < 1e-11, "surface normal must be normalized")
    dt = dot(v, n)
    discriminant = 1 - ni_over_nt^2 * (1 - dt^2)
    if discriminant > 0
        return ni_over_nt*(v - n*dt) - sqrt(discriminant)*n
    else
        return nothing
    end
end

"""
    schlick(cosine, ref_idx)

Schlick's approximation of specular reflection coefficient.
"""
function schlick(cosine, ref_idx)
    r0 = ((1 - ref_idx) / (1 + ref_idx))^2
    return r0 + (1 - r0) * (1 - cosine)^5
end

abstract type Material end

"""
Lambertian surface (ideal diffusive reflection),
specified by albedo (reflectance) per color channel.
"""
struct Lambertian <: Material
    albedo::Vector
end

"""
    scatter(lambertian, ray, rec)

Compute scattered ray and color attenuation factors for a ray hitting a lambertian surface.
"""
function scatter(lambertian::Lambertian, ray::Ray, rec::HitRecord)
    scattered = Ray(rec.point, rec.normal + random_in_unit_sphere())
    return (scattered, lambertian.albedo)
end

"""
Metal surface, specified by albedo (reflectance) per color channel
and fuzziness factor (scales random additive permutation of reflected ray).
"""
struct Metal <: Material
    # reflectance per color channel
    albedo::Vector
    # fuzziness factor
    fuzz::Real

    Metal(albedo, fuzz) = new(albedo, min(fuzz, 1))
end

"""
    scatter(metal, ray, rec)

Compute scattered ray and color attenuation factors for a ray hitting a metal surface.
"""
function scatter(metal::Metal, ray::Ray, rec::HitRecord)
    nraydir = unitvector(ray.direction)
    reflected = reflect(nraydir, rec.normal)
    scattered = Ray(rec.point, reflected + metal.fuzz*random_in_unit_sphere())
    if dot(scattered.direction, rec.normal) > 0
        return (scattered, metal.albedo)
    else
        return (nothing, metal.albedo)
    end
end

"""
Dielectric surface, specified by ratio of the indices of refraction.
"""
struct Dielectric <: Material
    ref_idx::Real
end

"""
    scatter(dielectric, ray, rec)

Compute scattered ray and color attenuation factors for a ray hitting a dielectric surface material.
"""
function scatter(dielectric::Dielectric, ray::Ray, rec::HitRecord)
    # normalized ray direction
    nraydir = unitvector(ray.direction)

    reflected = reflect(nraydir, rec.normal)

    cosine = dot(nraydir, rec.normal)
    if cosine > 0
        refracted = refract(nraydir, -rec.normal, dielectric.ref_idx)
    else
        refracted = refract(nraydir, rec.normal, 1 / dielectric.ref_idx)
        cosine = -cosine
    end

    if refracted != nothing
        reflect_prob = schlick(cosine, dielectric.ref_idx)
    else
        reflect_prob = 1
    end

    # randomly choose between reflection or refraction
    if rand() < reflect_prob
        return (Ray(rec.point, reflected), ones(3))
    else
        return (Ray(rec.point, refracted), ones(3))
    end
end


text:
function with_temp_dir(f::Function)
    original_directory = pwd()
    tmp_dir = mktempdir()
    atexit( () -> rm(tmp_dir; force = true, recursive = true) )
    cd(tmp_dir)
    result = f()
    cd(original_directory)
    rm(tmp_dir; force = true, recursive = true)
    return result
end

function with_cloned_repo(f::Function, repo_url::AbstractString)
    original_directory = pwd()
    result = with_temp_dir() do
        run(`git clone $(repo_url) MYCLONEDREPO`)
        cd("MYCLONEDREPO")
        return f()
    end
    cd(original_directory)
    return result
end

function with_branch(f::Function,
                     repo_url::AbstractString,
                     base_branch::AbstractString,
                     head_branch::AbstractString)
    original_directory = pwd()
    result = with_cloned_repo(repo_url) do
        run(`git checkout $(base_branch)`)
        run(`git branch $(head_branch)`)
        run(`git checkout $(head_branch)`)
        return f()
    end
    cd(original_directory)
    return result
end


text:
# Copyright 2017 Raytheon BBN Technologies
#
# 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

using Qlab

include("test_digitize.jl")
include("test_common_fits.jl")
include("test_SavitzkyGolay.jl")
include("test_tomography.jl")


text:
@testset "StrictlyLessThan" begin
    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, -1 <= x <= 5, Int)
    @variable(m, -1 <= y <= 5, Int)
    @variable(m, -5 <= z <= 6, Int)
    @constraint(m, x + y + z < 4)
    optimize!(m)
    com = CS.get_inner_model(m)

    constraint = com.constraints[1]
    @test CS.is_constraint_solved(constraint, constraint.fct, constraint.set, [1, 1, 1])
    @test !CS.is_constraint_solved(constraint, constraint.fct, constraint.set, [1, 1, 2])

    constr_indices = constraint.indices
    @test !CS.still_feasible(
        com,
        constraint,
        constraint.fct,
        constraint.set,
        constr_indices[3],
        6,
    )
    @test CS.still_feasible(
        com,
        constraint,
        constraint.fct,
        constraint.set,
        constr_indices[3],
        5,
    )

    @test CS.fix!(com, com.search_space[constr_indices[2]], 0)
    @test !CS.still_feasible(
        com,
        constraint,
        constraint.fct,
        constraint.set,
        constr_indices[3],
        5,
    )

    # need to create a backtrack_vec to reverse pruning
    dummy_backtrack_obj = CS.BacktrackObj(com)
    dummy_backtrack_obj.step_nr = 1
    push!(com.backtrack_vec, dummy_backtrack_obj)
    # reverse previous fix
    CS.reverse_pruning!(com, 1)
    com.c_backtrack_idx = 1

    # now setting it to 5 should be feasible
    @test CS.still_feasible(
        com,
        constraint,
        constraint.fct,
        constraint.set,
        constr_indices[3],
        5,
    )

    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, -1 <= x <= 5, Int)
    @variable(m, -1 <= y <= 5, Int)
    @variable(m, -5 <= z <= 6, Int)
    @constraint(m, x + y + z < 4)
    optimize!(m)
    com = CS.get_inner_model(m)
    constraint = com.constraints[1]
    constr_indices = constraint.indices

    @test CS.prune_constraint!(com, constraint, constraint.fct, constraint.set)
    @test sort(CS.values(com.search_space[1])) == -1:5
    @test sort(CS.values(com.search_space[2])) == -1:5
    @test sort(CS.values(com.search_space[3])) == -5:5 # 6 not possible

    @test CS.fix!(com, com.search_space[constr_indices[3]], 5)
    CS.changed!(com, constraint, constraint.fct, constraint.set)
    @test CS.prune_constraint!(com, constraint, constraint.fct, constraint.set)
    @test CS.value(com.search_space[1]) == -1
    @test CS.value(com.search_space[2]) == -1
    @test CS.value(com.search_space[3]) == 5
    @test CS.isfixed(com.search_space[1])
    @test CS.isfixed(com.search_space[2])
    @test CS.isfixed(com.search_space[3])

    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, -1 <= x <= 5, Int)
    @variable(m, -1 <= y <= 5, Int)
    @variable(m, -5 <= z <= 6, Int)
    @constraint(m, x + y + z < 3.999)
    optimize!(m)
    com = CS.get_inner_model(m)

    constraint = com.constraints[1]
    @test constraint.rhs ≈ 3
    @test constraint.strict_rhs ≈ 3.999
    @test CS.is_constraint_solved(constraint, constraint.fct, constraint.set, [1, 1, 1])
    @test !CS.is_constraint_solved(constraint, constraint.fct, constraint.set, [1, 1, 2])

    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, -1 <= x <= 5, Int)
    @variable(m, -1 <= y <= 5, Int)
    @variable(m, -5 <= z <= 6, Int)
    @constraint(m, x + y + z < 4.1)
    optimize!(m)
    com = CS.get_inner_model(m)

    constraint = com.constraints[1]
    @test constraint.rhs ≈ 4
    @test constraint.strict_rhs ≈ 4.1
    @test CS.is_constraint_solved(constraint, constraint.fct, constraint.set, [1, 1, 1])
    @test CS.is_constraint_solved(constraint, constraint.fct, constraint.set, [1, 1, 2])

    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, -1 <= x <= 5, Int)
    @variable(m, -1 <= y <= 5, Int)
    @variable(m, -5 <= z <= 6, Int)
    @constraint(m, x + 0.9y + 1.1z < 4.1)
    optimize!(m)
    com = CS.get_inner_model(m)

    constraint = com.constraints[1]
    @test constraint.strict_rhs ≈ 4.1
    @test CS.is_constraint_solved(constraint, constraint.fct, constraint.set, [1, 1, 1])
    @test !CS.is_constraint_solved(constraint, constraint.fct, constraint.set, [1, 1, 2])
end

@testset "StrictlyGreaterThan (bridged)" begin
    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, -5 <= x <= 1, Int)
    @variable(m, -1 <= y <= 5, Int)
    @variable(m, -6 <= z <= 6, Int)
    @constraint(m, -x + y + z > 4)
    optimize!(m)
    com = CS.get_inner_model(m)

    constraint = com.constraints[1]
    @test !CS.is_constraint_solved(constraint, constraint.fct, constraint.set, [-1, 1, 1])
    @test !CS.is_constraint_solved(constraint, constraint.fct, constraint.set, [-1, 1, 2])
    @test CS.is_constraint_solved(constraint, constraint.fct, constraint.set, [-1, 2, 2])

    constr_indices = constraint.indices
    @test !CS.still_feasible(
        com,
        constraint,
        constraint.fct,
        constraint.set,
        constr_indices[3],
        -6,
    )
    @test CS.still_feasible(
        com,
        constraint,
        constraint.fct,
        constraint.set,
        constr_indices[3],
        -5,
    )

    @test CS.fix!(com, com.search_space[constr_indices[2]], 0)
    @test !CS.still_feasible(
        com,
        constraint,
        constraint.fct,
        constraint.set,
        constr_indices[3],
        -5,
    )

    # need to create a backtrack_vec to reverse pruning
    dummy_backtrack_obj = CS.BacktrackObj(com)
    dummy_backtrack_obj.step_nr = 1
    push!(com.backtrack_vec, dummy_backtrack_obj)
    # reverse previous fix
    CS.reverse_pruning!(com, 1)
    com.c_backtrack_idx = 1

    # now setting it to -5 should be feasible
    @test CS.still_feasible(
        com,
        constraint,
        constraint.fct,
        constraint.set,
        constr_indices[3],
        -5,
    )

    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, -5 <= x <= 1, Int)
    @variable(m, -1 <= y <= 5, Int)
    @variable(m, -6 <= z <= 6, Int)
    @constraint(m, -x + y + z > 4)
    optimize!(m)
    com = CS.get_inner_model(m)
    constraint = com.constraints[1]
    constr_indices = constraint.indices

    @test CS.prune_constraint!(com, constraint, constraint.fct, constraint.set)
    @test sort(CS.values(com.search_space[1])) == -5:1
    @test sort(CS.values(com.search_space[2])) == -1:5
    @test sort(CS.values(com.search_space[3])) == -5:6 # -6 not possible

    @test CS.fix!(com, com.search_space[constr_indices[3]], -5)
    CS.changed!(com, constraint, constraint.fct, constraint.set)
    @test CS.prune_constraint!(com, constraint, constraint.fct, constraint.set)
    @test CS.value(com.search_space[1]) == -5
    @test CS.value(com.search_space[2]) == 5
    @test CS.value(com.search_space[3]) == -5
    @test CS.isfixed(com.search_space[1])
    @test CS.isfixed(com.search_space[2])
    @test CS.isfixed(com.search_space[3])

    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, -1 <= x <= 5, Int)
    @variable(m, -1 <= y <= 5, Int)
    @variable(m, -5 <= z <= 6, Int)
    @constraint(m, x + y + z > 3.999)
    optimize!(m)
    com = CS.get_inner_model(m)

    constraint = com.constraints[1]
    @test constraint.rhs ≈ -4
    @test constraint.strict_rhs ≈ -3.999
    @test !CS.is_constraint_solved(constraint, constraint.fct, constraint.set, [1, 1, 1])
    @test CS.is_constraint_solved(constraint, constraint.fct, constraint.set, [1, 1, 2])
    @test CS.is_constraint_solved(constraint, constraint.fct, constraint.set, [1, 1, 3])

    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, -1 <= x <= 5, Int)
    @variable(m, -1 <= y <= 5, Int)
    @variable(m, -5 <= z <= 6, Int)
    @constraint(m, x + y + z > 4.1)
    optimize!(m)
    com = CS.get_inner_model(m)

    constraint = com.constraints[1]
    @test constraint.rhs ≈ -5
    @test constraint.strict_rhs ≈ -4.1
    @test !CS.is_constraint_solved(constraint, constraint.fct, constraint.set, [1, 1, 1])
    @test !CS.is_constraint_solved(constraint, constraint.fct, constraint.set, [1, 1, 2])

    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, -1 <= x <= 5, Int)
    @variable(m, -1 <= y <= 5, Int)
    @variable(m, -5 <= z <= 6, Int)
    @constraint(m, x + 0.9y + 1.1z > 4.1)
    optimize!(m)
    com = CS.get_inner_model(m)

    constraint = com.constraints[1]
    @test constraint.strict_rhs ≈ -4.1
    @test !CS.is_constraint_solved(constraint, constraint.fct, constraint.set, [1, 1, 1])
    @test CS.is_constraint_solved(constraint, constraint.fct, constraint.set, [1, 1, 3])
end

@testset "StrictlyGreaterThan with LP solver" begin
    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)

    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true,
                "logging" => [], "lp_optimizer" => cbc_optimizer,))
    @variable(m, -1 <= x <= 5, Int)
    @variable(m, -1 <= y <= 5, Int)
    @variable(m, -5 <= z <= 6, Int)
    @constraint(m, x + 0.9y + 1.1z > 4.1)
    optimize!(m)
    com = CS.get_inner_model(m)

    constraint = com.constraints[1]
    @test constraint.rhs ≈ -4.2
    @test constraint.strict_rhs ≈ -4.1
    @test !CS.is_constraint_solved(constraint, constraint.fct, constraint.set, [1, 1, 1])
    @test CS.is_constraint_solved(constraint, constraint.fct, constraint.set, [1, 1, 3])
end

@testset "strictly less than is_constraint_violated test" begin
    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, -5 <= x[1:2] <= 5, Int)
    @constraint(m, sum(x) < 7)
    optimize!(m)
    com = CS.get_inner_model(m)

    constraint = com.constraints[1]

    variables = com.search_space
    @test CS.fix!(com, variables[constraint.indices[1]], 5; check_feasibility = false)
    @test CS.fix!(com, variables[constraint.indices[2]], 2; check_feasibility = false)
    @test CS.is_constraint_violated(com, constraint, constraint.fct, constraint.set)

    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, -5 <= x[1:2] <= 5, Int)
    @constraint(m, sum(x) < 4)
    optimize!(m)
    com = CS.get_inner_model(m)

    constraint = com.constraints[1]

    variables = com.search_space
    @test CS.fix!(com, variables[constraint.indices[1]], 5; check_feasibility = false)
    @test !CS.is_constraint_violated(com, constraint, constraint.fct, constraint.set)
end

@testset "constraint without variables" begin
    m = Model(optimizer_with_attributes(CS.Optimizer, "logging" => []))
    @variable(m, -5 <= x[1:2] <= 5, Int)
    @constraint(m, x[2] - x[1] > -x[1] + x[2])
    optimize!(m)
    @test JuMP.termination_status(m) == MOI.INFEASIBLE

    m = Model(optimizer_with_attributes(CS.Optimizer, "logging" => []))
    @variable(m, -5 <= x[1:2] <= 5, Int)
    @constraint(m, x[1] - x[2] + (-x[1]) + x[2] < 0.001)
    optimize!(m)
    @test JuMP.termination_status(m) == MOI.OPTIMAL
end

@testset "constraint where one variable has coefficient 0" begin
    m = Model(optimizer_with_attributes(CS.Optimizer, "logging" => []))
    @variable(m, -5 <= x[1:2] <= 5, Int)
    @constraint(m, dot([0,1.2], x) > 5)
    optimize!(m)
    @test JuMP.termination_status(m) == MOI.OPTIMAL
    @test JuMP.value(x[2]) ≈ 5

    m = Model(optimizer_with_attributes(CS.Optimizer, "logging" => []))
    @variable(m, -5 <= x[1:2] <= 5, Int)
    @constraint(m, dot([0,1.2], x) > 6)
    optimize!(m)
    @test JuMP.termination_status(m) == MOI.INFEASIBLE
end


text:
# equation 6.20 fromn GPML book.
# without using squared.
function fitpdfviaSDP(  y::Vector{T},
                        Kp::Matrix{T},
                        μ::T,
                        max_iterations::Int)::Vector{T} where T <: Real

    # parse.
    N = length(y)

    # set up.
    α = Convex.Variable(N)

    constraints = collect( α[i] >= 0 for i = 1:N )
    obj_expr = ( Convex.sumsquares(Kp*α-y) + μ*Convex.quadform(α, Kp) )
    problem = Convex.minimize(obj_expr, constraints)

    # optimize.
    #@time Convex.solve!(problem, SCS.SCSSolver(max_iters = max_iterations))
    @time Convex.solve!(problem, SCS.Optimizer(max_iters = max_iterations))

    # diagnostics.
    println("status: ", problem.status)

    return vec(α.value)
end

max_stars_repo_name: RoyCCWang/PatchMixtureKriging

max_issues_repo_name: RoyCCWang/PatchMixtureKriging


text:
# # How to Make a fuel plot

# See [How to set up plots](3.0_set_up_plots.md) to get started

# ### Define a system and solve for results

# use the [Power Simulations package](http://github.com/nrel/PowerSimulations.jl.git)

# to set a [Power Systems](http://github.com/nrel/PowerSystems.jl.git) system and solve a problem.

# ```julia
# using PowerGraphics
# using PowerSystems
# const PSY = PowerSystems
# const PG = PowerGraphics
# ```
# ### Make fuel plots of the results

# ```julia
# fuel_plot(results, system)
# ```

# ![this one](plots-3/Example_Fuel_Plot_Stack.png)
# ![this one](plots-3/Example_Fuel_Plot_Bar.png)

# ### Make fuel plots with a load line

# ```julia
# fuel_plot(results, system; load = true)
# ```
# ![this one](plots-3/Fuel_Plot_with_Load_Stack.png)
# ![this one](plots-3/Fuel_Plot_with_Load_Bar.png)

# ### Make fuel plots with a curtailment line

# ```julia
# fuel_plot(results, system; curtailment = true)
# ```
# ![this one](plots-3/Fuel_Plot_with_Curtailment_Stack.png)
# ![this one](plots-3/Fuel_Plot_with_Curtailment_Bar.png)

# ### Set different colors for the plots

# ```julia
# colors = [:pink :green :blue :magenta :black]
# fuel_plot(results, system; seriescolor = colors)
# ```

# ![this one](plots-3/Example_Fuel_Plot_with_Other_Colors_Stack.png)
# ![this one](plots-3/Example_Fuel_Plot_with_Other_Colors_Bar.png)

# ### Create a Stair Plot, instead of interpolating between values

# ```julia
# fuel_plot(results, system; stair = true)
# ```

# ![this one](plots-3/Fuel_Stair.png)
# ![this one](plots-3/Fuel_Bar.png)

# ### Set a title

# ```julia
# title = "Example of a Title"
# fuel_plot(results, system; title = title)
# ```

# ![this one](plots-3/Example_of_a_Title_Stack.png)
# ![this one](plots-3/Example_of_a_Title_Bar.png)

# ### For saving the plot with the PlotlyJS backend, you can set a different format for saving

# ```julia
# fuel_plot(results, system; save = path, format = "html")
# ```

# Default format for saving is png.
# Optional formats for saving include png, html, and svg.

# ## Alternatively, set up a dictionary of the desired sorting method using `make_fuel_dictionary()`

# For example, if the generators are to be sorted by a different method than the default fuels

# ```julia
# Categories = Dict{String, NamedTuple}(
#     "coal" => NamedTuple{(:primemover, :fuel)}, (PSY.PrimeMovers.ST, PSY.ThermalFuels.COAL),
#     "oil" =>  NamedTuple{(:primemover, :fuel)}, (PSY.PrimeMovers.ST, PSY.ThermalFuels.DISTILLATE_FUEL_OIL))
# generators = make_fuel_dictionary(results, system; Categories = Categories)
# fuel_plot(results, generators)
# ```