text:
# Test general queries
@testset "General Queries" begin
    # Setup the infinite model
    optimizer = () -> MOIU.MockOptimizer(MOIU.UniversalFallback(MOIU.Model{Float64}()),
                                         eval_objective_value=false)
    m = InfiniteModel(optimizer)
    @infinite_parameter(m, 0 <= par <= 1, supports = [0, 1])
    @infinite_variable(m, inf(par))
    @hold_variable(m, g)
    @objective(m, Min, g^2)
    @constraint(m, c1, 2 * inf * g <= 1)
    tm = transcription_model(m)
    JuMP.optimize!(m)
    # setup the results
    mockoptimizer = JuMP.backend(tm).optimizer.model
    MOI.set(mockoptimizer, MOI.TerminationStatus(), MOI.OPTIMAL)
    MOI.set(mockoptimizer, MOI.RawStatusString(), "solver specific string")
    MOI.set(mockoptimizer, MOI.ResultCount(), 2)
    MOI.set(mockoptimizer, MOI.PrimalStatus(1), MOI.FEASIBLE_POINT)
    MOI.set(mockoptimizer, MOI.DualStatus(1), MOI.FEASIBLE_POINT)
    MOI.set(mockoptimizer, MOI.PrimalStatus(2), MOI.FEASIBLE_POINT)
    MOI.set(mockoptimizer, MOI.DualStatus(2), MOI.FEASIBLE_POINT)
    MOI.set(mockoptimizer, MOI.SimplexIterations(), 4)
    MOI.set(mockoptimizer, MOI.BarrierIterations(), 7)
    MOI.set(mockoptimizer, MOI.NodeCount(), 2)
    MOI.set(mockoptimizer, MOI.SolveTime(), 0.42)
    # test termination_status
    @testset "JuMP.termination_status" begin
        @test termination_status(m) == MOI.OPTIMAL
    end
    # test primal_status
    @testset "JuMP.primal_status" begin
        @test primal_status(m) == MOI.FEASIBLE_POINT
        @test primal_status(m, result = 2) == MOI.FEASIBLE_POINT
    end
    # test dual_status
    @testset "JuMP.dual_status" begin
        @test dual_status(m) == MOI.FEASIBLE_POINT
        @test dual_status(m, result = 2) == MOI.FEASIBLE_POINT
    end
    # test solve_time
    @testset "JuMP.solve_time" begin
        @test solve_time(m) == 0.42
    end
    # test raw_status
    @testset "JuMP.raw_status" begin
        @test raw_status(m) == "solver specific string"
    end
    # test simplex_iterations
    @testset "JuMP.simplex_iterations" begin
        @test simplex_iterations(m) == 4
    end
    # test barrier_iterations
    @testset "JuMP.barrier_iterations" begin
        @test barrier_iterations(m) == 7
    end
    # test node_count
    @testset "JuMP.node_count" begin
        @test node_count(m) == 2
    end
end

# Test objective queries
@testset "Objective Queries" begin
    # Setup the infinite model
    optimizer = () -> MOIU.MockOptimizer(MOIU.UniversalFallback(MOIU.Model{Float64}()),
                                         eval_objective_value=false)
    m = InfiniteModel(optimizer)
    @infinite_parameter(m, 0 <= par <= 1, supports = [0, 1])
    @infinite_variable(m, inf(par))
    @hold_variable(m, g)
    @objective(m, Min, g^2)
    @constraint(m, c1, 2 * g <= 1)
    tm = transcription_model(m)
    JuMP.optimize!(m)
    # setup the results
    mockoptimizer = JuMP.backend(tm).optimizer.model
    MOI.set(mockoptimizer, MOI.TerminationStatus(), MOI.OPTIMAL)
    MOI.set(mockoptimizer, MOI.ObjectiveValue(), -1.0)
    MOI.set(mockoptimizer, MOI.DualObjectiveValue(), -2.0)
    MOI.set(mockoptimizer, MOI.ObjectiveBound(), 2.0)
    MOI.set(mockoptimizer, MOI.ResultCount(), 1)
    MOI.set(mockoptimizer, MOI.PrimalStatus(), MOI.FEASIBLE_POINT)
    MOI.set(mockoptimizer, MOI.DualStatus(), MOI.FEASIBLE_POINT)
    MOI.set(mockoptimizer, MOI.ConstraintDual(), JuMP.optimizer_index(c1), -2.0)
    @testset "JuMP.objective_bound" begin
        @test objective_bound(m) == 2.0
    end
    # test objective_bound
    @testset "JuMP.objective_value" begin
        @test objective_value(m) == -1.
    end
    # test dual_objective_bound
    @testset "JuMP.dual_objective_value" begin
        @test dual_objective_value(m) == -2.
    end
end

# Test variable queries
@testset "Variable Queries" begin
    # Setup the infinite model
    optimizer = () -> MOIU.MockOptimizer(MOIU.UniversalFallback(MOIU.Model{Float64}()),
                                         eval_objective_value=false)
    m = InfiniteModel(optimizer)
    @infinite_parameter(m, 0 <= par <= 1, supports = [0, 1],
                        derivative_method = OrthogonalCollocation(3))
    @infinite_parameter(m, 0 <= par2 <= 2, supports = [0, 2])
    @infinite_variable(m, inf(par))
    @infinite_variable(m, inf2(par2, par))
    @finite_parameter(m, fin, 42)
    @hold_variable(m, g)
    d1 = @deriv(inf, par)
    var = build_variable(error, inf2, Dict{Int, Float64}(2 => 0))
    rv = add_variable(m, var)
    @objective(m, Min, g^2)
    @constraint(m, c1, 2 * g <= 1)
    tm = transcription_model(m)
    JuMP.optimize!(m)
    inft = transcription_variable(inf, label = All)
    gt = transcription_variable(g)
    inf2t = transcription_variable(inf2, label = All)
    d1t = transcription_variable(d1, label = All)
    rvt = transcription_variable(rv, label = All)
    # setup the optimizer
    mockoptimizer = JuMP.backend(tm).optimizer.model
    MOI.set(mockoptimizer, MOI.TerminationStatus(), MOI.OPTIMAL)
    MOI.set(mockoptimizer, MOI.ResultCount(), 1)
    MOI.set(mockoptimizer, MOI.PrimalStatus(1), MOI.FEASIBLE_POINT)
    MOI.set(mockoptimizer, MOI.VariablePrimal(1), JuMP.optimizer_index(gt), 1.0)
    MOI.set(mockoptimizer, MOI.VariablePrimal(1), JuMP.optimizer_index(inft[1]), 2.0)
    MOI.set(mockoptimizer, MOI.VariablePrimal(1), JuMP.optimizer_index(inft[2]), 1.0)
    MOI.set(mockoptimizer, MOI.VariablePrimal(1), JuMP.optimizer_index(inft[3]), 2.0)
    MOI.set(mockoptimizer, MOI.VariablePrimal(1), JuMP.optimizer_index(rvt[1]), -2.0)
    MOI.set(mockoptimizer, MOI.VariablePrimal(1), JuMP.optimizer_index(rvt[2]), -1.0)
    MOI.set(mockoptimizer, MOI.VariablePrimal(1), JuMP.optimizer_index(d1t[1]), 2.0)
    MOI.set(mockoptimizer, MOI.VariablePrimal(1), JuMP.optimizer_index(d1t[2]), 1.0)
    MOI.set(mockoptimizer, MOI.VariablePrimal(1), JuMP.optimizer_index(d1t[3]), 2.0)
    # test has_values
    @testset "JuMP.has_values" begin
        @test has_values(m)
    end
    # test map_value
    @testset "map_value" begin
        @test InfiniteOpt.map_value(inf, Val(:TransData), 1, label = All) == [2., 1., 2.]
        @test InfiniteOpt.map_value(g, Val(:TransData), 1) == 1.
        @test InfiniteOpt.map_value(rv, Val(:TransData), 1, label = All) == [-2., -1.]
    end
    # test value
    @testset "JuMP.value" begin
        @test value(inf) == [2., 2.]
        @test value(inf, label = All) == [2., 1., 2.]
        @test value(inf, label = All, ndarray = true) == [2., 1., 2.]
        @test value(d1) == [2., 2.]
        @test value(d1, label = All) == [2., 1., 2.]
        @test value(g) == 1.
        @test value(rv) == [-2., -1.]
        @test value(par) == [0., 1.]
        @test value(par, label = All) == [0., 0.5, 1.]
        @test value(fin) == 42
        @test value(fin, label = All) == 42
    end
    # test map_optimizer_index
    @testset "map_optimizer_index" begin
        @test isa(InfiniteOpt.map_optimizer_index(g, Val(:TransData)), MOI.VariableIndex)
        @test isa(InfiniteOpt.map_optimizer_index(inf, Val(:TransData)), Vector{MOI.VariableIndex})
        @test isa(InfiniteOpt.map_optimizer_index(rv, Val(:TransData)), Vector{MOI.VariableIndex})
    end
    # test optimizer_index
    @testset "JuMP.optimizer_index" begin
        @test isa(optimizer_index(g), MOI.VariableIndex)
        @test isa(optimizer_index(inf, label = InternalLabel), Vector{MOI.VariableIndex})
        @test isa(optimizer_index(rv), Vector{MOI.VariableIndex})
    end
    # test dual
    @testset "JuMP.dual" begin
        @test_throws ErrorException dual(g)
    end
end

# Test expression/measure queries
@testset "Expression/Measure Queries" begin
    # Setup the infinite model
    optimizer = () -> MOIU.MockOptimizer(MOIU.UniversalFallback(MOIU.Model{Float64}()),
                                         eval_objective_value=false)
    m = InfiniteModel(optimizer)
    @infinite_parameter(m, 0 <= par <= 1, supports = [0, 1])
    @infinite_parameter(m, 0 <= par2 <= 2, supports = [0, 2])
    @infinite_variable(m, inf(par))
    @infinite_variable(m, inf2(par2, par))
    @hold_variable(m, g)
    meas1 = support_sum(2inf, par)
    meas2 = support_sum(inf2, par2)
    @objective(m, Min, g^2)
    @constraint(m, c1, 2 * g <= 1)
    tm = transcription_model(m)
    JuMP.optimize!(m)
    inft = transcription_variable(inf)
    gt = transcription_variable(g)
    inf2t = transcription_variable(inf2)
    # setup the optimizer
    mockoptimizer = JuMP.backend(tm).optimizer.model
    MOI.set(mockoptimizer, MOI.TerminationStatus(), MOI.OPTIMAL)
    MOI.set(mockoptimizer, MOI.ResultCount(), 1)
    MOI.set(mockoptimizer, MOI.PrimalStatus(1), MOI.FEASIBLE_POINT)
    MOI.set(mockoptimizer, MOI.VariablePrimal(1), JuMP.optimizer_index(gt), 1.0)
    MOI.set(mockoptimizer, MOI.VariablePrimal(1), JuMP.optimizer_index(inft[1]), 2.0)
    MOI.set(mockoptimizer, MOI.VariablePrimal(1), JuMP.optimizer_index(inft[2]), 0.0)
    MOI.set(mockoptimizer, MOI.VariablePrimal(1), JuMP.optimizer_index(inf2t[1]), 3.0)
    MOI.set(mockoptimizer, MOI.VariablePrimal(1), JuMP.optimizer_index(inf2t[2]), -1.0)
    MOI.set(mockoptimizer, MOI.VariablePrimal(1), JuMP.optimizer_index(inf2t[3]), -3.0)
    MOI.set(mockoptimizer, MOI.VariablePrimal(1), JuMP.optimizer_index(inf2t[4]), -2.0)
    # test has_values
    @testset "JuMP.has_values" begin
        @test has_values(m)
    end
    # test map_value
    @testset "map_value" begin
        @test InfiniteOpt.map_value(meas1, Val(:TransData), 1) == 4.
        @test InfiniteOpt.map_value(meas2, Val(:TransData), 1) == [0., -3.]
        @test InfiniteOpt.map_value(3g - 1, Val(:TransData), 1) == 2.
        @test InfiniteOpt.map_value(inf^2 + g, Val(:TransData), 1) == [5., 1.]
        @test InfiniteOpt.map_value(zero(AffExpr) + 1, Val(:TransData), 1) == 1.
    end
    # test value
    @testset "JuMP.value" begin
        @test value(meas1, label = All) == 4.
        @test value(meas2, label = UserDefined) == [0., -3.]
        @test value(3g - 1) == 2.
        @test value(inf^2 + g - 2) == [3., -1.]
        @test value(inf^2 + g - 2, ndarray = true) == [3., -1.]
        @test value(zero(JuMP.GenericAffExpr{Float64, GeneralVariableRef}) - 42) == -42.
    end
    # test dual
    @testset "JuMP.dual" begin
        @test_throws ErrorException dual(meas1)
    end
end

# Test constraint queries
@testset "Constraint Queries" begin
    # Setup the infinite model
    optimizer = () -> MOIU.MockOptimizer(MOIU.UniversalFallback(MOIU.Model{Float64}()),
                                         eval_objective_value=false)
    m = InfiniteModel(optimizer)
    @infinite_parameter(m, 0 <= par <= 1, supports = [0, 1])
    @infinite_variable(m, inf(par))
    @hold_variable(m, g)
    @constraint(m, c1, g <= 0)
    @constraint(m, c2, inf >= 0)
    @objective(m, Min, g^2)
    tm = transcription_model(m)
    JuMP.optimize!(m)
    inft = transcription_variable(inf)
    gt = transcription_variable(g)
    c1t = transcription_constraint(c1)
    c2t = transcription_constraint(c2)
    # setup optimizer info
    mockoptimizer = JuMP.backend(tm).optimizer.model
    MOI.set(mockoptimizer, MOI.TerminationStatus(), MOI.OPTIMAL)
    MOI.set(mockoptimizer, MOI.ResultCount(), 1)
    MOI.set(mockoptimizer, MOI.PrimalStatus(), MOI.FEASIBLE_POINT)
    MOI.set(mockoptimizer, MOI.DualStatus(), MOI.FEASIBLE_POINT)
    MOI.set(mockoptimizer, MOI.VariablePrimal(), JuMP.optimizer_index(inft[1]), -1.0)
    MOI.set(mockoptimizer, MOI.VariablePrimal(), JuMP.optimizer_index(inft[2]), 0.0)
    MOI.set(mockoptimizer, MOI.VariablePrimal(), JuMP.optimizer_index(gt), 1.0)
    MOI.set(mockoptimizer, MOI.ConstraintDual(), JuMP.optimizer_index(c1t), -1.0)
    MOI.set(mockoptimizer, MOI.ConstraintDual(), JuMP.optimizer_index(c2t[1]), 0.0)
    MOI.set(mockoptimizer, MOI.ConstraintDual(), JuMP.optimizer_index(c2t[2]), 1.0)
    # test map_value
    @testset "map_value" begin
        @test InfiniteOpt.map_value(c1, Val(:TransData), 1) == 1.
        @test InfiniteOpt.map_value(c2, Val(:TransData), 1) == [-1., 0.]
    end
    # test value
    @testset "JuMP.value" begin
        @test value(c1) == 1.
        @test value(c2, label = UserDefined) == [-1., 0.]
        @test value(c2, label = UserDefined, ndarray = true) == [-1., 0.]
    end
    # test map_optimizer_index
    @testset "map_optimizer_index" begin
        @test isa(InfiniteOpt.map_optimizer_index(c1, Val(:TransData)), MOI.ConstraintIndex)
        @test isa(InfiniteOpt.map_optimizer_index(c2, Val(:TransData)), Vector{<:MOI.ConstraintIndex})
    end
    # test optimizer_index
    @testset "JuMP.optimizer_index" begin
        @test isa(optimizer_index(c1), MOI.ConstraintIndex)
        @test isa(optimizer_index(c2, label = All), Vector{<:MOI.ConstraintIndex})
        @test isa(optimizer_index(c2, label = All, ndarray = true), Vector{<:MOI.ConstraintIndex})
    end
    # test has_values
    @testset "JuMP.has_duals" begin
        @test has_duals(m)
    end
    # test map_dual
    @testset "map_dual" begin
        @test InfiniteOpt.map_dual(c1, Val(:TransData), 1) == -1.
        @test InfiniteOpt.map_dual(c2, Val(:TransData), 1) == [0., 1.]
    end
    # test dual
    @testset "JuMP.dual" begin
        @test dual(c1) == -1.
        @test dual(c2, label = UserDefined) == [0., 1.]
        @test dual(c2, label = UserDefined, ndarray = true) == [0., 1.]
    end
    # test map_shadow_price
    @testset "map_shadow_price" begin
        @test InfiniteOpt.map_shadow_price(c1, Val(:TransData)) == -1.
        @test InfiniteOpt.map_shadow_price(c2, Val(:TransData)) == [-0., -1.]
    end
    # test shadow_price
    @testset "JuMP.shadow_price" begin
        @test shadow_price(c1) == -1.
        @test shadow_price(c1, ndarray = true) == [-1.]
        @test shadow_price(c2, label = PublicLabel) == [-0., -1.]
    end
end

# Test LP sensitivity methods
@testset "LP Sensitivities" begin
    # Setup the infinite model
    optimizer = () -> MOIU.MockOptimizer(MOIU.UniversalFallback(MOIU.Model{Float64}()),
                                         eval_objective_value=false)
    m = InfiniteModel(optimizer)
    @infinite_parameter(m, 0 <= par <= 1, supports = [0, 1])
    @infinite_variable(m, inf(par))
    @hold_variable(m, g)
    @constraint(m, c1, g <= 0)
    @constraint(m, c2, inf >= 0)
    @objective(m, Max, 2g)
    tm = transcription_model(m)
    optimize!(m)
    inft = transcription_variable(inf)
    gt = transcription_variable(g)
    c1t = transcription_constraint(c1)
    c2t = transcription_constraint(c2)
    # setup the optimizer info
    mockoptimizer = JuMP.backend(tm).optimizer.model;
    MOI.set(mockoptimizer, MOI.TerminationStatus(), MOI.OPTIMAL)
    MOI.set(mockoptimizer, MOI.ResultCount(), 1)
    MOI.set(mockoptimizer, MOI.PrimalStatus(), MOI.FEASIBLE_POINT)
    MOI.set(mockoptimizer, MOI.DualStatus(), MOI.FEASIBLE_POINT)
    MOI.set(mockoptimizer, MOI.ObjectiveValue(), 0.0)
    MOI.set(mockoptimizer, MOI.VariablePrimal(), JuMP.optimizer_index(gt), 0.0)
    MOI.set(mockoptimizer, MOI.VariablePrimal(), JuMP.optimizer_index(inft[1]), 0.0)
    MOI.set(mockoptimizer, MOI.VariablePrimal(), JuMP.optimizer_index(inft[2]), 0.0)
    MOI.set(mockoptimizer, MOI.ConstraintDual(), JuMP.optimizer_index(c1t), -2.0)
    MOI.set(mockoptimizer, MOI.ConstraintDual(), JuMP.optimizer_index(c2t[1]), 0.0)
    MOI.set(mockoptimizer, MOI.ConstraintDual(), JuMP.optimizer_index(c2t[2]), 0.0)
    MOI.set(mockoptimizer, MOI.ConstraintBasisStatus(), JuMP.optimizer_index(c1t), MOI.NONBASIC)
    MOI.set(mockoptimizer, MOI.ConstraintBasisStatus(), JuMP.optimizer_index(c2t[1]), MOI.NONBASIC)
    MOI.set(mockoptimizer, MOI.ConstraintBasisStatus(), JuMP.optimizer_index(c2t[2]), MOI.NONBASIC)
    # test map_lp_rhs_perturbation_range
    @testset "map_lp_rhs_perturbation_range" begin
        @test InfiniteOpt.map_lp_rhs_perturbation_range(c1, Val(:TransData), 1e-8) == (-Inf, Inf)
        @test InfiniteOpt.map_lp_rhs_perturbation_range(c2, Val(:TransData), 1e-8) == [(-Inf, Inf),
                                                                                (-Inf, Inf)]
    end
    # test lp_rhs_perturbation_range
    @testset "JuMP.lp_rhs_perturbation_range" begin
        @test lp_rhs_perturbation_range(c1) == (-Inf, Inf)
        @test lp_rhs_perturbation_range(c2, label = All) == [(-Inf, Inf), (-Inf, Inf)]
        @test lp_rhs_perturbation_range(c2, ndarray = true) == [(-Inf, Inf), (-Inf, Inf)]
    end
    # test map_lp_objective_perturbation_range
    @testset "map_lp_objective_perturbation_rangee" begin
        @test InfiniteOpt.map_lp_objective_perturbation_range(g, Val(:TransData), 1e-8) == (-2.0, Inf)
        @test InfiniteOpt.map_lp_objective_perturbation_range(inf, Val(:TransData), 1e-8) == [(-Inf, 0.0),
                                                                                (-Inf, 0.0)]
    end
    # test lp_objective_perturbation_range
    @testset "JuMP.lp_objective_perturbation_range" begin
        @test lp_objective_perturbation_range(g) == (-2.0, Inf)
        @test lp_objective_perturbation_range(inf, label = UserDefined) == [(-Inf, 0.0), (-Inf, 0.0)]
        @test lp_objective_perturbation_range(inf, ndarray = true) == [(-Inf, 0.0), (-Inf, 0.0)]
    end
end


text:
"""
    origin()

Reset the current matrix, and then set the 0/0 origin to the center of the drawing
(otherwise it will stay at the top left corner, the default).

You can refer to the 0/0 point as `O`. (O = `Point(0, 0)`),
"""
function origin()
    setmatrix([1.0, 0.0, 0.0, 1.0, 0.0, 0.0])
    Cairo.translate(get_current_cr(), current_width()/2.0, current_height()/2.0)
end

function origin(x, y)
    setmatrix([1.0, 0.0, 0.0, 1.0, 0.0, 0.0])
    Cairo.translate(get_current_cr(), x, y)
end

"""
    origin(pt:Point)

Reset the current matrix, then move the `0/0` position to `pt`.
"""
function origin(pt)
    setmatrix([1.0, 0.0, 0.0, 1.0, 0.0, 0.0])
    Cairo.translate(get_current_cr(), pt.x, pt.y)
end

"""
    rescale(x, from_min, from_max, to_min=0.0, to_max=1.0)

Convert `x` from one linear scale (`from_min` to `from_max`) to another
(`to_min` to `to_max`).

The scales can also be supplied in tuple form:

```
rescale(x, (from_min, from_max), (to_min, to_max))
```

```
using Luxor
julia> rescale(15, 0, 100, 0, 1)
0.15

julia> rescale(15, (0, 100), (0, 1))
0.15

julia> rescale(pi/20, 0, 2pi, 0, 1)
0.025

julia> rescale(pi/20, (0, 2pi), (0, 1))
0.025

julia> rescale(25, 0, 1, 0, 1.609344)
40.2336

julia> rescale(15, (0, 100), (1000, 0))
850.0

```

"""
rescale(x, from_min, from_max, to_min=0.0, to_max=1.0) =
    ((x - from_min) / (from_max - from_min)) * (to_max - to_min) + to_min
rescale(x, from::NTuple{2, Number}, to::NTuple{2, Number}) =
    ((x - from[1]) / (from[2] - from[1])) * (to[2] - to[1]) + to[1]

"""
    background(color)

Fill the canvas with a single color. Returns the (red, green, blue, alpha) values.

Examples:

    background("antiquewhite")
    background(1, 0.0, 1.0)
    background(1, 0.0, 1.0, .5)

If Colors.jl is installed:

    background(RGB(0, 1, 0))
    background(RGBA(0, 1, 0))
    background(RGBA(0, 1, 0, .5))
    background(Luv(20, -20, 30))

If you don't specify a background color for a PNG drawing, the background will
be transparent. You can set a partly or completely transparent background for
PNG files by passing a color with an alpha value, such as this 'transparent
black':

    background(RGBA(0, 0, 0, 0))

or

    background(0, 0, 0, 0)

Returns a tuple `(r, g, b, a)` of the color that was used to paint the background.
"""
function background(col::Colors.Colorant)
    gsave()
    setcolor(col)
    paint()
    r, g, b, a = get_current_redvalue(), get_current_greenvalue(), get_current_bluevalue(), get_current_alpha()
    grestore()
    return (r, g, b, a)
end

function background(col::T) where T <: AbstractString
    return background(parse(Colorant, col))
end

function background(r, g, b, a = 1)
    return background(RGBA(r, g, b, a))
end

"""
    setantialias(n)

Set the current antialiasing to a value between 0 and 6:

    antialias_default  = 0, the default antialiasing for the subsystem and target device
    antialias_none     = 1, use a bilevel alpha mask
    antialias_gray     = 2, use single-color antialiasing (using shades of gray for black text on a white background, for example)
    antialias_subpixel = 3, take advantage of the order of subpixel elements on devices such as LCD panels
    antialias_fast     = 4, perform some antialiasing but prefer speed over quality
    antialias_good     = 5, balance quality against performance
    antialias_best     = 6, render at the highest quality, sacrificing speed if necessary

This affects graphics, but not text, and it doesn't apply to all types of output file.
"""
setantialias(n) = Cairo.set_antialias(get_current_cr(), n)

"""
    newpath()

Create a new path. This is Cairo's `new_path()` function.
"""
newpath() = Cairo.new_path(get_current_cr())

"""
    newsubpath()

Add a new subpath to the current path. This is Cairo's `new_sub_path()`
function. It can be used for example to make holes in shapes.
"""
newsubpath() = Cairo.new_sub_path(get_current_cr())

"""
    closepath()

Close the current path. This is Cairo's `close_path()` function.
"""
closepath() = Cairo.close_path(get_current_cr())

"""
    strokepath()

Stroke the current path with the current line width, line join, line cap, dash,
and stroke scaling settings. The current path is then cleared.
"""
strokepath() = get_current_strokescale() ? Cairo.stroke_transformed(get_current_cr()) : Cairo.stroke(get_current_cr())

"""
    fillpath()

Fill the current path according to the current settings. The current path is then cleared.
"""
fillpath() = Cairo.fill(get_current_cr())

"""
    paint()

Paint the current clip region with the current settings.
"""
paint() = Cairo.paint(get_current_cr())

"""
    strokepreserve()

Stroke the current path with current line width, line join, line cap, dash, and
stroke scaling settings, but then keep the path current.
"""
strokepreserve()    = get_current_strokescale() ? Cairo.stroke_preserve_transformed(get_current_cr()) : Cairo.stroke_preserve(get_current_cr())

"""
    fillpreserve()

Fill the current path with current settings, but then keep the path current.
"""
fillpreserve()      = Cairo.fill_preserve(get_current_cr())

"""
    fillstroke()

Fill and stroke the current path.
"""
function fillstroke()
    fillpreserve()
    strokepath()
end

"""
    do_action(action)

This is usually called by other graphics functions. Actions for graphics commands include
`:fill`, `:stroke`, `:clip`, `:fillstroke`, `:fillpreserve`, `:strokepreserve`, `:none`, and
`:path`.
"""
function do_action(action)
    if action == :fill
        fillpath()
    elseif action == :stroke
        strokepath()
    elseif action == :clip
        clip()
    elseif action == :fillstroke
        fillstroke()
    elseif action == :fillpreserve
        fillpreserve()
    elseif action == :strokepreserve
        strokepreserve()
    elseif action == :none
    end
    return true
end

"""
    clip()

Establish a new clipping region by intersecting the current clipping region with the
current path and then clearing the current path.

An existing clipping region is enforced through and after a `gsave()`-`grestore()` block, but a clipping region set inside
a `gsave()`-`grestore()` block is lost after `grestore()`. [?]
"""
clip() = Cairo.clip(get_current_cr())

"""
    clippreserve()

Establish a new clipping region by intersecting the current clipping region with the current
path, but keep the current path.
"""
clippreserve() = Cairo.clip_preserve(get_current_cr())

"""
    clipreset()

Reset the clipping region to the current drawing's extent.
"""
clipreset() = Cairo.reset_clip(get_current_cr())

"""
    setline(n)

Set the line width, in points.
"""
setline(n) = Cairo.set_line_width(get_current_cr(), n)

"""
    setlinecap(s)

Set the line ends. `s` can be "butt" or `:butt` (the
default), "square" or `:square`, or "round" or `:round`.
"""
function setlinecap(str::String="butt")
    if str == "round"
        Cairo.set_line_cap(get_current_cr(), Cairo.CAIRO_LINE_CAP_ROUND)
    elseif str == "square"
        Cairo.set_line_cap(get_current_cr(), Cairo.CAIRO_LINE_CAP_SQUARE)
    else
        Cairo.set_line_cap(get_current_cr(), Cairo.CAIRO_LINE_CAP_BUTT)
    end
end

function setlinecap(sym::Symbol)
    if sym == :round
        Cairo.set_line_cap(get_current_cr(), Cairo.CAIRO_LINE_CAP_ROUND)
    elseif sym == :square
        Cairo.set_line_cap(get_current_cr(), Cairo.CAIRO_LINE_CAP_SQUARE)
    else
        Cairo.set_line_cap(get_current_cr(), Cairo.CAIRO_LINE_CAP_BUTT)
    end
end

"""
    setlinejoin("miter")
    setlinejoin("round")
    setlinejoin("bevel")

Set the line join style, or how to render the junction of two lines when stroking.
"""
function setlinejoin(str="miter")
    if str == "round"
        Cairo.set_line_join(get_current_cr(), Cairo.CAIRO_LINE_JOIN_ROUND)
    elseif str == "bevel"
        Cairo.set_line_join(get_current_cr(), Cairo.CAIRO_LINE_JOIN_BEVEL)
    else
        Cairo.set_line_join(get_current_cr(), Cairo.CAIRO_LINE_JOIN_MITER)
    end
end

"""
    setdash("dot")

Set the dash pattern to one of: "solid", "dotted", "dot", "dotdashed",
"longdashed", "shortdashed", "dash", "dashed", "dotdotdashed",
"dotdotdotdashed".

Use `setdash(dashes::Vector)` to specify the pattern numerically.
"""
function setdash(dashing::AbstractString)
    Cairo.set_line_type(get_current_cr(), dashing)
end

"""
    setdash(dashes::Vector, offset=0.0)

Set the dash pattern to the values in `dashes`. The first number is the length of the ink, the second the gap, and so on.

The `offset` specifies an offset into the pattern at which the stroke begins. So an offset of 10 means that the stroke starts at `dashes[1] + 10` into the pattern.

Or use `setdash("dot")` etc.
"""
function setdash(dashes::Vector, offset=0.0)
    # no negative dashes
    Cairo.set_dash(get_current_cr(), abs.(Float64.(dashes)), offset)
end

"""
    setstrokescale()

Return the current stroke scaling setting.
"""
setstrokescale() = get_current_strokescale()

"""
    setstrokescale(state::Bool)

Enable/disable stroke scaling for the current drawing.
"""
setstrokescale(state::Bool) = set_current_strokescale(state)

"""
    move(pt)

Move to a point.
"""
move(x, y)      = Cairo.move_to(get_current_cr(), x, y)
move(pt)        = move(pt.x, pt.y)

"""
    rmove(pt)

Move relative to current position by the `pt`'s x and y:
"""
rmove(x, y)     = Cairo.rel_move_to(get_current_cr(), x, y)
rmove(pt)       = rmove(pt.x, pt.y)

"""
    line(pt)

Draw a line from the current position to the `pt`.
"""
line(x, y)      = Cairo.line_to(get_current_cr(), x, y)
line(pt)        = line(pt.x, pt.y)

"""
    line(pt1::Point, pt2::Point, action=:none)

Make a line between two points, `pt1` and `pt2` and do an action.
"""
function line(pt1::Point, pt2::Point, action=:none)
    move(pt1)
    line(pt2)
    do_action(action)
end

"""
    rline(pt)

Draw a line relative to the current position to the `pt`.
"""
rline(x, y)     = Cairo.rel_line_to(get_current_cr(), x, y)
rline(pt)       = rline(pt.x, pt.y)

"""
    rule(pos, theta;
        boundingbox=BoundingBox(),
        vertices=false)

Draw a straight line through `pos` at an angle `theta` from the x axis.

By default, the line spans the entire drawing, but you can supply a BoundingBox
to change the extent of the line.

```
rule(O)       # draws an x axis
rule(O, pi/2) # draws a  y axis
```

The function:

```
rule(O, pi/2, boundingbox=BoundingBox()/2)
```

draws a line that spans a bounding box half the width and height of the drawing, and returns a Set of end points. If you just want the vertices and don't want to draw anything, use `vertices=true`.
"""
function rule(pos, theta=0.0;
        boundingbox=BoundingBox(),
        vertices=false)
    #TODO interaction with clipping regions needs work

    #TODO which boundingbox is providing the default???
    bbox       = box(boundingbox, vertices=true)
    topside    = bbox[2:3]
    rightside  = bbox[3:4]
    bottomside = vcat(bbox[4], bbox[1])
    leftside   = bbox[1:2]

    # ruled line could be as long as the diagonal so add a bit extra
    r = boxdiagonal(boundingbox) + 10
    rcosa = r * cos(theta)
    rsina = r * sin(theta)
    ruledline = (pos - (rcosa, rsina), pos + (rcosa, rsina))

    interpoints = Array{Point, 1}()

    # check for intersection with top of bounding box
    flag, ip = intersectionlines(ruledline[1], ruledline[2], topside[1], topside[2], crossingonly=true)
    if flag
        if !(ip.x > topside[2].x || ip.x < topside[1].x)
            push!(interpoints, ip)
        end
    end

    # check for right intersection
    flag, ip = intersectionlines(ruledline[1], ruledline[2], rightside[1], rightside[2], crossingonly=true)
    if flag
        if !(ip.y > rightside[2].y || ip.y < rightside[1].y)
            push!(interpoints, ip)
        end
    end

    # check for bottom intersection
    flag, ip = intersectionlines(ruledline[1], ruledline[2], bottomside[1], bottomside[2], crossingonly=true)
    if flag
        if !(ip.x < bottomside[2].x || ip.x > bottomside[1].x)
            push!(interpoints, ip)
        end
    end

    # check for left intersection
    flag, ip = intersectionlines(ruledline[1], ruledline[2], leftside[1], leftside[2], crossingonly=true)
    if flag
        if !(ip.y > leftside[1].y || ip.y < leftside[2].y)
            push!(interpoints, ip)
        end
    end

    # eliminate duplicates due to rounding errors
    ruledline = unique(interpoints)

    # finally draw the line if we have two points
    if vertices == false && length(ruledline) == 2
        line(ruledline..., :stroke)
    end
    return ruledline
end

saved_colors = Tuple{Float64, Float64, Float64, Float64}[]

# I originally used simple Cairo save() but the colors/opacity
# thing I've got going didn't save/restore properly, hence the stack
"""
    gsave()

Save the current color settings on the stack.
"""
function gsave()
    Cairo.save(get_current_cr())
    r, g, b, a = (get_current_redvalue(),
                  get_current_greenvalue(),
                  get_current_bluevalue(),
                  get_current_alpha()
                 )
    push!(saved_colors, (r, g, b, a))
    return (r, g, b, a)
end

"""
    grestore()

Replace the current graphics state with the one on top of the stack.
"""
function grestore()
    Cairo.restore(get_current_cr())
    try
        (r, g, b, a) =  pop!(saved_colors)
        set_current_redvalue(r)
        set_current_greenvalue(g)
        set_current_bluevalue(b)
        set_current_alpha(a)
    catch err
         println("$err Not enough colors on the stack to restore.")
    end
end

"""
    The `layer` macro is a shortcut for `gsave()` ... `grestore()`.
"""
macro layer(a)
    quote
        gsave()
        $(esc(a))
        grestore()
    end
end

"""
    scale(x, y)

Scale workspace by `x` and `y`.

Example:

```
scale(0.2, 0.3)
```

"""
scale(sx::Real, sy::Real) = Cairo.scale(get_current_cr(), sx, sy)

"""
    scale(f)

Scale workspace by `f` in both `x` and `y`.
"""
scale(f::Real) = Cairo.scale(get_current_cr(), f, f)

"""
    rotate(a::Float64)

Rotate workspace by `a` radians clockwise (from positive x-axis to positive y-axis).
"""
rotate(a) = Cairo.rotate(get_current_cr(), a)

"""
    translate(point)
    translate(x::Real, y::Real)

Translate the workspace to `x` and `y` or to `pt`.
"""
translate(tx::Real, ty::Real)        = Cairo.translate(get_current_cr(), tx, ty)
translate(pt::Point)     = translate(pt.x, pt.y)

"""
    getpath()

Get the current path and return a CairoPath object, which is an array of `element_type` and
`points` objects. With the results you can step through and examine each entry:

```
o = getpath()
for e in o
      if e.element_type == Cairo.CAIRO_PATH_MOVE_TO
          (x, y) = e.points
          move(x, y)
      elseif e.element_type == Cairo.CAIRO_PATH_LINE_TO
          (x, y) = e.points
          # straight lines
          line(x, y)
          strokepath()
          circle(x, y, 1, :stroke)
      elseif e.element_type == Cairo.CAIRO_PATH_CURVE_TO
          (x1, y1, x2, y2, x3, y3) = e.points
          # Bezier control lines
          circle(x1, y1, 1, :stroke)
          circle(x2, y2, 1, :stroke)
          circle(x3, y3, 1, :stroke)
          move(x, y)
          curve(x1, y1, x2, y2, x3, y3)
          strokepath()
          (x, y) = (x3, y3) # update current point
      elseif e.element_type == Cairo.CAIRO_PATH_CLOSE_PATH
          closepath()
      else
          error("unknown CairoPathEntry " * repr(e.element_type))
          error("unknown CairoPathEntry " * repr(e.points))
      end
  end
```
"""
getpath()      = Cairo.convert_cairo_path_data(Cairo.copy_path(get_current_cr()))

"""
    getpathflat()

Get the current path, like `getpath()` but flattened so that there are no Bèzier curves.

Returns a CairoPath which is an array of `element_type` and `points` objects.
"""
getpathflat()  = Cairo.convert_cairo_path_data(Cairo.copy_path_flat(get_current_cr()))

"""
    rulers()

Draw and label two rulers starting at `O`, the current 0/0, and continuing out
along the current positive x and y axes.
"""
function rulers()
    @layer begin
        n = 200
        w = 20
        setopacity(0.5)
        setline(0.25)
        sethue(1.0, 0.549, 0.0) # darkorange
        #x axis
        box(O, O + (n, -w), :fillstroke)
        #y axis
        box(O - (w, 0), O + (0, n), :fillstroke)
        sethue(0.722, 0.525, 0.043) # darkgoldenrod
        setopacity(1)
        [line(Point(x, 0), Point(x, -w/4), :stroke) for x in 0:10:n]
        [line(Point(-w/4, y), Point(0, y), :stroke) for y in 0:10:n]
        [line(Point(x, 0), Point(x, -w/6), :stroke) for x in 0:5:n]
        [line(Point(-w/6, y), Point(0, y), :stroke) for y in 0:5:n]
        fontsize(2)
        [text(string(x), Point(x, -w/3), halign=:right) for x in 10:10:n]
        @layer begin
                rotate(pi/2)
                [text(string(x), Point(x, w/3), halign=:right) for x in 10:10:n]
        end
        fontsize(15)
        text("X", O + (n-w/2, w), halign=:right, valign=:middle)
        text("Y", O + (-3w/2, n-w), halign=:right, valign=:middle, angle=pi/2)
        sethue(1.0, 1.0, 1.0)
        text("X", O + (w, -w/2), halign=:right, valign=:middle)
        text("Y", O + (-w/3, w/3), halign=:right, valign=:middle, angle=pi/2)
        #center
        circle(O, 2, :strokepreserve)
        setopacity(0.5)
        fillpath()
    end
    return true
end

ext: jl


text:
using POMDPs, POMDPToolbox
using MetaGraphs
using Distributions, StatsBase
using Roadnet_MDP
include("network_library.jl")

######### MDP common definitions
samenode(s1::roadnet_pursuer_state,s2::roadnet_pursuer_state) = s1.node == s2.node && s1.pnode == s1.pnode

POMDPs.n_states(mdp::roadnet_with_pursuer) = length(POMDPs.states(mdp)) # gives acces to number of states
POMDPs.n_actions(mdp::roadnet_with_pursuer) = length(action_set(mdp)) # gives access to number of actions
POMDPs.discount(mdp::roadnet_with_pursuer) = mdp.discount # the discount factor
POMDPs.isterminal(mdp::roadnet_with_pursuer, s::roadnet_pursuer_state) = s.node in mdp.exit_nodes || s.node==s.pnode # dictates when MDP solver will terminate
POMDPs.initial_state(mdp::roadnet_with_pursuer) = roadnet_with_pursuer(1,10) #could be a distribution, but we'll fix it for now

function POMDPs.states(mdp::roadnet_with_pursuer)
    # return all possible states
    s = Vector{roadnet_pursuer_state}(0)
    for v=1:nv(mdp.road_net), p=1:nv(mdp.road_net)
        push!(s,roadnet_pursuer_state(v,p))
    end
    return s
end

function POMDPs.state_index(mdp::roadnet_with_pursuer,state::roadnet_pursuer_state)
    idx = sub2ind((nv(mdp.road_net),2),state.node,state.pnode)
    return idx
end

function POMDPs.actions(mdp::roadnet_with_pursuer)
    return action_set(mdp)
end

function POMDPs.action_index(mdp::roadnet_with_pursuer,act::Symbol)
    @assert act in action_set(mdp)

    return parse(Int,match(r"[0-9]+$",string(act)).match)
end

function POMDPs.generate_s(mdp::roadnet_with_pursuer,state::roadnet_pursuer_state,action::Symbol,rng::AbstractRNG)
    # make sure this is a POMDP MetaGraph -- i.e. one we created with the special metadata
    @assert :POMDPgraph in keys(mdp.road_net.gprops)

    # generate next UGV state
    a = action
    n = state.node
    p = state.pnode

    ## pursuer move
    ## pursuer doesn't have uncertainty in movements, there is no selected action, the pursuer state just updates according to a specified method
    pursuer_neighbors_idx = neighbors(mdp.road_net,p)
    newPursuer_node = pursuer_movement(pursuer_neighbors_idx,p,method=:random)
    #  newPursuer_node = pursuer_neighbors_idx[1]::Int #just choose first neighbor for now

    ## UGV move
    state_neighbors_idx = neighbors(mdp.road_net,n)
    num_neighbors = length(state_neighbors_idx)

    target_idx = get_target(state_neighbors_idx,a)

    UGV_probability = fill(0.0,num_neighbors)

    # make transition probability
    if target_idx == -1
        # action not valid, stay put
        return roadnet_pursuer_state(state.node,newPursuer_node)
    else
        UGV_probability[target_idx] = mdp.tprob

        # spread remaining probability to other valid neighbors
        # i.e. spread among remaining neighbors num_neighbors-1
        remaining_prob = (1.0-mdp.tprob)/(num_neighbors-1)
        left_over = collect(1:num_neighbors)

        UGV_probability[left_over[left_over.!=target_idx]] = remaining_prob
    end

    newUGV_node = sample(state_neighbors_idx,Weights(UGV_probability))

    exit = is_exit(mdp,newUGV_node)

    return roadnet_pursuer_state(newUGV_node,newPursuer_node)

end

function convert_s(::Type{V} where V<: AbstractVector{Float64}, s::roadnet_pursuer_state, mdp::roadnet_with_pursuer)
#implement this function
end

function convert_s(::Type, ::AbstractArray, ::roadnet_with_pursuer)
#implement this function
end

function POMDPs.reward(mdp::roadnet_with_pursuer,state::roadnet_pursuer_state,action::Symbol,statep::roadnet_pursuer_state)
    #  println(mdp.road_net.gprops[:exit_nodes])
    if state.node in mdp.road_net.gprops[:exit_nodes]
        # reached exit
        #  println("#####escaped######")
        return mdp.reward_vals[state.node]
    end
    if statep.node == statep.pnode
        # we were aprehended
        #  println("####aprehended#####")
        #  println(mdp.road_net.gprops[:reward_dict][:caught])
        return mdp.road_net.gprops[:reward_dict][:caught]
    end

    r = 0.
    n = length(mdp.reward_states)
    # otherwise look up the reward value
    for i = 1:n
        if samenode(statep,mdp.reward_states[i])
            # if statep is a rewards state add the reward
            r += mdp.reward_vals[i]
            #  if r == 1000.
                #  println("escaped")
            #  elseif r == 2000.
                #  println("caught")
            #  else
                #  println("moved")
            #  end
        end
    end
    return r
end

#  function POMDPs.initial_state_distribution(p::roadnet_with_pursuer)
    #  return roadnet_pursuer_state(1,3)
#  end

##### helper functions
function is_exit(g::roadnet_with_pursuer,u::Int64)
    return g.road_net.vprops[u][:state_property] == :exit
end

function pursuer_movement(n::Array{Int64},p::Int64;method::Symbol=:random)::Int64
    if method == :random
        # random walk, allow pursuer to also stay put
        # TODO, pursuer and ugv can avoid collision if they move to each other's nodes simultaneously, maybe fix this?
        # TODO, pursuer can occupy the exit, possibly fix at some point
        n_copy = copy(n)
        push!(n_copy,p)
        q = sample(n_copy)
        #  println("p: $p, n: $n, n_copy: $n_copy, q: $q")
        return q
    elseif method == :first
        # choose the first neighbor
        return n[1]
    else
        println("Method not yet supported")
    end
end
function max_d(g::roadnet_with_pursuer)
    # make sure this is a POMDP MetaGraph -- i.e. one we created with the special metadata
    @assert :POMDPgraph in keys(g.road_net.gprops)

    return g.road_net.gprops[:net_stats].maximum_degree
end

function get_target_node(g::MetaGraph,n::Int64,a::Symbol)::Int64
    # make sure this is a POMDP MetaGraph -- i.e. one we created with the special metadata
    @assert :POMDPgraph in keys(g.gprops)

    node_neighbors = neighbors(g.graph,n)
    target = get_target(neighbors(g,n),a)

    return node_neighbors[target]
end

function get_target(neighbors::Vector{Int64},a::Symbol)::Int64
    idx = parse(Int,match(r"[0-9]+$",string(a)).match)
    if idx in 1:length(neighbors)
        return idx
    else
        return -1 # this node doesn't have that many neighbors
    end
end

function action_set(g::roadnet_with_pursuer;map::Bool=false)
    # make sure this is a POMDP MetaGraph -- i.e. one we created with the special metadata
    @assert :POMDPgraph in keys(g.road_net.gprops)

    if !map
        return [Symbol("n",i) for i in 1:max_d(g)]
    else
        return collect(1:max_d(g))
    end

end

# MCTS node labels
function MCTS.node_tag(s::roadnet_pursuer_state)
    if s.exit
        return "EXIT"
    elseif s.node == s.pnode
        return "CAUGHT"
    else
        return "[$(s.node),$(s.pnode)]"
    end
end

MCTS.node_tag(a::Symbol) = "[$a]"

max_stars_repo_name: COHRINT/FaMSeC

max_issues_repo_name: COHRINT/FaMSeC

max_forks_repo_name: COHRINT/FaMSeC


text:
using ChainRulesCore
using ChainRulesTestUtils
using FiniteDifferences
using ITensors
using ITensors.NDTensors
using ITensorsGrad
using OptimKit
using Random
using Test
using Zygote

using Zygote: ZygoteRuleConfig

#
# ITensor extensions
#

#
# For ITensor compatibility with FiniteDifferences
#

function FiniteDifferences.to_vec(A::ITensor)
  # TODO: generalize to sparse tensors
  # TODO: define `itensor([1.0])` as well
  # as `itensor([1.0], ())` to help with generic code.
  function vec_to_ITensor(x)
    return isempty(inds(A)) ? ITensor(x[]) : itensor(x, inds(A))
  end
  return vec(array(A)), vec_to_ITensor
end

function FiniteDifferences.to_vec(x::Index)
  return (Bool[], _ -> x)
end

function FiniteDifferences.to_vec(x::Tuple{Vararg{Index}})
  return (Bool[], _ -> x)
end

function FiniteDifferences.to_vec(x::Pair{<:Tuple{Vararg{Index}},<:Tuple{Vararg{Index}}})
  return (Bool[], _ -> x)
end

function FiniteDifferences.rand_tangent(rng::AbstractRNG, A::ITensor)
  # TODO: generalize to sparse tensors
  return isempty(inds(A)) ? ITensor(randn(eltype(A))) : randomITensor(eltype(A), inds(A))
end

function FiniteDifferences.rand_tangent(rng::AbstractRNG, A::Tensor)
  # TODO: generalize to sparse tensors
  return randomTensor(eltype(A), inds(A))
end

function FiniteDifferences.rand_tangent(rng::AbstractRNG, x::Index)
  return NoTangent()
end

function FiniteDifferences.rand_tangent(rng::AbstractRNG, x::Tuple{Vararg{Index}})
  return NoTangent()
end

function FiniteDifferences.rand_tangent(
  rng::AbstractRNG, x::Pair{<:Tuple{Vararg{Index}},<:Tuple{Vararg{Index}}}
)
  return NoTangent()
end

#
# For ITensor compatibility with ChainRulesTestUtils
#

function ChainRulesTestUtils.test_approx(
  actual::ITensor, expected::ITensor, msg=""; kwargs...
)
  ChainRulesTestUtils.@test_msg msg isapprox(actual, expected; kwargs...)
end

function ChainRulesTestUtils.test_approx(
  actual::ITensor, expected::Number, msg=""; kwargs...
)
  ChainRulesTestUtils.@test_msg msg isapprox(actual[], expected; kwargs...)
end

function ChainRulesTestUtils.test_approx(
  actual::Number, expected::ITensor, msg=""; kwargs...
)
  ChainRulesTestUtils.@test_msg msg isapprox(actual, expected[]; kwargs...)
end

@testset "ITensorsGrad.jl" begin
  @testset "Basic rrules" begin
    i = Index(2, "i")
    A = randomITensor(i', dag(i))
    Ac = randomITensor(ComplexF64, i', dag(i))
    B = randomITensor(i', dag(i))
    C = ITensor(3.4)

    test_rrule(getindex, ITensor(3.4); check_inferred=false)
    test_rrule(getindex, A, 1, 2; check_inferred=false)
    test_rrule(*, A', A; check_inferred=false)
    test_rrule(*, 3.2, A; check_inferred=false)
    test_rrule(*, A, 4.3; check_inferred=false)
    test_rrule(+, A, B; check_inferred=false)
    test_rrule(prime, A; check_inferred=false)
    test_rrule(prime, A, 2; check_inferred=false)
    test_rrule(prime, A; fkwargs=(; tags="i"), check_inferred=false)
    test_rrule(prime, A; fkwargs=(; tags="x"), check_inferred=false)
    test_rrule(replaceprime, A, 1 => 2; check_inferred=false)
    test_rrule(swapprime, A, 0 => 1; check_inferred=false)
    test_rrule(addtags, A, "x"; check_inferred=false)
    test_rrule(addtags, A, "x"; fkwargs=(; plev=1), check_inferred=false)
    test_rrule(removetags, A, "i"; check_inferred=false)
    test_rrule(replacetags, A, "i" => "j"; check_inferred=false)
    test_rrule(
      swaptags,
      randomITensor(Index(2, "i"), Index(2, "j")),
      "i" => "j";
      check_inferred=false,
    )
    test_rrule(replaceind, A, i' => sim(i); check_inferred=false)
    test_rrule(replaceinds, A, (i, i') => (sim(i), sim(i)); check_inferred=false)
    test_rrule(swapind, A, i', i; check_inferred=false)
    test_rrule(swapinds, A, (i',), (i,); check_inferred=false)
    test_rrule(itensor, randn(2, 2), i', i; check_inferred=false)
    test_rrule(ITensor, randn(2, 2), i', i; check_inferred=false)
    test_rrule(ITensor, 2.3; check_inferred=false)
    test_rrule(dag, A; check_inferred=false)
    test_rrule(permute, A, reverse(inds(A)); check_inferred=false)

    f = x -> sin(scalar(x)^3)
    args = (C,)
    test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
    f = x -> sin(x[]^3)
    args = (C,)
    test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
    f = adjoint
    args = (A,)
    test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
    f = (x, y) -> (x * y)[1, 1]
    args = (A', A)
    test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
    f = x -> prime(x, 2)[1, 1]
    args = (A,)
    test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
    f = x -> x'[1, 1]
    args = (A,)
    test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
    f = x -> addtags(x, "x")[1, 1]
    args = (A,)
    test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
    f = x -> (x' * x)[1, 1]
    args = (A,)
    test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
    f = x -> (prime(x) * x)[1, 1]
    args = (A,)
    test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
    f = x -> ((x'' * x') * x)[1, 1]
    args = (A,)
    test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
    f = x -> (x'' * (x' * x))[1, 1]
    args = (A,)
    test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
    f = (x, y, z) -> (x * y * z)[1, 1]
    args = (A'', A', A)
    test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
    f = x -> (x'' * x' * x)[1, 1]
    args = (A,)
    test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
    f = x -> (x''' * x'' * x' * x)[1, 1]
    args = (A,)
    test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
    f = x -> (x''' * x'' * x' * x)[1, 1]
    args = (A,)
    test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
    f = (x, y) -> (x + y)[1, 1]
    args = (A, B)
    test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
    f = x -> (x + x)[1, 1]
    args = (A,)
    test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
    f = x -> (2x)[1, 1]
    args = (A,)
    test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
    f = x -> (x + 2x)[1, 1]
    args = (A,)
    test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
    f = x -> (x + 2 * mapprime(x' * x, 2 => 1))[1, 1]
    args = (A,)
    test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
    f = (x, y) -> (x * y)[]
    args = (A, δ(dag(inds(A))))
    test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
    f = x -> (x * x)[]
    args = (A,)
    test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
    f = x -> (x * δ(dag(inds(x))))[]
    args = (A,)
    test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
    f = function (x)
      y = x' * x
      tr = δ(dag(inds(y)))
      return (y * tr)[]
    end
    args = (A,)
    test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
    f = function (x)
      y = x'' * x' * x
      tr = δ(dag(inds(y)))
      return (y * tr)[]
    end
    args = (A,)
    test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
    f = x -> (x^2 * δ((i', i)))[1, 1]
    args = (6.2,)
    test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
    f = x -> (x^2 * δ(i', i))[1, 1]
    args = (5.2,)
    test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
    f = x -> itensor([x^2 x; x^3 x^4], i', i)
    args = (2.54,)
    test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
    f = x -> ITensor([x^2 x; x^3 x^4], i', i)
    args = (2.1,)
    test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
    f = x -> ITensor(x)
    args = (2.12,)
    test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
    f = function (x)
      j = Index(2)
      T = itensor([x^2 sin(x); x^2 exp(-2x)], j', dag(j))
      return real((dag(T) * T)[])
    end
    args = (2.8,)
    test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
    args = (2.8 + 3.1im,)
    test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
    f = function f(x)
      j = Index(2)
      v = itensor([exp(-3.2x), cos(2x^2)], j)
      T = itensor([x^2 sin(x); x^2 exp(-2x)], j', dag(j))
      return real((dag(v') * T * v)[])
    end
    args = (2.8,)
    test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
    args = (2.8 + 3.1im,)
    test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
    f = function (x)
      j = Index(2)
      return real((x^3 * ITensor([sin(x) exp(-2x); 3x^3 x+x^2], j', dag(j)))[1, 1])
    end
    args = (3.4 + 2.3im,)
    test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
    f = x -> prime(permute(x, reverse(inds(x))))[1, 1]
    args = (A,)
    test_rrule(ZygoteRuleConfig(), f, args...; rrule_f=rrule_via_ad, check_inferred=false)
    f = x -> prime(x; plev=1)[1, 1]
    args = (A,)
    @test_throws ErrorException f'(args...)
  end
  @testset "Energy minimization" begin
    N = 6
    s = siteinds("S=1/2", N; conserve_qns=true)
    os = OpSum()
    for n in 1:(N - 1)
      os .+= 0.5, "S+", n, "S-", n + 1
      os .+= 0.5, "S-", n, "S+", n + 1
      os .+= "Sz", n, "Sz", n + 1
    end
    Hmpo = MPO(os, s)
    ψ₀mps = randomMPS(s, n -> isodd(n) ? "↑" : "↓")
    H = prod(Hmpo)
    ψ₀ = prod(ψ₀mps)
    # The Rayleigh quotient to minimize
    function E(H::ITensor, ψ::ITensor)
      ψdag = dag(ψ)
      return (ψdag' * H * ψ)[] / (ψdag * ψ)[]
    end
    E(ψ::ITensor) = E(H, ψ)
    ∇E(ψ::ITensor) = E'(ψ)
    fg(ψ::ITensor) = (E(ψ), ∇E(ψ))
    linesearch = HagerZhangLineSearch(;
      c₁=0.1, c₂=0.9, ϵ=1e-6, θ=1 / 2, γ=2 / 3, ρ=5.0, verbosity=0
    )
    algorithm = LBFGS(3; maxiter=30, gradtol=1e-8, linesearch=linesearch)
    ψ, fψ, gψ, numfg, normgradhistory = optimize(fg, ψ₀, algorithm)
    D, _ = eigen(H; ishermitian=true)
    @test E(H, ψ) < E(H, ψ₀)
    @test E(H, ψ) ≈ minimum(D)
  end
  @testset "Energy minimization (MPS)" begin
    N = 4
    χ = 4
    s = siteinds("S=1/2", N; conserve_qns=true)
    os = OpSum()
    for n in 1:(N - 1)
      os .+= 0.5, "S+", n, "S-", n + 1
      os .+= 0.5, "S-", n, "S+", n + 1
      os .+= "Sz", n, "Sz", n + 1
    end
    Hmpo = MPO(os, s)
    ψ₀mps = randomMPS(s, n -> isodd(n) ? "↑" : "↓"; linkdims=χ)
    H = ITensors.data(Hmpo)
    ψ₀ = ITensors.data(ψ₀mps)
    # The Rayleigh quotient to minimize
    function E(H::Vector{ITensor}, ψ::Vector{ITensor})
      N = length(ψ)
      ψdag = dag.(addtags.(ψ, "bra"; tags="Link"))
      ψ′dag = prime.(ψdag)
      e = ITensor(1.0)
      for n in 1:N
        e = e * ψ′dag[n] * H[n] * ψ[n]
      end
      norm = ITensor(1.0)
      for n in 1:N
        norm = norm * ψdag[n] * ψ[n]
      end
      return e[] / norm[]
    end
    E(ψ) = E(H, ψ)
    ∇E(ψ) = E'(ψ)
    fg(ψ) = (E(ψ), ∇E(ψ))
    linesearch = HagerZhangLineSearch(; c₁=0.1, c₂=0.9, ϵ=1e-6, θ=1 / 2, γ=2 / 3, ρ=5.0)
    algorithm = LBFGS(5; maxiter=100, gradtol=1e-8, linesearch=linesearch, verbosity=0)
    ψ, fψ, gψ, numfg, normgradhistory = optimize(fg, ψ₀, algorithm)
    sweeps = Sweeps(5)
    setmaxdim!(sweeps, χ)
    fψmps, ψmps = dmrg(Hmpo, ψ₀mps, sweeps; outputlevel=0)
    time_Eψ = @elapsed E(ψ)
    time_∇Eψ = @elapsed E'(ψ)
    @test E(H, ψ) ≈ inner(ψmps, Hmpo, ψmps) / inner(ψmps, ψmps)
  end
end

max_stars_repo_name: mtfishman/ITensorsGrad

max_issues_repo_name: mtfishman/ITensorsGrad

max_forks_repo_name: mtfishman/ITensorsGrad


text:
const D3=3
const POINT=15
const UNSET = 0

function GmshDiscreteModel(mshfile)

  if !isfile(mshfile)
    error("Msh file not found: $mshfile")
  end

  gmsh.initialize()
  gmsh.option.setNumber("General.Terminal", 1)
  gmsh.open(mshfile)

  D = _setup_dim(gmsh)

  node_to_coords = _setup_node_coords(gmsh,D)

  grid, cell_to_entity = _setup_grid(gmsh,D,node_to_coords)

  dim_to_gface_to_nodes, dim_gface_to_entity = _setup_faces(gmsh,D)

  push!(dim_to_gface_to_nodes,cells(grid))
  push!(dim_gface_to_entity,cell_to_entity)

  dim_to_group_to_entities = _setup_dim_to_group_to_entities(gmsh)

  dim_to_group_to_name = _setup_dim_to_group_to_name(gmsh)

  dim_to_offset = _setup_dim_to_offset(gmsh)

  nnodes = length(node_to_coords)
  node_to_label = _setup_node_to_label(gmsh,dim_to_offset,nnodes)

  gmsh.finalize()

  graph = FullGridGraph(grid)

  facelabels = _setup_face_labels(
    graph,
    dim_to_gface_to_nodes,
    dim_gface_to_entity,
    dim_to_offset,
    dim_to_group_to_entities,
    dim_to_group_to_name,
    node_to_label)

  grids = _setup_model_grids(graph,grid,D)

  DiscreteModelFromData( grids, graph, facelabels)

end

function _setup_node_to_label(gmsh,dim_to_offset,nnodes)
  node_to_label = zeros(Int,nnodes)
  entities = gmsh.model.getEntities()
  for (dim, entity) in entities
    nodes, _, _ = gmsh.model.mesh.getNodes(dim, entity)
    node_to_label[nodes] .= (entity + dim_to_offset[dim+1])
  end
  node_to_label
end

function _setup_grid(gmsh,D,node_to_coords)

  cell_to_nodes, nminD, cell_to_extrusion, cell_to_order =
    _setup_connectivity(gmsh,D)

  grid = UnstructuredGrid(
    node_to_coords,
    cell_to_nodes.data,
    cell_to_nodes.ptrs,
    cell_to_extrusion,
    cell_to_order)

  cell_to_entity = _setup_cell_to_entity(
    gmsh,D,length(cell_to_nodes),nminD)

  (grid, cell_to_entity)

end

function _setup_faces(gmsh,D)

  dim_to_gface_to_nodes = []
  dim_gface_to_entity = []

  for d in 0:(D-1)
    face_to_nodes, nmin, _, _ = _setup_connectivity(gmsh,d)
    face_to_entity = _setup_cell_to_entity(gmsh,d,length(face_to_nodes),nmin)
    push!(dim_to_gface_to_nodes,face_to_nodes)
    push!(dim_gface_to_entity,face_to_entity)
  end

  (dim_to_gface_to_nodes, dim_gface_to_entity)

end

function _setup_dim_to_group_to_entities(gmsh)
  dim_to_group_to_entities = Vector{Vector{Int}}[]
  for d in 0:D3
    pgs = gmsh.model.getPhysicalGroups(d)
    n = length(pgs)
    group_to_entities = Vector{Int}[]
    for pg in pgs
      es = gmsh.model.getEntitiesForPhysicalGroup(pg...)
      entities = Int[]
      for e in es
        push!(entities,e)
      end
      push!(group_to_entities,entities)
    end
    push!(dim_to_group_to_entities,group_to_entities)
  end
  dim_to_group_to_entities
end

function _setup_dim_to_group_to_name(gmsh)
  u = 1
  dim_to_group_to_name = Vector{String}[]
  for d in 0:D3
    pgs = gmsh.model.getPhysicalGroups(d)
    n = length(pgs)
    names = String[]
    for pg in pgs
      _name = gmsh.model.getPhysicalName(pg...)
      if _name == ""
        name = "untitled$u"
        u += 1
      else
        name = _name
      end
      push!(names,name)
    end
    push!(dim_to_group_to_name,names)
  end
  dim_to_group_to_name
end

function _setup_face_labels(
  graph,
  dim_to_gface_to_nodes,
  dim_gface_to_entity,
  dim_to_offset,
  dim_to_group_to_entities,
  dim_to_group_to_name,
  node_to_label)

  D = ndims(graph)

  dim_to_face_to_label = [node_to_label,]

  for d in 1:D
    z = fill(UNSET,length(connections(graph,d,0)))
    push!(dim_to_face_to_label,z)
  end

  _fill_dim_to_face_to_label!(
    dim_to_face_to_label,
    graph,
    dim_to_gface_to_nodes,
    dim_gface_to_entity,
    dim_to_offset)

  tag_to_name, tag_to_groups, _ = _setup_tag_to_name(dim_to_group_to_name)

  tag_to_labels = _setup_tag_to_labels(
    tag_to_groups,dim_to_group_to_entities,dim_to_offset)

  FaceLabels(dim_to_face_to_label,tag_to_labels,tag_to_name)

end

function _fill_dim_to_face_to_label!(
  dim_to_face_to_label,
  graph,
  dim_to_gface_to_nodes,
  dim_gface_to_entity,
  dim_to_offset)

  D = length(dim_to_face_to_label)-1
  d = D
  cell_to_entity =  dim_gface_to_entity[d+1]
  cell_to_label = dim_to_face_to_label[d+1]
  offset = dim_to_offset[d+1]
  _apply_offset_for_cells!(cell_to_label,cell_to_entity,offset)

  for d in 0:(D-1)

    gface_to_nodes = dim_to_gface_to_nodes[d+1]
    face_to_nodes = connections(graph,d,0)
    node_to_faces = connections(graph,0,d)
    gface_to_face = _setup_gface_to_face(
      face_to_nodes,
      node_to_faces,
      gface_to_nodes)
    face_to_label = dim_to_face_to_label[d+1]
    gface_to_entity = dim_gface_to_entity[d+1]
    offset = dim_to_offset[d+1]
    _apply_offset_for_faces!(face_to_label,gface_to_entity,gface_to_face,offset)

  end

  for d = 0:(D-1)
    for j in (d+1):D
      dface_to_jfaces = connections(graph,d,j)
      dface_to_label = dim_to_face_to_label[d+1]
      jface_to_label = dim_to_face_to_label[j+1]
      _fix_dface_to_label!(dface_to_label,jface_to_label,dface_to_jfaces)
    end
  end

end

function _fix_dface_to_label!(dface_to_label,jface_to_label,dface_to_jfaces)

  ndfaces = length(dface_to_label)
  @assert ndfaces == length(dface_to_jfaces)

  for dface in 1:ndfaces

    dlabel = dface_to_label[dface]
    if dlabel != UNSET
      continue
    end

    jfaces = dface_to_jfaces[dface]
    for jface in jfaces
      jlabel = jface_to_label[jface]
      if jlabel != UNSET
        dface_to_label[dface] = jlabel
        break
      end
    end

  end

end

function _setup_gface_to_face(
  face_to_nodes,
  node_to_faces,
  gface_to_nodes)

  gface_to_face = find_gface_to_face(
    face_to_nodes.data,
    face_to_nodes.ptrs,
    node_to_faces.data,
    node_to_faces.ptrs,
    gface_to_nodes.data,
    gface_to_nodes.ptrs)

  gface_to_face

end

function _apply_offset_for_cells!(cell_to_label,cell_to_entity,offset)
  ncells = length(cell_to_label)
  for cell in 1:ncells
    entity = cell_to_entity[cell]
    cell_to_label[cell] = entity+offset
  end
end

function _apply_offset_for_faces!(face_to_label,gface_to_entity,gface_to_face,offset)
  ngfaces = length(gface_to_face)
  for gface in 1:ngfaces
    entity = gface_to_entity[gface]
    face = gface_to_face[gface]
    face_to_label[face] = entity+offset
  end
end

function _setup_tag_to_labels(
  tag_to_groups,dim_to_group_to_entities,dim_to_offset)

  tag_to_labels = Vector{Int}[]
  for groups in tag_to_groups

    labels = Int[]
    for group in groups
      dim, id = group
      offset = dim_to_offset[dim]
      entities = dim_to_group_to_entities[dim][id]
      for entity in entities
        label = entity + offset
        push!(labels,label)
      end
    end
    push!(tag_to_labels,labels)

  end

  tag_to_labels

end

function _setup_tag_to_name(dim_to_group_to_name)
  tag_to_name = String[]
  tag_to_groups = Vector{Tuple{Int,Int}}[]
  name_to_tag = Dict{String,Int}()
  tag = 1
  for (dim, group_to_name) in enumerate(dim_to_group_to_name)
    for (id,name) in enumerate(group_to_name)
      group = (dim,id)
      if !haskey(name_to_tag,name)
        push!(tag_to_name,name)
        groups = [group,]
        push!(tag_to_groups,groups)
        name_to_tag[name] = tag
        tag += 1
      else
        _tag = name_to_tag[name]
        groups = tag_to_groups[_tag]
        push!(groups,group)
      end
    end
  end
  (tag_to_name, tag_to_groups, name_to_tag)
end

function _setup_dim_to_offset(gmsh)
  entities = gmsh.model.getEntities()
  dim_to_nentities = zeros(Int,D3+1)
  for e in entities
    d = e[1]
    id = e[2]
    _id = dim_to_nentities[d+1]
    dim_to_nentities[d+1] = max(id,_id)
  end
  dim_to_offset = zeros(Int,D3+1)
  for d=1:D3
    dim_to_offset[d+1] += dim_to_nentities[d-1+1] + dim_to_offset[d-1+1]
  end
  dim_to_offset
end

function _setup_dim(gmsh)
  entities = gmsh.model.getEntities()
  D = -1
  for e in entities
    D = max(D,e[1])
  end
  if D == -1
    gmsh.finalize()
    error("No entities in the msh file.")
  end
  D
end

function _setup_cell_to_entity(gmsh,d,ncells,nmin)

  cell_to_entity = fill(UNSET,ncells)
  entities = gmsh.model.getEntities(d)
  for e in entities
    _, elemTags, _ = gmsh.model.mesh.getElements(e[1], e[2])
    _fill_cell_to_entity!(cell_to_entity,elemTags,e[2],nmin-1)
  end

  cell_to_entity

end

function _fill_cell_to_entity!(cell_to_entity,elemTags,entity,o)

  for i_to_cell in elemTags
    for cell in i_to_cell
      cell_to_entity[cell-o] = entity
    end
  end

end

function _setup_model_grids(graph,grid,D)

  grids = Grid[]

  order = cellorders(grid).value
  extrusionD = celltypes(grid).value
  node_to_coords = points(grid)

  for d in 0:(D-1)

    grid_d = _setup_grid_from_graph(graph,d,node_to_coords,extrusionD,order)
    push!(grids,grid_d)

  end

  push!(grids,grid)

  grids

end

function _setup_grid_from_graph(graph,d,node_to_coords,extrusionD,order)

  if !all(extrusionD .== extrusionD[1])
    gmsh.finalize()
    error("Only lines, tris, quads, tets, and hexs allowed for the moment.")
  end

  extrusion = extrusionD[1:d]
  face_to_nodes = connections(graph,d,0)
  n = length(face_to_nodes)
  face_to_extrusion = ConstantCellValue(extrusion,n)
  face_to_order = ConstantCellValue(order,n)

  UnstructuredGrid(
    node_to_coords,
    face_to_nodes.data,
    face_to_nodes.ptrs,
    face_to_extrusion,
    face_to_order)

end

function _setup_node_coords(gmsh,D)

  nodeTags, coord, parametricCoord = gmsh.model.mesh.getNodes()

  nmin = minimum(nodeTags)
  nmax = maximum(nodeTags)
  nnodes = length(nodeTags)

  if !(nmax == nnodes && nmin == 1)
    gmsh.finalize()
    error("Only consecutive node tags allowed.")
  end

  node_to_coords = zeros(Point{D,Float64},nnodes)
  _fill_node_coords!(node_to_coords,nodeTags,coord,D)

  node_to_coords

end

function _fill_node_coords!(node_to_coords,nodeTags,coord,D)

  m = zero(MVector{D,Float64})

  for node in nodeTags
    for j in 1:D
      k = (node-1)*D3 + j
      xj = coord[k]
      m[j] = xj
    end
    node_to_coords[node] = m
  end

end

function _check_cell_tags(elemTags)
  nmin::Int = minimum([minimum(t) for t in elemTags])
  nmax::Int = maximum([maximum(t) for t in elemTags])
  ncells = sum([length(t) for t in elemTags])
  if !( (nmax-nmin+1) == ncells)
    gmsh.finalize()
    error("Only consecutive elem tags allowed.")
  end
  (ncells,nmin,nmax)
end

function _setup_connectivity(gmsh,d)

  elemTypes, elemTags, nodeTags = gmsh.model.mesh.getElements(d)

  if length(elemTypes) == 0

    ncells = 0
    ndata = 0
    nmin = 1
    cell_to_nodes_prts = zeros(Int,ncells+1)
    cell_to_nodes_data = zeros(Int,ndata)
    cell_to_nodes = CellVectorFromDataAndPtrs(
      cell_to_nodes_data,cell_to_nodes_prts)

    cell_to_extrusion = nothing
    cell_to_order = nothing

    nmin = 1
    return (cell_to_nodes, nmin, cell_to_extrusion, cell_to_order)

  end

  ncells, nmin, nmax = _check_cell_tags(elemTags)

  etype_to_nlnodes = _setup_etype_to_nlnodes(elemTypes,gmsh)

  ndata = sum([length(t) for t in nodeTags])

  cell_to_nodes_prts = zeros(Int,ncells+1)
  cell_to_nodes_data = zeros(Int,ndata)

  _fill_connectivity!(
    cell_to_nodes_data,
    cell_to_nodes_prts,
    nmin-1,
    etype_to_nlnodes,
    elemTypes,
    elemTags,
    nodeTags,
    d)

  cell_to_nodes = CellVectorFromDataAndPtrs(cell_to_nodes_data,cell_to_nodes_prts)

  if !(length(elemTypes)==1)
    gmsh.finalize()
    s = """
    Only one element type per dimension allowed for the moment.
    Dimension $d has $(length(elemTypes)) different element types
    """
    error(s)
  end

  etype::Int = elemTypes[1]
  name, dim, order::Int, numv, parv = gmsh.model.mesh.getElementProperties(etype)
  extrusion = _extrussion_from_etype(Val(etype))
  if extrusion == nothing
    gmsh.finalize()
    error("Unsupported element type. elemType: $etype")
  end

  if order == 0 && etype == POINT
    order = 1
  end

  cell_to_extrusion = ConstantCellValue(extrusion,ncells)
  cell_to_order = ConstantCellValue(order,ncells)

  (cell_to_nodes, nmin, cell_to_extrusion, cell_to_order)

end

const t = TET_AXIS
const h = HEX_AXIS
_extrussion_from_etype(::Val{1}) = (h,)
_extrussion_from_etype(::Val{2}) = (t,t)
_extrussion_from_etype(::Val{3}) = (h,h)
_extrussion_from_etype(::Val{4}) = (t,t,t)
_extrussion_from_etype(::Val{5}) = (h,h,h)
_extrussion_from_etype(::Val{POINT}) = ()
_extrussion_from_etype(v) = nothing

function  _fill_connectivity!(
    cell_to_nodes_data,
    cell_to_nodes_prts,
    o,
    etype_to_nlnodes,
    elemTypes,
    elemTags,
    nodeTags,
    d)

  for (j,etype) in enumerate(elemTypes)
    nlnodes = etype_to_nlnodes[etype]
    i_to_cell = elemTags[j]
    for cell in i_to_cell
      cell_to_nodes_prts[cell+1-o] = nlnodes
    end
  end

  length_to_ptrs!(cell_to_nodes_prts)

  for (j,etype) in enumerate(elemTypes)
    nlnodes = etype_to_nlnodes[etype]
    i_to_cell = elemTags[j]
    i_lnode_to_node = nodeTags[j]
    if (nlnodes == d+1)
      _orient_simplex_connectivities!(nlnodes,i_lnode_to_node)
    end
    for (i,cell) in enumerate(i_to_cell)
      a = cell_to_nodes_prts[cell-o]-1
      for lnode in 1:nlnodes
        node = i_lnode_to_node[(i-1)*nlnodes+lnode]
        cell_to_nodes_data[a+lnode] = node
      end
    end
  end

end

function _setup_etype_to_nlnodes(elemTypes,gmsh)
  netypes = maximum(elemTypes)
  etype_to_nlnodes = zeros(Int,netypes)
  for etype in elemTypes
    name, dim, order, numv, parv =
      gmsh.model.mesh.getElementProperties(etype)
    etype_to_nlnodes[etype] = numv
  end
  etype_to_nlnodes
end

# TODO to be moved to UnstructuredGrids (begin)

using UnstructuredGrids.Kernels: max_cells_arround_vertex
using UnstructuredGrids.Kernels: _fill_cells_around_scratch!
using UnstructuredGrids.Kernels: _set_intersection!

function find_gface_to_face(
  face_to_nodes_data,
  face_to_nodes_ptrs,
  node_to_faces_data::AbstractVector{T},
  node_to_faces_ptrs,
  gface_to_nodes_data,
  gface_to_nodes_ptrs) where T

  ngfaces = length(gface_to_nodes_ptrs) - 1
  gface_to_face = zeros(T,ngfaces)
  n = max_cells_arround_vertex(node_to_faces_ptrs)
  faces_around = fill(UNSET,n)
  faces_around_scratch = fill(UNSET,n)

  _fill_gface_to_face!(
    gface_to_face,
    face_to_nodes_data,
    face_to_nodes_ptrs,
    node_to_faces_data,
    node_to_faces_ptrs,
    gface_to_nodes_data,
    gface_to_nodes_ptrs,
    faces_around,
    faces_around_scratch)

  gface_to_face

end

function  _fill_gface_to_face!(
  gface_to_face,
  face_to_nodes_data,
  face_to_nodes_ptrs,
  node_to_faces_data,
  node_to_faces_ptrs,
  gface_to_nodes_data,
  gface_to_nodes_ptrs,
  faces_around,
  faces_around_scratch)

  ngfaces = length(gface_to_nodes_ptrs) - 1

  nfaces_around = UNSET
  nfaces_around_scratch = UNSET

  for gface in 1:ngfaces

    a = gface_to_nodes_ptrs[gface]-1
    b = gface_to_nodes_ptrs[gface+1]
    nlnodes = b-(a+1)

    for lnode in 1:nlnodes
      node = gface_to_nodes_data[lnode+a]
      if lnode == 1
        nfaces_around = _fill_cells_around_scratch!(
          faces_around,
          node,
          node_to_faces_data,
          node_to_faces_ptrs)
      else
        nfaces_around_scratch = _fill_cells_around_scratch!(
          faces_around_scratch,
          node,
          node_to_faces_data,
          node_to_faces_ptrs)
        _set_intersection!(
          faces_around,faces_around_scratch,
          nfaces_around,nfaces_around_scratch)
      end
    end

    for face in faces_around
      if face != UNSET
        gface_to_face[gface] = face
        break
      end
    end

  end

end

function _orient_simplex_connectivities!(nlnodes,i_lnode_to_node)
  aux = zeros(eltype(i_lnode_to_node),nlnodes)
  offset = nlnodes-1
  for i in 1:nlnodes:length(i_lnode_to_node)
    aux = i_lnode_to_node[i:i+offset]
    sort!(aux)
    i_lnode_to_node[i:i+offset] = aux
  end
end

# TODO to be moved to UnstructuredGrids (end)


text:
# Tests for distances
function t_dbscan()

X = [randn(10,200)*2+rand(10,200) rand(10,200)+10 -5+randn(10,400)]

Test.@test try

	f1(x,y) = Distances.pairwise(Distances.Euclidean(),x,y)
	f2(dist,agg)= (p,S) -> agg([dist(p, j4pr.getobs(S,i)) for i in j4pr.nobs(S)])

	W = [	j4pr.dist(X)+j4pr.dbscan(0.5, 10, f1), # this version of dbscan expects distance matrix input
      		j4pr.dbscan(0.7, f2(Distances.euclidean, mean), min_neighbors=2, min_cluster_size=5)
	]

	Xd = j4pr.datacell(X)

	for w in W
		for dtr in [X,Xd], dts in [X,Xd]
			wt = dtr |> w
			dts |> wt
		end
	end
	true
catch
	false
end

end


text:
module Forward

include("pass.jl")

include("lib/real.jl")

end


text:
using BenchmarkTools
using PWS

system = PWS.chemotaxis_system()
algorithm = DirectMCEstimate(10_000)
dtimes = collect(range(0.0, 2.0, length=51)[2:end])
initial = generate_configuration(system)

cond_ensemble = PWS.ConditionalEnsemble(system)
marg_ensemble = PWS.MarginalEnsemble(system)

@benchmark PWS.simulate(algorithm, initial, cond_ensemble, dtimes)
@benchmark PWS.simulate(algorithm, PWS.marginal_configuration(initial), marg_ensemble, dtimes)


text:
using DSGE, ModelConstructors, Distributed
using DataFrames, OrderedCollections, Dates, Plots
GR.inline("pdf")

####################################
# To use:
# Just run in the Julia REPL
# include("decompose_forecast.jl")
####################################

use_parallel = false # set to true if you want to use parallel workers
n_workers    = 10
make_plots   = true  # set to true if you want to produce forecast decomposition plots

############
# Estimate
############

# Initialize a new instance of an An Schorfheide model,
# a canonical New Keynesian DSGE model augmented for
# Bayesian estimation.
m = AnSchorfheide()

# See "test_smc.jl" for more information on setting up SMC
# We choose settings so the estimation finishes faster,
# so the particles will not approximate the actual posterior well.
m <= Setting(:data_vintage, "200101", true, "vint", "Data vintage")
m <= Setting(:date_forecast_start, DSGE.quartertodate("2020-Q1"))
m <= Setting(:test_decomp_fcast, "true", true, "test_decomp_fcast",
             "Tag for estimation output file name indicating this file is for testing a forecast decomposition")
m <= Setting(:use_parallel_workers, use_parallel)
m <= Setting(:n_particles, 800)
m <= Setting(:n_Φ, 50) # number of stages in SMC
m <= Setting(:λ, 2.0) # the curvature of the tempering schedule
m <= Setting(:adaptive_tempering_target_smc, false) # Setting a fixed tempering schedule
m <= Setting(:reoptimize, false)
m <= Setting(:sampling_method, :SMC)
forecast_string = "new" # string to distinguish the forecast for the new model

# Load data and create fake old data
df = load_data(m)
df_old = deepcopy(df[1:end - 4, :]) # subset data into a pseudo "old" and "new" data so old data is 1 year old

# Set up our "fake" old model
m_old = deepcopy(m)
m_old <= Setting(:data_vintage, "190101", true, "vint", "Data vintage") # fake data vintage, 1 year old
m_old <= Setting(:date_forecast_start, DSGE.quartertodate("2019-Q1"))
forecast_string_old = "old" # forecast string for the old forecast

# Estimate
if use_parallel
    myprocs = addprocs(n_workers)
    @everywhere using DSGE, OrderedCollections
end

estimate(m, df_to_matrix(m, df), verbose = :low, run_csminwel = false)
estimate(m_old, df_to_matrix(m_old, df_old), verbose = :low, run_csminwel = false)

if use_parallel
    rmprocs(myprocs)
end

#####################
# Decompose forecast
#####################
# We now produce a forecast decomposition, which decomposes the difference
# in the forecasts of the old model/data and the new model/data into the effect from
# data revisions, news (e.g. new data), and parameter re-estimation.
# To save on time, we only compute the modal forecast decomposition,
# and we will not use the DSGE's indivdual shocks for the decomposition.

# Run forecast
DSGE.decompose_forecast(m, m_old, df, df_old, :mode, :none, :none, [:obs, :pseudo],
                        check = true, forecast_string_new = forecast_string,
                        forecast_string_old = forecast_string_old)
DSGE.decomposition_means(m, m_old, :mode, :none, :none, [:obs, :pseudo],
                         verbose = :none, forecast_string_new = forecast_string,
                         forecast_string_old = forecast_string_old)

if make_plots
    # Set up plots
    kwargs = Dict{Symbol, Any}()
    kwargs[:tick_size]      = 1
    kwargs[:hist_label]     = ""
    kwargs[:forecast_label] = ""
    kwargs[:start_date]     = DSGE.iterate_quarters(date_mainsample_end(m), -4)
    kwargs[:end_date]       = DSGE.iterate_quarters(date_mainsample_end(m), 8)  # 2 years out

    # Make plots
    DSGE.plot_forecast_decomposition(m, m_old, collect(keys(m.observables)),
                                     :obs, :mode, :none, :none;
                                     individual_shocks = false,
                                     forecast_string_new = forecast_string,
                                     forecast_string_old = forecast_string_old,
                                     kwargs...)
    DSGE.plot_forecast_decomposition(m, m_old, collect(keys(m.pseudo_observables)),
                                     :pseudo, :mode, :none, :none;
                                     individual_shocks = false,
                                     forecast_string_new = forecast_string,
                                     forecast_string_old = forecast_string_old,
                                     kwargs...)
end

nothing


text:
sizes = [10, 20, 30]
if Threads.nthreads() > 1
    threaded = true
else
    threaded = false
end
for T in [Float64, Float32]
    @info "" T sizes threaded
    benchmark_result = runbench(
        T;
        sizes = sizes,
        threaded = threaded,
    )
    @test benchmark_result isa BLASBenchmarksCPU.BenchmarkResult
    @test benchmark_result_type(benchmark_result) === T
    df = benchmark_result_df(benchmark_result)
    @test df isa BLASBenchmarksCPU.DataFrame
    plot_directory = mktempdir()
    BLASBenchmarksCPU.plot(
        benchmark_result;
        plot_directory = plot_directory,
        displayplot = false
    )
    BLASBenchmarksCPU.plot(
        benchmark_result;
        plot_directory = plot_directory,
        measure = :median,
        displayplot = false
    )
    BLASBenchmarksCPU.plot(
        benchmark_result;
        plot_directory = plot_directory,
        measure = :mean,
        displayplot = false
    )
    BLASBenchmarksCPU.plot(
        benchmark_result;
        plot_directory = plot_directory,
        measure = :maximum,
        displayplot = false
    )
end