text:
using BernsteinEllipses
using BernsteinEllipses: rsmo, radius, semimajor, semiminor

using Test

const BitsFloats = (Float32,Float64)
const Floats = (BitsFloats..., BigFloat)
const Reals = (Int, Floats...)
rnc(reals::Type) = (reals,complex.(reals))
rnc(reals::Tuple) = (reals...,complex.(reals)...)

@testset "rsmo" begin
    for T in rnc(Reals)
        @inferred rsmo(one(T))
    end

    @test rsmo(π) ≈ rsmo(big(π))
    @test rsmo(π*im) ≈ rsmo(big(π)*im)

    # Because integer types do not have a signed zero, (-im)^2 gets mapped
    # to the second quadrant (-1+0im) instead of the third (-1-im).
    # Make sure rsmo takes care of this case.
    @test rsmo(-im) == rsmo(-1.0im)

    x = -2.0; @test rsmo(complex(x,0.0)) == rsmo(complex(x,-0.0))
    x = -1.0; @test rsmo(complex(x,0.0)) == rsmo(complex(x,-0.0))
    x =  1.0; @test rsmo(complex(x,0.0)) == rsmo(complex(x,-0.0))
    x =  2.0; @test rsmo(complex(x,0.0)) == rsmo(complex(x,-0.0))

    y = -2.0; @test rsmo(complex(0.0,y)) == rsmo(complex(-0.0,y))
    y = -1.0; @test rsmo(complex(0.0,y)) == rsmo(complex(-0.0,y))
    y =  0.0; @test rsmo(complex(0.0,y)) == rsmo(complex(-0.0,y))
    y =  1.0; @test rsmo(complex(0.0,y)) == rsmo(complex(-0.0,y))
    y =  2.0; @test rsmo(complex(0.0,y)) == rsmo(complex(-0.0,y))

    x = -0.5; @test rsmo(complex(x,0.0)) == -rsmo(complex(x,-0.0))
    x =  0.0; @test rsmo(complex(x,0.0)) == -rsmo(complex(x,-0.0))
    x =  0.5; @test rsmo(complex(x,0.0)) == -rsmo(complex(x,-0.0))
end

@testset "jouk" begin
    for T in rnc(Reals)
        @inferred jouk(one(T))
        @inferred ijouk(one(T))
        @inferred ijouk(one(T),halfplane=Val(false),branch=Val(true))
        @inferred ijouk(one(T),halfplane=Val(false),branch=Val(false))
        @inferred ijouk(one(T),halfplane=Val(true),branch=Val(true))
        @inferred ijouk(one(T),halfplane=Val(true),branch=Val(false))
    end

    @test jouk(ijouk(π   )) ≈ π
    @test jouk(ijouk(π*im)) ≈ π*im
    @test jouk(ijouk(big(π)   )) ≈ π
    @test jouk(ijouk(big(π)*im)) ≈ π*im
    @test jouk(ijouk(1/π ,halfplane=Val(true))) ≈ 1/π
    @test jouk(ijouk(π*im,halfplane=Val(true))) ≈ π*im
    @test jouk(ijouk(1/big(π) ,halfplane=Val(true))) ≈ 1/big(π)
    @test jouk(ijouk(big(π)*im,halfplane=Val(true))) ≈ π*im

    for x in LinRange(-1,1,11)
        @test ijouk(complex(x,+0.0)) ≈ 1/ijouk(complex(x,-0.0))
        if x != 0
            @test ijouk(complex(1/x,+0.0), halfplane=Val(true)) ≈ 1/ijouk(complex(1/x,-0.0), halfplane=Val(true))
        end
    end

    @test abs(ijouk(π*im)) >= 1
    @test abs(ijouk(π*im, branch=Val(false))) <= 1
    @test imag(ijouk(π*im, halfplane=Val(true))) >= 0
    @test imag(ijouk(π*im, halfplane=Val(true), branch=Val(false))) <= 0
end

@testset "radius & friends" begin
    for T in rnc(Reals)
        @inferred radius(one(T))
        @inferred semimajor(one(T))
        @inferred semiminor(one(T))
        @inferred semiminor(one(T),one(real(T)))
    end

    @test radius(1) ≈ 1
    @test radius(1+1/4) ≈ 2
    @test radius(1+1/big(4)) ≈ big(2)
    @test semimajor(π) ≈ π
    @test semimajor(big(π)) ≈ big(π)
    @test semiminor(π*im) ≈ π
    @test semiminor(big(π)*im) ≈ big(π)
    @test semiminor(ComplexF64(π,ℯ),π) ≈ ℯ
    @test semiminor(Complex{BigFloat}(π,ℯ),big(π)) ≈ big(ℯ)
end


text:
module UnitfulBatsrus

import Unitful
import Unitful: q, c, μ0, ϵ0, k, me, mp
using Unitful: @unit, Unitlike
export @bu_str
export getunit, getunits

# lengths
@unit R_bu  "Rₑ"     EarthRadii        (6378)*Unitful.km        false
@unit Rg_bu  "Rg"    GanymedeRadii     (2634)*Unitful.km        false # Ganymede's radius
@unit Rm_bu  "Rm"    MercuryRadii      (2444)*Unitful.km        false # Mercury's radius

# masses
@unit me_bu "mₑ"     ElectronMass      me                       false
@unit mp_bu "mₚ"     ProtonMass        mp                       false
@unit u_bu  "u"      AMU               1*Unitful.u              false

# charges
@unit q_bu  "q"      ElementaryCharge  q                        false
@unit k_bu  "k"      BoltzmannConstant k                        false

# densities
@unit amucc_bu "amu/cc" AtomicDensity  1*Unitful.u/Unitful.cm^3 false

# velocities
@unit V_bu  "V"      Speed             1*Unitful.km/Unitful.s   false

# current densities
@unit ampm2_bu  "μA/m²" CurrentDensity  1*Unitful.μA/Unitful.m^2 false

# Others
@unit tm2_bu "nT/m²" MagneticFieldDivergence 1*Unitful.nT/Unitful.m^2 false
@unit vm2_bu "V/m²"  ElectricFieldDivergence 1*Unitful.V/Unitful.m^2  false

include("batsrusmacro.jl")

# Some gymnastics required here because if we precompile, we cannot add to
# Unitful.basefactors at compile time and expect the changes to persist to runtime.
const localunits = Unitful.basefactors
function __init__()
   merge!(Unitful.basefactors, localunits)
   Unitful.register(UnitfulBatsrus)
end

function getunit(data, var)

   # Batsrus has a bug in the 2D cuts of 3D runs: it always outputs the 3
   # coordinate units in the headline. To temporarily deal with it, I plus
   # the index by 1.
   var_ = findfirst(x->x==lowercase(var), lowercase.(data.head.variables)) + 1
   isnothing(var_) && error("$(var) not found in file header variables!")
   if data.head.headline == "normalized variables"
      var_unit = nothing
   else
      var_unit_strs = split(data.head.headline)

      if var_unit_strs[var_] == "R"
         var_unit = bu"R"
      elseif var_unit_strs[var_] == "Mp/cc"
         var_unit = bu"amucc"
      elseif var_unit_strs[var_] == "uA/m2"
         var_unit = bu"ampm2"
      elseif var_unit_strs[var_] == "V/m2"
         var_unit = bu"vm2"
      else
         var_unit = Unitful.uparse(var_unit_strs[var_])
      end
   end

   var_unit
end

function getunits(data)

   var_unit_strs = split(data.head.headline)
   var_units = [] # needs to be improved!
   for var_unit_str in var_unit_strs
      if var_unit_str == "R"
         var_unit = bu"R"
      elseif var_unit_str == "Mp/cc"
         var_unit = bu"amucc"
      elseif var_unit_str == "uA/m2"
         var_unit = bu"ampm2"
      elseif var_unit_str == "V/m2"
         var_unit = bu"vm2"
      else
         var_unit = UnitfulBatsrus.Unitful.uparse(var_unit_str)
      end
      push!(var_units, var_unit)
   end
   var_units
end

end

ext: jl


text:
immutable MassCG
    mass::Float64
    cg::XYZ
    inertia::Mat3x3
end
zero(::Type{MassCG}) = MassCG(0.0,zero(XYZ),zero(Mat3x3))
MassCG() = zero(MassCG)

function parallel_axis_inertia(mcg::MassCG,new_cg::XYZ)
    vec = new_cg - mcg.cg
    return mcg.inertia + mcg.mass*(sumabs2(vec)*eye(Mat3x3) - vec*vec')
end

function +(a::MassCG,b::MassCG)
    mass = a.mass + b.mass
    if mass == 0
        cg = zero(XYZ)
    else
        cg = (a.mass*a.cg + b.mass*b.cg)/mass
    end
    inertia = parallel_axis_inertia(a,cg) + parallel_axis_inertia(b,cg)
    MassCG(mass,cg,inertia)
end

function -(a::MassCG,b::MassCG)
    mass = a.mass - b.mass
    if mass == 0
        cg = zero(XYZ)
    else
        cg = (a.mass*a.cg - b.mass*b.cg)/mass
    end
    inertia = parallel_axis_inertia(a,cg) - parallel_axis_inertia(b,cg)
    MassCG(mass,cg,inertia)
end

MassCG(card::CAERO1,model::NastranModel) = MassCG()
MassCG(card::SPLINE1,model::NastranModel) = MassCG()
MassCG(card::SPLINE2,model::NastranModel) = MassCG()
MassCG(card::CELAS1,model::NastranModel) = MassCG()
MassCG(card::CELAS2,model::NastranModel) = MassCG()
MassCG(card::CBUSH,model::NastranModel) = MassCG()
MassCG(card::RBE2,model::NastranModel) = MassCG()
MassCG(card::RBE3,model::NastranModel) = MassCG()

function MassCG(card::CBEAM,model::NastranModel)
    prop = model.properties[card.prop_id]
    mat = model.materials[prop.mat_id]
    if prop.nonstructural_mass != 0
        throw(ExceptionError("Nonstructural Mass not supported on CBEAM"))
    end
    if mat.rho != 0
        throw(ExceptionError("Structural Mass not supported on CBEAM"))
    end
    MassCG()
end
function MassCG(card::CROD,model::NastranModel)
    prop = model.properties[card.prop_id]
    mat = model.materials[prop.mat_id]
    if prop.nonstructural_mass != 0
        throw(ExceptionError("Nonstructural Mass not supported on CROD"))
    end
    if mat.rho != 0
        throw(ExceptionError("Structural Mass not supported on CROD"))
    end
    MassCG()
end
function MassCG(card::CTRIAR,model::NastranModel)
    prop = model.properties[card.prop_id]
    mat = model.materials[prop.mat1_id]
    if prop.nonstructural_mass != 0
        throw(ExceptionError("Nonstructural Mass not supported on CTRIAR"))
    end
    if mat.rho != 0
        throw(ExceptionError("Structural Mass not supported on CTRIAR"))
    end
    MassCG()
end

function MassCG(card::CQUADR,model::NastranModel)
    prop = model.properties[card.prop_id]
    mat = model.materials[prop.mat1_id]
    if prop.nonstructural_mass != 0
        throw(ExceptionError("Nonstructural Mass not supported on CQUADR"))
    end
    if mat.rho != 0
        throw(ExceptionError("Structural Mass not supported on CQUADR"))
    end
    MassCG()
end

function inertia_from_tetra(points::Vector{XYZ})
    I11 = 0.0
    I22 = 0.0
    I33 = 0.0
    I21 = 0.0
    I31 = 0.0
    I32 = 0.0
    for i in 1:4
        for j in i:4
            I11 += points[i][2]*points[j][2] + points[i][3]*points[j][3]
            I22 += points[i][1]*points[j][1] + points[i][3]*points[j][3]
            I33 += points[i][1]*points[j][1] + points[i][2]*points[j][2]
        end
    end
    for i in 1:4
        for j in 1:4
            I21 -= points[i][1]*points[j][2]
            I31 -= points[i][1]*points[j][3]
            I32 -= points[i][2]*points[j][3]
        end
        I21 -= points[i][1]*points[i][2]
        I31 -= points[i][1]*points[i][3]
        I32 -= points[i][2]*points[i][3]
    end
    I21 /= 2.0
    I31 /= 2.0
    I32 /= 2.0
    Mat3x3([
         I11 I21 I31;
         I21 I22 I32;
         I31 I32 I33;
         ])/10.0
end

# http://docsdrive.com/pdfs/sciencepublications/jmssp/2005/8-11.pdf
function masscg_from_tetra(points::Vector{XYZ})
    cg = sum(points)/4
    vol = dot(cross(points[2] - points[1],points[3] - points[2]),points[4]-points[1])/6
    delta = map(p-> p-cg,points)
    I = inertia_from_tetra(delta)*vol
    return MassCG(vol,cg,I)
end

function masscg_from_pyramid(center::XYZ,points::XYZ...)
    face_center = sum(points)/length(points)
    mcg = MassCG()
    for i in 1:length(points)
        j = (i % length(points)) + 1
        mcg += masscg_from_tetra([points[i],points[j],face_center,center])
    end
    mcg
end

function MassCG(card::CHEXA,model::NastranModel)
    xyz1 = get_global_xyz(model,card.point1)
    xyz2 = get_global_xyz(model,card.point2)
    xyz3 = get_global_xyz(model,card.point3)
    xyz4 = get_global_xyz(model,card.point4)
    xyz5 = get_global_xyz(model,card.point5)
    xyz6 = get_global_xyz(model,card.point6)
    xyz7 = get_global_xyz(model,card.point7)
    xyz8 = get_global_xyz(model,card.point8)
    xyz9 = card.point9 == 0 ? (xyz1 + xyz2)/2 : get_global_xyz(model,card.point9)
    xyz10 = card.point10 == 0 ? (xyz2 + xyz3)/2 : get_global_xyz(model,card.point10)
    xyz11 = card.point11 == 0 ? (xyz3 + xyz4)/2 : get_global_xyz(model,card.point11)
    xyz12 = card.point12 == 0 ? (xyz4 + xyz1)/2 : get_global_xyz(model,card.point12)
    xyz13 = card.point13 == 0 ? (xyz1 + xyz5)/2 : get_global_xyz(model,card.point13)
    xyz14 = card.point14 == 0 ? (xyz2 + xyz6)/2 : get_global_xyz(model,card.point14)
    xyz15 = card.point15 == 0 ? (xyz3 + xyz7)/2 : get_global_xyz(model,card.point15)
    xyz16 = card.point16 == 0 ? (xyz4 + xyz8)/2 : get_global_xyz(model,card.point16)
    xyz17 = card.point17 == 0 ? (xyz5 + xyz6)/2 : get_global_xyz(model,card.point17)
    xyz18 = card.point18 == 0 ? (xyz6 + xyz7)/2 : get_global_xyz(model,card.point18)
    xyz19 = card.point19 == 0 ? (xyz7 + xyz8)/2 : get_global_xyz(model,card.point19)
    xyz20 = card.point20 == 0 ? (xyz8 + xyz5)/2 : get_global_xyz(model,card.point20)

    center = (xyz1 + xyz2 + xyz3 + xyz4 + xyz5 + xyz6 + xyz7 + xyz8)/8
    mcg = MassCG()
    mcg += masscg_from_pyramid(center,xyz1,xyz9,xyz2,xyz10,xyz3,xyz11,xyz4,xyz12)
    mcg += masscg_from_pyramid(center,xyz20,xyz8,xyz19,xyz7,xyz18,xyz6,xyz17,xyz5)
    mcg += masscg_from_pyramid(center,xyz1,xyz13,xyz5,xyz17,xyz6,xyz14,xyz2,xyz9)
    mcg += masscg_from_pyramid(center,xyz2,xyz14,xyz6,xyz18,xyz7,xyz15,xyz3,xyz10)
    mcg += masscg_from_pyramid(center,xyz3,xyz15,xyz7,xyz19,xyz8,xyz16,xyz4,xyz11)
    mcg += masscg_from_pyramid(center,xyz4,xyz16,xyz8,xyz20,xyz5,xyz13,xyz1,xyz12)
    mcg
end

function MassCG(card::CONM2,model::NastranModel)
    inertia = Mat3x3([
                        card.I11 card.I21 card.I31;
                        card.I21 card.I22 card.I32;
                        card.I31 card.I32 card.I33;
                        ])
    if card.csys_id == -1
        MassCG(card.mass,XYZ(card.x,card.y,card.z),inertia)
    else
        vec = XYZ(card.x,card.y,card.z)
        grid = get_global_xyz(model,card.grid_id)
        if card.csys_id == 0
            loc = grid + vec
        else
            csys = CoordinateSystem(model,card.csys_id)
            loc = grid + rotate_to_global(csys,vec)
            inertia = rotate_to_global(csys,inertia)
        end
        MassCG(card.mass,loc,inertia)
    end
end

MassCG(model::NastranModel) = mapreduce(card -> MassCG(card,model), +, values(model.elements))
MassCG(models::Vector{NastranModel}) = mapreduce(model -> MassCG(model), +, MassCG(), models)

MassCG(filename::AbstractString) = MassCG(NastranModel(filename))
MassCG(filenames::Vector{AbstractString}) =
    mapreduce(filename->MassCG(NastranModel(filename)),+,MassCG(),filenames)


text:
using ..TopOpt.TopOptProblems: StiffnessTopOptProblem, Metadata

get_fixities_node_set_name(i) = "fixed_u$(i)"

@params struct TrussProblem{xdim,T,N,M} <: StiffnessTopOptProblem{xdim,T}
    truss_grid::TrussGrid{xdim,T,N,M} # ground truss mesh
    materials::Vector{TrussFEAMaterial{T}}
    ch::ConstraintHandler{<:DofHandler{xdim,<:Ferrite.Cell{xdim,N,M},T},T}
    force::Dict{Int,SVector{xdim,T}}
    black::AbstractVector
    white::AbstractVector
    varind::AbstractVector{Int} # variable dof => free dof, based on black & white
    metadata::Metadata
end
# - `force_dof`: dof number at which the force is applied

gettrussgrid(sp::TrussProblem) = sp.truss_grid
getFerritegrid(sp::TrussProblem) = sp.truss_grid.grid

TopOpt.TopOptProblems.getE(sp::TrussProblem) = [m.E for m in sp.materials]
TopOpt.TopOptProblems.getν(sp::TrussProblem) = [m.ν for m in sp.materials]
getA(sp::TrussProblem) = [cs.A for cs in sp.truss_grid.crosssecs]
Ferrite.getnnodes(problem::StiffnessTopOptProblem) = Ferrite.getnnodes(getdh(problem).grid)

function TrussProblem(
    ::Type{Val{CellType}},
    node_points::Dict{iT,SVector{xdim,T}},
    elements::Dict{iT,Tuple{iT,iT}},
    loads::Dict{iT,SVector{xdim,T}},
    supports::Dict{iT,SVector{xdim,fT}},
    mats=TrussFEAMaterial{T}(1.0, 0.3),
    crosssecs=TrussFEACrossSec{T}(1.0),
) where {xdim,T,iT,fT,CellType}
    # unify number type
    # _T = promote_type(eltype(sizes), typeof(mats), typeof(ν), typeof(force))
    # if _T <: Integer
    #     T = Float64
    # else
    #     T = _T
    # end
    if CellType === :Linear
        truss_grid = TrussGrid(node_points, elements, supports; crosssecs)
        geom_order = 1
    else
        @assert false "Other cell type not implemented"
    end
    # reference domain dimension for a line element
    ξdim = 1
    ncells = getncells(truss_grid)

    if mats isa Vector
        @assert length(mats) == ncells
        mats = convert(Vector{TrussFEAMaterial{T}}, mats)
    elseif mats isa TrussFEAMaterial
        mats = [convert(TrussFEAMaterial{T}, mats) for i in 1:ncells]
    else
        error("Invalid mats: $(mats)")
    end

    # * load nodeset
    # the grid node ordering coincides with the input node_points
    if haskey(truss_grid.grid.nodesets, "load")
        pop!(truss_grid.grid.nodesets, "load")
    end
    load_nodesets = Set{Int}()
    for (k, _) in loads
        push!(load_nodesets, k)
    end
    addnodeset!(truss_grid.grid, "load", load_nodesets)

    # * support nodeset
    for i in 1:xdim
        if haskey(truss_grid.grid.nodesets, get_fixities_node_set_name(i))
            pop!(truss_grid.grid.nodesets, get_fixities_node_set_name(i))
        end
        support_nodesets = Set{Int}()
        for (nodeidx, condition) in supports
            if condition[i]
                push!(support_nodesets, nodeidx)
            end
        end
        addnodeset!(truss_grid.grid, get_fixities_node_set_name(i), support_nodesets)
    end

    # Create displacement field u
    dh = DofHandler(truss_grid.grid)
    if CellType === :Linear
        # truss linear
        # interpolation_space
        ip = Lagrange{ξdim,RefCube,geom_order}()
        push!(dh, :u, xdim, ip)
    else
        # TODO truss 2-order
        @assert false "not implemented"
        # ip = Lagrange{2, RefCube, 2}()
        # push!(dh, :u, xdim, ip)
    end
    close!(dh)

    ch = ConstraintHandler(dh)
    for i in 1:xdim
        dbc = Dirichlet(
            :u,
            getnodeset(truss_grid.grid, get_fixities_node_set_name(i)),
            (x, t) -> zeros(T, 1),
            [i],
        )
        add!(ch, dbc)
    end
    close!(ch)

    # update the DBC to current time
    t = T(0)
    update!(ch, t)

    metadata = Metadata(dh)

    # fnode = Tuple(getnodeset(truss_grid.grid, "load"))[1]
    # node_dofs = metadata.node_dofs
    # force_dof = node_dofs[2, fnode]

    black, white = find_black_and_white(dh)
    varind = find_varind(black, white)

    return TrussProblem(truss_grid, mats, ch, loads, black, white, varind, metadata)
end

function Base.show(io::Base.IO, mime::MIME"text/plain", sp::TrussProblem)
    println(io, "TrussProblem:")
    print(io, "    ")
    Base.show(io, mime, sp.truss_grid)
    println(io, "    E: $(sp.E)")
    println(io, "    point loads: $(length(sp.force))")
    return println(io, "    active vars: $(sum(sp.varind .!= 0))")
end

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

TopOpt.TopOptProblems.nnodespercell(p::TrussProblem) = nnodespercell(p.truss_grid)

"""
    getcloaddict(TrussProblem{xdim,T})

Get a dict (node_idx => force vector) for concentrated loads
"""
function TopOpt.TopOptProblems.getcloaddict(p::TrussProblem{xdim,T}) where {xdim,T}
    return p.force
end

function default_quad_order(::TrussProblem)
    return 1
end


text:
@testset "acquisition" begin
model = GPE([1.], [2.], MeanZero(), SEIso(1., 0.))
ac = BO.MaxMean()
opt = BO.nlopt_setup(ac, model, [-5.], [5.],
                     merge(BO.defaultoptions(typeof(model), typeof(ac)),
                           (maxtime = 3., ftol_abs = eps())))
@test opt.maxeval == 2000
@test opt.maxtime == 3.
@test opt.ftol_abs == eps()

maxf, maxx = BO.acquire_max(opt, [-5.], [5.], 10)
@test maxx ≈ [1.]
end


text:
using Plots, Printf, LinearAlgebra

const do_visu = parse(Bool, ENV["DO_VIZ"])
const do_save = parse(Bool, ENV["DO_SAVE"])
const nx = parse(Int, ENV["NX"])
const ny = parse(Int, ENV["NY"])

# enable plotting by default
# if !@isdefined do_visu; do_visu = false end

# macros to avoid array allocation
macro qHx(ix,iy)  esc(:( -(0.5*(H[$ix,$iy+1]+H[$ix+1,$iy+1]))*(0.5*(H[$ix,$iy+1]+H[$ix+1,$iy+1]))*(0.5*(H[$ix,$iy+1]+H[$ix+1,$iy+1])) * (H[$ix+1,$iy+1]-H[$ix,$iy+1])*_dx )) end
macro qHy(ix,iy)  esc(:( -(0.5*(H[$ix+1,$iy]+H[$ix+1,$iy+1]))*(0.5*(H[$ix+1,$iy]+H[$ix+1,$iy+1]))*(0.5*(H[$ix+1,$iy]+H[$ix+1,$iy+1])) * (H[$ix+1,$iy+1]-H[$ix+1,$iy])*_dy )) end
macro dtau(ix,iy) esc(:(  (1.0/(min_dxy2 / (H[$ix+1,$iy+1]*H[$ix+1,$iy+1]*H[$ix+1,$iy+1]) / 4.1) + _dt)^-1  )) end

function compute_update!(H2, dHdtau, H, Hold, _dt, damp, min_dxy2, _dx, _dy)
    Threads.@threads for iy=1:size(dHdtau,2)
    # for iy=1:size(dHdtau,2)
        for ix=1:size(dHdtau,1)
            dHdtau[ix,iy] = -(H[ix+1, iy+1] - Hold[ix+1, iy+1])*_dt +
                             (-(@qHx(ix+1,iy)-@qHx(ix,iy))*_dx -(@qHy(ix,iy+1)-@qHy(ix,iy))*_dy) +
                             damp*dHdtau[ix,iy]                        # damped rate of change
            H2[ix+1,iy+1] = H[ix+1,iy+1] + @dtau(ix,iy)*dHdtau[ix,iy]  # update rule, sets the BC as H[1]=H[end]=0
        end
    end
    return
end

function compute_residual!(ResH, H, Hold, _dt, _dx, _dy)
    Threads.@threads for iy=1:size(ResH,2)
    # for iy=1:size(ResH,2)
        for ix=1:size(ResH,1)
            ResH[ix,iy] = -(H[ix+1, iy+1] - Hold[ix+1, iy+1])*_dt +
                           (-(@qHx(ix+1,iy)-@qHx(ix,iy))*_dx -(@qHy(ix,iy+1)-@qHy(ix,iy))*_dy)
        end
    end
    return
end

function assign!(Hold, H)
    Threads.@threads for iy=1:size(H,2)
    # for iy=1:size(H,2)
        for ix=1:size(H,1)
            Hold[ix,iy] = H[ix,iy]
        end
    end
    return
end

@views function diffusion_2D_damp_perf_loop_fun(; do_visu=true, save_fig=false)
    # Physics
    lx, ly = 10.0, 10.0   # domain size
    ttot   = 1.0          # total simulation time
    dt     = 0.2          # physical time step
    # Numerics
    # nx, ny = 512, 512     # numerical grid resolution
    nout   = 100          # check error every nout
    tol    = 1e-6         # tolerance
    itMax  = 1e5          # max number of iterations
    damp   = 1-35/nx      # damping (this is a tuning parameter, dependent on e.g. grid resolution)
    # Derived numerics
    dx, dy = lx/nx, ly/ny # grid size
    xc, yc = LinRange(dx/2, lx-dx/2, nx), LinRange(dy/2, ly-dy/2, ny)
    # Array allocation
    ResH   = zeros(nx-2, ny-2) # normal grid, without boundary points
    dHdtau = zeros(nx-2, ny-2) # normal grid, without boundary points
    # Initial condition
    H      = exp.(.-(xc.-lx/2).^2 .-(yc.-ly/2)'.^2)
    Hold   = copy(H)
    H2     = copy(H)
    _dx, _dy, _dt = 1.0/dx, 1.0/dy, 1.0/dt
    min_dxy2 = min(dx,dy)^2
    t = 0.0; it = 0; ittot = 0; t_tic = 0.0; niter = 0
    # Physical time loop
    while t<ttot
        iter = 0; err = 2*tol
        # Picard-type iteration
        while err>tol && iter<itMax
            if (it==1 && iter==0) t_tic = Base.time(); niter = 0 end
            compute_update!(H2, dHdtau, H, Hold, _dt, damp, min_dxy2, _dx, _dy)
            H, H2 = H2, H  # pointer swap
            if iter % nout == 0
                compute_residual!(ResH, H, Hold, _dt, _dx, _dy)
                err = norm(ResH)/length(ResH)
            end
            iter += 1; niter += 1
        end
        ittot += iter; it += 1; t += dt
        assign!(Hold, H)
    end
    t_toc = Base.time() - t_tic
    A_eff = (2*2+1)/1e9*nx*ny*sizeof(Float64)  # Effective main memory access per iteration [GB]
    t_it  = t_toc/niter                        # Execution time per iteration [s]
    T_eff = A_eff/t_it                         # Effective memory throughput [GB/s]
    @printf("Time = %1.3f sec, T_eff = %1.2f GB/s (niter = %d)\n", t_toc, round(T_eff, sigdigits=2), niter)
    # Visualize
    if do_visu
        fontsize = 12
        opts = (aspect_ratio=1, yaxis=font(fontsize, "Courier"), xaxis=font(fontsize, "Courier"),
                ticks=nothing, framestyle=:box, titlefontsize=fontsize, titlefont="Courier", colorbar_title="",
                xlabel="Lx", ylabel="Ly", xlims=(xc[1],xc[end]), ylims=(yc[1],yc[end]), clims=(0.,1.))
        display(heatmap(xc, yc, H'; c=:davos, title="damped diffusion (nt=$it, iters=$ittot)", opts...))
        if save_fig savefig("diff2D_damp.png") end
    end
    if do_save open("./output/out_diffusion_2D_damp_perf_loop_fun.txt","a") do io; println(io, "$(nx) $(ny) $(t_toc) $(T_eff)") end end
    return
end

diffusion_2D_damp_perf_loop_fun(; do_visu=do_visu)


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

# compliance table for move up

# initialise data relevant to move up
function initCompTable!(sim::Simulation, compTable::CompTable)
	# shorthand names:
	@unpack numAmbs, numStations = sim
	ctd = sim.moveUpData.compTableData

	# set compliance table
	checkCompTable(compTable, sim) # not checking station capacities
	ctd.compTable = compTable

	ctd.compTableStationSlots = [vcat([[stationIndex for j = 1:stationSlotCount] for (stationIndex, stationSlotCount) in enumerate(ctd.compTable[i,:])]...) for i = 1:numAmbs] # compTableStationSlots[i] gives the indices of stations as many times as the number of ambulances required at the station for compTable[i,:]

	# set array sizes
	ctd.ambMovable = Vector{Bool}(undef, numAmbs)
end
function initCompTable!(sim::Simulation, nestedCompTable::NestedCompTable)
	compTable = unnestCompTable(nestedCompTable, sim.numStations)
	initCompTable!(sim, compTable)
end
function initCompTable!(sim::Simulation, compTableFilename::String)
	compTable = readCompTableFile(compTableFilename)
	initCompTable!(sim, compTable)
end

function compTableMoveUp(sim::Simulation)
	@assert(sim.moveUpData.useMoveUp)

	# shorthand:
	@unpack ambulances, stations, numAmbs, numStations = sim
	ctd = sim.moveUpData.compTableData
	@unpack compTable, ambMovable = ctd

	# get movable ambulances (movableAmbs)
	for i = 1:numAmbs
		ambMovable[i] = isAmbMovable(ambulances[i])
	end
	movableAmbs = ambulances[ambMovable]
	numMovableAmbs = length(movableAmbs)

	if numMovableAmbs == 0
		return moveUpNull()
	end

	# calculate travel time for each movable ambulance to reach every station that requires >= 1 amb
	ambToStationTimes = zeros(Float, numMovableAmbs, numStations) # ambToStationTimes[i,j] = time for movable ambulance i to travel to station j
	moveUpStations = findall(!iszero, compTable[numMovableAmbs,:]) # indices of stations that require >= 1 amb
	for i = 1:numMovableAmbs
		ambToStationTimes[i,moveUpStations] = ambMoveUpTravelTimes!(sim, movableAmbs[i]; stations = stations[moveUpStations])
	end

	# solve as an integer program
	# ambStationIndices = solveCompTableIP(compTable[numMovableAmbs,:], ambToStationTimes)

	# solve as an assignment problem
	ambStationIndices = solveCompTableAssignmentProblem(ctd.compTableStationSlots[numMovableAmbs], ambToStationTimes)

	ambStations = stations[ambStationIndices]

	if checkMode
		# check that compTable is followed
		compTableRow = compTable[numMovableAmbs,:]
		for j in ambStationIndices
			compTableRow[j] -= 1
		end
		@assert(iszero(compTableRow))
	end

	return movableAmbs, ambStations
end

function solveCompTableIP(compTableRow::Vector{Int}, ambToStationCost::Array{Float,2})
	# solve compliance table problem as an IP, minimising total cost
	# compTableRow[j] gives the number of ambulances to allocate to station j
	# ambToStationCost[i,j] is the cost of assigning movable ambulance i to station j

	# shorthand:
	a = numMovableAmbs = sum(compTableRow)
	s = numStations = length(compTableRow)

	# using JuMP
	model = Model()
	set_optimizer(model, with_optimizer(GLPK.Optimizer)) # solve speed not tested
	# set_optimizer(model, with_optimizer(Cbc.Optimizer)) # not sure how to mute output from cbc

	@variables(model, begin
		(x[i=1:a,j=1:s], Bin) # x[i,j] = 1 if ambulance: movableAmbs[i] should be moved to station: stations[j]
		# maxTravelCost # if minimising max travel cost of all ambs, instead of total
	end)

	@constraints(model, begin
		(followCompTable[j=1:s], sum(x[i,j] for i=1:a) == compTableRow[j])
		(ambAtOneLocation[i=1:a], sum(x[i,j] for j=1:s) == 1) # each ambulance must be assigned to one station
		# maxTravelCostBound[i=1:a], sum(x[i,j] * ambToStationCost[i,j] for j=1:s) <= maxTravelCost # max travel cost of all ambulances
	end)

	@expressions(model, begin
		totalAmbTravelCost, sum(x[i,j] * ambToStationCost[i,j] for i=1:a, j=1:s)
	end)

	# # testing: giving back fake results, for testing runtime without solving IP model
	# if true
		# return moveUpNull()
	# end

	# solve
	@objective(model, Min, totalAmbTravelCost) # or maxTravelCost
	optimize!(model)
	@assert(termination_status(model) == MOI.OPTIMAL)
	xValue = JuMP.value.(x) # JuMP and LightXML both export value()

	# extract solution
	sol = convert(Array{Bool,2}, round.(xValue))
	I = findall(sol)
	(ambIndices, stationIndices) = (getindex.(I, 1), getindex.(I, 2))
	ambStationIndices = zeros(Int,numMovableAmbs) # ambStationIndices[i] gives index of the station that movable ambulance i should be assigned to
	ambStationIndices[ambIndices] = stationIndices

	# if checkMode
		# # check that followCompTable constraint was met
		# for j = 1:numStations
			# @assert(sum(sol[:,j]) == compTableRow[j])
		# end
		# # check that ambAtOneLocation constraint was met
		# for i = 1:numMovableAmbs
			# @assert(sum(sol[i,:]) == 1)
		# end
	# end

	return ambStationIndices
end

# solve the compliance table problem with the Hungarian algorithm (for the assignment problem),
# this can only be used if the objective is to minimise the sum of individual ambulance redeployment costs
function solveCompTableAssignmentProblem(stationSlots::Vector{Int}, ambToStationCost::Array{Float,2}; results::Dict = Dict())
	# stationSlots gives the indices of stations as many times as the number of ambulances required at the station
	# ambToStationCost[i,j] is the cost of assigning movable ambulance i to station j
	# formulated as an assignment problem and solved with the Hungarian algorithm

	(numMovableAmbs, numStations) = size(ambToStationCost)
	@assert(length(stationSlots) == numMovableAmbs)
	@assert(all(stationIndex -> 1 <= stationIndex <= numStations, stationSlots))

	# formulate assignment problem and solve
	# need to create as many copies of station j according to values in stationSlots
	weights = ambToStationCost[:,stationSlots] # size(weights) is (numMovableAmbs, numMovableAmbs)
	matching = Hungarian.munkres(weights) # returns a sparse matrix
	# (assignment, cost) = hungarian(weights) # slightly slower than Hungarian.munkres(weights); it allows for dummy nodes but I don't need this functionality

	# extract solution
	I = findall(matching .== Hungarian.STAR)
	(ambIndices, stationSlotIndices) = (getindex.(I, 1), getindex.(I, 2))
	ambStationIndices = zeros(Int,numMovableAmbs) # ambStationIndices[i] gives index of the station that movable ambulance i should be assigned to
	ambStationIndices[ambIndices] = stationSlots[stationSlotIndices]

	# if checkMode
		# check that compliance table is followed
		# @assert(Set(ambStationIndices) == Set(stationSlots)) # this may be slow
	# end

	results[:objVal] = sum(weights[I])

	return ambStationIndices
end

# check that compliance table is valid
function checkCompTable(compTable::CompTable;
	numAmbs::Int = nullIndex, numStations::Int = nullIndex,
	stationCapacities::Union{Vector{Int},Nothing} = nothing)

	(m,n) = size(compTable)
	@assert(all(v -> v >= 0, compTable)) # need non-negative integers
	for i = 1:m
		@assert(sum(compTable[i,:]) == i) # row sum should match row index (= number of movable ambulances)
	end
	@assert(numAmbs == nullIndex || numAmbs == m)
	@assert(numStations == nullIndex || numStations == n)
	if stationCapacities != nothing
		@assert(length(stationCapacities) == n)
		for j = 1:n
			@assert(all(v -> v <= stationCapacities[j], compTable[:,j]))
			# @assert(all(compTable[:,j] .<= stationCapacities[j]))
		end
	end
	return true
end
function checkCompTable(compTable::CompTable, sim::Simulation)
	return checkCompTable(compTable; numAmbs = sim.numAmbs, numStations = sim.numStations, stationCapacities = [s.capacity for s in sim.stations])
end

function isCompTableNested(compTable::CompTable)
	checkCompTable(compTable)
	for i = 2:size(compTable,1)
		if count(!iszero, compTable[i,:] - compTable[i-1,:]) != 1 # should only have one change between rows
			return false
		end
	end
	return true
end

# take a compliance table and return the nested form, if possible
function nestCompTable(compTable::CompTable)::NestedCompTable
	@assert(isCompTableNested(compTable)) # check if nesting is possible
	(m,n) = size(compTable)
	nestedCompTable = NestedCompTable(undef,m)
	nestedCompTable[1] = findfirst(!iszero, compTable[1,:])
	for i = 2:m
		nestedCompTable[i] = findfirst(!iszero, compTable[i,:] - compTable[i-1,:])
	end
	return nestedCompTable
end

# convert nested comp table to normal compliance table (list to array)
function unnestCompTable(nestedCompTable::NestedCompTable, numStations::Int)::CompTable
	m = length(nestedCompTable)
	n = numStations
	compTable = CompTable(undef,m,n)
	compTable[:] .= 0
	for i = 1:m
		compTable[i:m, nestedCompTable[i]] .+= 1
	end
	return compTable
end

# returns a randomly generated nested comp table
function makeRandNestedCompTable(numAmbs::Int, numStations::Int;
	stationCapacities::Union{Vector{Int},Nothing} = nothing,
	rng::AbstractRNG = GLOBAL_RNG)::NestedCompTable

	@assert(numStations > 0)
	if stationCapacities == nothing
		return rand(rng, 1:numStations, numAmbs)
	else
		@assert(numStations == length(stationCapacities))
		@assert(numAmbs <= sum(stationCapacities))
		remainingCapacity = copy(stationCapacities)
		unfilledStations = Set(findall(!iszero, remainingCapacity))
		nestedCompTable = NestedCompTable(undef,numAmbs)
		nestedCompTable[:] .= 0
		for i = 1:numAmbs
			j = rand(rng, unfilledStations) # station index
			nestedCompTable[i] = j
			remainingCapacity[j] -= 1
			if remainingCapacity[j] == 0
				delete!(unfilledStations, j)
			end
		end
		return nestedCompTable
	end
end


text:
#type alias for string or symbol
StringOrSymbol = Union{String,Symbol}
@fenicsclass Mesh  #https://fenicsproject.org/olddocs/dolfin/1.5.0/python/programmers-reference/cpp/mesh/Mesh.html
#are converted automatically by PyCall

# Get the cell orientations set.
cell_orientations(mesh::Mesh) = fenicspycall(mesh, :cell_orientations)
#returns cell connectivity
cells(mesh::Mesh) = fenicspycall(mesh, :cells)
#Compute minimum cell diameter.
hmin(mesh::Mesh) = fenicspycall(mesh, :hmin)
#Compute maximum cell diameter.
hmax(mesh::Mesh) = fenicspycall(mesh, :hmax)
init(mesh::Mesh) = fenicspycall(mesh, :init)
init(mesh::Mesh, dim::Int) = fenicspycall(mesh, :init, dim) # version with dims
init_global(mesh::Mesh) = fenicspycall(mesh,:init_global)
#returns coordinates of all vertices
coordinates(mesh::Mesh) = fenicspycall(mesh,:coordinates)

data(mesh::Mesh) = fenicspycall(mesh,:data)
#Get  mesh (sub)domains
domains(mesh::Mesh) = fenicspycall(mesh, :domains)
topology(mesh::Mesh) = fenicspycall(mesh,:topology)

#get mesh geometry
geometry(mesh::Mesh) = fenicspycall(mesh,:geometry)
#returns number of cells
num_cells(mesh::Mesh) = fenicspycall(mesh, :num_cells)
#returns number of edges
num_edges(mesh::Mesh) = fenicspycall(mesh, :num_edges)
#Get number of entities of given topological dimension.
num_entities(mesh::Mesh, dim::Int) = fenicspycall(mesh, :num_entities, dim)
#Get number of faces in mesh.
num_faces(mesh::Mesh) = fenicspycall(mesh, :num_faces)
#Get number of facets in mesh.
num_facets(mesh::Mesh) = fenicspycall(mesh, :num_facets)
#Get number of vertices in mesh.
num_vertices(mesh::Mesh) = fenicspycall(mesh, :num_vertices)
#hash(mesh::Mesh) = fenicspycall(Mesh, :hash)
bounding_box_tree(mesh::Mesh) = fenicspycall(mesh,:bounding_box_tree) #this object is a pyobject
#Compute maximum cell inradius.
rmax(mesh::Mesh)=fenicspycall(mesh, :rmax)
#Compute minimum cell inradius.
rmin(mesh::Mesh)=fenicspycall(mesh, :rmin)
#Get number of local entities of given topological dimension.
size(mesh::Mesh, dim::Int) = fenicspycall(mesh, :size, dim) # version with dims
#Returns the ufl cell of the mesh.
ufl_cell(mesh::Mesh)=fenicspycall(mesh, :ufl_cell)
#Returns the ufl Domain corresponding to the mesh.
ufl_domain(mesh::Mesh)=fenicspycall(mesh, :ufl_domain)
#Returns an id that UFL can use to decide if two objects are the same.
ufl_id(mesh::Mesh)=fenicspycall(mesh, :ufl_id)

# Return symbolic cell diameter
CellDiameter(mesh::Mesh) = Expression(fenics.CellDiameter(mesh.pyobject))
# Return symbolic cell normal
CellNormal(mesh::Mesh) = Expression(fenics.CellNormal(mesh.pyobject))
# Return symbolic cell volume
CellVolume(mesh::Mesh) = Expression(fenics.CellVolume(mesh.pyobject))

export cell_orientations,cells,hmin , hmax, init, init_global, coordinates, data,
domains, geometry,topology, num_cells,num_edges,num_entities,num_faces,num_facets,num_vertices, bounding_box_tree,
rmax, rmin, size, ufl_cell , ufl_domain, ufl_id, CellDiameter, CellNormal, CellVolume

# This constant is initialized in __init__
export CellType

"""
Mesh(path::StringOrSymbol) \n
Creates a Mesh based on a specified filename(path)
"""
Mesh(path::StringOrSymbol) = Mesh(fenics.Mesh(path))

"""
Mesh(object::Mesh) \n
Creates a copy of a mesh
"""
Mesh(object::Mesh) = Mesh(fenics.Mesh(object.pyobject))

"""
UnitTriangleMesh() \n
A mesh consisting of a single triangle with vertices at \n
(0, 0) (1, 0) (0, 1)
"""
UnitTriangleMesh() = Mesh(fenics.cpp.generation.UnitTriangleMesh.create())

"""
UnitTetrahedronMesh() \n
A mesh consisting of a single tetrahedron with vertices at \n
(0, 0, 0) (1, 0, 0) (0, 1, 0) (0, 0, 1)
"""
UnitTetrahedronMesh() = Mesh(fenics.cpp.generation.UnitTetrahedronMesh.create())

"""
UnitSquareMesh(nx::Int, ny::Int, diagonal::StringOrSymbol="right" ) \n

Triangular/quadrilateral mesh of the 2D unit square [0,1] x [0,1]. \n
Given the number of cells (nx, ny) in each direction, the total number of triangles \n
will be 2*nx*ny and the total number of vertices will be (nx + 1)*(ny + 1) \n
diagonal ("left", "right", "right//left", "left//right", or "crossed") indicates the direction of the diagonals.
"""
UnitSquareMesh(nx::Int, ny::Int, diagonal::StringOrSymbol="right") = Mesh(fenics.UnitSquareMesh(nx, ny, diagonal))
UnitSquareMesh(nx::Int, ny::Int, cellType::PyObject) = Mesh(fenics.UnitSquareMesh.create(nx, ny, cellType))

function UnitQuadMesh(nx::Int,ny::Int)
    println("Deprecated in FEniCS v.2018, remove in .7 Julia")
end

"""
UnitIntervalMesh(nx::Int) \n

A mesh of the unit interval (0, 1) with a given number of cells (nx) in the axial direction. \n
The total number of intervals will be nx and the total number of vertices will be (nx + 1).
"""
UnitIntervalMesh(nx::Int) = Mesh(fenics.UnitIntervalMesh(nx))

"""
UnitCubeMesh(nx::Int, ny::Int, nz::Int) \n

Tetrahedral/hexahedral mesh of the 3D unit cube [0,1] x [0,1] x [0,1]. \n
Given the number of cells (nx, ny, nz) in each direction, the total number of \n
tetrahedra will be 6*nx*ny*nz and the total number of vertices will be (nx + 1)*(ny + 1)*(nz + 1).

"""
UnitCubeMesh(nx::Int, ny::Int, nz::Int) = Mesh(fenics.UnitCubeMesh(nx,ny,nz))
UnitCubeMesh(nx::Int, ny::Int, nz::Int, cellType::PyObject) = Mesh(fenics.UnitCubeMesh.create(nx, ny, nz, cellType))
"""
BoxMesh(p0, p1, nx::Int, ny::Int, nz::Int) \n

Tetrahedral mesh of the 3D rectangular prism spanned by two points p0 and p1. \n
Given the number of cells (nx, ny, nz) in each direction, the total number of \n
tetrahedra will be 6*nx*ny*nz and the total number of vertices will be (nx + 1)*(ny + 1)*(nz + 1).
"""
BoxMesh(p0, p1, nx::Int, ny::Int, nz::Int)= Mesh(fenics.BoxMesh(p0,p1,nx,ny,nz))
BoxMesh(p::NTuple{2, PyObject}, n::NTuple{3, Int}, cellType::PyObject) = Mesh(fenics.BoxMesh.create(p, n, cellType))

"""
RectangleMesh(p0,p1,nx::Int,ny::Int,diagdir::StringOrSymbol="right") \n
Triangular mesh of the 2D rectangle spanned by two points p0 and p1. \n
Given the number of cells (nx, ny) in each direction, the total number \n
of triangles will be 2*nx*ny and the total number of vertices will be (nx + 1)*(ny + 1) \n
diagdir ("left", "right", "right/left", "left/right", or "crossed") indicates the direction of the diagonals.
"""
RectangleMesh(p0,p1,nx::Int,ny::Int,diagdir::StringOrSymbol="right") = Mesh(fenics.RectangleMesh(p0,p1,nx,ny,diagdir))
RectangleMesh(p::NTuple{2, PyObject}, n::NTuple{2, Int}, cellType::PyObject) = Mesh(fenics.RectangleMesh.create(p, n, cellType))

"""
BoundaryMesh(mesh::Mesh,type_boundary::StringOrSymbol="exterior",order=true) \n

A BoundaryMesh  is a mesh over the boundary of some given mesh. \n
The cells of the boundary mesh (facets of the original mesh) are oriented to \n
produce outward pointing normals relative to the original mesh. \n
The type_boundary can be "exterior", "interior" or "local". "exterior" is the globally \n
external boundary, "interior" is the inter-process mesh and "local" is the boundary \n
of the local (this process) mesh. \n
order:(bool) Optional argument which can be used to control whether or not the \n
boundary mesh should be ordered according to the UFC ordering convention. \n
If set to false, the boundary mesh will be ordered with right-oriented facets \n
(outward-pointing unit normals). The default value is true.
"""
BoundaryMesh(mesh::Mesh,type_boundary::StringOrSymbol="exterior",order=true) = Mesh(fenics.BoundaryMesh(mesh.pyobject, type_boundary, order))

export UnitTriangleMesh, UnitTetrahedronMesh, UnitSquareMesh, UnitQuadMesh,
UnitIntervalMesh, UnitCubeMesh, BoxMesh, RectangleMesh, Mesh, BoundaryMesh

function pyUnitTriangleMesh()
  fenics.cpp.generation.UnitTriangleMesh.create()
end

function pyUnitTetrahedronMesh()
  fenics.cpp.generation.UnitTetrahedronMesh.create()
end

function pyUnitCubeMesh(nx::Int, ny::Int, nz::Int)
  pycall(fenics.UnitCubeMesh::PyObject,PyObject::Type,nx,ny,nz)
end

function pyBoxMesh(p0, p1, nx::Int, ny::Int, nz::Int) # look at array types to declare p0,p1
  pycall(fenics.BoxMesh::PyObject,PyObject::Type,p0,p1,nx,ny,nz)
end

function pyRectangleMesh(p0,p1,nx::Int,ny::Int,diagdir::StringOrSymbol="right")
  pycall(fenics.RectangleMesh::PyObject,PyObject::Type,p0,p1,nx,ny,diagdir)
end

"""
For the diagdir, the possible options can be found below (these indicate the direction of the diagonals)
  (“left”, “right”, “right/left”, “left/right”, or “crossed”).
"""
function pyUnitSquareMesh(nx::Int,ny::Int,diagdir::StringOrSymbol="right")
  pycall(fenics.UnitSquareMesh::PyObject,PyObject::Type,nx,ny,diagdir)
end

function pyUnitQuadMesh(nx::Int,ny::Int)
  pycall(fenics.UnitSquareMesh::PyObject,PyObject::Type,nx,ny)
end #https://fenicsproject.org/olddocs/dolfin/2016.2.0/python/programmers-reference/cpp/mesh/UnitQuadMesh.html
#states that the UnitQuadMesh code is experimental. Nevertheless I plan to add it , and maybe remove it at the final
#iteration

function pyUnitIntervalMesh(nx::Int)
  pycall(fenics.UnitIntervalMesh::PyObject,PyObject::Type,nx)
end

function pyMesh(path::StringOrSymbol)
  pycall(fenics.Mesh::PyObject,PyObject::Type,path)
end

function Point(point::Union{Vector,Tuple})
  pycall(fenics.Point::PyObject,PyObject::Type,point)
end

export pyUnitTriangleMesh, pyUnitTetrahedronMesh, pyUnitSquareMesh, pyUnitQuadMesh,
pyUnitIntervalMesh, pyUnitCubeMesh, pyBoxMesh, pyRectangleMesh,pyMesh, Point

@fenicsclass Cell #https://fenicsproject.org/docs/dolfin/1.5.0/python/programmers-reference/cpp/mesh/Cell.html
Cell(mesh::MeshImpl, i::Int) = Cell(fenics.Cell(mesh.pyobject, i))

# Get coordinates of cell vertices
get_vertex_coordinates(cell::Cell) = fenicspycall(cell, :get_vertex_coordinates)
# Compute greatest distance between between two vertices of a cell
h(cell::Cell) = fenicspycall(cell, :h)
# Compute midpoint of a cell
midpoint(cell::Cell) = fenicspycall(cell, :midpoint)
# Compute volume of cell
volume(cell::Cell) = fenicspycall(cell, :volume)

export Cell, get_vertex_coordinates, h, midpoint, volume


text:
using JuMP: REPLMode, IJuliaMode

# Test string creation
@testset "String Creators" begin
    # initialize model and attributes
    m = InfiniteModel()
    @infinite_parameter(m, 0 <= par1 <= 1)
    @infinite_parameter(m, pars[1:2] in MvNormal([1, 1], 1))
    @infinite_parameter(m, pars2[1:2] in [0, 2])
    @infinite_parameter(m, pars3[1:2] in [0, 1], independent = true)
    @infinite_variable(m, x(par1))
    @infinite_variable(m, z(pars))
    @infinite_variable(m, inf(pars, par1, pars3))
    @hold_variable(m, y)
    d1 = @deriv(x, par1)
    d2 = @deriv(inf, pars[1], par1)
    d3 = @deriv(z, pars[2])
    set_name(d3, "d3")
    @objective(m, Min, 2 + y)
    @constraint(m, c1, x + y - 2 <= 0)
    ac1 = @constraint(m, x + y - 2 <= 0)
    @constraint(m, c2, y^2 - 3 == 0)
    ac2 = @constraint(m, y^2 - 3 == 0)
    @BDconstraint(m, c3(par1 in [0, 0.5]), x == 5)
    ac3 = @BDconstraint(m, (par1 in [0, 0.5]), x == 5)
    # test _infopt_math_symbol (REPL)
    @testset "_infopt_math_symbol (REPL)" begin
        if Sys.iswindows()
            @test InfiniteOpt._infopt_math_symbol(REPLMode, :intersect) == "and"
            @test InfiniteOpt._infopt_math_symbol(REPLMode, :partial) == "d"
            @test InfiniteOpt._infopt_math_symbol(REPLMode, :integral) == "integral"
            @test InfiniteOpt._infopt_math_symbol(REPLMode, :expect) == "E"
        else
            @test InfiniteOpt._infopt_math_symbol(REPLMode, :intersect) == "∩"
            @test InfiniteOpt._infopt_math_symbol(REPLMode, :partial) == "∂"
            @test InfiniteOpt._infopt_math_symbol(REPLMode, :integral) == "∫"
            @test InfiniteOpt._infopt_math_symbol(REPLMode, :expect) == "𝔼"
        end
        @test InfiniteOpt._infopt_math_symbol(REPLMode, :times) == "*"
        @test InfiniteOpt._infopt_math_symbol(REPLMode, :prop) == "~"
    end
    # test _infopt_math_symbol (IJulia)
    @testset "_infopt_math_symbol (IJulia)" begin
        @test InfiniteOpt._infopt_math_symbol(IJuliaMode, :intersect) == "\\cap"
        @test InfiniteOpt._infopt_math_symbol(IJuliaMode, :eq) == "="
        @test InfiniteOpt._infopt_math_symbol(IJuliaMode, :prop) == "\\sim"
        @test InfiniteOpt._infopt_math_symbol(IJuliaMode, :partial) == "\\partial"
        @test InfiniteOpt._infopt_math_symbol(IJuliaMode, :open_rng) == "\\left["
        @test InfiniteOpt._infopt_math_symbol(IJuliaMode, :close_rng) == "\\right]"
        @test InfiniteOpt._infopt_math_symbol(IJuliaMode, :integral) == "\\int"
        @test InfiniteOpt._infopt_math_symbol(IJuliaMode, :expect) == "\\mathbb{E}"
    end
    # test _plural
    @testset "_plural" begin
        @test InfiniteOpt._plural(1) == ""
        @test InfiniteOpt._plural(2) == "s"
    end
    # test set_string (IntervalSet)
    @testset "set_string (IntervalSet)" begin
        # test simple case
        set = IntervalSet(0, 1)
        @test InfiniteOpt.set_string(REPLMode, set) == "[0, 1]"
        @test InfiniteOpt.set_string(IJuliaMode, set) == "[0, 1]"
        # test rounding case
        set = IntervalSet(-0, 1)
        @test InfiniteOpt.set_string(REPLMode, set) == "[0, 1]"
        @test InfiniteOpt.set_string(IJuliaMode, set) == "[0, 1]"
        # test decimal case
        set = IntervalSet(0.1, 1.3)
        @test InfiniteOpt.set_string(REPLMode, set) == "[0.1, 1.3]"
        @test InfiniteOpt.set_string(IJuliaMode, set) == "[0.1, 1.3]"
    end
    # test set_string (DistributionSet)
    @testset "set_string (DistributionSet)" begin
        # test univariate set
        set = UniDistributionSet(Uniform())
        @test InfiniteOpt.set_string(REPLMode, set) == "Uniform{Float64}(a=0.0, b=1.0)"
        @test InfiniteOpt.set_string(IJuliaMode, set) == "Uniform{Float64}(a=0.0, b=1.0)"
        # test mulivariate set
        set = MultiDistributionSet(MvNormal([1], 1))
        str = "IsoNormal(\ndim: 1\nμ: [1.0]\nΣ: [1.0]\n)\n"
        @test InfiniteOpt.set_string(REPLMode, set) == str
        @test InfiniteOpt.set_string(IJuliaMode, set) == str
    end
    # test set_string (CollectionSet)
    @testset "set_string (CollectionSet)" begin
        set = CollectionSet([IntervalSet(0, 1), IntervalSet(0, 0.1)])
        str = "CollectionSet with 2 sets:\n [0, 1]\n [0, 0.1]"
        @test InfiniteOpt.set_string(REPLMode, set) == str
        @test InfiniteOpt.set_string(IJuliaMode, set) == str
    end
    # test set_string (Fallback)
    @testset "set_string (Fallback)" begin
        set = BadSet()
        @test InfiniteOpt.set_string(REPLMode, set) == "BadSet()"
        @test InfiniteOpt.set_string(IJuliaMode, set) == "BadSet()"
    end
    # test in_set_string (IntervalSet)
    @testset "JuMP.in_set_string (Interval)" begin
        # test simple case
        set = IntervalSet(0, 1)
        str = JuMP._math_symbol(REPLMode, :in) * " [0, 1]"
        @test in_set_string(REPLMode, set) == str
        str = JuMP._math_symbol(IJuliaMode, :in) * " [0, 1]"
        @test in_set_string(IJuliaMode, set) == str
        # test rounding case
        set = IntervalSet(-0, 1)
        str = JuMP._math_symbol(REPLMode, :in) * " [0, 1]"
        @test in_set_string(REPLMode, set) == str
        str = JuMP._math_symbol(IJuliaMode, :in) * " [0, 1]"
        @test in_set_string(IJuliaMode, set) == str
        # test decimal case
        set = IntervalSet(0.1, 1.3)
        str = JuMP._math_symbol(REPLMode, :in) * " [0.1, 1.3]"
        @test in_set_string(REPLMode, set) == str
        str = JuMP._math_symbol(IJuliaMode, :in) * " [0.1, 1.3]"
        @test in_set_string(IJuliaMode, set) == str
        # test finite case
        set = IntervalSet(0.1, 0.1)
        str = JuMP._math_symbol(REPLMode, :eq) * " 0.1"
        @test in_set_string(REPLMode, set) == str
        str = JuMP._math_symbol(IJuliaMode, :eq) * " 0.1"
        @test in_set_string(IJuliaMode, set) == str
    end
    # test in_set_string (Distribution)
    @testset "JuMP.in_set_string (Distribution)" begin
        # test univariate set
        set = UniDistributionSet(Uniform())
        str = InfiniteOpt._infopt_math_symbol(REPLMode, :prop) * " Uniform"
        @test in_set_string(REPLMode, set) == str
        str = InfiniteOpt._infopt_math_symbol(IJuliaMode, :prop) * " Uniform"
        @test in_set_string(IJuliaMode, set) == str
        # test mulivariate set
        set = MultiDistributionSet(MvNormal([1], 1))
        str = InfiniteOpt._infopt_math_symbol(REPLMode, :prop) * " MvNormal(dim: (1))"
        @test in_set_string(REPLMode, set) == str
        str = InfiniteOpt._infopt_math_symbol(IJuliaMode, :prop) * " MvNormal(dim: (1))"
        @test in_set_string(IJuliaMode, set) == str
        # test matrix set
        set = MultiDistributionSet(MatrixBeta(2, 2, 2))
        str = InfiniteOpt._infopt_math_symbol(REPLMode, :prop) * " MatrixBeta(dims: (2, 2))"
        @test in_set_string(REPLMode, set) == str
        str = InfiniteOpt._infopt_math_symbol(IJuliaMode, :prop) * " MatrixBeta(dims: (2, 2))"
        @test in_set_string(IJuliaMode, set) == str
    end
    # test in_set_string (Fallback)
    @testset "JuMP.in_set_string (Fallback)" begin
        set = BadSet()
        in1 = JuMP._math_symbol(REPLMode, :in)
        in2 = JuMP._math_symbol(IJuliaMode, :in)
        @test in_set_string(REPLMode, set) == in1 * " BadSet()"
        @test in_set_string(IJuliaMode, set) == in2 * " BadSet()"
    end
    # test in_set_string (IntervalSet with Bounds)
    @testset "JuMP.in_set_string (IntervalSet with Bounds)" begin
        # test in bounds
        bounds = ParameterBounds((par1 => IntervalSet(0, 0),))
        set = IntervalSet(0, 1)
        str = JuMP._math_symbol(REPLMode, :eq) * " 0"
        @test in_set_string(REPLMode, par1, set, bounds) == str
        str = JuMP._math_symbol(IJuliaMode, :eq) * " 0"
        @test in_set_string(IJuliaMode, par1, set, bounds) == str
        # test not in bounds
        str = JuMP._math_symbol(REPLMode, :in) * " [0, 1]"
        @test in_set_string(REPLMode, pars[1], set, bounds) == str
        str = JuMP._math_symbol(IJuliaMode, :in) * " [0, 1]"
        @test in_set_string(IJuliaMode, pars[1], set, bounds) == str
    end
    # test in_set_string (InfiniteScalarSet with Bounds)
    @testset "JuMP.in_set_string (InfiniteScalarSet with Bounds)" begin
        # test in bounds
        bounds = ParameterBounds((par1 => IntervalSet(0, 0),))
        set = UniDistributionSet(Uniform())
        str = JuMP._math_symbol(REPLMode, :eq) * " 0"
        @test in_set_string(REPLMode, par1, set, bounds) == str
        str = JuMP._math_symbol(IJuliaMode, :eq) * " 0"
        @test in_set_string(IJuliaMode, par1, set, bounds) == str
        # test in bounds and not equality
        bounds = ParameterBounds((par1 => IntervalSet(0, 1),))
        set = UniDistributionSet(Uniform())
        str = InfiniteOpt._infopt_math_symbol(REPLMode, :prop) * " Uniform " *
              InfiniteOpt._infopt_math_symbol(REPLMode, :intersect) * " [0, 1]"
        @test in_set_string(REPLMode, par1, set, bounds) == str
        str = InfiniteOpt._infopt_math_symbol(IJuliaMode, :prop) * " Uniform " *
              InfiniteOpt._infopt_math_symbol(IJuliaMode, :intersect) * " [0, 1]"
        @test in_set_string(IJuliaMode, par1, set, bounds) == str
        # test not in bounds
        str = InfiniteOpt._infopt_math_symbol(REPLMode, :prop) * " Uniform"
        @test in_set_string(REPLMode, pars[1], set, bounds) == str
        str = InfiniteOpt._infopt_math_symbol(IJuliaMode, :prop) * " Uniform"
        @test in_set_string(IJuliaMode, pars[1], set, bounds) == str
    end
    # test measure_data_string with 1-D DiscreteMeasureData/FunctionalDiscreteMeasureData
    @testset "measure_data_string (1-D)" begin
        # test with bounds
        data = FunctionalDiscreteMeasureData(par1, ones, 0, All, NoGenerativeSupports(),
                                             default_weight, 0, 1, false)
        str = "par1 " * JuMP._math_symbol(REPLMode, :in) * " [0, 1]"
        @test InfiniteOpt.measure_data_string(REPLMode, data) == str
        str = "par1 " * JuMP._math_symbol(IJuliaMode, :in) * " [0, 1]"
        @test InfiniteOpt.measure_data_string(IJuliaMode, data) == str
        # test without bounds
        data = FunctionalDiscreteMeasureData(par1, ones, 0, All, NoGenerativeSupports(),
                                             default_weight, NaN, NaN, false)
        @test InfiniteOpt.measure_data_string(REPLMode, data) == "par1"
        @test InfiniteOpt.measure_data_string(IJuliaMode, data) == "par1"
    end
    # test measure_data_string with Multi-D DiscreteMeasureData/FunctionalDiscreteMeasureData
    @testset "measure_data_string (Multi-D)" begin
        # test with homogenous bounds
        data = FunctionalDiscreteMeasureData(pars2, ones, 0, All, default_weight, [0, 0], [1, 1], false)
        str = "pars2 " * JuMP._math_symbol(REPLMode, :in) * " [0, 1]^2"
        @test InfiniteOpt.measure_data_string(REPLMode, data) == str
        str = "pars2 " * JuMP._math_symbol(IJuliaMode, :in) * " [0, 1]^2"
        @test InfiniteOpt.measure_data_string(IJuliaMode, data) == str
        # test heterogeneous bounds
        data = FunctionalDiscreteMeasureData(pars2, ones, 0, All, default_weight, [0, 0], [0.5, 1], false)
        str = "pars2[1] " * JuMP._math_symbol(REPLMode, :in) * " [0, 0.5], " *
              "pars2[2] " * JuMP._math_symbol(REPLMode, :in) * " [0, 1]"
        @test InfiniteOpt.measure_data_string(REPLMode, data) == str
        str = "pars2_{1} " * JuMP._math_symbol(IJuliaMode, :in) * " [0, 0.5], " *
              "pars2_{2} " * JuMP._math_symbol(IJuliaMode, :in) * " [0, 1]"
        @test InfiniteOpt.measure_data_string(IJuliaMode, data) == str
        # test no bounds with homogenous names
        data = FunctionalDiscreteMeasureData(pars2, ones, 0, All, default_weight, [NaN, NaN], [NaN, NaN], false)
        @test InfiniteOpt.measure_data_string(REPLMode, data) == "pars2"
        @test InfiniteOpt.measure_data_string(IJuliaMode, data) == "pars2"
        # test no bounds with heterogeneous names
        data = FunctionalDiscreteMeasureData([par1, pars[1]], ones, 0, All, default_weight, [NaN, NaN], [NaN, NaN], false)
        @test InfiniteOpt.measure_data_string(REPLMode, data) == "[par1, pars[1]]"
        @test InfiniteOpt.measure_data_string(IJuliaMode, data) == "[par1, pars[1]]"
    end
    # test _get_root_parameter_name
    @testset "_get_root_parameter_name" begin
        # test single
        data = TestData(par1, 0, 1)
        @test InfiniteOpt._get_root_parameter_name(data) == "par1"
        @test InfiniteOpt._get_root_parameter_name(data) == "par1"
        # test with homogenous names
        data = TestData(pars2, 0, 1)
        @test InfiniteOpt._get_root_parameter_name(data) == "pars2"
        @test InfiniteOpt._get_root_parameter_name(data) == "pars2"
        # test with heterogeneous names
        data = TestData([par1, pars[1]], 0, 1)
        @test InfiniteOpt._get_root_parameter_name(data) == "[par1, pars[1]]"
        @test InfiniteOpt._get_root_parameter_name(data) == "[par1, pars[1]]"
    end
    # test measure_data_string (Fallback)
    @testset "measure_data_string (Fallback)" begin
        # test single
        data = TestData(par1, 0, 1)
        @test InfiniteOpt.measure_data_string(REPLMode, data) == "par1"
        @test InfiniteOpt.measure_data_string(IJuliaMode, data) == "par1"
        # test with homogenous names
        data = TestData(pars2, 0, 1)
        @test InfiniteOpt.measure_data_string(REPLMode, data) == "pars2"
        @test InfiniteOpt.measure_data_string(IJuliaMode, data) == "pars2"
        # test with heterogeneous names
        data = TestData([par1, pars[1]], 0, 1)
        @test InfiniteOpt.measure_data_string(REPLMode, data) == "[par1, pars[1]]"
        @test InfiniteOpt.measure_data_string(IJuliaMode, data) == "[par1, pars[1]]"
    end
    # test variabel_string (MeasureRef)
    @testset "variable_string (MeasureRef)" begin
        # test non measure toolbox measures
        data = FunctionalDiscreteMeasureData(pars2, ones, 0, All, default_weight, [0, 0], [1, 1], false)
        meas = dispatch_variable_ref(measure(y, data, name = "test"))
        str = "test{pars2 " * JuMP._math_symbol(REPLMode, :in) * " [0, 1]^2}[y]"
        @test InfiniteOpt.variable_string(REPLMode, meas) == str
        str = "\\text{test}_{pars2 " * JuMP._math_symbol(IJuliaMode, :in) * " [0, 1]^2}\\left[y\\right]"
        @test InfiniteOpt.variable_string(IJuliaMode, meas) == str
        # test measure toolbox special cases for integrals and expectations
        @infinite_parameter(m, t in [0, 1])
        meas = dispatch_variable_ref(expect(y, t))
        str = "\\mathbb{E}_{t}\\left[y\\right]"
        @test InfiniteOpt.variable_string(IJuliaMode, meas) == str
        str = InfiniteOpt._infopt_math_symbol(REPLMode, :expect) * "{t}[y]"
        @test InfiniteOpt.variable_string(REPLMode, meas) == str
        meas = dispatch_variable_ref(integral(y, t))
        str = "\\int_{t " * JuMP._math_symbol(IJuliaMode, :in) * " [0, 1]}ydt"
        @test InfiniteOpt.variable_string(IJuliaMode, meas) == str
        int = InfiniteOpt._infopt_math_symbol(REPLMode, :integral)
        str = int * "{t " * JuMP._math_symbol(REPLMode, :in) * " [0, 1]}[y]"
        @test InfiniteOpt.variable_string(REPLMode, meas) == str
        meas = dispatch_variable_ref(integral(y, par1))
        str = "\\int_{par1 " * JuMP._math_symbol(IJuliaMode, :in) * " [0, 1]}yd(par1)"
        @test InfiniteOpt.variable_string(IJuliaMode, meas) == str
    end
    # test _get_base_name
    @testset "_get_base_name (REPL)" begin
        @test InfiniteOpt._get_base_name(REPLMode, pars[1]) == "pars[1]"
        bad_ref = HoldVariableRef(m, HoldVariableIndex(-1))
        @test InfiniteOpt._get_base_name(REPLMode, bad_ref) == "noname"
    end
    # test _get_base_name
    @testset "_get_base_name (IJulia)" begin
        @test InfiniteOpt._get_base_name(IJuliaMode, pars[1]) == "pars_{1}"
        bad_ref = HoldVariableRef(m, HoldVariableIndex(-1))
        @test InfiniteOpt._get_base_name(IJuliaMode, bad_ref) == "noname"
    end
    # test variable_string (InfiniteVariableRef)
    @testset "variable_string (InfiniteVariableRef)" begin
        # test bad ref
        bad_ref = InfiniteVariableRef(m, InfiniteVariableIndex(-1))
        @test InfiniteOpt.variable_string(REPLMode, bad_ref) == "noname"
        @test InfiniteOpt.variable_string(IJuliaMode, bad_ref) == "noname"
        # test complex normal one
        dvref = dispatch_variable_ref(inf)
        @test InfiniteOpt.variable_string(REPLMode, dvref) == "inf(pars, par1, pars3)"
        @test InfiniteOpt.variable_string(IJuliaMode, dvref) == "inf(pars, par1, pars3)"
        # test mixed types
        @infinite_variable(m, inf2([par1, pars3[2]]))
        dvref = dispatch_variable_ref(inf2)
        @test InfiniteOpt.variable_string(REPLMode, dvref) == "inf2([par1, pars3[2]])"
        @test InfiniteOpt.variable_string(IJuliaMode, dvref) == "inf2([par1, pars3[2]])"
    end
    # test variable_string (ParameterFunctionRef)
    @testset "variable_string (ParameterFunctionRef)" begin
        # test bad ref
        bad_ref = ParameterFunctionRef(m, ParameterFunctionIndex(-1))
        @test InfiniteOpt.variable_string(REPLMode, bad_ref) == "noname"
        @test InfiniteOpt.variable_string(IJuliaMode, bad_ref) == "noname"
        # test normal one
        dvref = dispatch_variable_ref(parameter_function(sin, par1))
        @test InfiniteOpt.variable_string(REPLMode, dvref) == "sin(par1)"
        @test InfiniteOpt.variable_string(IJuliaMode, dvref) == "sin(par1)"
    end
    # test variable_string (DerivativeRef)
    @testset "variable_string (DerivativeRef)" begin
        # test bad ref
        bad_ref = DerivativeRef(m, DerivativeIndex(-1))
        @test InfiniteOpt.variable_string(REPLMode, bad_ref) == "noname"
        @test InfiniteOpt.variable_string(IJuliaMode, bad_ref) == "noname"
        # test normal one
        dref = dispatch_variable_ref(d1)
        d_re = InfiniteOpt._infopt_math_symbol(REPLMode, :partial)
        @test InfiniteOpt.variable_string(REPLMode, dref) == "$d_re/$(d_re)par1[x(par1)]"
        @test InfiniteOpt.variable_string(IJuliaMode, dref) == "\\frac{\\partial}{\\partial par1}\\left[x(par1)\\right]"
        # test nested one
        dref = dispatch_variable_ref(d2)
        expected = "$d_re/$(d_re)par1[$d_re/$(d_re)pars[1][inf(pars, par1, pars3)]]"
        @test InfiniteOpt.variable_string(REPLMode, dref) == expected
        expected = "\\frac{\\partial}{\\partial par1}\\left[\\frac{\\partial}{\\partial pars_{1}}\\left[inf(pars, par1, pars3)\\right]\\right]"
        @test InfiniteOpt.variable_string(IJuliaMode, dref) == expected
        # test has explicit name
        dref = dispatch_variable_ref(d3)
        @test InfiniteOpt.variable_string(REPLMode, dref) == "d3(pars)"
        @test InfiniteOpt.variable_string(IJuliaMode, dref) == "d3(pars)"
    end
    # _make_str_value (Number)
    @testset "_make_str_value (Number)" begin
        @test InfiniteOpt._make_str_value(1.0) == "1"
    end
    # _make_str_value (Array)
    @testset "_make_str_value (Array)" begin
        # test single
        @test InfiniteOpt._make_str_value([1.1]) == "1.1"
        # test short array
        values = [1., 2., 3.]
        @test InfiniteOpt._make_str_value(values) == "[1, 2, 3]"
        # test long array
        values = [1., 2., 3., 4., 5., 6.]
        @test InfiniteOpt._make_str_value(values) == "[1, ..., 6]"
    end
    # test variable_string (PointVariableRef)
    @testset "variable_string (PointVariableRef)" begin
        # test bad ref
        bad_ref = PointVariableRef(m, PointVariableIndex(-1))
        @test InfiniteOpt.variable_string(REPLMode, bad_ref) == "noname"
        @test InfiniteOpt.variable_string(IJuliaMode, bad_ref) == "noname"
        # test short one
        dvref = dispatch_variable_ref(@point_variable(m, x(0)))
        @test InfiniteOpt.variable_string(REPLMode, dvref) == "x(0)"
        @test InfiniteOpt.variable_string(IJuliaMode, dvref) == "x(0)"
        # test complex one
        dvref = dispatch_variable_ref(@point_variable(m, inf([0, 0], 0, [0, 0])))
        @test InfiniteOpt.variable_string(REPLMode, dvref) == "inf([0, 0], 0, [0, 0])"
        @test InfiniteOpt.variable_string(IJuliaMode, dvref) == "inf([0, 0], 0, [0, 0])"
        # test with alias name
        dvref = dispatch_variable_ref(@point_variable(m, z([0, 0]), z0))
        @test InfiniteOpt.variable_string(REPLMode, dvref) == "z0"
        @test InfiniteOpt.variable_string(IJuliaMode, dvref) == "z0"
        # test with derivative
        dvref = dispatch_variable_ref(@point_variable(m, d1(0)))
        d_re = InfiniteOpt._infopt_math_symbol(REPLMode, :partial)
        @test InfiniteOpt.variable_string(REPLMode, dvref) == "$(d_re)/$(d_re)par1[x(par1)](0)"
        # test named derivative
        dvref = dispatch_variable_ref(@point_variable(m, d3([0, 0])))
        @test InfiniteOpt.variable_string(REPLMode, dvref) == "d3([0, 0])"
    end
    # test variable_string (ReducedVariableRef)
    @testset "variable_string (ReducedVariableRef)" begin
        # test bad ref
        bad_ref = ReducedVariableRef(m, ReducedVariableIndex(-1))
        @test InfiniteOpt.variable_string(REPLMode, bad_ref) == "noname"
        @test InfiniteOpt.variable_string(IJuliaMode, bad_ref) == "noname"
        # test short one
        eval_supps = Dict{Int, Float64}(1 => 0, 2 => 0, 4=>0)
        var = build_variable(error, inf, eval_supps, check = false)
        dvref = dispatch_variable_ref(add_variable(m, var))
        @test InfiniteOpt.variable_string(REPLMode, dvref) == "inf([0, 0], par1, [0, pars3[2]])"
        @test InfiniteOpt.variable_string(IJuliaMode, dvref) == "inf([0, 0], par1, [0, pars3[2]])"
        # test named derivative
        eval_supps = Dict{Int, Float64}(1 => 0)
        var = build_variable(error, d3, eval_supps, check = false)
        dvref = dispatch_variable_ref(add_variable(m, var))
        @test InfiniteOpt.variable_string(REPLMode, dvref) == "d3([0, pars[2]])"
        # unnamed derivative
        eval_supps = Dict{Int, Float64}(1 => 0)
        var = build_variable(error, d1, eval_supps, check = false)
        dvref = dispatch_variable_ref(add_variable(m, var))
        d_re = InfiniteOpt._infopt_math_symbol(REPLMode, :partial)
        @test InfiniteOpt.variable_string(REPLMode, dvref) == "$(d_re)/$(d_re)par1[x(par1)](0)"
    end
    # test variable_string (Fallback)
    @testset "variable_string (Fallback)" begin
        @test InfiniteOpt.variable_string(REPLMode, y) == "y"
        @test InfiniteOpt.variable_string(IJuliaMode, y) == "y"
    end
    # test the function_string extensions
    @testset "JuMP.function_string" begin
        @test JuMP.function_string(REPLMode, dispatch_variable_ref(y)) == "y"
        @test JuMP.function_string(IJuliaMode, dispatch_variable_ref(inf)) == "inf(pars, par1, pars3)"
        @test JuMP.function_string(REPLMode, y) == "y"
        @test JuMP.function_string(IJuliaMode, inf) == "inf(pars, par1, pars3)"
        d_re = InfiniteOpt._infopt_math_symbol(REPLMode, :partial)
        @test JuMP.function_string(REPLMode, d1) == "$d_re/$(d_re)par1[x(par1)]"
    end
    # test bound_string
    @testset "bound_string" begin
        # test with single bound
        bounds = ParameterBounds(Dict(par1 => IntervalSet(0.5, 0.7)))
        str = "par1 " * JuMP._math_symbol(REPLMode, :in) * " [0.5, 0.7]"
        @test InfiniteOpt.bound_string(REPLMode, bounds) == str
        str = "par1 " *  JuMP._math_symbol(IJuliaMode, :in) * " [0.5, 0.7]"
        @test InfiniteOpt.bound_string(IJuliaMode, bounds) == str
    end
    # test constraint_string (Finite constraint)
    @testset "JuMP.constraint_string (Finite)" begin
        # test named
        str = "c2 : y² " * JuMP._math_symbol(REPLMode, :eq) * " 3.0"
        @test constraint_string(REPLMode, c2) == str
        str =  "c2 : \$ y^2 " * JuMP._math_symbol(IJuliaMode, :eq) * " 3.0 \$"
        @test constraint_string(IJuliaMode, c2) == str
        # test unnamed
        str = "y² " * JuMP._math_symbol(REPLMode, :eq) * " 3.0"
        @test constraint_string(REPLMode, ac2) == str
        str =  "\$ y^2 " * JuMP._math_symbol(IJuliaMode, :eq) * " 3.0 \$"
        @test constraint_string(IJuliaMode, ac2) == str
        # test named in math mode
        str = "c2 : y² " * JuMP._math_symbol(REPLMode, :eq) * " 3.0"
        @test constraint_string(REPLMode, c2, in_math_mode = true) == str
        str =  "y^2 " * JuMP._math_symbol(IJuliaMode, :eq) * " 3.0"
        @test constraint_string(IJuliaMode, c2, in_math_mode = true) == str
    end
    # test _param_domain_string (IndependentParameter)
    @testset "_param_domain_string (IndependentParameter)" begin
        bounds = ParameterBounds((par1 => IntervalSet(0, 0),))
        idx = index(par1)
        str = "par1 " * JuMP._math_symbol(REPLMode, :eq) * " 0"
        @test InfiniteOpt._param_domain_string(REPLMode, m, idx, bounds) == str
        str = "par1 " * JuMP._math_symbol(IJuliaMode, :eq) * " 0"
        @test InfiniteOpt._param_domain_string(IJuliaMode, m, idx, bounds) == str
    end
    # test _param_domain_string (DependentParameters)
    @testset "_param_domain_string (DependentParameters)" begin
        # Collection set
        bounds = ParameterBounds((pars2[1] => IntervalSet(0, 1),))
        idx = index(pars2[1]).object_index
        str = "pars2[1] " * JuMP._math_symbol(REPLMode, :in) * " [0, 1], " *
              "pars2[2] " * JuMP._math_symbol(REPLMode, :in) * " [0, 2]"
        @test InfiniteOpt._param_domain_string(REPLMode, m, idx, bounds) == str
        str = "pars2_{1} " * JuMP._math_symbol(IJuliaMode, :in) * " [0, 1], " *
              "pars2_{2} " * JuMP._math_symbol(IJuliaMode, :in) * " [0, 2]"
        @test InfiniteOpt._param_domain_string(IJuliaMode, m, idx, bounds) == str
        # other set with equalities
        bounds = ParameterBounds((pars[1] => IntervalSet(0, 0),
                                  pars[2] => IntervalSet(1, 1)))
        idx = index(pars[1]).object_index
        str = InfiniteOpt.bound_string(REPLMode, bounds)
        str2 = string(split(str, ", ")[2], ", ", split(str, ", ")[1])
        @test InfiniteOpt._param_domain_string(REPLMode, m, idx, bounds) in [str, str2]
        str = InfiniteOpt.bound_string(IJuliaMode, bounds)
        str2 = string(split(str, ", ")[2], ", ", split(str, ", ")[1])
        @test InfiniteOpt._param_domain_string(IJuliaMode, m, idx, bounds) in [str, str2]
        # other set without equalities and including in the bounds
        bounds = ParameterBounds((pars[1] => IntervalSet(0, 1),))
        str = "pars " * InfiniteOpt._infopt_math_symbol(REPLMode, :prop) *
              " MvNormal(dim: (2)) " *
              InfiniteOpt._infopt_math_symbol(REPLMode, :intersect) *
              " (pars[1] " * JuMP._math_symbol(REPLMode, :in) * " [0, 1])"
        @test InfiniteOpt._param_domain_string(REPLMode, m, idx, bounds) == str
        str = "pars " * InfiniteOpt._infopt_math_symbol(IJuliaMode, :prop) *
              " MvNormal(dim: (2)) " *
              InfiniteOpt._infopt_math_symbol(IJuliaMode, :intersect) *
              " (pars_{1} " * JuMP._math_symbol(IJuliaMode, :in) * " [0, 1])"
        @test InfiniteOpt._param_domain_string(IJuliaMode, m, idx, bounds) == str
        # other set without equalities and not included in bounds
        bounds = ParameterBounds((par1 => IntervalSet(0, 1),))
        str = "pars " * InfiniteOpt._infopt_math_symbol(REPLMode, :prop) *
              " MvNormal(dim: (2))"
        @test InfiniteOpt._param_domain_string(REPLMode, m, idx, bounds) == str
        str = "pars " * InfiniteOpt._infopt_math_symbol(IJuliaMode, :prop) *
              " MvNormal(dim: (2))"
        @test InfiniteOpt._param_domain_string(IJuliaMode, m, idx, bounds) == str
    end
    # test constraint_string (infinite constraint)
    @testset "JuMP.constraint_string (Infinite)" begin
        # test c1 with name
        str = "c1 : x(par1) + y " * JuMP._math_symbol(REPLMode, :leq) * " 2.0, " *
              JuMP._math_symbol(REPLMode, :for_all) * " par1 " *
              JuMP._math_symbol(REPLMode, :in) * " [0, 1]"
        @test constraint_string(REPLMode, c1) == str
        str = "c1 : \$ x(par1) + y " * JuMP._math_symbol(IJuliaMode, :leq) * " 2.0, " *
              JuMP._math_symbol(IJuliaMode, :for_all) * " par1 " *
              JuMP._math_symbol(IJuliaMode, :in) * " [0, 1] \$"
        @test constraint_string(IJuliaMode, c1) == str
        # test c1 without name
        str = "x(par1) + y " * JuMP._math_symbol(REPLMode, :leq) * " 2.0, " *
              JuMP._math_symbol(REPLMode, :for_all) * " par1 " *
              JuMP._math_symbol(REPLMode, :in) * " [0, 1]"
        @test constraint_string(REPLMode, ac1) == str
        str = "\$ x(par1) + y " * JuMP._math_symbol(IJuliaMode, :leq) * " 2.0, " *
              JuMP._math_symbol(IJuliaMode, :for_all) * " par1 " *
              JuMP._math_symbol(IJuliaMode, :in) * " [0, 1] \$"
        @test constraint_string(IJuliaMode, ac1) == str
        # test c1 with name and in_math_mode
        str = "c1 : x(par1) + y " * JuMP._math_symbol(REPLMode, :leq) * " 2.0, " *
              JuMP._math_symbol(REPLMode, :for_all) * " par1 " *
              JuMP._math_symbol(REPLMode, :in) * " [0, 1]"
        @test constraint_string(REPLMode, c1, in_math_mode = true) == str
        str = "x(par1) + y " * JuMP._math_symbol(IJuliaMode, :leq) * " 2.0, " *
              JuMP._math_symbol(IJuliaMode, :for_all) * " par1 " *
              JuMP._math_symbol(IJuliaMode, :in) * " [0, 1]"
        @test constraint_string(IJuliaMode, c1, in_math_mode = true) == str
        # test c3 with name
        str = "c3 : x(par1) " * JuMP._math_symbol(REPLMode, :eq) * " 5.0, " *
              JuMP._math_symbol(REPLMode, :for_all) * " par1 " *
              JuMP._math_symbol(REPLMode, :in) * " [0, 0.5]"
        @test constraint_string(REPLMode, c3) == str
        str =  "c3 : \$ x(par1) " * JuMP._math_symbol(IJuliaMode, :eq) * " 5.0, " *
               JuMP._math_symbol(IJuliaMode, :for_all) * " par1 " *
               JuMP._math_symbol(IJuliaMode, :in) * " [0, 0.5] \$"
        @test constraint_string(IJuliaMode, c3) == str
        # test c3 without name
        str = "x(par1) " * JuMP._math_symbol(REPLMode, :eq) * " 5.0, " *
              JuMP._math_symbol(REPLMode, :for_all) * " par1 " *
              JuMP._math_symbol(REPLMode, :in) * " [0, 0.5]"
        @test constraint_string(REPLMode, ac3) == str
        str =  "\$ x(par1) " * JuMP._math_symbol(IJuliaMode, :eq) * " 5.0, " *
               JuMP._math_symbol(IJuliaMode, :for_all) * " par1 " *
               JuMP._math_symbol(IJuliaMode, :in) * " [0, 0.5] \$"
        @test constraint_string(IJuliaMode, ac3) == str
    end
    # test constraints_string
    @testset "JuMP.constraints_string" begin
        # test REPLMode
        strings = Vector{String}(undef, 6)
        strings[1] = "c1 : x(par1) + y " * JuMP._math_symbol(REPLMode, :leq) *
                     " 2.0, " * JuMP._math_symbol(REPLMode, :for_all) * " par1 " *
                     JuMP._math_symbol(REPLMode, :in) * " [0, 1]"
        strings[2] = "x(par1) + y " * JuMP._math_symbol(REPLMode, :leq) *
                     " 2.0, " * JuMP._math_symbol(REPLMode, :for_all) * " par1 " *
                     JuMP._math_symbol(REPLMode, :in) * " [0, 1]"
        strings[3] = "c2 : y² " * JuMP._math_symbol(REPLMode, :eq) * " 3.0"
        strings[4] = "y² " * JuMP._math_symbol(REPLMode, :eq) * " 3.0"
        strings[5] = "c3 : x(par1) " * JuMP._math_symbol(REPLMode, :eq) * " 5.0, " *
                     JuMP._math_symbol(REPLMode, :for_all) * " par1 " *
                     JuMP._math_symbol(REPLMode, :in) * " [0, 0.5]"
        strings[6] = "x(par1) " * JuMP._math_symbol(REPLMode, :eq) * " 5.0, " *
                     JuMP._math_symbol(REPLMode, :for_all) * " par1 " *
                     JuMP._math_symbol(REPLMode, :in) * " [0, 0.5]"
        @test constraints_string(REPLMode, m) == strings
        # test IJuliaMode
        strings = Vector{String}(undef, 6)
        strings[1] = "x(par1) + y " * JuMP._math_symbol(IJuliaMode, :leq) *
                     " 2.0, " * JuMP._math_symbol(IJuliaMode, :for_all) * " par1 " *
                     JuMP._math_symbol(IJuliaMode, :in) * " [0, 1]"
        strings[2] = "x(par1) + y " * JuMP._math_symbol(IJuliaMode, :leq) *
                     " 2.0, " * JuMP._math_symbol(IJuliaMode, :for_all) * " par1 " *
                     JuMP._math_symbol(IJuliaMode, :in) * " [0, 1]"
        strings[3] = "y^2 " * JuMP._math_symbol(IJuliaMode, :eq) * " 3.0"
        strings[4] = "y^2 " * JuMP._math_symbol(IJuliaMode, :eq) * " 3.0"
        strings[5] = "x(par1) " * JuMP._math_symbol(IJuliaMode, :eq) * " 5.0, " *
                     JuMP._math_symbol(IJuliaMode, :for_all) * " par1 " *
                     JuMP._math_symbol(IJuliaMode, :in) * " [0, 0.5]"
        strings[6] = "x(par1) " * JuMP._math_symbol(IJuliaMode, :eq) * " 5.0, " *
                     JuMP._math_symbol(IJuliaMode, :for_all) * " par1 " *
                     JuMP._math_symbol(IJuliaMode, :in) * " [0, 0.5]"
        @test constraints_string(IJuliaMode, m) == strings
    end
    # test objective_function_string
    @testset "JuMP.objective_function_string" begin
        @test objective_function_string(REPLMode, m) == "y + 2"
        @test objective_function_string(IJuliaMode, m) == "y + 2"
    end
end

# Helper function to test IO methods work correctly
function show_test(mode, obj, exp_str; repl=:both)
    if mode == REPLMode
        repl != :show  && @test sprint(print, obj) == exp_str
        repl != :print && @test sprint(show,  obj) == exp_str
    else
        @test sprint(show, "text/latex", obj) == exp_str
    end
end

# Make another helper function for other io methods
io_test(f::Function, exp_str::String, args...) = begin
    io = IOBuffer()
    f(io, args...)
    @test String(take!(io)) == exp_str
end

# test infinite model show
@testset "Show InfiniteModel" begin
    # initialize models
    m = InfiniteModel()
    @infinite_parameter(m, 0 <= par1 <= 1)
    @infinite_parameter(m, pars[1:2] in MvNormal([1, 1], 1))
    @infinite_variable(m, x(par1))
    @infinite_variable(m, z(pars))
    @hold_variable(m, y)
    @objective(m, Min, 2 + y)
    @constraint(m, c1, x + y -2 <= 0)
    bounds = ParameterBounds(Dict(par1 => IntervalSet(0.1, 1)))
    @BDconstraint(m, c3(par1 in [0, 0.5]), x <= 5)
    mockoptimizer = () -> MOIU.MockOptimizer(MOIU.UniversalFallback(MOIU.Model{Float64}()),
                                             eval_objective_value=false)
    # test Base.show (IntervalSet in REPL)
    @testset "Base.show (REPL IntervalSet)" begin
        show_test(REPLMode, IntervalSet(0, 1), "[0, 1]")
    end
    # test Base.show (IntervalSet in IJulia)
    @testset "Base.show (IJulia IntervalSet)" begin
        show_test(IJuliaMode, IntervalSet(0, 1), "[0, 1]")
    end
    # test Base.show (DistributionSet in REPL)
    @testset "Base.show (REPL DistributionSet)" begin
        show_test(REPLMode, UniDistributionSet(Uniform()), string(Uniform()))
    end
    # test Base.show (DistributionSet in IJulia)
    @testset "Base.show (IJulia DistributionSet)" begin
        show_test(IJuliaMode, UniDistributionSet(Uniform()), string(Uniform()))
    end
    # test Base.show (CollectionSet in REPL)
    @testset "Base.show (REPL CollectionSet)" begin
        show_test(REPLMode, CollectionSet([IntervalSet(0, 0)]),
                  "CollectionSet with 1 set:\n [0, 0]")
    end
    # test Base.show (CollectionSet in IJulia)
    @testset "Base.show (IJulia CollectionSet)" begin
        show_test(IJuliaMode, CollectionSet([IntervalSet(0, 0), IntervalSet(0, 2)]),
                  "CollectionSet with 2 sets:\n [0, 0]\n [0, 2]")
    end
    # test Base.show (ParameterBounds in REPL)
    @testset "Base.show (REPL ParameterBounds)" begin
        str = "Subdomain bounds (1): par1 " * JuMP._math_symbol(REPLMode, :in) *
              " [0.1, 1]"
        show_test(REPLMode, bounds, str)
    end
    # test Base.show (ParameterBounds in IJulia)
    @testset "Base.show (IJulia ParameterBounds)" begin
        str = "Subdomain bounds (1): par1 " * JuMP._math_symbol(IJuliaMode, :in) *
              " [0.1, 1]"
        show_test(IJuliaMode, bounds, str)
    end
    # test Base.show (GeneralVariableRef in IJulia)
    @testset "Base.show (IJulia GeneralVariableRef)" begin
        show_test(IJuliaMode, y, "\$\$ y \$\$")
    end
    # test Base.show (GeneralVariableRef in REPL)
    @testset "Base.show (REPL GeneralVariableRef)" begin
        show_test(REPLMode, y, "y")
    end
    # test Base.show (constraint in REPL)
    @testset "Base.show (REPL Constraint)" begin
        # test normal
        str = "c1 : x(par1) + y " * JuMP._math_symbol(REPLMode, :leq) *
              " 2.0, " * JuMP._math_symbol(REPLMode, :for_all) * " par1 " *
              JuMP._math_symbol(REPLMode, :in) * " [0, 1]"
        show_test(REPLMode, c1, str)
        # test bounded
        str = "c3 : x(par1) " * JuMP._math_symbol(REPLMode, :leq) * " 5.0, " *
              JuMP._math_symbol(REPLMode, :for_all) * " par1 " *
              JuMP._math_symbol(REPLMode, :in) * " [0, 0.5]"
        show_test(REPLMode, c3, str)
    end
    # test Base.show (constraint in IJulia)
    @testset "Base.show (IJulia Constraint)" begin
        # test normal
        str = "c1 : \$ x(par1) + y " * JuMP._math_symbol(IJuliaMode, :leq) *
              " 2.0, " * JuMP._math_symbol(IJuliaMode, :for_all) * " par1 " *
              JuMP._math_symbol(IJuliaMode, :in) * " [0, 1] \$"
        show_test(IJuliaMode, c1, str)
        # test bounded
        str =  "c3 : \$ x(par1) " * JuMP._math_symbol(IJuliaMode, :leq) * " 5.0, " *
               JuMP._math_symbol(IJuliaMode, :for_all) * " par1 " *
               JuMP._math_symbol(IJuliaMode, :in) * " [0, 0.5] \$"
        show_test(IJuliaMode, c3, str)
    end
    # test show_backend_summary
    @testset "JuMP.show_backend_summary" begin
        # test without optimizer
        str = "Optimizer model backend information: \nModel mode: AUTOMATIC\n" *
              "CachingOptimizer state: NO_OPTIMIZER\nSolver name: No optimizer" *
              " attached."
        io_test(show_backend_summary, str, m)
        # test with optimizer
        set_optimizer(optimizer_model(m), mockoptimizer)
        str = "Optimizer model backend information: \nModel mode: AUTOMATIC\n" *
              "CachingOptimizer state: EMPTY_OPTIMIZER\nSolver name: Mock"
        io_test(show_backend_summary, str, m)
    end
    # test show_objective_function_summary
    @testset "JuMP.show_objective_function_summary" begin
        str = "Objective function type: GenericAffExpr{Float64,GeneralVariableRef}\n"
        io_test(show_objective_function_summary, str, m)
    end
    # test show_constraints_summary
    @testset "JuMP.show_constraints_summary" begin
        # test the main function
        str = "`GenericAffExpr{Float64,GeneralVariableRef}`-in-`MathOptInter" *
              "face.LessThan{Float64}`: 2 constraints\n"
        io_test(show_constraints_summary, str, m)
    end
    # test show_objective_function_summary
    @testset "Base.show (InfiniteModel)" begin
        # test minimization
        str = "An InfiniteOpt Model\nMinimization problem with:\nFinite " *
              "Parameters: 0\nInfinite Parameters: 3\nVariables: 3" *
              "\nDerivatives: 0\nMeasures: 0" *
              "\nObjective function type: GenericAffExpr{Float64,General" *
              "VariableRef}\n`GenericAffExpr{Float64,GeneralVariableRef}`-in-" *
              "`MathOptInterface.LessThan{Float64}`: 2 constraints" *
              "\nNames registered in the model: c1, c3, par1, " *
              "pars, x, y, z\nOptimizer model backend information: \nModel " *
              "mode: AUTOMATIC\nCachingOptimizer state: EMPTY_OPTIMIZER\n" *
              "Solver name: Mock"
        show_test(REPLMode, m, str, repl=:show)
        # test maximization
        set_objective_sense(m, MOI.MAX_SENSE)
        str = "An InfiniteOpt Model\nMaximization problem with:\nFinite " *
              "Parameters: 0\nInfinite Parameters: 3\nVariables: 3" *
              "\nDerivatives: 0\nMeasures: 0" *
              "\nObjective function type: GenericAffExpr{Float64,General" *
              "VariableRef}\n`GenericAffExpr{Float64,GeneralVariableRef}`-in-" *
              "`MathOptInterface.LessThan{Float64}`: 2 constraints" *
              "\nNames registered in the model: c1, c3, par1, " *
              "pars, x, y, z\nOptimizer model backend information: \nModel " *
              "mode: AUTOMATIC\nCachingOptimizer state: EMPTY_OPTIMIZER\n" *
              "Solver name: Mock"
        show_test(REPLMode, m, str, repl=:show)
        # test feasibility
        set_objective_sense(m, MOI.FEASIBILITY_SENSE)
        str = "An InfiniteOpt Model\nFeasibility problem with:\nFinite " *
              "Parameters: 0\nInfinite Parameters: 3\nVariables: 3" *
              "\nDerivatives: 0\nMeasures: 0" *
              "\n`GenericAffExpr{Float64,GeneralVariableRef}`-in-`MathOpt" *
              "Interface.LessThan{Float64}`: 2 constraints" *
              "\nNames registered in the model: c1, c3, par1, " *
              "pars, x, y, z\nOptimizer model backend information: \nModel " *
              "mode: AUTOMATIC\nCachingOptimizer state: EMPTY_OPTIMIZER\n" *
              "Solver name: Mock"
        show_test(REPLMode, m, str, repl=:show)
    end
end


text:
function packit(ti)
    tbname, tbio = mktemp(cleanup = false)
    mktempdir(SystemPath) do path
        rowi = NamedTuple(ti)
        tbl = merge(rowi, (resfile = basename(rowi.resfile), poi_videofile = basename(rowi.poi_videofile), calib_videofile = basename(rowi.calib_videofile)))
        CSV.write(joinpath(path, "csvfile.csv"), [tbl])
        cp(rowi.resfile, joinpath(path, basename(rowi.resfile)))
        map(unique((rowi.poi_videofile, rowi.calib_videofile))) do file
            cp(file, joinpath(path, basename(file)))
        end
        Tar.create(string(path), tbio)
        close(tbio)
        @info "an error has occurred! please send me this file:" tbname
    end
end
function debugging(ti, ex::Exception)
    packit(ti)
    throw(ex)
end
debugging(csvfile, i::Int; delim = nothing) = packit(loadcsv(csvfile, delim)[i])

function debug(tbname)
    Memoization.empty_all_caches!();
    tmp = pwd()
    for file in readdir(tmp, join = true)
        if last(splitext(file)) ≠ ".toml"
            rm(file, force = true, recursive = true)
        end
    end
    tbname = download(tbname)
    files = Tar.extract(tbname)
    map(readdir(files, join = true)) do file
        mv(file, joinpath(tmp, basename(file)))
    end
    process_csv("csvfile.csv")
end


text:
# A Using Case:
# TODO: multiple circuit mode
"""
yao> nbits(4)

yao> begin
    2 => X
    3 => C, 1=>C, 4=>X
end

yao> 2 => C, 5=>X

yao> download_dense()

yao> download_coo()

yao> download_tensornetwork()

yao> download_tensornetwork()

yao> reset()
"""

const WEBCONFIG = Dict{Symbol, Any}(
    :maxn_dense => 10,
    :maxn_sparse => 15,
    :PORT => 8080,
    :LOCALIP => string(Sockets.getipaddr()),
)

mutable struct YaoWSState
    block::ChainBlock
    parseinfo::YaoBlocks.ParseInfo
end
YaoWSState() = YaoWSState(chain(-1), YaoBlocks.ParseInfo(-1, ""))

function yaoweb_parse(str)
      Meta.parse(str)
end

yaoweb_trans(ex, state) = @match ex begin
    :(nbits($n)) => :(setnbits!($n, $state))

    :(download_tensornetwork()) => :(download(Val(:tensornetwork), $(state.block)))
    :(download_dense()) => :(download(Val(:dense), $(state.block)))
    :(download_coo()) => :(download(Val(:coo), $(state.block)))
    :(download_script()) => :(download(Val(:yaoscript), $(state.block)))
    _ => begin
        blk = YaoBlocks.parse_ex(ex, state.parseinfo)
        :(push!($(state.block), $blk))
    end
end

function setnbits!(n::Int, state::YaoWSState)
    state.block = chain(n)
    state.parseinfo.nbit = Int(n)
    "circuit reseted, qubit number is set to $n"
end

function download(::Val{:dense}, b::AbstractBlock{N}) where N
    nmax = WEBCONFIG[:maxn_dense]
    N > nmax && return "Error: circuit size too large to get dense matrix, should be <= $nmax"
    io = IOBuffer()
    writedlm(io, Matrix(b))
    data = String(take!(io))
    return "download::yao-dense.dat::$data"
end

function download(::Val{:coo}, b::AbstractBlock{N}) where N
    nmax = WEBCONFIG[:maxn_sparse]
    N > nmax && return "Error: circuit size too large to get sparse matrix, should be <= $nmax"
    data = join([join((i,j,v), "  ") for (i,j,v) in zip(LuxurySparse.findnz(mat(b))...)], "\n")
    return "download::yao-sparse-coo.dat::$data"
end

function download(::Val{:yaoscript}, b::AbstractBlock)
    data = YaoBlocks.dump_gate(b)
    return "download::yao-script.yao::$data"
end

function download(::Val{:tensornetwork}, b::AbstractBlock)
    data = graph2strings(opennetwork(simplify_blocktypes(b)))
    data = join(data, "::")
    return "download::yao-tensornetwork::$data"
end

function yaoweb_handler(str::AbstractString, state)
    eval(yaoweb_trans(yaoweb_parse(str), state))
end


text:
using Downloads: download

depspath(args...) = joinpath(@__DIR__, args...)

for path in ["anova", "regression", "mcmc", "nonlinear", "univariate"]
    rm(depspath(path), force=true, recursive=true)
    mkpath(depspath(path))
end

anova = ["SiRstv", "SmLs01", "SmLs02", "SmLs03", "AtmWtAg", "SmLs04", "SmLs05", "SmLs06",
         "SmLs07", "SmLs08", "SmLs09"]

lls = ["Norris", "Pontius", "NoInt1", "NoInt2", "Filip", "Longley", "Wampler1",
       "Wampler2", "Wampler3", "Wampler4", "Wampler5"]

mcmc = ["mcmc01", "mcmc02", "mcmc03", "mcmc04", "mcmc05", "mcmc06"]

nls = ["Misra1a", "Chwirut2", "Chwirut1", "Lanczos3", "Gauss1", "Gauss2", "DanWood",
        "Misra1b", "Kirby2", "Hahn1", "Nelson", "MGH17", "Lanczos1", "Lanczos2", "Gauss3",
        "Misra1c", "Misra1d", "Roszman1", "ENSO", "MGH09", "Thurber", "BoxBOD", "Rat42", "MGH10",
        "Eckerle4", "Rat43", "Bennett5"]

univ = ["PiDigits", "Lottery", "Lew", "Mavro", "Michelso", "NumAcc1", "NumAcc2",
        "NumAcc3", "NumAcc4"]

for file in map(x -> x * ".dat", anova)
    download("https://www.itl.nist.gov/div898/strd/anova/$file", depspath("anova", file))
end
for file in map(x -> x * ".dat", lls)
    download("https://www.itl.nist.gov/div898/strd/lls/data/LINKS/DATA/$file", depspath("regression", file))
end
for file in map(x -> x * ".dat", mcmc)
    download("https://www.itl.nist.gov/div898/strd/mcmc/$file", depspath("mcmc", file))
end
for file in map(x -> x * ".dat", nls)
    download("https://www.itl.nist.gov/div898/strd/nls/data/LINKS/DATA/$file", depspath("nonlinear", file))
end
for file in map(x -> x * ".dat", univ)
    download("https://www.itl.nist.gov/div898/strd/univ/data/$file", depspath("univariate", file))
end