text:
using VisualCircuits, Compose

begin
    notQ1 = Symbol("~Q1")
    notQ2 = Symbol("~Q2")

    # JK flip flop
    IC7476(x, y) = IC(x, y, "7476", [
        :K1     :Q1     notQ1   :GND    :K2     :Q2     notQ2   :J2;
        :CLK1   :PRE1   :CLR1   :J1     :VCC    :CLK2   :PRE2   :CLR2
    ])
    # 4x 2 input and
    IC7408(x, y) = IC(x, y, "7408", [
        :VCC    :B4     :A4     :Y4     :B3     :A3     :Y3;
        :A1     :B1     :Y1     :A2     :B2     :Y2     :GND
    ])

    width, height = 11, 7.3
    jk1 = IC7476(3.5, 4.8)
    jk2 = IC7476(8.5, 4.8)
    and1 = IC7408(5.25, 1.4)

    # edges
    le = 0
    re = width

    # guide (not displayed)
    box = compose(
        context(),
        rectangle(0, 0, width, height),
        stroke("cyan"), fill(nothing)
    )

    # Wires
    y1 = 0.2
    y2 = height - 0.2
    power = compose(
        context(),
        connect((le, y1), (re, y1)),
        connect(and1[:VCC], (and1[:VCC][1], y1)),

        connect((le, y2), (re, y2)),
        connect(jk2[:VCC], (jk2[:VCC][1], y2)),
        connect(jk1[:VCC], (jk1[:VCC][1], y2)),

        connect(jk1[:PRE1], (jk1[:PRE1][1], y2)),
        connect(jk1[:PRE2], (jk1[:PRE2][1], y2)),
        connect(jk2[:PRE1], (jk2[:PRE1][1], y2)),
        connect(jk2[:PRE2], (jk2[:PRE2][1], y2)),

        connect(jk1[:J1], (jk1[:J1][1], y2)),
        connect(jk1[:K1], (jk1[:K1][1], y1)),

        dots(
            (jk1[:VCC][1], y2), (jk2[:VCC][1], y2),
            (and1[:VCC][1], y1),
            (jk1[:PRE1][1], y2), (jk1[:PRE2][1], y2),
            (jk2[:PRE1][1], y2), (jk2[:PRE2][1], y2),
            (jk1[:J1][1], y2), (jk1[:K1][1], y1),
            color = "darkred"
        ),
        stroke("darkred")
    )

    y = 2.9
    ground = compose(
        context(),
        connect((le, y), (re, y)),
        connect(jk1[:GND], (jk1[:GND][1], y)),
        connect(jk2[:GND], (jk2[:GND][1], y)),
        connect(and1[:GND], (and1[:GND][1], y)),
        dots(
            (jk1[:GND][1], y), (jk2[:GND][1], y),
            (and1[:GND][1], y), color = "darkblue"
        ),
        stroke("darkblue")
    )

    x = 1
    y = 6.3
    clear = compose(
        context(),
        connect(jk1[:CLR1], (jk1[:CLR1][1], y)),
        connect(jk1[:CLR2], (jk1[:CLR2][1], y)),
        connect(jk2[:CLR1], (jk2[:CLR1][1], y)),
        connect(jk2[:CLR2], (jk2[:CLR2][1], y), (x, y)),
        dots(
            (jk1[:CLR1][1], y), (jk1[:CLR2][1], y),
            (jk2[:CLR1][1], y), color = "purple"
        ),
        stroke("purple")
    )

    x = 1
    y = 6.7
    clock = compose(
        context(),
        connect(jk1[:CLK1], (jk1[:CLK1][1], y)),
        connect(jk1[:CLK2], (jk1[:CLK2][1], y)),
        connect(jk2[:CLK1], (jk2[:CLK1][1], y)),
        connect(jk2[:CLK2], (jk2[:CLK2][1], y), (x, y)),
        dots(
            (jk1[:CLK1][1], y), (jk1[:CLK2][1], y),
            (jk2[:CLK1][1], y), color = "darkcyan"),
        stroke("darkcyan")
    )

    # Q1 -> and, J2, K2
    y = 3.7
    logic1 = compose(
        context(),
        connect(jk1[:Q1], (jk1[:Q1][1], y), (jk1[:J2][1], y), jk1[:J2]),
        connect(jk1[:K2], and1[:A1]),
        connect(jk1[:Q2], and1[:B1]),
        dots((jk1[:K2][1], y)),
        stroke("black")
    )

    x = 5.8
    y1 = 2.5
    y2 = 3.5
    y3 = 3.7
    y9 = 5.9
    logic2 = compose(
        context(),
        connect(and1[:Y1], (and1[:Y1][1], y1), (and1[:B2][1], y1), and1[:B2]),
        connect(and1[:B2], (x, y9), (jk2[:J1][1], y9), jk2[:J1]),
        connect((x, y3), (jk2[:K1][1], y3), jk2[:K1]),

        # connect(and1[:Y1], (and1[:Y1][1], y1), (and1[:A2][1], y1), and1[:A2]),
        # connect(and1[:Y1], (and1[:Y1][1], y2), (jk2[:K1][1], y2), jk2[:K1]),
        # connect((x, y2), (x, y9), (jk2[:J1][1], y9), jk2[:J1]),

        # connect(and1[:Y1], (and1[:Y1][1], y1), (and1[:A2][1], y1), and1[:A2]),
        # connect(and1[:A2], (and1[:A2][1], y2), (x, y2), (x, y9), (jk2[:J1][1], y9), jk2[:J1]),
        # connect((x, y3), (jk2[:K1][1], y3), jk2[:K1]),
        dots((and1[:B2][1], y1), (x, y3)),
        stroke("black")
    )

    y1 = 2.5
    y2 = 3.3
    # y3 = 3.9
    logic3 = compose(
        context(),
        connect(and1[:A2], (and1[:A2][1], y2), (jk2[:Q1][1], y2), jk2[:Q1]),
        connect(and1[:Y2], (and1[:Y2][1], y1), (jk2[:J2][1], y1), jk2[:J2]),
        connect(jk2[:K2], (jk2[:K2][1], y1)),
        dots((jk2[:K2][1], y1)),
        stroke("black")
    )

    x1 = 3.0
    x2 = 7.6
    y1 = 3.3
    y2 = 2.5
    y3 = 3.7
    y9 = 1.5
    outputs = compose(
        context(),
        connect(jk1[:Q1], (jk1[:Q1][1], y9)),
        connect((jk1[:Q2][1], y2), (x1, y2), (x1, y9)),
        connect((jk2[:Q1][1], y1), (x2, y1), (x2, y9)),
        connect(jk2[:Q2], (jk2[:Q2][1], y9)),
        dots(
            (jk1[:Q1][1], y9), (x1, y9),
            (x2, y9), (jk2[:Q2][1], y9),
            color = "orange"
        ),
        dots((jk1[:Q2][1], y2), (jk2[:Q1][1], y1), (jk1[:Q1][1], y3)),
        stroke("orange")
    )

    lables = compose(
        context(),
        text(0.9, 6.3, "CLR", hright, vcenter),
        text(0.9, 6.7, "CLK", hright, vcenter),
        text(2.2, 1.3, "OUT1", hcenter, vbottom),
        text(3.0, 1.3, "OUT2", hcenter, vbottom),
        text(7.6, 1.3, "OUT3", hcenter, vbottom),
        text(9.2, 1.3, "OUT4", hcenter, vbottom),
        fontsize(20pt), font("Helvetica-Bold"), linewidth(0mm)
    )

    circuit = compose(
        context(units = UnitBox(0, 0, width, height)),
        # box,
        render(jk1), render(jk2), render(and1),

        power,
        ground,
        clear,
        #inputs,
        clock,
        logic1,
        logic2,
        logic3,
        outputs,
        lables,

        stroke("black"), linewidth(2mm)
    )

    circuit |> SVG("counter.svg", autosize(width, height)...)
    # circuit |> PDF("counter.pdf", autosize(width, height)...)
end


text:
function create_space(config::Config)
    base_space = LightGraphs.barabasi_albert(config.n_nodes, config.n0, config.k)
    return Agents.GraphSpace(base_space)
end

function initialize_model(config::Config, space::Agents.GraphSpace)
    properties = get_model_properties(config)
    model = Agents.AgentBasedModel(Citizen, space; properties=properties, scheduler=random_activation)
    populate_space!(model, config)
    distribute_littering!(model, config)
    model.trash_aggregated = sum(model.trash_in_spaces)
    model.new_trash = deepcopy(model.trash_aggregated)
    return model
end

function get_model_properties(config::Config)
    return Dict(
        :agentcount => config.agentcount,
        :cost_matrix => [2 2; 3 1],
        :clock => 0,
        :garbage_collection => config.garbage_collection,
        :trash_in_spaces => zeros(Int64, config.n_nodes),
        :trash_aggregated => 0,
        :new_trash => 0
    )
end

function populate_space!(model::Agents.ABM, config::Config)
    ids_initial_litterers = StatsBase.sample(1:config.agentcount, config.initial_litterers, replace=false)
    for id in 1:config.agentcount
        pos = StatsBase.sample(1:config.n_nodes)
        if !(id in ids_initial_litterers)
            Agents.add_agent_pos!(Citizen(id, pos, false), model)
        else
            Agents.add_agent_pos!(Citizen(id, pos, true), model)
        end
    end
    return model
end

function distribute_littering!(model::Agents.ABM, config::Config)
    for i in 1:config.initial_trash
        space_to_trash = StatsBase.sample(1:config.n_nodes)
        model.trash_in_spaces[space_to_trash] += 1
    end
    return model
end

ext: jl


text:
@with_kw mutable struct BatchSolver <: Solver
    agent::PolicyParams
    S::AbstractSpace # State space
    max_steps::Int = 100 # Maximum number of steps per episode
    𝒟_train # training data
    param_optimizers::Dict{Any, TrainingParams} = Dict() # Training parameters for the parameters
    a_opt::TrainingParams # Training parameters for the actor
    c_opt::Union{Nothing, TrainingParams} = nothing # Training parameters for the discriminator
    target_fn = nothing # the target function for value-based methods
    target_update = (π⁻, π; kwargs...) -> polyak_average!(π⁻, π, 0.005f0) # Function for updating the target network
    𝒫::NamedTuple = (;) # Parameters of the algorithm
    log::Union{Nothing, LoggerParams} = nothing # The logging parameters
    required_columns = Symbol[] # Extra columns to sample
    epoch = 0 # Number of epochs of training
end

function POMDPs.solve(𝒮::BatchSolver, mdp)
    γ = Float32(discount(mdp))
    # Sampler for logging performance
    s = Sampler(mdp, 𝒮.agent, S=𝒮.S, max_steps=𝒮.max_steps, required_columns=𝒮.required_columns)
    isnothing(𝒮.log.sampler) && (𝒮.log.sampler = s)

    # Log initial performance
    log(𝒮.log, 𝒮.epoch)

    # Loop over the number of epochs
    infos = []
    grad_steps = 0
    for 𝒮.epoch=𝒮.epoch:𝒮.epoch + 𝒮.a_opt.epochs
        minibatch_infos = [] # stores the info from each minibatch

        # Shuffle the experience buffer
        shuffle!(𝒮.𝒟_train)

        # Call train for each minibatch
        batches = partition(1:length(𝒮.𝒟_train), 𝒮.a_opt.batch_size)
        for batch in batches
            mb = minibatch(𝒮.𝒟_train, batch)
            info = Dict()

            # Train parameters
            for (θs, p_opt) in 𝒮.param_optimizers
                train!(θs, (;kwargs...) -> p_opt.loss(𝒮.agent.π, 𝒮.𝒫, mb; kwargs...), p_opt, info=info)
            end

            # Compute target
            y = !isnothing(𝒮.target_fn) ? 𝒮.target_fn(𝒮.agent.π⁻, 𝒮.𝒫, mb, γ) : nothing

            # Optionally train the critic
            if !isnothing(𝒮.c_opt)
                train!(critic(𝒮.agent.π), (;kwargs...)->𝒮.c_opt.loss(𝒮.agent.π, 𝒮.𝒫, mb, y; kwargs...), 𝒮.c_opt, info=info)

                if !isnothing(y)
                    𝒮.target_update(𝒮.agent.π⁻, 𝒮.agent.π)
                end
            end

            # Train the actor
            train!(actor(𝒮.agent.π), (;kwargs...)->𝒮.a_opt.loss(𝒮.agent.π, 𝒮.𝒫, mb; kwargs...), 𝒮.a_opt, info=info)

            grad_steps += 1
            log(𝒮.log, grad_steps, info)

            push!(minibatch_infos, info)
        end
        push!(infos, aggregate_info(minibatch_infos))

        # Early stopping
        𝒮.a_opt.early_stopping(infos) && break
    end

    𝒮.agent.π
end


text:
import UnicodePlots, Plots
using  Printf, LoopVectorization, LinearAlgebra, SparseArrays, MAT, Base.Threads

include("CreateMeshFCFV.jl")
include("VisuFCFV.jl")
include("DiscretisationFCFV_Stokes.jl")
include("SolversFCFV_Stokes.jl")
include("EvalAnalDani_v2.jl")
include("DiscretisationFCFV_Stokes_o2.jl")
BLAS.set_num_threads(4)
# include("MarkerRoutines.jl")

#--------------------------------------------------------------------#

mutable struct Markers
    x         ::  Array{Float64,1}
    y         ::  Array{Float64,1}
    phase     ::  Array{Float64,1}
    cellx     ::  Array{Int64,1}#Vector{CartesianIndex{2}}
    celly     ::  Array{Int64,1}
    nmark     ::  Int64
    nmark_max ::  Int64
end

@views function LocateMarkers(p,dx,dy,xc,yc,xmin,xmax,ymin,ymax)
    # Find marker cell indices
    @threads for k=1:p.nmark
        if (p.x[k]<xmin || p.x[k]>xmax || p.y[k]<ymin || p.y[k]>ymax)
            p.phase[k] = -1
        end
        if p.phase[k]>=0
            dstx         = p.x[k] - xc[1]
            i            = ceil(Int, dstx / dx + 0.5)
            dsty         = p.y[k] - yc[1]
            j            = ceil(Int, dsty / dy + 0.5)
            p.cellx[k]   = i
            p.celly[k]   = j
        end
    end
end

@views function Markers2Cells2(p,phase,xc,yc,dx,dy,ncx,ncy,prop,avg)
    weight      = zeros(Float64, (ncx, ncy))
    phase_th    = [similar(phase) for _ = 1:nthreads()] # per thread
    weight_th   = [similar(weight) for _ = 1:nthreads()] # per thread
    @threads for tid=1:nthreads()
        fill!(phase_th[tid] , 0)
        fill!(weight_th[tid], 0)
    end
    chunks = Iterators.partition(1:p.nmark, p.nmark ÷ nthreads())
    @sync for chunk in chunks
        Threads.@spawn begin
            tid = threadid()
            # fill!(phase_th[tid], 0)  # DON'T
            # fill!(weight_th[tid], 0)
            for k in chunk
                if p.phase[k]>=0
                # Get the indices:
                i = p.cellx[k]
                j = p.celly[k]
                # Relative distances
                dxm = 2.0 * abs(xc[i] - p.x[k])
                dym = 2.0 * abs(yc[j] - p.y[k])
                # Increment cell counts
                area = (1.0 - dxm / dx) * (1.0 - dym / dy)
                val  =  prop[Int64(p.phase[k])]
                if avg==0 phase_th[tid][i,  j] += val       * area end
                if avg==1 phase_th[tid][i,  j] += (1.0/val) * area end
                if avg==2 phase_th[tid][i,  j] += log(val) * area end
                weight_th[tid][i, j] += area
                end
            end
        end
    end
    phase  .= reduce(+, phase_th)
    weight .= reduce(+, weight_th)
    phase ./= weight
    if avg==1
        phase .= 1.0 ./ phase
    end
    if avg==2
        phase .= exp.(phase)
    end
    return
end

function SetMarkers!( p, R )
@tturbo for k=1:p.nmark
    in         = (p.x[k]^2 + p.y[k]^2) < R^2
    p.phase[k] = (in==0)* 1.0 + (in==1)*2.0
end
end

#--------------------------------------------------------------------#

function SetUpProblem!(mesh, P, Vx, Vy, Sxx, Syy, Sxy, VxDir, VyDir, SxxNeu, SyyNeu, SxyNeu, SyxNeu, sx, sy, rc, eta, gbar)
    # Evaluate T analytic on cell faces
    etam, etac = eta[1], eta[2]
    gr  = 0;                        # Simple shear: gr=1, er=0
    er  = -1;                       # Strain rate
    @tturbo for ifac=1:mesh.nf
        x, y                           = mesh.xf[ifac], mesh.yf[ifac]
        in                             = sqrt(x^2.0 + y^2.0)<=rc
        pa                             = SolutionFields_p( etam, etac, rc, er, gr, x, y, in)
        vxa, vya                       = SolutionFields_v( etam, etac, rc, er, gr, x, y, in)
        VxDir[ifac] = vxa
        VyDir[ifac] = vya
        # Pseudo-tractions
        dvxdxa, dvxdya, dvydxa, dvydya = SolutionFields_dv(etam, etac, rc, er, gr, x, y, in)
        SxxNeu[ifac] = - pa +  etam*dvxdxa
        SyyNeu[ifac] = - pa +  etam*dvydya
        SxyNeu[ifac] =         etam*dvxdya
        SyxNeu[ifac] =         etam*dvydxa
    end
    # Evaluate T analytic on barycentres
    @tturbo for iel=1:mesh.nel
        x, y                           = mesh.xc[iel], mesh.yc[iel]
        in                             = sqrt(x^2.0 + y^2.0)<=rc
        pa                             = SolutionFields_p( etam, etac, rc, er, gr, x, y, in)
        vxa, vya                       = SolutionFields_v( etam, etac, rc, er, gr, x, y, in)
        dvxdxa, dvxdya, dvydxa, dvydya = SolutionFields_dv(etam, etac, rc, er, gr, x, y, in)
        P[iel]       = pa
        Vx[iel]      = vxa
        Vy[iel]      = vya
        etaa         = (in==0) * etam + (in==1) * etac
        Sxx[iel]     = -pa + 2.0*etaa*dvxdxa
        Syy[iel]     = -pa + 2.0*etaa*dvydya
        Sxy[iel]     = etaa*(dvxdya+dvydxa)
        sx[iel]      = 0.0
        sy[iel]      = 0.0
        out          = mesh.phase[iel] == 1.0
        mesh.ke[iel] = (out==1) * 1.0*etam + (out!=1) * 1.0*etac
    end

    # Compute jump condition
    @tturbo for iel=1:mesh.nel
        x, y                           = mesh.xc[iel], mesh.yc[iel]
        for ifac=1:mesh.nf_el
            # Face
            nodei  = mesh.e2f[iel,ifac]
            xF     = mesh.xf[nodei]
            yF     = mesh.yf[nodei]
            nodei  = mesh.e2f[iel,ifac]
            bc     = mesh.bc[nodei]
            dAi    = mesh.dA[iel,ifac]
            ni_x   = mesh.n_x[iel,ifac]
            ni_y   = mesh.n_y[iel,ifac]
            phase1 = Int64(mesh.phase[iel])
            in     = (phase1==1) * 0 + (phase1==2) * 1
            pre1                           = SolutionFields_p( etam, etac, rc, er, gr, x, y, in)
            dVxdx1, dVxdy1, dVydx1, dVydy1 = SolutionFields_dv(etam, etac, rc, er, gr, x, y, in)
            eta_face1 = (in==0) * etam + (in==1) * etac
            tL_x  = ni_x*(-pre1 + eta_face1*dVxdx1) + ni_y*eta_face1*dVxdy1
            tL_y  = ni_y*(-pre1 + eta_face1*dVydy1) + ni_x*eta_face1*dVydx1
            # From element 2
            ineigh = (mesh.e2e[iel,ifac]>0) * mesh.e2e[iel,ifac] + (mesh.e2e[iel,ifac]<1) * iel
            phase2 = Int64(mesh.phase[ineigh])
            in     = (phase2==1) * 0 + (phase2==2) * 1
            pre2                           = SolutionFields_p( etam, etac, rc, er, gr, x, y, in)
            dVxdx2, dVxdy2, dVydx2, dVydy2 = SolutionFields_dv(etam, etac, rc, er, gr, x, y, in)
            eta_face2 = (in==0) * etam + (in==1) * etac
            tR_x = ni_x*(-pre2 + eta_face2*dVxdx2) + ni_y*eta_face2*dVxdy2
            tR_y = ni_y*(-pre2 + eta_face2*dVydy2) + ni_x*eta_face2*dVydx2
            gbar[iel,ifac,1] = 0.5 * (tL_x - tR_x)
            gbar[iel,ifac,2] = 0.5 * (tL_y - tR_y)
        end
    end
    return
end

#--------------------------------------------------------------------#

function ComputeError( mesh, Vxe, Vye, Txxe, Tyye, Txye, Pe, rc, eta )
    etam, etac = eta[1], eta[2]
    gr  = 0;                        # Simple shear: gr=1, er=0
    er  = -1;                       # Strain rate
    eVx  = zeros(mesh.nel)
    eVy  = zeros(mesh.nel)
    eTxx = zeros(mesh.nel)
    eTyy = zeros(mesh.nel)
    eTxy = zeros(mesh.nel)
    eP   = zeros(mesh.nel)
    Vxa  = zeros(mesh.nel)
    Vya  = zeros(mesh.nel)
    Txxa = zeros(mesh.nel)
    Tyya = zeros(mesh.nel)
    Txya = zeros(mesh.nel)
    Pa   = zeros(mesh.nel)
    Tiia = zeros(mesh.nel)
    Va   = zeros(mesh.nel)
    eV   = zeros(mesh.nel)
    eTii = zeros(mesh.nel)
    @tturbo for iel=1:mesh.nel
        x, y                           = mesh.xc[iel], mesh.yc[iel]
        in                             = sqrt(x^2.0 + y^2.0)<=rc
        pa                             = SolutionFields_p( etam, etac, rc, er, gr, x, y, in)
        vxa, vya                       = SolutionFields_v( etam, etac, rc, er, gr, x, y, in)
        dvxdxa, dvxdya, dvydxa, dvydya = SolutionFields_dv(etam, etac, rc, er, gr, x, y, in)
        Pa[iel], Vxa[iel], Vya[iel]    = pa, vxa, vya
        Txxa[iel] = 2.0*mesh.ke[iel]*dvxdxa
        Tyya[iel] = 2.0*mesh.ke[iel]*dvydya
        Txya[iel] = 1.0*mesh.ke[iel]*(dvxdya + dvydxa)
        eVx[iel]  = Vxe[iel]  - Vxa[iel]
        eVy[iel]  = Vye[iel]  - Vya[iel]
        eTxx[iel] = Txxe[iel] - Txxa[iel]
        eTyy[iel] = Tyye[iel] - Tyya[iel]
        eTxy[iel] = Txye[iel] - Txya[iel]
        eP[iel]   = Pe[iel]   - Pa[iel]
        Va[iel]   = sqrt(Vxa[iel]^2 + Vya[iel]^2)
        Ve        = sqrt(Vxe[iel]^2 + Vye[iel]^2)
        Tiia[iel] = sqrt(0.5*(Txxa[iel]^2 + Tyya[iel]^2) + Txya[iel]^2)
        Tiie      = sqrt(0.5*(Txxe[iel]^2 + Tyye[iel]^2) + Txye[iel]^2)
        eV[iel]   = Ve   - Va[iel]
        eTii[iel] = Tiie - Tiia[iel]
    end
    errVx  = norm(eVx) /norm(Vxa)
    errVy  = norm(eVy) /norm(Vya)
    errTxx = norm(eTxx)/norm(Txxa)
    errTyy = norm(eTyy)/norm(Tyya)
    errTxy = norm(eTxy)/norm(Txya)
    errP   = norm(eP)  /norm(Pa)
    errV   = norm(eV)  /norm(Va)
    errTii = norm(eTii)/norm(Tiia)
    return errVx, errVy, errTxx, errTyy, errTxy, errP, errV, errTii, Txxa, Tyya, Txya
end

#--------------------------------------------------------------------#

function StabParam(tau, dA, Vol, mesh_type, coeff)
    # if mesh_type=="Quadrangles";        taui = coeff*tau  end
    if mesh_type=="Quadrangles";        taui = coeff*tau/dA  end
    if mesh_type=="UnstructTriangles";  taui = coeff*tau*dA  end
    return taui
end

#--------------------------------------------------------------------#

@views function main(n, mesh_type, order)

    println("\n******** FCFV STOKES ********")

    # Create sides of mesh
    xmin, xmax = -3.0, 3.0
    ymin, ymax = -3.0, 3.0
    # n          = 1
    nx, ny     = n, n
    solver     = 1
    R          = 1.0
    inclusion  = 1
    eta        = [1.0 100.0]
    # mesh_type  = "Quadrangles"
    # mesh_type  = "UnstructTriangles"
    BC         = [2; 1; 1; 1] # S E N W --- 1: Dirichlet / 2: Neumann
    o2         = order-1
    nmx        = 4            # 2 marker per cell in x
    nmy        = 4            # 2 marker per cell in y
    # Generate mesh
    if mesh_type=="Quadrangles"
        # tau  = 1.0
        # if o2==0 tau  = 5e0 end
        # if o2==1 tau  = 1e1 end
        if o2==0 tau  = 6e-2 end
        if o2==1 tau  = 6e-2 end
        mesh = MakeQuadMesh( nx, ny, xmin, xmax, ymin, ymax, inclusion, R, BC )
    elseif mesh_type=="UnstructTriangles"
        if o2==0 tau  = 1.0 end
        if o2==1 tau  = 1e4 end
        mesh = MakeTriangleMesh( nx, ny, xmin, xmax, ymin, ymax, inclusion, R, BC )
        # area = 1.0 # area factor: SETTING REPRODUCE THE RESULTS OF MATLAB CODE USING TRIANGLE
        # ninc = 29  # number of points that mesh the inclusion: SETTING REPRODUCE THE RESULTS OF MATLAB CODE USING TRIANGLE
        # mesh = MakeTriangleMesh( nx, ny, xmin, xmax, ymin, ymax, inclusion, R, BC, area, ninc )
    end
    println("Number of elements: ", mesh.nel, " --- Number of dofs: ", 2*mesh.nf+mesh.nel, " --- tau0 =  ", tau)

    # Initialise markers
    ncx, ncy  = nx, ny
    nmark0    = ncx*ncy*nmx*nmy; # total initial number of marker in grid
    dx,  dy   = (xmax-xmin)/ncx, (ymax-ymin)/ncy
    dxm, dym  = dx/nmx, dy/nmy
    xm1d      =  LinRange(xmin+dxm/2, xmax-dxm/2, ncx*nmx)
    ym1d      =  LinRange(ymin+dym/2, ymax-dym/2, ncy*nmy)
    (xmi,ymi) = ([x for x=xm1d,y=ym1d], [y for x=xm1d,y=ym1d])
    xc        =  LinRange(xmin+dx/2, xmax-dx/2, ncx)
    yc        =  LinRange(ymin+dy/2, ymax-dy/2, ncy)
    # Over allocate markers
    nmark_max = 1*nmark0;
    phm    = zeros(Float64, nmark_max)
    xm     = zeros(Float64, nmark_max)
    ym     = zeros(Float64, nmark_max)
    cellxm = zeros(Int64,   nmark_max)
    cellym = zeros(Int64,   nmark_max)
    xm[1:nmark0]     = vec(xmi)
    ym[1:nmark0]     = vec(ymi)
    phm[1:nmark0]    = zeros(Float64, size(xmi))
    cellxm[1:nmark0] = zeros(Int64,   size(xmi)) #zeros(CartesianIndex{2}, size(xm))
    cellym[1:nmark0] = zeros(Int64,   size(xmi))
    p      = Markers( xm, ym, phm, cellxm, cellym, nmark0, nmark_max )

    # Define phase
    SetMarkers!( p, R )

    # Update cell info on markers
    LocateMarkers(p,dx,dy,xc,yc,xmin,xmax,ymin,ymax)

    # Source term and BCs etc...
    Pa     = zeros(mesh.nel)
    Vxa    = zeros(mesh.nel)
    Vya    = zeros(mesh.nel)
    Sxxa   = zeros(mesh.nel)
    Syya   = zeros(mesh.nel)
    Sxya   = zeros(mesh.nel)
    sex    = zeros(mesh.nel)
    sey    = zeros(mesh.nel)
    VxDir  = zeros(mesh.nf)
    VyDir  = zeros(mesh.nf)
    SxxNeu = zeros(mesh.nf)
    SyyNeu = zeros(mesh.nf)
    SxyNeu = zeros(mesh.nf)
    SyxNeu = zeros(mesh.nf)
    gbar   = zeros(mesh.nel,mesh.nf_el,2)
    println("Model configuration :")
    @time SetUpProblem!(mesh, Pa, Vxa, Vya, Sxxa, Syya, Sxya, VxDir, VyDir, SxxNeu, SyyNeu, SxyNeu, SyxNeu, sex, sey, R, eta, gbar)

    k2d = ones(ncx,ncy)
    @time Markers2Cells2(p,k2d,xc,yc,dx,dy,ncx,ncy,eta,1)
    if mesh_type=="Quadrangles"
         mesh.ke .= k2d[:]
    end

    # Compute some mesh vectors
    println("Compute FCFV vectors:")
    # @time ae, be, ze = ComputeFCFV(mesh, sex, sey, VxDir, VyDir, SxxNeu, SyyNeu, SxyNeu, SyxNeu, tau)
    @time ae, be, be_o2, ze, pe, mei, pe, rjx, rjy = ComputeFCFV_o2(mesh, sex, sey, VxDir, VyDir, SxxNeu, SyyNeu, SxyNeu, SyxNeu, tau, o2)

    # Assemble element matrices and RHS
    println("Compute element matrices:")
    # @time Kuu_v, fu_v, Kup_v, fp = ElementAssemblyLoop(mesh, ae, be, ze, VxDir, VyDir, SxxNeu, SyyNeu, SxyNeu, SyxNeu, gbar, tau)
    @time Kuu_v, fu_v, Kup_v, fp = ElementAssemblyLoop_o2(mesh, ae, be, be_o2, ze, mei, pe, rjx, rjy, VxDir, VyDir, SxxNeu, SyyNeu, SxyNeu, SyxNeu, gbar, tau, o2)

    # Assemble triplets and sparse
    println("Assemble triplets and sparse:")
    @time Kuu, fu, Kup = CreateTripletsSparse(mesh, Kuu_v, fu_v, Kup_v)
    # display(UnicodePlots.spy(Kuu))
    # display(UnicodePlots.spy(Kup))

    # Solve for hybrid variable
    println("Linear solve:")
    @time Vxh, Vyh, Pe = StokesSolvers(mesh, Kuu, Kup, fu, fp, solver)

    # # Reconstruct element values
    println("Compute element values:")
    # @time Vxe, Vye, Txxe, Tyye, Txye = ComputeElementValues(mesh, Vxh, Vyh, Pe, ae, be, ze, VxDir, VyDir, tau)
    @time Vxe, Vye, Txxe, Tyye, Txye = ComputeElementValues_o2(mesh, Vxh, Vyh, Pe, ae, be, be_o2, ze, rjx, rjy, mei, VxDir, VyDir, tau, o2)

    # # Compute discretisation errors
    err_Vx, err_Vy, err_Txx, err_Tyy, err_Txy, err_P, err_V, err_Tii, Txxa, Tyya, Txya = ComputeError( mesh, Vxe, Vye, Txxe, Tyye, Txye, Pe, R, eta )
    @printf("Error in Vx : %2.2e\n", err_Vx )
    @printf("Error in Vy : %2.2e\n", err_Vy )
    @printf("Error in Txx: %2.2e\n", err_Txx)
    @printf("Error in Tyy: %2.2e\n", err_Tyy)
    @printf("Error in Txy: %2.2e\n", err_Txy)
    @printf("Error in P  : %2.2e\n", err_P  )

    Perr = abs.(Pa.-Pe)
    Verr = sqrt.( (Vxe.-Vxa).^2 .+ (Vye.-Vya).^2 )
    println("L_inf P error: ", maximum(Perr), " --- L_inf V error: ", maximum(Verr))

    # Visualise
    # println("Visualisation:")
    # PlotMakie(mesh, v, xmin, xmax, ymin, ymax; cmap = :viridis, min_v = minimum(v), max_v = maximum(v))
    # @time PlotMakie( mesh, Vxe, xmin, xmax, ymin, ymax, :jet1, minimum(Vxe), maximum(Vxe) )
    # @time PlotMakie( mesh, Verr, xmin, xmax, ymin, ymax, :jet1, minimum(Verr), maximum(Verr) )
    @time PlotMakie( mesh, Pe, xmin, xmax, ymin, ymax, :jet1, minimum(Pa), maximum(Pa) )
    # @time PlotMakie( mesh, Perr, xmin, xmax, ymin, ymax, :jet1, minimum(Perr), maximum(Perr) )
    # @time PlotMakie( mesh, Txxe, xmin, xmax, ymin, ymax, :jet1, -6.0, 2.0 )
    # @time PlotMakie( mesh, (mesh.ke), xmin, xmax, ymin, ymax, :jet1 )
    # @time PlotMakie( mesh, mesh.phase, xmin, xmax, ymin, ymax, :jet1)

    ndof = 2*mesh.nf+mesh.nel
    return ndof, err_Vx, err_Vy, err_Txx, err_Tyy, err_Txy, err_P, err_V, err_Tii#, maximum(Perr), maximum(Verr)
end

#################### ORDER 1
order = 1

N             = 30 .* [1; 2; 3; 4;  5; 6; 7; 8; 9; 10; 11; 12; 13; 14 ] #
println(N)
mesh_type  = "Quadrangles"
eV_quad    = zeros(size(N))
eP_quad    = zeros(size(N))
eTau_quad  = zeros(size(N))
t_quad     = zeros(size(N))
ndof_quad  = zeros(size(N))
for k=1:length(N)
    t_quad[k]    = @elapsed ndof, err_Vx, err_Vy, err_Txx, err_Tyy, err_Txy, err_P, err_V, err_Tii = main( N[k], mesh_type, order )
    eV_quad[k]   = err_V
    eP_quad[k]   = err_P
    eTau_quad[k] = err_Tii
    ndof_quad[k] = ndof
end

mesh_type  = "UnstructTriangles"
eV_tri     = zeros(size(N))
eP_tri     = zeros(size(N))
eTau_tri   = zeros(size(N))
t_tri      = zeros(size(N))
ndof_tri   = zeros(size(N))
for k=1:length(N)
    t_tri[k]     = @elapsed ndof, err_Vx, err_Vy, err_Txx, err_Tyy, err_Txy, err_P, err_V, err_Tii  = main( N[k], mesh_type, order )
    eV_tri[k]    = err_V
    eP_tri[k]    = err_P
    eTau_tri[k]  = err_Tii
    ndof_tri[k]  = ndof
end

#################### ORDER 2
order = 2

mesh_type     = "Quadrangles"
eV_quad_o2    = zeros(size(N))
eP_quad_o2    = zeros(size(N))
eTau_quad_o2  = zeros(size(N))
t_quad_o2     = zeros(size(N))
ndof_quad_o2  = zeros(size(N))
for k=1:length(N)
    t_quad_o2[k]    = @elapsed ndof, err_Vx, err_Vy, err_Txx, err_Tyy, err_Txy, err_P, err_V, err_Tii = main( N[k], mesh_type, order )
    eV_quad_o2[k]   = err_V
    eP_quad_o2[k]   = err_P
    eTau_quad_o2[k] = err_Tii
    ndof_quad_o2[k] = ndof
end

mesh_type     = "UnstructTriangles"
eV_tri_o2     = zeros(size(N))
eP_tri_o2     = zeros(size(N))
eTau_tri_o2   = zeros(size(N))
t_tri_o2      = zeros(size(N))
ndof_tri_o2   = zeros(size(N))
for k=1:length(N)
    t_tri_o2[k]     = @elapsed ndof, err_Vx, err_Vy, err_Txx, err_Tyy, err_Txy, err_P, err_V, err_Tii  = main( N[k], mesh_type, order )
    eV_tri_o2[k]    = err_V
    eP_tri_o2[k]    = err_P
    eTau_tri_o2[k]  = err_Tii
    ndof_tri_o2[k]  = ndof
end

#######################################

p = Plots.plot(  log10.(1.0 ./ N) , log10.(eV_quad),   markershape=:rect,      color=:blue,                         label="Quads V O1"                          )
p = Plots.plot!( log10.(1.0 ./ N) , log10.(eP_quad),   markershape=:rect,      color=:blue,      linestyle = :dot,  label="Quads P O1"                          )
p = Plots.plot!( log10.(1.0 ./ N) , log10.(eTau_quad), markershape=:rect,      color=:blue,      linestyle = :dash, label="Quads Tau O1"                        )
p = Plots.plot!( log10.(1.0 ./ N) , log10.(eV_tri),    markershape=:dtriangle, color=:blue,                         label="Triangles V O1"                     )
p = Plots.plot!( log10.(1.0 ./ N) , log10.(eP_tri),    markershape=:dtriangle, color=:blue,      linestyle = :dot,  label="Triangles P O1"                     )
p = Plots.plot!( log10.(1.0 ./ N) , log10.(eTau_tri),  markershape=:dtriangle, color=:blue,      linestyle = :dash, label="Triangles Tau O1")#, legend=:bottomright, xlabel = "log_10(h_x)", ylabel = "log_10(err_T)" )

p = Plots.plot!(  log10.(1.0 ./ N) , log10.(eV_quad_o2),   markershape=:rect,      color=:red,                         label="Quads V O2"                          )
p = Plots.plot!( log10.(1.0 ./ N) , log10.(eP_quad_o2),   markershape=:rect,      color=:red,      linestyle = :dot,  label="Quads P O2"                          )
p = Plots.plot!( log10.(1.0 ./ N) , log10.(eTau_quad_o2), markershape=:rect,      color=:red,      linestyle = :dash, label="Quads Tau O2"                        )
p = Plots.plot!( log10.(1.0 ./ N) , log10.(eV_tri_o2),    markershape=:dtriangle, color=:red,                         label="Triangles V O2"                     )
p = Plots.plot!( log10.(1.0 ./ N) , log10.(eP_tri_o2),    markershape=:dtriangle, color=:red,      linestyle = :dot,  label="Triangles P O2"                     )
p = Plots.plot!( log10.(1.0 ./ N) , log10.(eTau_tri_o2),  markershape=:dtriangle, color=:red,      linestyle = :dash, label="Triangles Tau O2", legend=:outertopright, xlabel = "log_10(h_x)", ylabel = "log_10(err_T)" )
order1 = [2e-4, 1e-4]
order2 = [4e-4, 1e-4]
n      = [35, 70]
p = Plots.plot!( log10.(1.0 ./ n) , log10.(order1), color=:black, label="Order 1")
p = Plots.plot!( log10.(1.0 ./ n) , log10.(order2), color=:black, label="Order 2", linestyle = :dash)
p = Plots.annotate!(log10.(1.0 ./ N[1]), log10(order1[1]), "O1", :black)
p = Plots.annotate!(log10.(1.0 ./ N[1]), log10(order2[1]), "O2", :black)

# p = Plots.plot(  ndof_quad[2:end], t_quad[2:end], markershape=:rect,      label="Quads"                          )
# p = Plots.plot!( ndof_tri[2:end],  t_tri[2:end],  markershape=:dtriangle, label="Triangles", legend=:bottomright, xlabel = "ndof", ylabel = "time" )
display(p)

# n = 2
# tau = 1.0
# main(n, tau)

# # n    = collect(1:1:16)
# # n    = collect(17:1:20) # p2
# # tau  = collect(4:1:25)
# n    = collect(1:1:20)
# tau  = collect(1:1:4)
# resu = zeros(length(n), length(tau))
# resp = zeros(length(n), length(tau))

# for in = 1:length(n)
#     for it = 1:length(tau)
#         rp, ru = main(n[in], tau[it])
#         resu[in,it] = ru
#         resp[in,it] = rp
#     end
# end

# # p2 = Plots.heatmap(n, tau, resp', c=:jet1 )
# # display(Plots.plot(p2))

# file = matopen(string(@__DIR__,"/results/MaxPerr_p3.mat"), "w" )
# write(file, "n",        n )
# write(file, "tau",    tau )
# write(file, "resu",  resu )
# write(file, "resp",  resp )
# close(file)


text:
@testset "MonteCarlo" begin
    mc = MonteCarlo(1000)

    @test isa(mc, MonteCarlo)
    @test mc.n == 1000
end

@testset "QuasiMonteCarlo" begin

    # Make sure that all QMC Types implement qmc_samples
    for qmc in subtypes(AbstractQuasiMonteCarlo)
        @test hasmethod(qmc_samples, (qmc, Integer))
    end

    @testset "SobolSampling" begin
        sobol = SobolSampling(1000)

        @test isa(sobol, SobolSampling)
        @test sobol.n == 1000

        @testset "sample" begin
            inputs = [RandomVariable.(Uniform(), [:a, :b]); Parameter(1, :c)]

            samples = sample(inputs, SobolSampling(4))

            @test isapprox(samples.a, [0.375, 0.875, 0.625, 0.125])
            @test isapprox(samples.b, [0.375, 0.875, 0.125, 0.625])
            @test samples.c == [1.0, 1.0, 1.0, 1.0]
        end
    end

    @testset "HaltonSampling" begin
        halton = HaltonSampling(1000)

        @test isa(halton, HaltonSampling)
        @test halton.n == 1000

        @testset "sample" begin
            inputs = [RandomVariable.(Uniform(), [:a, :b]); Parameter(1, :c)]

            samples = sample(inputs, HaltonSampling(4))

            @test isapprox(samples.a, [0.5, 0.25, 0.75, 0.125])
            @test isapprox(samples.b, [1 / 3, 2 / 3, 1 / 9, 4 / 9])
            @test samples.c == [1.0, 1.0, 1.0, 1.0]
        end
    end
end

@testset "LatticeRuleSampling" begin
    lattice = LatticeRuleSampling(1000)
    @test isa(lattice, LatticeRuleSampling)
    @test lattice.n == 1000
    @testset "sample" begin
        inputs = [RandomVariable.(Uniform(), [:a, :b]); Parameter(1, :c)]

        Random.seed!(1)
        s = QuasiMonteCarlo.sample(4, [0, 0], [1, 1], LatticeRuleSample())
        Random.seed!(1)
        samples = sample(inputs, LatticeRuleSampling(4))
        Random.seed!()

        @test isapprox(samples.a, s[1, :])
        @test isapprox(samples.b, s[2, :])
        @test samples.c == [1.0, 1.0, 1.0, 1.0]
    end
end

max_stars_repo_name: FriesischScott/julia-uq

max_issues_repo_name: FriesischScott/julia-uq

max_forks_repo_name: FriesischScott/julia-uq


text:
include("../src/utils/model.jl")
include("../src/solvers/descent.jl")
include("../src/solvers/annealing.jl")
include("../src/solvers/lp_timing.jl")
include("../src/solvers/neighbour_operators.jl")

# Création d'une struct encapsulant la spécification d'un test LP
struct LpSolverSpec
    algo::Union{Symbol,Nothing}
    external_mip_solver::Union{Symbol,Nothing}
    lp_timing_solver::Union{DataType,Nothing}
    backend_model::Union{DataType,Nothing}
end

function test_one_lp_descent(spec)
    ln2("Création sol avec algo ")
    lg2(spec.algo, ":", spec.external_mip_solver)
    lg2(" (", spec.lp_timing_solver, "... ")
    Args.set(:external_mip_solver, spec.external_mip_solver)
    sol = Solution(inst, algo=spec.algo)
    model = sol.solver.model
    initial_sort!(sol, presort=:shuffle)

    @test isa(sol.solver, spec.lp_timing_solver)
    # On vérifie par exemple que le nom du solver externe "Clp" effectif
    # correspondant bien au symbole :clp du solveur demandé
    @test Symbol(lowercase(solver_name(model))) == spec.external_mip_solver
    lg2("ok (avec sol.cost=$(sol.cost))\n")

    # ===========
    lg2("Création DescentSolver nb_cons_reject_max=$(Args.get(:itermax))")
    lg2(" (cost=700<=850.0?)... ")
    sv = DescentSolver(inst)
    # println("\nconstruction faite")
    sv.do_save_bestsol = false
    sv.nb_cons_reject_max = Args.get(:itermax) # Devra être écrasé
    sv.durationmax = 1                         # Devra être écrasé
    # solve(sv, startsol=sol, nb_cons_reject_max=Args.get(:itermax))
    # Args.set(:level, 3)
    solve(sv, small_shift!, startsol=sol, nb_cons_reject_max=200, durationmax=3)
    # println("\nrésolution faite")
    @test sv.durationmax == 3
    @test sv.nb_cons_reject_max == 200
    @test isa(sv.bestsol.solver, spec.lp_timing_solver)
    # @test sv.bestsol.cost <= 850.0
    # @test sv.bestsol.cost <= 900.0
    # @test sv.bestsol.cost >= 700.0
    @test 700 <= sv.bestsol.cost <= 1_000.0
    lg1("ok ($(sv.bestsol.cost))\n")
end

function test_one_lp_annealing(spec)
    ln2("Création sol avec algo ")
    lg2(spec.algo, ":", spec.external_mip_solver)
    lg2(" (", spec.lp_timing_solver, "... ")
    Args.set(:external_mip_solver, spec.external_mip_solver)
    sol = Solution(inst, algo=spec.algo)
    model = sol.solver.model
    initial_sort!(sol, presort=:shuffle)

    @test isa(sol.solver, spec.lp_timing_solver)
    # On vérifie par exemple que le nom du solver externe "Clp" effectif
    # correspondant bien au symbole :clp du solveur demandé
    @test Symbol(lowercase(solver_name(model))) == spec.external_mip_solver
    lg2("ok (avec sol.cost=$(sol.cost))\n")

    # ===========
    lg2("Création AnnealingSolver nb_cons_reject_max=$(Args.get(:itermax))")
    lg2(" (cost=700<=850.0?)... ")
    sv = AnnealingSolver(inst, startsol=sol, nb_cons_reject_max=200,)
    solve(sv, small_shift!, durationmax=3)

    @test sv.durationmax == 3
    @test sv.nb_cons_reject_max == 200
    @test isa(sv.bestsol.solver, spec.lp_timing_solver)
    # @test sv.bestsol.cost <= 850.0
    # @test sv.bestsol.cost <= 900.0
    # @test sv.bestsol.cost >= 700.0
    @test 700 <= sv.bestsol.cost <= 1_000.0
    lg1("ok ($(sv.bestsol.cost))\n")
end


text:
# only need to be run once to install packages
#Pkg.clone("https://github.com/JuliaData/SplitApplyCombine.jl.git")
#Pkg.clone("https://github.com/xiaodaigh/FastGroupBy.jl.git")

using FastGroupBy, PooledArrays, Compat, BenchmarkTools

const N = Int(2e9/8)
const K = UInt(100)

const id3 = rand(Int32(1):Int32(round(N/K)), N)
const id6 = rand(Int64(1):Int64(round(N/K)), N)
const v1 =  rand(Int32(1):Int32(5), N)

@belapsed sumby(id6[1:2],v1[1:2])
@btime sumby(id6,v1)

# generate string ids
function randstrarray1(pool, N)
    K = length(pool)
    PooledArray(PooledArrays.RefArray(rand(1:K, N)), pool)
end

pool1 = [@sprintf "id%010d" k for k in 1:(N/K)]
id3 = randstrarray1(pool1, N)

# treat it as Pooledarray
@time sumby(id3, v1)

# treat by as strings and use dictionary method
id3_str = rand(pool1, N)
@time sumby_dict(id3_str, v1)

# parallelized sum
@time addprocs() # create Julia workers
@time using FastGroupBy
@everywhere using FastGroupBy
@everywhere using SplitApplyCombine
@time psumby(id6,v1) # 35 seconds


text:
const DYNAMIC_DSL_TRACE = Symbol("@trace")

function arg_to_ast(arg::Argument)
    ast = esc(arg.name)
    if (arg.default != nothing)
        default = something(arg.default)
        ast = Expr(:kw, ast, esc(default))
    end
    ast
end

function escape_default(arg)
    (arg.default == nothing ? nothing :
        Expr(:call, :Some, esc(something(arg.default))))
end

function make_dynamic_gen_function(name, args, body, return_type, annotations)
    escaped_args = map(arg_to_ast, args)
    gf_args = [esc(state), escaped_args...]

    julia_fn_name = gensym(name)
    julia_fn_defn = Expr(:function,
        Expr(:call, esc(julia_fn_name), gf_args...),
        esc(body))
    arg_types = map((arg) -> esc(arg.typ), args)
    arg_defaults = map(escape_default, args)
    has_argument_grads = map(
        (arg) -> (DSL_ARG_GRAD_ANNOTATION in arg.annotations), args)
    accepts_output_grad = DSL_RET_GRAD_ANNOTATION in annotations

    quote
        # first define the underlying Julia function
        $julia_fn_defn

        # now wrap it in a DynamicDSLFunction value
        Core.@__doc__ $(esc(name)) = DynamicDSLFunction(
            Type[$(arg_types...)],
            Union{Some{Any},Nothing}[$(arg_defaults...)],
            $(esc(julia_fn_name)),
            $has_argument_grads,
            $(esc(return_type)),
            $accepts_output_grad)
    end
end


text:
# Test core DispatchVariableRef extensions
@testset "Core Data Accessers" begin
    # Setup data
    m = InfiniteModel()
    obj_idx = DependentParametersIndex(1)
    idx = DependentParameterIndex(obj_idx, 1)
    domain = CollectionDomain([IntervalDomain(0, 1), IntervalDomain(0, 1)])
    methods = [InfiniteOpt.DefaultDerivativeMethod for i = 1:2]
    params = DependentParameters(domain, Dict{Vector{Float64}, Set{DataType}}(), 5, methods)
    object = MultiParameterData(params, 1, 1:2, ["p1", "p2"])
    pref = DependentParameterRef(m, idx)
    gvref = GeneralVariableRef(m, 1, DependentParameterIndex, 1)
    bad_idx = DependentParameterIndex(DependentParametersIndex(-1), 2)
    bad_pref = DependentParameterRef(m, bad_idx)
    # test dispatch_variable_ref
    @testset "dispatch_variable_ref" begin
        @test isequal(dispatch_variable_ref(m, idx), pref)
        @test isequal(dispatch_variable_ref(gvref), pref)
    end
    # test _add_data_object
    @testset "_add_data_object" begin
        @test InfiniteOpt._add_data_object(m, object) == obj_idx
        @test InfiniteOpt._param_object_indices(m)[end] == obj_idx
    end
    # test _data_dictionary
    @testset "_data_dictionary" begin
        @test InfiniteOpt._data_dictionary(pref) === m.dependent_params
        @test InfiniteOpt._data_dictionary(gvref) === m.dependent_params
        @test InfiniteOpt._data_dictionary(m, DependentParameters) === m.dependent_params
    end
    # test _data_object
    @testset "_data_object" begin
        @test InfiniteOpt._data_object(pref) == object
        @test InfiniteOpt._data_object(gvref) == object
        @test_throws ErrorException InfiniteOpt._data_object(bad_pref)
    end
    # test _core_variable_object
    @testset "_core_variable_object" begin
        @test InfiniteOpt._core_variable_object(pref) == params
        @test InfiniteOpt._core_variable_object(gvref) == params
    end
    # test _num_parameters
    @testset "_num_parameters" begin
        @test InfiniteOpt._num_parameters(pref) == 2
    end
    # test _delete_data_object
    @testset "_delete_data_object" begin
        @test is_valid(m, pref)
        @test InfiniteOpt._delete_data_object(pref) isa Nothing
        @test !is_valid(m, pref)
    end
end

# Test Definition
@testset "Definition" begin
    # Setup data
    m = InfiniteModel();
    sdomain1 = IntervalDomain(0, 1)
    sdomain2 = UniDistributionDomain(Uniform())
    domain1 = CollectionDomain([sdomain1, sdomain1])
    domain2 = MultiDistributionDomain(MvNormal(Diagonal(map(abs2, ones(2)))))
    domain3 = MultiDistributionDomain(MatrixBeta(2, 2, 2))
    domain4 = CollectionDomain([sdomain1, sdomain2])
    inds1 = IC.indices(ones(2))
    inds2 = IC.indices(ones(2, 2))
    inds3 = IC.indices(JuMPC.SparseAxisArray(Dict((i,) => domain2 for i in 1:2)))
    method = InfiniteOpt.DefaultDerivativeMethod
    # test _make_array_domain
    @testset "_make_array_domain" begin
        # multi-dimensional distributions
        @test InfiniteOpt._make_array_domain(error, [domain2, domain2], inds1) == domain2
        @test InfiniteOpt._make_array_domain(error, [domain3 for i in 1:4], inds2) == domain3
        @test_throws ErrorException InfiniteOpt._make_array_domain(error, [domain3], inds1)
        @test_throws ErrorException InfiniteOpt._make_array_domain(error, [domain2, domain2], inds3)
        @test_throws ErrorException InfiniteOpt._make_array_domain(error, [domain2, domain3], inds1)
        # collection domains
        @test InfiniteOpt._make_array_domain(error, [domain4, domain4], inds1) == domain4
        @test_throws ErrorException InfiniteOpt._make_array_domain(error, [domain4, domain1], inds1)
        @test_throws ErrorException InfiniteOpt._make_array_domain(error, [domain4, domain4], inds2)
        @test_throws ErrorException InfiniteOpt._make_array_domain(error, [domain4, domain4], inds3)
        # InfiniteArrayDomains
        @test InfiniteOpt._make_array_domain(error, [BadArrayDomain() for i in 1:2], inds1) == BadArrayDomain()
        @test_throws ErrorException InfiniteOpt._make_array_domain(error, [domain1, domain2], inds1)
        # scalar domains
        @test InfiniteOpt._make_array_domain(error, [sdomain1, sdomain1], inds1).domains == domain1.domains
        @test InfiniteOpt._make_array_domain(error, [sdomain1, sdomain2], inds1).domains == domain4.domains
        # fallback
        @test_throws ErrorException InfiniteOpt._make_array_domain(error, domain2, inds1)
    end
    # test _process_supports
    @testset "_process_supports" begin
        # single support
        supps = Dict([0., 0.] => Set([UserDefined]))
        @test InfiniteOpt._process_supports(error, [0, 0], domain1, 2) == supps
        @test_throws ErrorException InfiniteOpt._process_supports(error, [2, 2], domain1, 2)
        # multiple supports
        supps = Dict([0., 0.] => Set([UserDefined]), [1., 1.] => Set([UserDefined]))
        @test InfiniteOpt._process_supports(error, [[0, 1], [0, 1]], domain1, 2) == supps
        @test_throws ErrorException InfiniteOpt._process_supports(error, [[0, 0], [1]], domain1, 2)
        @test_throws ErrorException InfiniteOpt._process_supports(error, [[0, 2], [0, 2]], domain1, 2)
    end
    # test _process_derivative_methods
    @testset "_process_derivative_methods" begin
        same = [method, method]
        domains = [domain1, domain1]
        @test InfiniteOpt._process_derivative_methods(error, same, domains) == same
        mix = [method, TestMethod()]
        @test InfiniteOpt._process_derivative_methods(error, mix, domains) == mix
        @test InfiniteOpt._process_derivative_methods(error, [], domains) == same
        @test_throws ErrorException InfiniteOpt._process_derivative_methods(error, [1, 2], domains)
    end
    # test _build_parameters
    @testset "_build_parameters" begin
        # test errors
        @test_throws ErrorException InfiniteOpt._build_parameters(error, [domain1, domain1], inds1, bad = 42)
        @test_throws ErrorException InfiniteOpt._build_parameters(error, [domain1, domain2], inds1)
        @test_throws ErrorException InfiniteOpt._build_parameters(error, [domain2, domain2], inds2)
        @test_throws ErrorException InfiniteOpt._build_parameters(error, [domain1, domain1], inds1, supports = [[0, 0], [2]])
        @test_throws ErrorException InfiniteOpt._build_parameters(error, [domain1, domain1], inds1, derivative_method = [2, 2])
        # test vector
        @test InfiniteOpt._build_parameters(error, [domain1, domain1], inds1).domain == domain1
        @test InfiniteOpt._build_parameters(error, [domain1, domain1], inds1).supports == Dict{Vector{Float64}, Set{DataType}}()
        @test InfiniteOpt._build_parameters(error, [domain1, domain1],  inds1, supports = [0, 0]).supports == Dict([0., 0.] => Set([UserDefined]))
        @test InfiniteOpt._build_parameters(error, [domain1, domain1], inds1).sig_digits isa Int
        @test InfiniteOpt._build_parameters(error, [domain1, domain1], inds1).derivative_methods == [method, method]
        # test array
        @test InfiniteOpt._build_parameters(error, [domain3 for i in 1:4], inds2).domain == domain3
        @test InfiniteOpt._build_parameters(error, [domain3 for i in 1:4], inds2, supports = [[0, 1] for i in 1:4]).supports == Dict(zeros(4) => Set([UserDefined]), ones(4) => Set([UserDefined]))
        @test length(InfiniteOpt._build_parameters(error, [domain3 for i in 1:4], inds2, num_supports = 2).supports) == 2
        @test InfiniteOpt._build_parameters(error, [domain3 for i in 1:4], inds2, sig_digits = 3).sig_digits == 3
        @test InfiniteOpt._build_parameters(error, [domain3 for i in 1:4], inds2, derivative_method = [TestMethod() for i in 1:4]).derivative_methods isa Vector{TestMethod}
    end
    # test add_parameters
    @testset "add_parameters" begin
        # test default
        prefs = [GeneralVariableRef(m, 1, DependentParameterIndex, i) for i in 1:2]
        params = InfiniteOpt._build_parameters(error, [domain1, domain1], inds1)
        @test isequal(add_parameters(m, params, ["", ""]), prefs)
        @test InfiniteOpt._param_object_indices(m) == [index(prefs[1]).object_index]
        @test InfiniteOpt._last_param_num(m) == 2
        @test name.(prefs) == ["", ""]
        # test vector build
        prefs = [GeneralVariableRef(m, 2, DependentParameterIndex, i) for i in 1:2]
        params = InfiniteOpt._build_parameters(error, [domain2, domain2], inds1)
        @test isequal(add_parameters(m, params, ["p1", "p2"]), prefs)
        @test InfiniteOpt._param_object_indices(m)[2] == index(prefs[1]).object_index
        @test InfiniteOpt._last_param_num(m) == 4
        @test name.(prefs) == ["p1", "p2"]
        # test array build
        prefs = [GeneralVariableRef(m, 3, DependentParameterIndex, i) for i in 1:4]
        params = InfiniteOpt._build_parameters(error, [domain3 for i in 1:4], inds2)
        @test isequal(add_parameters(m, params, ["p$i" for i in 1:4]), prefs)
        @test InfiniteOpt._param_object_indices(m)[3] == index(prefs[1]).object_index
        @test InfiniteOpt._last_param_num(m) == 8
        @test name.(prefs) == ["p$i" for i in 1:4]
        # test name error
        @test_throws ErrorException add_parameters(m, params, ["p"])
    end
end

# test macro definition
@testset "Macro Definition" begin
    # Setup data
    m = InfiniteModel();
    dist1 = Uniform()
    dist2 = MvNormal(Diagonal(map(abs2, ones(2))))
    dist3 = MatrixBeta(2, 2, 2)
    sdomain1 = IntervalDomain(0, 1)
    sdomain2 = UniDistributionDomain(dist1)
    domain1 = CollectionDomain([sdomain1, sdomain1])
    domain2 = MultiDistributionDomain(dist2)
    domain3 = MultiDistributionDomain(dist3)
    domain4 = CollectionDomain([sdomain1, sdomain2])
    domain5 = CollectionDomain([sdomain2, sdomain2])
    # test @infinite_parameter
    @testset "@infinite_parameter" begin
        # test macro errors
        @test_macro_throws ErrorException @infinite_parameter(m)
        @test_macro_throws ErrorException @infinite_parameter(m, p[1:2], 42)
        @test_macro_throws ErrorException @infinite_parameter(m, param)
        @test_macro_throws ErrorException @infinite_parameter(m, "bob"[1:2])
        @test_macro_throws ErrorException @infinite_parameter(m, [1:2] in [0, 1],
                                                              independent = a)
        @test_macro_throws ErrorException @infinite_parameter(m, [1:2] in [0, 1],
                                                              distribution = dist1)
        @test_macro_throws ErrorException @infinite_parameter(m, [i = 1:2] in [0, 1],
                                                              num_supports = i)
        # test domain errors
        @test_macro_throws ErrorException @infinite_parameter(m, p[1:2] in 42)
        @test_macro_throws ErrorException @infinite_parameter(m, p[1:2] in dist1)
        @test_macro_throws ErrorException @infinite_parameter(m, p[1:2] ~ [0, 1])
        # test build errors
        @test_macro_throws ErrorException @infinite_parameter(m, p[1:2] ~ dist1, bob = 42)
        @test_macro_throws ErrorException @infinite_parameter(m, p[i = 1:2] in [domain1, domain2][i])
        @test_macro_throws ErrorException @infinite_parameter(m, p[1:2] ~ dist3)
        @test_macro_throws ErrorException @infinite_parameter(m, p[i = 1:3] in domain1)
        @test_macro_throws ErrorException @infinite_parameter(m, p[1:2] in domain1, supports = 4)
        @test_macro_throws ErrorException @infinite_parameter(m, p[i = 1:2] in domain1, supports = [[1, 0], 1][i])
        # test simple explict build
        prefs = [GeneralVariableRef(m, 1, DependentParameterIndex, i) for i in 1:2]
        @test isequal(@infinite_parameter(m, a[1:2] ~ dist1, num_supports = 10), prefs)
        @test name.(prefs) == ["a[1]", "a[2]"]
        @test length(InfiniteOpt._core_variable_object(prefs[1]).supports) == 10
        @test WeightedSample in first(InfiniteOpt._core_variable_object(prefs[1]).supports)[2]
        @test InfiniteOpt._core_variable_object(prefs[1]).domain.domains == domain5.domains
        # test another explicit build
        prefs = [GeneralVariableRef(m, 2, DependentParameterIndex, i) for i in 1:2]
        expected = JuMPC.DenseAxisArray(prefs, 3:4)
        @test isequal(@infinite_parameter(m, b[3:4] in domain2, supports = 0), expected)
        @test name.(prefs) == ["b[3]", "b[4]"]
        @test InfiniteOpt._core_variable_object(prefs[1]).supports == Dict{Vector{Float64}, Set{DataType}}(zeros(2) => Set([UserDefined]))
        # test explicit build with some args
        prefs = [GeneralVariableRef(m, 3, DependentParameterIndex, i) for i in 1:2]
        expected = convert(JuMPC.SparseAxisArray, prefs)
        @test all(isequal.(@infinite_parameter(m, c[1:2] in sdomain1, supports = 0,
                                  base_name = "z",
                                  container = SparseAxisArray), expected))
        @test name.(prefs) == ["z[1]", "z[2]"]
        # test explicit build again
        prefs = [GeneralVariableRef(m, 4, DependentParameterIndex, i) for i in 1:2]
        @test isequal(@infinite_parameter(m, d[1:2] in [0, 1], num_supports = 10), prefs)
        @test name.(prefs) == ["d[1]", "d[2]"]
        @test length(InfiniteOpt._core_variable_object(prefs[1]).supports) == 10
        @test UniformGrid in first(InfiniteOpt._core_variable_object(prefs[1]).supports)[2]
        @test InfiniteOpt._core_variable_object(prefs[1]).domain.domains == domain1.domains
        # test test anonymous
        prefs = [GeneralVariableRef(m, 5, DependentParameterIndex, i) for i in 1:4]
        prefs = reshape(prefs, (2, 2))
        @test isequal(@infinite_parameter(m, [1:2, 1:2], distribution = dist3), prefs)
        @test name.(prefs) == ["" ""; "" ""]
        @test isempty(InfiniteOpt._core_variable_object(prefs[1]).supports)
        # test anonymous with domain keyword
        prefs = [GeneralVariableRef(m, 6, DependentParameterIndex, i) for i in 1:2]
        @test isequal(@infinite_parameter(m, [1:2], domain = sdomain2,
                                  derivative_method = TestMethod()), prefs)
        # test anonymous with dist keyword and base_name
        prefs = [GeneralVariableRef(m, 7, DependentParameterIndex, i) for i in 1:2]
        @test isequal(@infinite_parameter(m, [1:2],
                                  distribution = dist1, base_name = "zz"), prefs)
        @test name.(prefs) == ["zz[1]", "zz[2]"]
    end
end

# test naming stuff
@testset "Naming" begin
    # setup data
    m = InfiniteModel();
    gvrefs = @infinite_parameter(m, a[1:2] in [0, 1])
    prefs = dispatch_variable_ref.(gvrefs)
    bad_idx = DependentParameterIndex(DependentParametersIndex(-1), 2)
    bad_pref = DependentParameterRef(m, bad_idx)
    # test _param_index
    @testset "_param_index" begin
        @test InfiniteOpt._param_index(prefs[1]) == 1
        @test InfiniteOpt._param_index(prefs[2]) == 2
    end
    # test _update_param_name_dict
    @testset "_update_param_name_dict" begin
        m.name_to_param = Dict{String, AbstractInfOptIndex}()
        @test InfiniteOpt._update_param_name_dict(m, m.dependent_params) isa Nothing
        @test m.name_to_param["a[1]"] == DependentParameterIndex(DependentParametersIndex(1), 1)
        @test m.name_to_param["a[2]"] == DependentParameterIndex(DependentParametersIndex(1), 2)
        m.name_to_param = nothing
    end
    # parameter_by_name
    @testset "parameter_by_name" begin
        @test isequal(parameter_by_name(m, "a[1]"), gvrefs[1])
        @test isequal(parameter_by_name(m, "a[2]"), gvrefs[2])
        @test isa(parameter_by_name(m, "a[3]"), Nothing)
        @infinite_parameter(m, b[2:3] in [0, 1], base_name = "a")
        @test_throws ErrorException parameter_by_name(m, "a[2]")
    end
    # test name
    @testset "JuMP.name" begin
        @test name(prefs[1]) == "a[1]"
        @test name(prefs[2]) == "a[2]"
        @test name(gvrefs[1]) == "a[1]"
        @test name(gvrefs[2]) == "a[2]"
        @test name(bad_pref) == ""
    end
    # test set_name
    @testset "JuMP.set_name" begin
        @test set_name(prefs[1], "joe") isa Nothing
        @test name(prefs[1]) == "joe"
        @test set_name(prefs[2], "joe") isa Nothing
        @test name(prefs[2]) == "joe"
        @test set_name(gvrefs[2], "joe2") isa Nothing
        @test name(gvrefs[2]) == "joe2"
    end
    # test is_valid
    @testset "JuMP.is_valid" begin
        @test is_valid(m, prefs[1])
        @test !is_valid(InfiniteModel(), prefs[1])
        @test is_valid(m, gvrefs[1])
    end
end

# test dependency functions
@testset "Parameter Dependencies" begin
    # setup data
    m = InfiniteModel();
    gvrefs = @infinite_parameter(m, a[1:2] in [0, 1])
    prefs = dispatch_variable_ref.(gvrefs)
    data = InfiniteOpt._data_object(first(prefs))
    bad_idx = DependentParameterIndex(DependentParametersIndex(-1), 2)
    bad_pref = DependentParameterRef(m, bad_idx)
    # test _infinite_variable_dependencies
    @testset "_infinite_variable_dependencies" begin
        @test InfiniteOpt._infinite_variable_dependencies(prefs[1]) == data.infinite_var_indices
        @test InfiniteOpt._infinite_variable_dependencies(prefs[2]) == data.infinite_var_indices
        @test InfiniteOpt._infinite_variable_dependencies(gvrefs[1]) == data.infinite_var_indices
        @test_throws ErrorException InfiniteOpt._infinite_variable_dependencies(bad_pref)
    end
    # test _parameter_function_dependencies
    @testset "_parameter_function_dependencies" begin
        @test InfiniteOpt._parameter_function_dependencies(prefs[1]) == data.parameter_func_indices
        @test InfiniteOpt._parameter_function_dependencies(prefs[2]) == data.parameter_func_indices
        @test InfiniteOpt._parameter_function_dependencies(gvrefs[1]) == data.parameter_func_indices
        @test_throws ErrorException InfiniteOpt._parameter_function_dependencies(bad_pref)
    end
    # test _measure_dependencies
    @testset "_measure_dependencies" begin
        @test InfiniteOpt._measure_dependencies(prefs[1]) == data.measure_indices[1]
        @test InfiniteOpt._measure_dependencies(prefs[2]) == data.measure_indices[2]
        @test InfiniteOpt._measure_dependencies(gvrefs[1]) == data.measure_indices[1]
        @test_throws ErrorException InfiniteOpt._measure_dependencies(bad_pref)
    end
    # test _constraint_dependencies
    @testset "_constraint_dependencies" begin
        @test InfiniteOpt._constraint_dependencies(prefs[1]) == data.constraint_indices[1]
        @test InfiniteOpt._constraint_dependencies(prefs[2]) == data.constraint_indices[2]
        @test InfiniteOpt._constraint_dependencies(gvrefs[1]) == data.constraint_indices[1]
        @test_throws ErrorException InfiniteOpt._constraint_dependencies(bad_pref)
    end
    # test _derivative_dependencies
    @testset "_derivative_dependencies" begin
        @test InfiniteOpt._derivative_dependencies(prefs[1]) == data.derivative_indices[1]
        @test InfiniteOpt._derivative_dependencies(prefs[2]) == data.derivative_indices[2]
        @test InfiniteOpt._derivative_dependencies(gvrefs[1]) == data.derivative_indices[1]
        @test_throws ErrorException InfiniteOpt._derivative_dependencies(bad_pref)
    end
    # test used_by_infinite_variable
    @testset "used_by_infinite_variable" begin
        # test not used
        @test !used_by_infinite_variable(prefs[1])
        @test !used_by_infinite_variable(prefs[2])
        @test !used_by_infinite_variable(gvrefs[1])
        # test used
        push!(data.infinite_var_indices, InfiniteVariableIndex(1))
        @test used_by_infinite_variable(prefs[1])
        # undo changes
        empty!(data.infinite_var_indices)
    end
    # test used_by_parameter_function
    @testset "used_by_parameter_function" begin
        # test not used
        @test !used_by_parameter_function(prefs[1])
        @test !used_by_parameter_function(prefs[2])
        @test !used_by_parameter_function(gvrefs[1])
        # test used
        push!(data.parameter_func_indices, ParameterFunctionIndex(1))
        @test used_by_parameter_function(prefs[1])
        # undo changes
        empty!(data.parameter_func_indices)
    end
    # test used_by_derivative
    @testset "used_by_derivative" begin
        # test not used
        @test !used_by_derivative(prefs[1])
        @test !used_by_derivative(prefs[2])
        @test !used_by_derivative(gvrefs[1])
        # test used
        push!(data.derivative_indices[1], DerivativeIndex(1))
        @test used_by_derivative(prefs[1])
        # undo changes
        empty!(data.derivative_indices[1])
    end
    # test used_by_measure
    @testset "used_by_measure" begin
        # test not used
        @test !used_by_measure(prefs[1])
        @test !used_by_measure(prefs[2])
        @test !used_by_measure(gvrefs[1])
        # test used
        push!(data.measure_indices[1], MeasureIndex(1))
        @test used_by_measure(prefs[1])
        # undo changes
        empty!(data.measure_indices[1])
    end
    # test used_by_constraint
    @testset "used_by_constraint" begin
        # test not used
        @test !used_by_constraint(prefs[1])
        @test !used_by_constraint(prefs[2])
        @test !used_by_constraint(gvrefs[1])
        # test used
        push!(data.constraint_indices[1], InfOptConstraintIndex(1))
        @test used_by_constraint(prefs[1])
        # undo changes
        empty!(data.constraint_indices[1])
    end
    # test used_by_objective
    @testset "used_by_objective" begin
        @test !used_by_objective(prefs[1])
        @test !used_by_objective(prefs[2])
        @test !used_by_objective(gvrefs[1])
    end
    # test is_used
    @testset "is_used" begin
        # test not used
        @test !is_used(prefs[1])
        @test !is_used(prefs[2])
        @test !is_used(gvrefs[1])
        # test used
        push!(data.infinite_var_indices, InfiniteVariableIndex(1))
        @test is_used(prefs[1])
    end
end

# test Parameter Object Methods
@testset "Object Methods" begin
    # setup data
    m = InfiniteModel();
    gvrefs = @infinite_parameter(m, a[1:2] in [0, 1])
    prefs = dispatch_variable_ref.(gvrefs)
    data = InfiniteOpt._data_object(first(prefs))
    domain = CollectionDomain([IntervalDomain(0, 2), IntervalDomain(0, 2)])
    methods = [InfiniteOpt.DefaultDerivativeMethod for i = 1:2]
    params = DependentParameters(domain, Dict{Vector{Float64}, Set{DataType}}(), 10, methods)
    bad_idx = DependentParameterIndex(DependentParametersIndex(-1), 2)
    bad_pref = DependentParameterRef(m, bad_idx)
    # test _parameter_number
    @testset "_parameter_number" begin
        @test InfiniteOpt._parameter_number(prefs[1]) == 1
        @test InfiniteOpt._parameter_number(prefs[2]) == 2
        @test InfiniteOpt._parameter_number(gvrefs[1]) == 1
        @test_throws ErrorException InfiniteOpt._parameter_number(bad_pref)
    end
    # test _parameter_numbers
    @testset "_parameter_numbers" begin
        @test InfiniteOpt._parameter_numbers(prefs[1]) == [1]
        @test InfiniteOpt._parameter_numbers(prefs[2]) == [2]
        @test InfiniteOpt._parameter_numbers(gvrefs[1]) == [1]
    end
    # test _object_number
    @testset "_object_number" begin
        @test InfiniteOpt._object_number(prefs[1]) == 1
        @test InfiniteOpt._object_number(prefs[2]) == 1
        @test InfiniteOpt._object_number(gvrefs[1]) == 1
    end
    # test _object_numbers
    @testset "_object_numbers" begin
        @test InfiniteOpt._object_numbers(prefs[1]) == [1]
        @test InfiniteOpt._object_numbers(prefs[2]) == [1]
        @test InfiniteOpt._object_numbers(gvrefs[1]) == [1]
    end
    # test _adaptive_data_update
    @testset "_adaptive_data_update" begin
        # test with same data
        @test InfiniteOpt._adaptive_data_update(prefs[1], params, data) isa Nothing
        # test with different data
        domain = MultiDistributionDomain(MvNormal([0, 0], [1, 1]))
        ps = DependentParameters(domain, Dict{Vector{Float64}, Set{DataType}}(), 10, methods)
        @test InfiniteOpt._adaptive_data_update(prefs[2], ps, data) isa Nothing
        @test InfiniteOpt._core_variable_object(prefs[2]) == ps
    end
    # test _set_core_variable_object
    @testset "_set_core_variable_object" begin
        # test with different type
        @test InfiniteOpt._set_core_variable_object(prefs[1], params) isa Nothing
        @test InfiniteOpt._core_variable_object(prefs[2]) == params
        # test with same type
        @test InfiniteOpt._set_core_variable_object(prefs[2], params) isa Nothing
        @test InfiniteOpt._core_variable_object(prefs[1]) == params
    end
end

# test Derivative Method
@testset "Derivative Methods" begin
    # setup the data
    m = InfiniteModel();
    gvrefs = @infinite_parameter(m, a[1:2] in [0, 1], derivative_method = TestMethod())
    prefs = dispatch_variable_ref.(gvrefs)
    push!(InfiniteOpt._constraint_dependencies(prefs[1]), InfOptConstraintIndex(1))
    # test has_generative_supports
    @testset "has_generative_supports" begin
        @test !has_generative_supports(prefs[1])
        @test !has_generative_supports(gvrefs[2])
    end
    # test add_generative_supports
    @testset "add_generative_supports" begin
        @test add_generative_supports(prefs[1]) isa Nothing
        @test add_generative_supports(gvrefs[2]) isa Nothing
    end
    # test _set_has_generative_supports
    @testset "_set_has_generative_supports" begin
        @test InfiniteOpt._set_has_generative_supports(prefs[1], true) isa Nothing
        @test InfiniteOpt._set_has_generative_supports(gvrefs[2], false) isa Nothing
    end
    # test has_derivative_constraints
    @testset "has_derivative_constraints" begin
        @test !has_derivative_constraints(prefs[1])
        InfiniteOpt._data_object(prefs[1]).has_deriv_constrs[1] = true
        @test has_derivative_constraints(gvrefs[1])
    end
    # test _set_has_derivative_constraints
    @testset "_set_has_derivative_constraints" begin
        @test InfiniteOpt._set_has_derivative_constraints(prefs[1], true) isa Nothing
        @test has_derivative_constraints(prefs[1])
        @test InfiniteOpt._set_has_derivative_constraints(gvrefs[1], false) isa Nothing
        @test !has_derivative_constraints(prefs[1])
    end
    # test _derivative_methods
    @testset "_derivative_methods" begin
        @test InfiniteOpt._derivative_methods(prefs[1]) isa Vector{TestMethod}
    end
    # test derivative_method
    @testset "derivative_method" begin
        @test derivative_method(prefs[1]) isa TestMethod
        @test derivative_method(gvrefs[2]) isa TestMethod
    end
    # test _adaptive_method_update
    @testset "_adaptive_method_update" begin
        # same type
        params = InfiniteOpt._core_variable_object(prefs[1])
        @test InfiniteOpt._adaptive_method_update(prefs[1], params, TestMethod()) isa Nothing
        @test derivative_method(prefs[1]) isa TestMethod
        # different type
        params = InfiniteOpt._core_variable_object(prefs[1])
        method = InfiniteOpt.DefaultDerivativeMethod
        @test InfiniteOpt._adaptive_method_update(prefs[2], params, method) isa Nothing
        @test InfiniteOpt._derivative_methods(prefs[1]) == [TestMethod(), method]
    end
    # test set_derivative_method
    @testset "set_derivative_method" begin
        # same type
        method = InfiniteOpt.DefaultDerivativeMethod
        @test set_derivative_method(prefs[1], method) isa Nothing
        @test derivative_method(prefs[1]) == method
        @test InfiniteOpt._derivative_methods(prefs[1]) == [method, method]
        # different type
        @test set_derivative_method(prefs[2], TestMethod()) isa Nothing
        @test InfiniteOpt._derivative_methods(prefs[1]) == [method, TestMethod()]
        # test error
        @test_throws ErrorException set_derivative_method(gvrefs[2], TestGenMethod())
    end
    # test set_all_derivative_methods
    @testset "set_all_derivative_methods" begin
        @infinite_parameter(m, pref in [0, 1])
        @test set_all_derivative_methods(m, TestMethod()) isa Nothing
        @test InfiniteOpt._derivative_methods(prefs[1]) == [TestMethod(), TestMethod()]
        @test derivative_method(pref) == TestMethod()
        @test_throws ErrorException set_all_derivative_methods(m, TestGenMethod())
    end
end

# test Infinite Domain Methods
@testset "Infinite Domain Methods" begin
    # setup data
    m = InfiniteModel();
    gvrefs1 = @infinite_parameter(m, a[1:2] in [0, 1], num_supports = 2)
    prefs1 = dispatch_variable_ref.(gvrefs1)
    gvrefs2 = @infinite_parameter(m, b[1:2, 1:2] ~ MatrixBeta(2, 2, 2))
    prefs2 = dispatch_variable_ref.(gvrefs2)
    push!(InfiniteOpt._constraint_dependencies(prefs1[1]), InfOptConstraintIndex(1))
    # test _parameter_domain (raw domain)
    @testset "_parameter_domain (Raw Domain)" begin
        @test InfiniteOpt._parameter_domain(prefs1[1]) isa CollectionDomain
        @test InfiniteOpt._parameter_domain(prefs2[1, 1]) isa MultiDistributionDomain
    end
    # test _parameter_domain (CollectionDomain)
    @testset "_parameter_domain (CollectionDomain)" begin
        domain = InfiniteOpt._parameter_domain(prefs1[1])
        @test InfiniteOpt._parameter_domain(domain, prefs1[1]) == IntervalDomain(0, 1)
    end
    # test _parameter_domain (Fallback)
    @testset "_parameter_domain (Fallback)" begin
        domain = InfiniteOpt._parameter_domain(prefs2[1])
        @test_throws ErrorException InfiniteOpt._parameter_domain(domain, prefs2[1])
    end
    # test infinite_domain (single)
    @testset "infinite_domain (Single)" begin
        @test infinite_domain(prefs1[1]) == IntervalDomain(0, 1)
        @test_throws ErrorException infinite_domain(prefs2[1])
        @test infinite_domain(gvrefs1[1]) == IntervalDomain(0, 1)
    end
    # test _check_complete_param_array
    @testset "_check_complete_param_array" begin
        @test InfiniteOpt._check_complete_param_array(prefs1) isa Nothing
        @test_throws ErrorException InfiniteOpt._check_complete_param_array(prefs2[:, 1])
    end
    # test infinite_domain (Array)
    @testset "infinite_domain (Array)" begin
        @test infinite_domain(prefs1) isa CollectionDomain
        @test_throws ErrorException infinite_domain(prefs2[:, 1])
        @test infinite_domain(prefs2) isa MultiDistributionDomain
        @test infinite_domain(gvrefs1) isa CollectionDomain
    end
    # test _update_parameter_domain
    @testset "_update_parameter_domain" begin
        old_domain = infinite_domain(prefs1)
        new_domain = CollectionDomain([IntervalDomain(0, 1), IntervalDomain(0, 2)])
        @test InfiniteOpt._update_parameter_domain(prefs1[1], new_domain) isa Nothing
        @test num_supports(prefs1) == 0
        @test !optimizer_model_ready(m)
        @test infinite_domain(prefs1) == new_domain
        @test InfiniteOpt._update_parameter_domain(prefs1[1], old_domain) isa Nothing
    end
    # test set_infinite_domain (Single)
    @testset "set_infinite_domain (Single)" begin
        # test normal
        @test set_infinite_domain(prefs1[2], IntervalDomain(-1, 1)) isa Nothing
        @test infinite_domain(prefs1[2]) == IntervalDomain(-1, 1)
        @test set_infinite_domain(gvrefs1[2], IntervalDomain(0, 1)) isa Nothing
        @test set_infinite_domain(prefs1[1], UniDistributionDomain(Normal())) isa Nothing
        @test infinite_domain(prefs1[1]) isa UniDistributionDomain
        @test set_infinite_domain(prefs1[1], IntervalDomain(0, 1)) isa Nothing
        # test errors
        @test_throws ErrorException set_infinite_domain(prefs2[1], IntervalDomain(0, 2))
        push!(InfiniteOpt._measure_dependencies(prefs1[1]), MeasureIndex(1))
        @test_throws ErrorException set_infinite_domain(prefs1[1], IntervalDomain(0, 2))
        empty!(InfiniteOpt._measure_dependencies(prefs1[1]))
    end
    # test set_infinite_domain (Array)
    @testset "set_infinite_domain (Array)" begin
        # test normal
        old_domain = infinite_domain(prefs1)
        new_domain = CollectionDomain([IntervalDomain(0, 1), IntervalDomain(0, 2)])
        @test set_infinite_domain(prefs1, new_domain) isa Nothing
        @test infinite_domain(prefs1) == new_domain
        @test set_infinite_domain(gvrefs1, old_domain) isa Nothing
        # test errors
        @test_throws ErrorException set_infinite_domain(prefs2[:, 1], new_domain)
        push!(InfiniteOpt._measure_dependencies(prefs1[1]), MeasureIndex(1))
        @test_throws ErrorException set_infinite_domain(prefs1, new_domain)
        empty!(InfiniteOpt._measure_dependencies(prefs1[1]))
    end
    # test JuMP.has_lower_bound
    @testset "JuMP.has_lower_bound" begin
        @test has_lower_bound(prefs1[1])
        @test has_lower_bound(gvrefs1[2])
        @test !has_lower_bound(prefs2[2])
    end
    # test JuMP.lower_bound
    @testset "JuMP.lower_bound" begin
        @test lower_bound(prefs1[1]) == 0
        @test lower_bound(gvrefs1[2]) == 0
        @test_throws ErrorException lower_bound(prefs2[2])
    end
    # test JuMP.set_lower_bound
    @testset "JuMP.set_lower_bound" begin
        @test set_lower_bound(prefs1[1], -1) isa Nothing
        @test lower_bound(prefs1[1]) == -1
        @test set_lower_bound(gvrefs1[2], 0) isa Nothing
        @test lower_bound(prefs1[2]) == 0
        @test_throws ErrorException set_lower_bound(prefs2[2], -3)
    end
    # test JuMP.has_upper_bound
    @testset "JuMP.has_upper_bound" begin
        @test has_upper_bound(prefs1[1])
        @test has_upper_bound(gvrefs1[2])
        @test !has_upper_bound(prefs2[2])
    end
    # test JuMP.upper_bound
    @testset "JuMP.upper_bound" begin
        @test upper_bound(prefs1[1]) == 1
        @test upper_bound(gvrefs1[2]) == 1
        @test_throws ErrorException upper_bound(prefs2[2])
    end
    # test JuMP.set_upper_bound
    @testset "JuMP.set_upper_bound" begin
        @test set_upper_bound(prefs1[1], 42) isa Nothing
        @test upper_bound(prefs1[1]) == 42
        @test set_upper_bound(gvrefs1[2], 0) isa Nothing
        @test upper_bound(prefs1[2]) == 0
        @test_throws ErrorException set_upper_bound(prefs2[2], 42)
    end
end

# test Support Methods
@testset "Support Methods" begin
    # Setup data
    m = InfiniteModel();
    gvrefs1 = @infinite_parameter(m, a[1:2] in [0, 1], num_supports = 2)
    prefs1 = dispatch_variable_ref.(gvrefs1)
    gvrefs2 = @infinite_parameter(m, b[1:2, 1:2] ~ MatrixBeta(2, 2, 2))
    prefs2 = dispatch_variable_ref.(gvrefs2)
    gvrefs3 = @infinite_parameter(m, c[1:2] in [0, 1])
    gvrefs4 = @infinite_parameter(m, d[1:2] in [0, 1], supports = 0)
    push!(InfiniteOpt._constraint_dependencies(prefs1[1]), InfOptConstraintIndex(1))
    # test has_internal_supports
    @testset "has_internal_supports" begin
        @test !has_internal_supports(prefs1[1])
        @test !has_internal_supports(gvrefs4[2])
    end
    # test _set_has_internal_supports
    @testset "_set_has_internal_supports" begin
        @test InfiniteOpt._set_has_internal_supports(prefs1[1], true) isa Nothing
        @test has_internal_supports(prefs1[1])
        @test InfiniteOpt._set_has_internal_supports(gvrefs1[1], false) isa Nothing
        @test !has_internal_supports(prefs1[1])
    end
    # test _parameter_supports
    @testset "_parameter_supports" begin
        @test sort(collect(keys(InfiniteOpt._parameter_supports(prefs1[1])))) == [[0., 0.], [1., 1.]]
        @test isempty(InfiniteOpt._parameter_supports(prefs2[2]))
    end
    # test significant_digits
    @testset "significant_digits" begin
        @test significant_digits(prefs1[1]) == InfiniteOpt.DefaultSigDigits
        @test significant_digits(gvrefs1[2]) == InfiniteOpt.DefaultSigDigits
    end
    # test num_supports (Single)
    @testset "num_supports (Single)" begin
        # test default
        @test num_supports(prefs1[1]) == 2
        @test num_supports(gvrefs1[2]) == 2
        @test num_supports(prefs2[2]) == 0
        # test label
        @test num_supports(prefs1[1], label = MCSample) == 0
        @test num_supports(prefs1[1], label = UniformGrid) == 2
        @test num_supports(gvrefs1[2], label = MCSample) == 0
        @test num_supports(prefs1[1], label = All) == 2
        @test num_supports(prefs1[1], label = SampleLabel) == 0
    end
    # test num_supports (Array)
    @testset "num_supports (Array)" begin
        # test default
        @test num_supports(prefs1) == 2
        @test num_supports(gvrefs1) == 2
        @test num_supports(prefs2) == 0
        # test label
        @test num_supports(prefs1, label = MCSample) == 0
        @test num_supports(prefs1, label = UniformGrid) == 2
        @test num_supports(gvrefs1, label = MCSample) == 0
        # test error
        @test_throws ErrorException num_supports(prefs2[:, 1])
    end
    # test has_supports (Single)
    @testset "has_supports (Single)" begin
        @test has_supports(prefs1[1])
        @test has_supports(gvrefs1[2])
        @test !has_supports(prefs2[2])
    end
    # test has_supports (Array)
    @testset "has_supports (Array)" begin
        @test has_supports(prefs1)
        @test has_supports(gvrefs1)
        @test !has_supports(prefs2)
        # test error
        @test_throws ErrorException has_supports(prefs2[:, 1])
    end
    # test supports (Single)
    @testset "supports (Single)" begin
        # test default
        @test sort!(supports(prefs1[1])) == Float64[0, 1]
        @test sort!(supports(gvrefs1[2])) == Float64[0, 1]
        @test supports(prefs2[2]) == Float64[]
        # test label
        @test supports(prefs1[1], label = MCSample) == Float64[]
        @test sort!(supports(prefs1[1], label = UniformGrid)) == Float64[0, 1]
        @test supports(gvrefs1[2], label = MCSample) == Float64[]
        @test supports(gvrefs4[1], label = All) == [0]
    end
    # test supports (AbstractArray)
    @testset "supports (AbstractArray)" begin
        # test default
        new_prefs1 = convert(JuMPC.SparseAxisArray, prefs1)
        new_gvrefs1 = convert(JuMPC.SparseAxisArray, gvrefs1)
        expected = convert.(JuMPC.SparseAxisArray, [zeros(2), ones(2)])
        @test sort!(supports(new_prefs1)) == expected
        @test sort!(supports(new_gvrefs1)) == expected
        @test supports(prefs2) == []
        # test label
        @test supports(new_prefs1, label = MCSample) == []
        @test sort!(supports(new_prefs1, label = UniformGrid)) == expected
        @test supports(new_gvrefs1, label = MCSample) == []
        @test supports(gvrefs4, label = All) == zeros(2, 1)
        @test supports(gvrefs4, label = AbstractSupportLabel) == zeros(2, 1)
    end
    # test supports (Vector)
    @testset "supports (Vector)" begin
        # test default
        @test sortcols(supports(prefs1)) == [0 1; 0 1]
        @test sortcols(supports(gvrefs1)) == [0 1; 0 1]
        @test supports(gvrefs3) == zeros(2, 0)
        @test supports(gvrefs4) == zeros(2, 1)
        # test label
        @test supports(prefs1, label = MCSample) == zeros(2, 0)
        @test sortcols(supports(prefs1, label = UniformGrid)) == [0 1; 0 1]
        @test supports(gvrefs1, label = MCSample) == zeros(2, 0)
        @test supports(gvrefs1, label = InternalLabel) == zeros(2, 0)
        # test error
        @test_throws ErrorException supports(prefs2[:, 1])
    end
    # test _update_parameter_supports
    @testset "_update_parameter_supports" begin
        old_supports = supports(prefs1)
        @test InfiniteOpt._update_parameter_supports(prefs1, ones(Int, 2, 3),
                                                     UserDefined) isa Nothing
        @test supports(prefs1) == ones(Float64, 2, 1)
        @test !optimizer_model_ready(m)
        @test InfiniteOpt._update_parameter_supports(prefs1, old_supports,
                                                     UniformGrid) isa Nothing
    end
    # test _make_support_matrix
    @testset "_make_support_matrix" begin
        # test errors
        supps = [ones(3, 2)]
        inds = IC.indices(prefs2)
        @test_throws ErrorException InfiniteOpt._make_support_matrix(supps, inds)
        supps = [ones(2, 2), ones(3, 2)]
        @test_throws ErrorException InfiniteOpt._make_support_matrix(supps, inds)
        # test normal
        supps = [zeros(2, 2), ones(2, 2)]
        @test InfiniteOpt._make_support_matrix(supps, inds) == hcat(zeros(4), ones(4))
    end
    # test set_supports (Single)
    @testset "set_supports (Single)" begin
        @test_throws ErrorException set_supports(prefs1[1], [0, 1])
        @test_throws ErrorException set_supports(gvrefs1[1], [0, 1], label = Mixture)
    end
    # test set_supports (Array)
    @testset "set_supports (Array)" begin
        # test errors
        @test_throws ErrorException set_supports(prefs2[:, 1], ones(2, 2))
        @test_throws ErrorException set_supports(prefs1, ones(2, 2))
        @test_throws ErrorException set_supports(prefs1, ones(2, 2) * 2, force = true)
        # test default
        old_supports = supports(prefs1)
        @test set_supports(prefs1, ones(2, 3), force = true) isa Nothing
        @test supports(prefs1) == ones(Float64, 2, 1)
        @test collect(values(InfiniteOpt._parameter_supports(prefs1[1]))) == [Set([UserDefined])]
        # test keywords
        @test set_supports(gvrefs1, old_supports, force = true,
                           label = UniformGrid) isa Nothing
        @test sortcols(supports(prefs1)) == sortcols(old_supports)
        @test collect(values(InfiniteOpt._parameter_supports(prefs1[1]))) == [Set([UniformGrid]) for i = 1:2]
        # test setting with internal supports
        @test set_supports(gvrefs1, old_supports, force = true,
                            label = InternalLabel) isa Nothing
        @test sortcols(supports(prefs1, label = All)) == sortcols(old_supports)
        @test has_internal_supports(prefs1[1])
        @test set_supports(gvrefs1, old_supports, force = true,
                           label = UniformGrid) isa Nothing
        @test !has_internal_supports(prefs1[1])
    end
    # test set_supports (Vector of Supports)
    @testset "set_supports (Vector of Supports)" begin
        # default
        supps = [ones(2, 2) for i in 1:3]
        @test set_supports(prefs2, supps) isa Nothing
        @test supports(prefs2) == [ones(Float64, 2, 2)]
        @test collect(values(InfiniteOpt._parameter_supports(prefs2[1]))) == [Set([UserDefined])]
        # test keywords
        # supps = [Float64[] for i in CartesianIndices(prefs2)]
        # @test set_supports(gvrefs2, supps, force = true, label = Mixture) isa Nothing
        # @test supports(prefs2) == [ones(Float64, 0) for i in CartesianIndices(prefs2)]
        empty!(InfiniteOpt._parameter_supports(prefs2[1]))
    end
    # test add_supports (Single)
    @testset "add_supports (Single)" begin
        @test_throws ErrorException add_supports(prefs1[1], [0, 1])
        @test_throws ErrorException add_supports(gvrefs1[1], [0, 1], label = Mixture)
    end
    # test add_supports (Vector)
    @testset "add_supports (Vector)" begin
        # test errors
        @test_throws ErrorException add_supports(prefs2[:, 1], ones(2, 2))
        @test_throws ErrorException add_supports(prefs1, ones(2, 2) * 2)
        # test default
        @test add_supports(prefs1, 0.1 * ones(2, 2)) isa Nothing
        @test sortcols(supports(prefs1)) == Float64[0 0.1 1; 0 0.1 1]
        @test InfiniteOpt._parameter_supports(prefs1[1])[0.1 * ones(2)] == Set([UserDefined])
        # test keywords
        @test add_supports(gvrefs1, zeros(2, 1), check = false,
                           label = InternalLabel) isa Nothing
        @test sortcols(supports(prefs1)) == Float64[0 0.1 1; 0 0.1 1]
        @test InfiniteOpt._parameter_supports(prefs1[1])[zeros(2)] == Set([UniformGrid, InternalLabel])
        # test warning
        InfiniteOpt._data_object(prefs1[1]).has_deriv_constrs[1] = true
        warn = "Support/method changes will invalidate existing derivative evaluation " *
               "constraints that have been added to the InfiniteModel. Thus, " *
               "these are being deleted."
        @test_logs (:warn, warn) add_supports(prefs1, 0.2 * ones(2, 2)) isa Nothing
        InfiniteOpt._data_object(prefs1[1]).has_deriv_constrs[1] = false
    end
    # test add_supports (AbstractArray)
    @testset "add_supports (AbstractArray)" begin
        # default
        supps = [ones(2, 2) for i = 1:3]
        @test add_supports(prefs2, supps) isa Nothing
        @test supports(prefs2) == [ones(Float64, 2, 2)]
        @test collect(values(InfiniteOpt._parameter_supports(prefs2[1]))) == [Set([UserDefined])]
        # test keywords
        supps = [ones(2, 2) * 0.5]
        @test add_supports(gvrefs2, supps, check = false, label = MCSample) isa Nothing
        @test sort!(supports(prefs2), by = first) == [0.5 * ones(2, 2), ones(2, 2)]
        @test Set([MCSample]) in values(InfiniteOpt._parameter_supports(prefs2[1]))
    end
    # test delete_supports (Single)
    @testset "delete_supports (Single)" begin
        @test_throws ErrorException delete_supports(prefs1[1])
        @test_throws ErrorException delete_supports(gvrefs1[1])
    end
    # test delete_supports (AbstractArray)
    @testset "delete_supports (AbstractArray)" begin
        # test errors
        @test_throws ErrorException delete_supports(prefs2[:, 1])
        push!(InfiniteOpt._measure_dependencies(prefs1[1]), MeasureIndex(1))
        @test_throws ErrorException delete_supports(prefs1)
        empty!(InfiniteOpt._measure_dependencies(prefs1[1]))
        # test label deletion
        @test delete_supports(prefs1, label = InternalLabel) isa Nothing
        @test !has_internal_supports(prefs1[1])
        # normal
        @test delete_supports(prefs1) isa Nothing
        @test supports(prefs1) == zeros(Float64, 2, 0)
        @test delete_supports(gvrefs2) isa Nothing
        @test supports(prefs2) == []
    end
end

# test Support Filling
@testset "Support Filling" begin
    # Setup data
    m = InfiniteModel();
    gvrefs1 = @infinite_parameter(m, a[1:2] in [0, 1], num_supports = 2)
    prefs1 = dispatch_variable_ref.(gvrefs1)
    gvrefs2 = @infinite_parameter(m, b[1:2, 1:2] ~ MatrixBeta(2, 2, 2))
    prefs2 = dispatch_variable_ref.(gvrefs2)
    gvref = @infinite_parameter(m, c in [0, 1])
    pref = dispatch_variable_ref(gvref)
    # test generate_and_add_supports!
    @testset "generate_and_add_supports!" begin
        # old_supports = supports(prefs1)
        domain = infinite_domain(prefs1)
        @test generate_and_add_supports!(prefs1, domain, num_supports = 2) isa Nothing
        @test sortcols(supports(prefs1)) == Float64[0 1; 0 1]
        expected = [Set([UniformGrid]) for i = 1:2]
        @test collect(values(InfiniteOpt._parameter_supports(prefs1[1]))) == expected
        # @test set_supports(prefs1, old_supports, force = true) isa Nothing
    end
    # test fill_in_supports! (Single)
    @testset "fill_in_supports! (Single)" begin
        @test_throws ErrorException fill_in_supports!(prefs1[1])
        @test_throws ErrorException fill_in_supports!(gvrefs1[1])
    end
    # test fill_in_supports! (AbstractArray)
    @testset "fill_in_supports! (AbstractArray)" begin
        # test error
        @test_throws ErrorException fill_in_supports!(prefs2[1, :])
        # test default
        @test fill_in_supports!(prefs1) isa Nothing
        @test num_supports(prefs1) == 8
        @test fill_in_supports!(gvrefs2) isa Nothing
        @test num_supports(prefs2) == 10
        # test keywords
        @test fill_in_supports!(prefs1, num_supports = 2) isa Nothing
        @test num_supports(prefs1) == 8
        @test fill_in_supports!(gvrefs2, modify = false,
                                num_supports = 20) isa Nothing
        @test num_supports(prefs2) == 10
        # delete additions
        @test delete_supports(prefs1) isa Nothing
        @test delete_supports(prefs2) isa Nothing
    end
    # test fill_in_supports! (InfiniteModel)
    @testset "fill_in_supports! (InfiniteModel)" begin
        # test default
        @test fill_in_supports!(m) isa Nothing
        @test num_supports(prefs1) == 10
        @test num_supports(prefs2) == 10
        @test num_supports(pref) == 10
        # test with keywords
        @test fill_in_supports!(m, num_supports = 12, modify = false) isa Nothing
        @test num_supports(prefs1) == 10
        @test num_supports(prefs2) == 10
        @test num_supports(pref) == 10
    end
end

# Test parameter counting and listing methods
@testset "General Queries" begin
    # Setup data
    m = InfiniteModel();
    prefs1 = @infinite_parameter(m, a[1:2] in [0, 1], num_supports = 2)
    prefs2 = @infinite_parameter(m, b[1:2, 1:2] ~ MatrixBeta(2, 2, 2))
    pref = @infinite_parameter(m, c in [0, 1])
    fpref = @finite_parameter(m, d == 10)
    # test num_parameters (Default)
    @testset "num_parameters (Default)" begin
        @test num_parameters(m) == 8
    end
    # test num_parameters (Specific Scalar)
    @testset "num_parameters (Specific Scalar)" begin
        @test num_parameters(m, IndependentParameter) == 1
        @test num_parameters(m, FiniteParameter) == 1
    end
    # test num_parameters (ScalarParameter)
    @testset "num_parameters (ScalarParameter)" begin
        @test num_parameters(m, ScalarParameter) == 2
    end
    # test num_parameters (DependentParameters)
    @testset "num_parameters (DependentParameters)" begin
        @test num_parameters(m, DependentParameters) == 6
    end
    # test num_parameters (InfiniteParameter)
    @testset "num_parameters (InfiniteParameter)" begin
        @test num_parameters(m, InfiniteParameter) == 7
    end
    # test all_parameters (Default)
    @testset "all_parameters (Default)" begin
        @test isequal(all_parameters(m), [pref; prefs1; [prefs2...]; fpref])
    end
    # test all_parameters (Specific Scalar)
    @testset "all_parameters (Specific Scalar)" begin
        @test isequal(all_parameters(m, IndependentParameter), [pref])
        @test isequal(all_parameters(m, FiniteParameter), [fpref])
    end
    # test all_parameters (ScalarParameter)
    @testset "all_parameters (ScalarParameter)" begin
        @test isequal(all_parameters(m, ScalarParameter), [pref, fpref])
    end
    # test all_parameters (DependentParameters)
    @testset "all_parameters (DependentParameters)" begin
        @test isequal(all_parameters(m, DependentParameters), [prefs1; [prefs2...]])
    end
    # test all_parameters (InfiniteParameter)
    @testset "all_parameters (InfiniteParameter)" begin
        @test isequal(all_parameters(m, InfiniteParameter), [prefs1; [prefs2...]; pref])
    end
end


text:
using Fst

#=
a = Wfst()
add_arc!(a, "0", "1", "a", "b", 0.1)
add_arc!(a, "1", "0", "a", "b", 0.2)
add_arc!(a, "1", "2", "b", "b", 0.3)
add_arc!(a, "1", "3", "b", "b", 0.4)
add_arc!(a, "2", "3", "a", "b", 0.5)
add_arc!(a, "3", "3", "a", "a", 0.6)
add_initial_state!(a, "0")
add_final_state!(a, "3")
a.final_weights["3"] = 0.7
println(wfst2dot(a))
create_pdf(a, "a.pdf")

b = Wfst()
add_arc!(b, "0", "1", "b", "b", 0.1)
add_arc!(b, "1", "1", "b", "a", 0.2)
add_arc!(b, "1", "2", "a", "b", 0.3)
add_arc!(b, "1", "3", "a", "b", 0.4)
add_arc!(b, "2", "3", "b", "a", 0.5)
add_initial_state!(b, "0")
add_final_state!(b, "3")
b.final_weights["3"] = 0.6
create_pdf(b, "b.pdf")

c = compose(a, b)
create_pdf(c, "c.pdf")
=#

# Tests for composition involving <epsilon> transitions
#t1 = Wfst()
#add_arc!(t1, 0, 1, "a", "a", 1.0)
#add_arc!(t1, 1, 2, "b", "<eps>", 1.0)
#add_arc!(t1, 2, 3, "c", "<eps>", 1.0)
#add_arc!(t1, 3, 4, "d", "d", 1.0)
#add_initial_state!(t1, 0, 1.0)
#add_final_state!(t1, 4, 1.0)
#create_pdf(t1, "t1.pdf")
#
#t2 = Wfst()
#add_arc!(t2, 0, 1, "a", "d", 1.0)
#add_arc!(t2, 1, 2, "<eps>", "e", 1.0)
#add_arc!(t2, 2, 3, "d", "a", 1.0)
#add_initial_state!(t2, 0, 1.0)
#add_final_state!(t2, 3, 1.0)
#create_pdf(t2, "t2.pdf")
#
#c = compose(t1, t2)
#create_pdf(c, "c.pdf")
#
#c_eps = compose_epsilon(t1,t2)
#create_pdf(c_eps, "c_eps.pdf")

a = read_wfst(readall("fsts/compose_test1/a.txt"), Probability_semiring)
b = read_wfst(readall("fsts/compose_test1/b.txt"), Probability_semiring)
c = compose_basic(a, b)
create_pdf(c, "fsts/compose_test1/c-Probability.pdf")

a = read_wfst(readall("fsts/compose_test1/a.txt"), Tropical_semiring)
b = read_wfst(readall("fsts/compose_test1/b.txt"), Tropical_semiring)
c = compose_basic(a, b)
create_pdf(c, "fsts/compose_test1/c-Tropical.pdf")

a = read_wfst(readall("fsts/compose_test1/a.txt"), Log_semiring)
b = read_wfst(readall("fsts/compose_test1/b.txt"), Log_semiring)
c = compose_basic(a, b)
create_pdf(c, "fsts/compose_test1/c-Log.pdf")


text:
export array_qudiff, prepare_init_state, bval, aval

"""
    Based on : arxiv.org/abs/1010.2745v2

    x' = Ax + b

    * A - input matrix.
    * b - input vector.
    * x - inital vector
    * N - dimension of b (as a power of 2).
    * h - step size.
    * tspan - time span.

    * array_qudiff(N_t,N,h,A) - generates matrix for k-step solver
    * prepare_init_state(b,x,h,N_t) - generates inital states

    LDEMSAlgHHL
    * step - step for multistep method
    * α - coefficients for xₙ
    * β - coefficent for xₙ'
"""

function bval(g::Function,alg::LDEMSAlgHHL,t,h)
    b = zero(g(1))
    for i in 1:(alg.step)
        b += alg.β[i]*g(t-(i-1)*h)
    end
    return b
end

function aval(g::Function,alg::LDEMSAlgHHL,t,h)
    sz, = size(g(1))
    A = Array{ComplexF64}(undef,sz,(alg.step + 1)*sz)
    i_mat = Matrix{Float64}(I, size(g(1)))
    A[1:sz,sz*(alg.step) + 1:sz*(alg.step + 1)] = i_mat
    for i in 1:alg.step
        A[1:sz,sz*(i - 1) + 1: sz*i] = -1*(alg.α[alg.step - i + 1]*i_mat + h*alg.β[alg.step - i + 1]*g(t - (alg.step - i)*h))
    end
    return A
end

function prepare_init_state(g::Function,alg::LDEMSAlgHHL,tspan::NTuple{2, Float64},x::Vector,h::Float64)
    N_t = round(Int, (tspan[2] - tspan[1])/h + 1) #number of time steps
    N = nextpow(2,2*N_t + 1) # To ensure we have a power of 2 dimension for matrix
    sz, = size(g(1))
    init_state = zeros(ComplexF64,2*(N)*sz)
    #inital value
    init_state[1:sz] = x
    for i in 2:N_t
        b = bval(alg,h*(i - 1) + tspan[1],h) do t g(t) end
        init_state[Int(sz*(i - 1) + 1):Int(sz*(i))] = h*b
    end
    return init_state
end

function array_qudiff(g::Function,alg::LDEMSAlgHHL,tspan::NTuple{2, Float64},h::Float64)
    sz, = size(g(1))
    i_mat = Matrix{Float64}(I, size(g(1)))
    N_t = round(Int, (tspan[2] - tspan[1])/h + 1) #number of time steps
    N = nextpow(2,2*N_t + 1) # To ensure we have a power of 2 dimension for matrix
    A_ = zeros(ComplexF64, N*sz, N*sz)
    # Generates First two rows
    @inbounds A_[1:sz, 1:sz] = i_mat
    @inbounds A_[sz + 1:2*sz, 1:sz] = -1*(i_mat + h*g(tspan[1]))
    @inbounds A_[sz + 1:2*sz,sz+1:sz*2] = i_mat
    #Generates additional rows based on k - step
    for i in 3:alg.step
        @inbounds A_[sz*(i - 1) + 1:sz*i, sz*(i - 3) + 1:sz*i] = aval(QuAB2(),(i-2)*h + tspan[1],h) do t g(t) end
    end
    for i in alg.step + 1:N_t
        @inbounds A_[sz*(i - 1) + 1:sz*(i), sz*(i - alg.step - 1) + 1:sz*i] = aval(alg,(i - 2)*h + tspan[1],h) do t g(t) end
    end
    #Generates half mirroring matrix
    for i in N_t + 1:N
        @inbounds A_[sz*(i - 1) + 1:sz*(i), sz*(i - 2) + 1:sz*(i - 1)] = -1*i_mat
        @inbounds A_[sz*(i - 1) + 1:sz*(i), sz*(i - 1) + 1:sz*i] = i_mat
    end
    A_ = [zero(A_) A_;A_' zero(A_)]
    return A_
end

function DiffEqBase.solve(prob::QuLDEProblem{uType,tType,isinplace, F, P}, alg::LDEMSAlgHHL; dt = (prob.tspan[2]-prob.tspan[1])/100, kwargs...) where {uType,tType,isinplace, F, P}
    A = prob.A
    b = prob.b
    tspan = prob.tspan
    x = prob.u0
    nreg = alg.nreg

    matx = array_qudiff(alg, tspan, dt) do t A(t) end
    initstate = prepare_init_state(alg, tspan, x, dt) do t b(t) end
    λ = maximum(eigvals(matx))
    C_value = minimum(eigvals(matx) .|> abs)*0.01;
    matx = 1/(λ*2)*matx
    initstate = initstate*1/(2*λ) |> normalize!
    res = hhlsolve(matx,initstate, nreg, C_value)
    res = res/λ
    return res
end;


text:
using Test
using Revise
using Pkg;
Pkg.activate(".")

using dlgo
using dlgo.Agent

function main()
    board_size = 5
    game = new_game(board_size)
    bot = RandomBot()

    while !is_over(game)
        println("\u1b[2J")
        print_board(game)

        if game.next_player == black
            print("-- ")
            human_move = readline(stdin)
            point = point_from_coords(strip(human_move))
            while !isvoid(game.board[point])
                println(game.board[point])
                human_move = readline(stdin)
                point = point_from_coords(strip(human_move))
            end
            move = play(point)
        else
            move = select_move(bot, game)
        end
        print_move(game.next_player, move)
        game = apply_move(game, move)
    end
end

main()


text:
# This file is a part of SimilaritySearch.jl

using SimilaritySearch
using Distances
using LinearAlgebra
using Test

function test_seq(db, dist::PreMetric, ksearch, valid_lower=1e-3)
    seq = ExhaustiveSearch(dist, db)

    for i in rand(1:length(db), 100)
        res = search(seq, db[i], ksearch)
        @test minimum(res) < valid_lower
    end
end

@testset "indexing vectors with ExhaustiveSearch" begin
    # NOTE: The following algorithms are complex enough to say we are testing it doesn't have syntax errors, a more grained test functions are required
    ksearch = 3
    db = create_vectors(1000, 4)

    for (recall_lower_bound, dist) in [
        (1.0, L2Distance()), # 1.0 -> metric, < 1.0 if dist is not a metric
        (1.0, L1Distance()),
        (1.0, LInftyDistance()),
        (0.1, SqL2Distance()),
        (1.0, LpDistance(3.0)),
        (0.1, LpDistance(0.5)),
        (1.0, AngleDistance()),
        (1.0, CosineDistance())
    ]
        test_seq(db, dist, ksearch)
    end
end

@testset "indexing sequences with ExhaustiveSearch" begin
    # NOTE: The following algorithms are complex enough to say we are testing it doesn't have syntax errors, a more grained test functions are required
    ksearch = 10
    db = [create_sequence(5, false) for i in 1:1000]
    # metric distances should achieve recall=1 (perhaps lesser because of numerical inestability)
    for dist in [
        CommonPrefixDissimilarity(),
        LevenshteinDistance(),
        LcsDistance(),
        StringHammingDistance()
    ]
        test_seq(db, dist, ksearch)
    end
end

@testset "indexing sets with ExhaustiveSearch" begin
    # NOTE: The following algorithms are complex enough to say we are testing it doesn't have syntax errors, a more grained test functions are required
    ksearch = 10
    db = [create_sequence(5, true) for i in 1:1000]
    # metric distances should achieve recall=1 (perhaps lesser because of numerical inestability)
    for dist in [
        JaccardDistance(),
        DiceDistance(),
        IntersectionDissimilarity()
    ]
        test_seq(db, dist, ksearch)
    end
end

@testset "Normalized Cosine and Normalized Angle distances" begin
    # cosine and angle distance
    ksearch = 10
    db = create_vectors(1000, 4, true)

    test_seq(db, NormalizedAngleDistance(), ksearch)
    test_seq(db, NormalizedAngleDistance(), ksearch)
    test_seq(db, NormalizedCosineDistance(), ksearch)
    test_seq(db, NormalizedCosineDistance(), ksearch)
end

@testset "Binary hamming distance" begin
    n = 1000 # number of items in the dataset
    ksearch = 10
    db = [rand(UInt32, 8) for i in 1:n]

    test_seq(db, BinaryHammingDistance(), ksearch)
end