text:
function Base.show(io::IO, container::OptimizationContainer)
    show(io, get_jump_model(container))
end

function Base.show(io::IO, ::MIME"text/plain", input::Union{ServiceModel, DeviceModel})
    _show_method(io, input, :auto)
end

function Base.show(io::IO, ::MIME"text/html", input::Union{ServiceModel, DeviceModel})
    _show_method(io, input, :html; standalone=false, tf=PrettyTables.tf_html_simple)
end

function _show_method(
    io::IO,
    model::Union{ServiceModel, DeviceModel},
    backend::Symbol;
    kwargs...,
)
    println(io)
    header = ["Device Type", "Formulation", "Slacks"]

    table = Matrix{String}(undef, 1, length(header))
    table[1, 1] = string(get_component_type(model))
    table[1, 2] = string(get_formulation(model))
    table[1, 3] = string(model.use_slacks)

    PrettyTables.pretty_table(
        io,
        table;
        header=header,
        backend=backend,
        title="Device Model",
        alignment=:l,
        kwargs...,
    )

    if !isempty(model.attributes)
        println(io)
        header = ["Name", "Value"]

        table = Matrix{String}(undef, length(model.attributes), length(header))
        for (ix, (k, v)) in enumerate(model.attributes)
            table[ix, 1] = string(k)
            table[ix, 2] = string(v)
        end

        PrettyTables.pretty_table(
            io,
            table;
            header=header,
            backend=backend,
            title="Attributes",
            alignment=:l,
            kwargs...,
        )
    end

    if !isempty(model.time_series_names)
        println(io)
        header = ["Parameter Name", "Time Series Name"]

        table = Matrix{String}(undef, length(model.time_series_names), length(header))
        for (ix, (k, v)) in enumerate(model.time_series_names)
            table[ix, 1] = string(k)
            table[ix, 2] = string(v)
        end

        PrettyTables.pretty_table(
            io,
            table;
            header=header,
            backend=backend,
            title="Time Series Names",
            alignment=:l,
            kwargs...,
        )
    end

    if !isempty(model.duals)
        println(io)
        table = string.(model.duals)

        PrettyTables.pretty_table(
            io,
            table;
            header=header,
            backend=backend,
            title="Duals",
            alignment=:l,
            kwargs...,
        )
    end

    if !isempty(model.feedforwards)
        println(io)
        table = string.(model.feedforwards)

        PrettyTables.pretty_table(
            io,
            table;
            header=header,
            backend=backend,
            title="Feedforwards",
            alignment=:l,
            kwargs...,
        )
    else
        println(io)
        print(io, "No FeedForwards Assigned")
    end
end

function Base.show(io::IO, ::MIME"text/plain", input::NetworkModel)
    _show_method(io, input, :auto)
end

function Base.show(io::IO, ::MIME"text/html", input::NetworkModel)
    _show_method(io, input, :html; standalone=false, tf=PrettyTables.tf_html_simple)
end

function _show_method(io::IO, network_model::NetworkModel, backend::Symbol; kwargs...)
    table = [
        "Network Model" string(get_network_formulation(network_model))
        "Slacks" get_use_slacks(network_model)
        "PTDF" !isnothing(get_PTDF(network_model))
        "Duals" join(string.(get_duals(network_model)), " ")
    ]

    PrettyTables.pretty_table(
        io,
        table;
        backend=backend,
        header=["Field", "Value"],
        title="Network Model",
        alignment=:l,
        kwargs...,
    )
    return
end

function Base.show(io::IO, ::MIME"text/plain", input::OperationModel)
    _show_method(io, input, :auto)
end

function Base.show(io::IO, ::MIME"text/html", input::OperationModel)
    _show_method(io, input, :html; standalone=false, tf=PrettyTables.tf_html_simple)
end

function _show_method(io::IO, model::OperationModel, backend::Symbol; kwargs...)
    _show_method(io, model.template, backend; kwargs...)
end

function Base.show(io::IO, ::MIME"text/plain", input::ProblemTemplate)
    _show_method(io, input, :auto)
end

function Base.show(io::IO, ::MIME"text/html", input::ProblemTemplate)
    _show_method(io, input, :html; standalone=false, tf=PrettyTables.tf_html_simple)
end

function _show_method(io::IO, template::ProblemTemplate, backend::Symbol; kwargs...)
    table = [
        "Network Model" string(get_network_formulation(template.network_model))
        "Slacks" get_use_slacks(template.network_model)
        "PTDF" !isnothing(get_PTDF(template.network_model))
        "Duals" join(string.(get_duals(template.network_model)), " ")
    ]

    PrettyTables.pretty_table(
        io,
        table;
        backend=backend,
        noheader=true,
        title="Network Model",
        alignment=:l,
        kwargs...,
    )

    println(io)
    header = ["Device Type", "Formulation", "Slacks"]

    table = Matrix{String}(undef, length(template.devices), length(header))
    for (ix, model) in enumerate(values(template.devices))
        table[ix, 1] = string(get_component_type(model))
        table[ix, 2] = string(get_formulation(model))
        table[ix, 3] = string(model.use_slacks)
    end

    PrettyTables.pretty_table(io, table; header=header, title="Device Models", alignment=:l)

    if !isempty(template.branches)
        println(io)
        header = ["Branch Type", "Formulation", "Slacks"]

        table = Matrix{String}(undef, length(template.branches), length(header))
        for (ix, model) in enumerate(values(template.branches))
            table[ix, 1] = string(get_component_type(model))
            table[ix, 2] = string(get_formulation(model))
            table[ix, 3] = string(model.use_slacks)
        end

        PrettyTables.pretty_table(
            io,
            table;
            header=header,
            backend=backend,
            title="Branch Models",
            alignment=:l,
            kwargs...,
        )
    end

    if !isempty(template.services)
        println(io)
        if isempty(first(keys(template.services))[1])
            header = ["Service Type", "Formulation", "Slacks", "Aggregated Model"]
        else
            header = ["Name", "Service Type", "Formulation", "Slacks", "Aggregated Model"]
        end

        table = Matrix{String}(undef, length(template.services), length(header))
        for (ix, (key, model)) in enumerate(template.services)
            if isempty(key[1])
                table[ix, 1] = string(get_component_type(model))
                table[ix, 2] = string(get_formulation(model))
                table[ix, 3] = string(model.use_slacks)
                table[ix, 4] =
                    string(get(model.attributes, "aggregated_service_model", "false"))
            else
                table[ix, 1] = key[1]
                table[ix, 2] = string(get_component_type(model))
                table[ix, 3] = string(get_formulation(model))
                table[ix, 4] = string(model.use_slacks)
                table[ix, 5] =
                    string(get(model.attributes, "aggregated_service_model", "false"))
            end
        end

        PrettyTables.pretty_table(
            io,
            table;
            header=header,
            backend=backend,
            title="Service Models",
            alignment=:l,
            kwargs...,
        )
    end
    return
end

function Base.show(io::IO, ::MIME"text/plain", input::SimulationModels)
    _show_method(io, input, :auto)
end

function Base.show(io::IO, ::MIME"text/html", input::SimulationModels)
    _show_method(io, input, :html; standalone=false, tf=PrettyTables.tf_html_simple)
end

_get_model_type(::DecisionModel{T}) where {T <: DecisionProblem} = T
_get_model_type(::EmulationModel{T}) where {T <: DecisionProblem} = T

function _show_method(io::IO, sim_models::SimulationModels, backend::Symbol; kwargs...)
    println(io)
    header = ["Model Name", "Model Type", "Status", "Output Directory"]

    table = Matrix{Any}(undef, length(sim_models.decision_models), length(header))
    for (ix, model) in enumerate(sim_models.decision_models)
        table[ix, 1] = string(get_name(model))
        table[ix, 2] = string(_get_model_type(model))
        table[ix, 3] = string(get_status(model))
        table[ix, 4] = get_output_dir(model)
    end

    PrettyTables.pretty_table(
        io,
        table;
        header=header,
        backend=backend,
        title="Decision Models",
        alignment=:l,
        kwargs...,
    )

    if !isnothing(sim_models.emulation_model)
        println(io)
        table = Matrix{Any}(undef, 1, length(header))
        table[1, 1] = string(get_name(sim_models.emulation_model))
        table[1, 2] = string(_get_model_type(sim_models.emulation_model))
        table[1, 3] = string(get_status(sim_models.emulation_model))
        table[1, 4] = get_output_dir(sim_models.emulation_model)

        PrettyTables.pretty_table(
            io,
            table;
            header=header,
            backend=backend,
            title="Emulator Models",
            alignment=:l,
            kwargs...,
        )
    else
        println(io)
        println(io, "No Emulator Model Specified")
    end
end

function Base.show(io::IO, ::MIME"text/plain", input::SimulationSequence)
    _show_method(io, input, :auto)
end

function Base.show(io::IO, ::MIME"text/html", input::SimulationSequence)
    _show_method(io, input, :html; standalone=false, tf=PrettyTables.tf_html_simple)
end

function _show_method(io::IO, sequence::SimulationSequence, backend::Symbol; kwargs...)
    println(io)
    table = [
        "Simulation Step Interval" Dates.Hour(get_step_resolution(sequence))
        "Number of Problems" length(sequence.executions_by_model)
    ]

    PrettyTables.pretty_table(
        io,
        table;
        backend=backend,
        noheader=true,
        title="Simulation Sequence",
        alignment=:l,
        kwargs...,
    )

    println(io)
    header = ["Model Name", "Horizon", "Interval", "Executions Per Step"]

    table = Matrix{Any}(undef, length(sequence.executions_by_model), length(header))
    for (ix, (model, executions)) in enumerate(sequence.executions_by_model)
        table[ix, 1] = string(model)
        table[ix, 2] = sequence.horizons[model]
        table[ix, 3] = Dates.Minute(sequence.intervals[model])
        table[ix, 4] = executions
    end

    PrettyTables.pretty_table(
        io,
        table;
        header=header,
        backend=backend,
        title="Simulation Problems",
        alignment=:l,
    )

    if !isempty(sequence.feedforwards)
        println(io)
        header = ["Model Name", "Feed Forward Type"]
        table = Matrix{Any}(undef, length(sequence.feedforwards), length(header))
        for (ix, (k, ff)) in enumerate(sequence.feedforwards)
            table[ix, 1] = k
            table[ix, 2] = join(string.(typeof.(ff)), " ")
        end
        PrettyTables.pretty_table(
            io,
            table;
            header=header,
            backend=backend,
            title="Feedforwards",
            alignment=:l,
            kwargs...,
        )
    end
end

function Base.show(io::IO, ::MIME"text/plain", input::Simulation)
    _show_method(io, input, :auto)
end

function Base.show(io::IO, ::MIME"text/html", input::Simulation)
    _show_method(io, input, :html; standalone=false, tf=PrettyTables.tf_html_simple)
end

function _get_initial_time_for_show(sim::Simulation)
    ini_time = get_initial_time(sim)
    if isnothing(ini_time)
        return "Unset Initial Time"
    else
        return string(ini_time)
    end
end

function _get_build_status_for_show(sim::Simulation)
    internal = sim.internal
    if isnothing(internal)
        return "EMPTY"
    else
        return string(internal.build_status)
    end
end

function _get_run_status_for_show(sim::Simulation)
    internal = sim.internal
    if isnothing(internal)
        return "NOT_READY"
    else
        return string(internal.status)
    end
end

function _show_method(io::IO, sim::Simulation, backend::Symbol; kwargs...)
    table = [
        "Simulation Name" get_name(sim)
        "Build Status" _get_build_status_for_show(sim)
        "Run Status" _get_run_status_for_show(sim)
        "Initial Time" _get_initial_time_for_show(sim)
        "Steps" get_steps(sim)
    ]

    PrettyTables.pretty_table(
        io,
        table;
        backend=backend,
        noheader=true,
        title="Simulation",
        alignment=:l,
        kwargs...,
    )

    _show_method(io, sim.models, backend; kwargs...)
    _show_method(io, sim.sequence, backend; kwargs...)
end

function Base.show(io::IO, ::MIME"text/plain", input::SimulationResults)
    _show_method(io, input, :auto)
end

function Base.show(io::IO, ::MIME"text/html", input::SimulationResults)
    _show_method(io, input, :html; standalone=false, tf=PrettyTables.tf_html_simple)
end

function _show_method(io::IO, results::SimulationResults, backend::Symbol; kwargs...)
    header = ["Problem Name", "Initial Time", "Resolution", "Last Solution Timestamp"]

    table = Matrix{Any}(undef, length(results.decision_problem_results), length(header))
    for (ix, (key, result)) in enumerate(results.decision_problem_results)
        table[ix, 1] = key
        table[ix, 2] = result.timestamps.start
        table[ix, 3] = Dates.Minute(result.timestamps.step)
        table[ix, 4] = result.timestamps.stop
    end
    println(io)
    PrettyTables.pretty_table(
        io,
        table;
        header=header,
        backend=backend,
        title="Decision Problem Results",
        alignment=:l,
    )

    println(io)
    table = [
        "Name" results.emulation_problem_results.problem
        "Resolution" Dates.Minute(results.emulation_problem_results.resolution)
        "Number of steps" length(results.emulation_problem_results.timestamps)
    ]
    PrettyTables.pretty_table(
        io,
        table;
        noheader=true,
        backend=backend,
        title="Emulator Results",
        alignment=:l,
        kwargs...,
    )
end

ProblemResultsTypes = Union{ProblemResults, SimulationProblemResults}
function Base.show(io::IO, ::MIME"text/plain", input::ProblemResultsTypes)
    _show_method(io, input, :auto)
end

function Base.show(io::IO, ::MIME"text/html", input::ProblemResultsTypes)
    _show_method(io, input, :html; standalone=false, tf=PrettyTables.tf_html_simple)
end

function _show_method(
    io::IO,
    results::T,
    backend::Symbol;
    kwargs...,
) where {T <: ProblemResultsTypes}
    timestamps = get_timestamps(results)

    if backend == :html
        println(io, "
 Start: $(timestamps.start)
")
        println(io, "
 End: $(timestamps.stop)
")
        println(io, "
 Resolution: $(Dates.Minute(timestamps.step))
")
    else
        println(io, "Start: $(timestamps.start)")
        println(io, "End: $(timestamps.stop)")
        println(io, "Resolution: $(Dates.Minute(timestamps.step))")
    end

    values = Dict{String, Vector{String}}(
        "Variables" => list_variable_names(results),
        "Auxiliary variables" => list_aux_variable_names(results),
        "Duals" => list_dual_names(results),
        "Expressions" => list_expression_names(results),
        "Parameters" => list_parameter_names(results),
    )

    if hasfield(T, :problem)
        name = results.problem
    else
        name = "PowerSimulations"
    end

    for (k, val) in values
        if !isempty(val)
            println(io)
            PrettyTables.pretty_table(
                io,
                val;
                noheader=true,
                backend=backend,
                title="$name Problem $k Results",
                alignment=:l,
                kwargs...,
            )
        end
    end
end

function Base.show(io::IO, ::MIME"text/plain", bounds::ConstraintBounds)
    println(io, "ConstraintBounds:")
    println(io, "Constraint Coefficient")
    show(io, MIME"text/plain"(), bounds.coefficient)
    println(io, "Constraint RHS")
    show(io, MIME"text/plain"(), bounds.rhs)
end

function Base.show(io::IO, ::MIME"text/plain", bounds::VariableBounds)
    println(io, "VariableBounds:")
    show(io, MIME"text/plain"(), bounds.bounds)
end

function Base.show(io::IO, ::MIME"text/plain", bounds::NumericalBounds)
    println(io, rpad("  Minimum", 20), "Maximum")
    println(io, rpad("  $(bounds.min)", 20), "$(bounds.max)")
end


text:
module DummyPackage

    version() = "v1"

end # module

ext: jl


text:
export getEdgeIntegralOfPolygonalChain,getStraightLineCurrentIntegral

#using Base.BLAS

"""
        function jInv.Mesh.getEdgeIntegralOfPolygonalChain

        s = getEdgeIntegralPolygonalChain(mesh,polygon)
        s = getEdgeIntegralPolygonalChain(mesh,polygon,normalize)

        Compute the integral of a piecewise linear edge grid basis projected onto
        the edges of a polygonal chain. This function can be used to evaluate the
        source term of a line current carried by the polygonal chain in electromagnetic
        modelling.

        The piecewise linear edge grid basis consists of the functions
          phix_i (i = 1 ... numXEdges)
          phiy_i (i = 1 ... numYEdges)
          phiz_i (i = 1 ... numZEdges)
        where phix_i, phiy_i, phiz_i = 1 at the i-th x/y/z-directed edge
          and phix_i, phiy_i, phiz_i = 0 at all other edges.

        INPUT
        mesh ...... Tensor mesh
        polygon ... vertices of polygonal chain as numVertices x 3 array
        normalize . divide line integral by length of integration path (boolean, optional)

        OUTPUT
        s ......... source vector

        For a closed current loop, specify polygon such that
        polygon[1,:] == polygon[end,:]
"""
function getEdgeIntegralOfPolygonalChain(M::TensorMesh3D,polygon::Array{Float64,2};
                                         normalize=false)

#Reformat input
nx, ny, nz =  M.n[1],  M.n[2],  M.n[3]
x0, y0, z0 = M.x0[1], M.x0[2], M.x0[3]
hx, hy, hz = M.h1, M.h2, M.h3
px, py, pz = polygon[:,1], polygon[:,2], polygon[:,3]

# nodal grid
x,y,z = getNodalAxes(M)
# = [x0; x0 + cumsum(hx)]
# y  = [y0; y0 + cumsum(hy)]
# z  = [z0; z0 + cumsum(hz)]

# discrete edge function
sx = zeros(nx  , ny+1, nz+1)
sy = zeros(nx+1, ny  , nz+1)
sz = zeros(nx+1, ny+1, nz  )

# number of line segments
np = length(px) - 1

# check that all polygon vertices are inside the mesh
for ip = 1:np+1
  ax = px[ip]
  ay = py[ip]
  az = pz[ip]
  ix = findlast((ax .>= x[1:nx-1]) .& (ax .<= x[2:nx]))
  iy = findlast((ay .>= y[1:ny-1]) .& (ay .<= y[2:ny]))
  iz = findlast((az .>= z[1:nz-1]) .& (az .<= z[2:nz]))
  if (ix < 1) | (iy < 1) | (iz < 1)
    msg = @sprintf("Polygon vertex (%d,%d,%d) is outside the mesh \n", ax, ay, az)
    error(msg)
  end
end

# integrate each line segment
for ip = 1:np

    # start and end vertices
    ax = px[ip]
    ay = py[ip]
    az = pz[ip]
    bx = px[ip+1]
    by = py[ip+1]
    bz = pz[ip+1]

    # find intersection with mesh planes
    dx = bx - ax
    dy = by - ay
    dz = bz - az
    d  = sqrt(dx^2 + dy^2 + dz^2)
    tol = d * eps(Float64)
    if abs(dx) > tol
      tx = (x .- ax) / dx
      tx = tx[(tx .>=0) .& (tx .<= 1)]
    else
      tx = Float64[]
    end
    if abs(dy) > tol
      ty = (y .- ay) / dy
      ty = ty[(ty .>=0) .& (ty .<= 1)]
    else
      ty = Float64[]
    end
    if abs(dz) > tol
      tz = (z .- az) / dz
      tz = tz[(tz .>=0) .& (tz .<= 1)]
    else
      tz = Float64[]
    end

    t  = sort(unique([-0.0;0.0;tx;ty;tz;1.0]))[2:end]
    nq = length(t) - 1
    tc = 0.5 * (t[1:nq] + t[2:nq+1])

    for iq = 1:nq

        cx = ax + tc[iq] * dx
        cy = ay + tc[iq] * dy
        cz = az + tc[iq] * dz

        # locate cell id
        ix = findlast((cx .>= x[1:nx-1]) .& (cx .<= x[2:nx]))
        iy = findlast((cy .>= y[1:ny-1]) .& (cy .<= y[2:ny]))
        iz = findlast((cz .>= z[1:nz-1]) .& (cz .<= z[2:nz]))

        # local coordinates
        hxloc = hx[ix]
        hyloc = hy[iy]
        hzloc = hz[iz]
        axloc = ax + t[iq]   * dx - x[ix]
        ayloc = ay + t[iq]   * dy - y[iy]
        azloc = az + t[iq]   * dz - z[iz]
        bxloc = ax + t[iq+1] * dx - x[ix]
        byloc = ay + t[iq+1] * dy - y[iy]
        bzloc = az + t[iq+1] * dz - z[iz]
        # integrate
        sxloc,syloc,szloc = getStraightLineCurrentIntegral(hxloc,hyloc,hzloc,axloc,ayloc,azloc,bxloc,byloc,bzloc)

        # integrate
        sx[ix,[iy,iy+1],[iz,iz+1]] = sx[ix,[iy,iy+1],[iz,iz+1]] + reshape(sxloc, (2, 2))
        sy[[ix,ix+1],iy,[iz,iz+1]] = sy[[ix,ix+1],iy,[iz,iz+1]] + reshape(syloc, (2, 2))
        sz[[ix,ix+1],[iy,iy+1],iz] = sz[[ix,ix+1],[iy,iy+1],iz] + reshape(szloc, (2, 2))
    end
end

s = [sx[:]; sy[:]; sz[:]]

# normalize
if normalize
	if all(polygon[1,:] .== polygon[np+1,:])

		# closed polygon: divide by enclosed area
		a  = 0.0
		px = polygon[2:np,1] .- polygon[1,1]
		py = polygon[2:np,2] .- polygon[1,2]
		pz = polygon[2:np,3] .- polygon[1,3]
		for ip = 1:np-2
			cx = py[ip] * pz[ip+1] - pz[ip] * py[ip+1]
			cy = pz[ip] * px[ip+1] - px[ip] * pz[ip+1]
			cz = px[ip] * py[ip+1] - py[ip] * px[ip+1]
			a += sqrt(cx * cx + cy * cy + cz * cz)
		end
		a *= 0.5

	else
		# open polygon: divide by length
		a = 0.0
		for ip = 1:np
			dx = polygon[ip+1,1] - polygon[ip,1]
			dy = polygon[ip+1,2] - polygon[ip,2]
			dz = polygon[ip+1,3] - polygon[ip,3]
			a += sqrt(dx * dx + dy * dy + dz * dz)
		end

	end

	s /= a

end

return s

end

"""
        function jInv.Mesh.getStraightLineCurrentIntegral

        [sx,sy,sz] = getStraightLineCurrentIntegral(hx,hy,hz,ax,ay,az,bx,by,bz)

        Compute integral int(W . J dx^3) in brick of size hx x hy x hz
        where W denotes the 12 local bilinear edge basis functions
        and where J prescribes a unit line current
        between points (ax,ay,az) and (bx,by,bz).
"""
function getStraightLineCurrentIntegral(hx,hy,hz,ax,ay,az,bx,by,bz)

  # length of line segment
  lx = bx - ax
  ly = by - ay
  lz = bz - az
  l  = sqrt(lx^2 + ly^2 + lz^2)

  if l == 0
    sx = zeros(4,1)
    sy = zeros(4,1)
    sz = zeros(4,1)
    return
  end

  # linear interpolation between a and b
  x(t)=ax + t * lx  #0 <= t <= 1
  y(t)=ay + t * ly
  z(t)=az + t * lz

  # edge basis functions
  wx(t)=[
    (1. - y(t) / hy) * (1. - z(t) / hz)
    (    y(t) / hy) * (1. - z(t) / hz)
    (1. - y(t) / hy) * (    z(t) / hz)
    (    y(t) / hy) * (    z(t) / hz)]
  wy(t) = [
    (1. - x(t) / hx) * (1. - z(t) / hz)
    (    x(t) / hx) * (1. - z(t) / hz)
    (1. - x(t) / hx) * (    z(t) / hz)
    (    x(t) / hx) * (    z(t) / hz)]
  wz(t) = [
    (1. - x(t) / hx) * (1. - y(t) / hy)
    (    x(t) / hx) * (1. - y(t) / hy)
    (1. - x(t) / hx) * (    y(t) / hy)
    (    x(t) / hx) * (    y(t) / hy)]

  # integration using Simpson's rule
  sx = (wx(0.) + 4. * wx(0.5) + wx(1.)) * (lx / 6.)
  sy = (wy(0.) + 4. * wy(0.5) + wy(1.)) * (ly / 6.)
  sz = (wz(0.) + 4. * wz(0.5) + wz(1.)) * (lz / 6.)

  return sx, sy, sz
end


text:
##############################################################################
##
## Constants used as NA patterns in R.
##
##############################################################################

"""
    Non-standard addition to NaN bit pattern to discriminate NA from NaN.
    0x000007a2 == UInt32(1954)
    (I assume 1954 is the year of Ross's birth or something like that.)
"""
const R_NA_FLOAT64_LOW = 0x000007a2

## Note that float NA are defined as UInt64 to workaround the win32 ABI
## issue (see JuliaStats/RData.jl#5 and JuliaLang/julia#17195).
if ENDIAN_BOM == 0x01020304
    const R_NA_FLOAT64 = (R_NA_FLOAT64_LOW % UInt64 << 32) | (Base.exponent_mask(Float64) >> 32)
else
    const R_NA_FLOAT64 = Base.exponent_mask(Float64) | R_NA_FLOAT64_LOW
end

# Some .rda files might use R_NA_FLOAT64 = 0x7ff80000000007a2,
# so the proper check for NA: isnan(x) & lowword(x) == R_NA_FLOAT64_LOW
isna_float64(x::UInt64) =
    ((x & Base.exponent_mask(Float64)) == Base.exponent_mask(Float64)) &&
    ((x % UInt32) == R_NA_FLOAT64_LOW)

const R_NA_INT32 = typemin(Int32)
const R_NA_STRING = "NA"

const LONG_VECTOR_SUPPORT = (Sys.WORD_SIZE > 32) # disable long vectors support on 32-bit machines

if LONG_VECTOR_SUPPORT
    typealias RVecLength Int64
else
    typealias RVecLength Int
end

typealias RString UTF8String     # default String container for R string
typealias Hash Dict{RString, Any}

const emptyhash = Hash()
const emptyhashkey = RString("\0")


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

using Base.Test
using FEMBasis

@testset "NSeg interpolate" begin
    basis = NSeg()
    @test size(basis) == (1, 2)
    @test length(basis) == 2
    N = zeros(1, 2)
    eval_basis!(basis, N, (0.0, ))
    @test isapprox(N, [0.5 0.5])
end

@testset "NSurf interpolate" begin
    basis = NSurf()
    @test size(basis) == (2, 4)
    @test length(basis) == 4
    N = zeros(1, 4)
    eval_basis!(basis, N, (0.0, 0.0))
    @test isapprox(N, [0.25 0.25 0.25 0.25])
end

@testset "NSolid interpolate" begin
    basis = NSolid()
    @test size(basis) == (3, 8)
    @test length(basis) == 8
    N = zeros(1, 8)
    eval_basis!(basis, N, (0.0, 0.0, 0.0))
    @test isapprox(N, [0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125])
end


text:
module TestCollect

using Folds
using Folds.Testing: eval_collection_data
using Transducers: PreferParallel
using Test

raw_testdata_executors = """
PreferParallel()
SequentialEx()
ThreadedEx()
ThreadedEx(basesize = 1)
DistributedEx()
DistributedEx(basesize = 1)
DistributedEx(threads_basesize = 1)
DistributedEx(threads_basesize = 1, basesize = 1)
"""

raw_testdata_iterables = """
(x^2 for x in 1:10)
(x^2 for x in 1:10 if isodd(x))
(y for x in 1:10 if isodd(x) for y in 1:x^2)
(x^2 for x in reshape(1:6, 2, 3))
((x, y) for x in 1:3, y in 1:4)                                         # nodist
((x, y) for x in 1:3 for y in 1:4)
((x, y) for x in 1:3, y in 1:4 if isodd(x))                             # nodist
"""

testdata_executors = map(split(raw_testdata_executors, "\n", keepempty = false)) do x
    x => Base.include_string((@__MODULE__), x)
end

tests = eval_collection_data(raw_testdata_iterables)

@testset "$elabel" for (elabel, ex) in testdata_executors
    @testset "$(t.label)" for t in tests
        :nodist ∈ t.tags && ex isa DistributedEx && continue
        itr = t.data
        @test Folds.collect(itr, ex)::Array == collect(itr)
    end
end

end  # module


text:
"""
Provides orbit propagators for the two-body problem,
and the n-body problem.
"""
module Propagators

using Reexport
@reexport using ..AstrodynamicsCore

include("../Misc/DocStringExtensions.jl")
include("../Misc/UnitfulAliases.jl")

using ..NBody
using ..TwoBody
using ..ThreeBody

using DifferentialEquations
using LinearAlgebra: norm, normalize, cross, ×, dot, ⋅
using ComponentArrays
using StaticArrays: SVector, @SVector, SMatrix, @SMatrix

export  Trajectory,
        propagate,
        RestrictedTwoBodyTic!,
        RestrictedBiasedTwoBodyTic!,
        RestrictedThreeBodyTic!,
        NBodyTic!,
        show

include("Trajectory.jl")
include("PropagateTwoBody.jl")
include("PropagateThreeBody.jl")
include("PropagateNBody.jl")

end


text:
# Test operator methods
@testset "Internal Methods" begin
    # initialize model and references
    m = InfiniteModel()
    @infinite_parameter(m, 0 <= par <= 1)
    @infinite_variable(m, inf(par))
    @point_variable(m, inf(0), pt)
    @hold_variable(m, hold)
    meas = MeasureRef(m, 1)
    m.meas_to_name[1] = "meas"
    # test _var_type_parser
    @testset "_var_type_parser" begin
        # test same types
        @test InfiniteOpt._var_type_parser(MeasureRef, MeasureRef) == MeasureRef
        @test InfiniteOpt._var_type_parser(ParameterRef,
                                           ParameterRef) == ParameterRef
        @test InfiniteOpt._var_type_parser(InfiniteVariableRef,
                                     InfiniteVariableRef) == InfiniteVariableRef
        @test InfiniteOpt._var_type_parser(PointVariableRef,
                                           PointVariableRef) == PointVariableRef
        @test InfiniteOpt._var_type_parser(HoldVariableRef,
                                         HoldVariableRef) == HoldVariableRef
        @test InfiniteOpt._var_type_parser(ReducedInfiniteVariableRef,
                       ReducedInfiniteVariableRef) == ReducedInfiniteVariableRef
        # test finite types
        @test InfiniteOpt._var_type_parser(PointVariableRef,
                                         HoldVariableRef) == FiniteVariableRef
        @test InfiniteOpt._var_type_parser(HoldVariableRef,
                                         FiniteVariableRef) == FiniteVariableRef
        # test measure finite types
        @test InfiniteOpt._var_type_parser(PointVariableRef,
                                         MeasureRef) == MeasureFiniteVariableRef
        @test InfiniteOpt._var_type_parser(HoldVariableRef,
                                         MeasureRef) == MeasureFiniteVariableRef
        # test other combos
        @test InfiniteOpt._var_type_parser(InfiniteVariableRef,
                                               MeasureRef) == GeneralVariableRef
        @test InfiniteOpt._var_type_parser(ParameterRef,
                                           MeasureRef) == GeneralVariableRef
        @test InfiniteOpt._var_type_parser(InfiniteVariableRef,
                                           ParameterRef) == GeneralVariableRef
        @test InfiniteOpt._var_type_parser(InfiniteVariableRef,
                                         PointVariableRef) == GeneralVariableRef
        @test InfiniteOpt._var_type_parser(ParameterRef,
                                         PointVariableRef) == GeneralVariableRef
        @test InfiniteOpt._var_type_parser(InfiniteVariableRef,
                                        HoldVariableRef) == GeneralVariableRef
        @test InfiniteOpt._var_type_parser(ParameterRef,
                                        HoldVariableRef) == GeneralVariableRef
    end
    # test extension to add_to_expression! (mixed variable quadratic addition)
    @testset "JuMP.add_to_expression! (Mixed)" begin
        # make expression
        quad = pt^2 + 3hold + 2
        # test with same types
        @test isa(add_to_expression!(copy(quad), 1., hold, pt),
                  GenericQuadExpr{Float64, FiniteVariableRef})
        pair = UnorderedPair{FiniteVariableRef}(hold, pt)
        @test add_to_expression!(copy(quad), 1., hold, pt).terms[pair] == 1
        # test with different types
        @test isa(add_to_expression!(copy(quad), 2., hold, inf),
                  GenericQuadExpr{Float64, GeneralVariableRef})
        pair = UnorderedPair{GeneralVariableRef}(hold, inf)
        @test add_to_expression!(copy(quad), 2., hold, inf).terms[pair] == 2
    end
    # test extension to add_to_expression! (variable + number quadratic addition)
    @testset "JuMP.add_to_expression! (Number1)" begin
        # make expression
        quad = pt^2 + 3hold + pt + 3
        # test with same types
        @test isa(add_to_expression!(copy(quad), 1., 2, pt),
                  GenericQuadExpr{Float64, FiniteVariableRef})
        @test add_to_expression!(copy(quad), 1., 2, pt).aff.terms[pt] == 3
        # test with different types
        @test isa(add_to_expression!(copy(quad), 2., 3, inf),
                  GenericQuadExpr{Float64, GeneralVariableRef})
        @test add_to_expression!(copy(quad), 2., 3, inf).aff.terms[inf] == 6
    end
    # test extension to add_to_expression! (number + variable quadratic addition)
    @testset "JuMP.add_to_expression! (Number2)" begin
        # make expression
        quad = pt^2 + 3hold + pt + 3
        # test with same types
        @test isa(add_to_expression!(copy(quad), 1., pt, 2),
                  GenericQuadExpr{Float64, FiniteVariableRef})
        @test add_to_expression!(copy(quad), 1., pt, 2).aff.terms[pt] == 3
        # test with different types
        @test isa(add_to_expression!(copy(quad), 2., inf, 2),
                  GenericQuadExpr{Float64, GeneralVariableRef})
        @test add_to_expression!(copy(quad), 2., inf, 2).aff.terms[inf] == 4
    end
    # test extension to add_to_expression! (number + number quadratic addition)
    @testset "JuMP.add_to_expression! (Number3)" begin
        # make expression
        quad = pt^2 + 3hold + pt + 3
        # test
        @test isa(add_to_expression!(copy(quad), 1., 2, 2),
                  GenericQuadExpr{Float64, FiniteVariableRef})
        @test add_to_expression!(copy(quad), 1., 2, 2).aff.constant == 7
        # another test
        @test isa(add_to_expression!(copy(quad), 2., -3, 2),
                  GenericQuadExpr{Float64, FiniteVariableRef})
        @test add_to_expression!(copy(quad), 2., -3, 2).aff.constant == -9
    end
end

# Test operations between 2 different variables
@testset "Variable--Variable" begin
    # initialize model and references
    m = InfiniteModel()
    @infinite_parameter(m, 0 <= par <= 1)
    @infinite_variable(m, inf(par))
    @point_variable(m, inf(0), pt)
    @hold_variable(m, hold)
    meas = MeasureRef(m, 1)
    m.meas_to_name[1] = "meas"
    # test addition
    @testset "Base.:+" begin
        # test infinite + measure
        @test isa(inf + meas, GenericAffExpr{Float64, GeneralVariableRef})
        @test (inf + meas).constant == 0.0
        @test (inf + meas).terms[meas] == 1
        @test (inf + meas).terms[inf] == 1
        # test point + hold
        @test isa(pt + hold, GenericAffExpr{Float64, FiniteVariableRef})
        @test (pt + hold).terms[hold] == 1
        # test other combos
        @test isa(pt + meas, GenericAffExpr{Float64, MeasureFiniteVariableRef})
        @test isa(par + pt, GenericAffExpr{Float64, GeneralVariableRef})
        @test isa(par + inf, GenericAffExpr{Float64, GeneralVariableRef})
        @test isa(par + meas, GenericAffExpr{Float64, GeneralVariableRef})
        # test same
        @test isa(meas + meas, GenericAffExpr{Float64, MeasureRef})
        @test (meas + meas).terms[meas] == 2
        @test isa(inf + inf, GenericAffExpr{Float64, InfiniteVariableRef})
        @test isa(pt + pt, GenericAffExpr{Float64, PointVariableRef})
    end
    # test subtraction
    @testset "Base.:-" begin
        # test hold - measure
        @test isa(hold - meas, GenericAffExpr{Float64, MeasureFiniteVariableRef})
        @test (hold - meas).constant == 0.0
        @test (hold - meas).terms[meas] == -1
        @test (hold - meas).terms[hold] == 1
        # test point - hold
        @test isa(pt - hold, GenericAffExpr{Float64, FiniteVariableRef})
        @test (pt - hold).terms[hold] == -1
        # test other combos
        @test isa(pt - meas, GenericAffExpr{Float64, MeasureFiniteVariableRef})
        @test isa(par - pt, GenericAffExpr{Float64, GeneralVariableRef})
        @test isa(par - inf, GenericAffExpr{Float64, GeneralVariableRef})
        @test isa(par - meas, GenericAffExpr{Float64, GeneralVariableRef})
        # test same
        @test isa(meas - meas, GenericAffExpr{Float64, MeasureRef})
        @test (meas - meas) == zero(GenericAffExpr{Float64, MeasureRef})
        @test isa(inf - inf, GenericAffExpr{Float64, InfiniteVariableRef})
        @test isa(pt - pt, GenericAffExpr{Float64, PointVariableRef})
    end
    # test multiplication
    @testset "Base.:*" begin
        # test point * measure
        @test isa(pt * meas, GenericQuadExpr{Float64, MeasureFiniteVariableRef})
        @test (pt * meas).aff.constant == 0.0
        pair = UnorderedPair{MeasureFiniteVariableRef}(pt, meas)
        @test (pt * meas).terms[pair] == 1
        # test point * hold
        @test isa(pt * hold, GenericQuadExpr{Float64, FiniteVariableRef})
        pair = UnorderedPair{FiniteVariableRef}(pt, hold)
        @test (pt * hold).terms[pair] == 1
        # test other combos
        @test isa(pt * meas, GenericQuadExpr{Float64, MeasureFiniteVariableRef})
        @test isa(par * pt, GenericQuadExpr{Float64, GeneralVariableRef})
        @test isa(par * inf, GenericQuadExpr{Float64, GeneralVariableRef})
        @test isa(par * meas, GenericQuadExpr{Float64, GeneralVariableRef})
        # test same
        @test isa(meas * meas, GenericQuadExpr{Float64, MeasureRef})
        pair = UnorderedPair{MeasureRef}(meas, meas)
        @test (meas * meas).terms[pair] == 1
        @test isa(inf * inf, GenericQuadExpr{Float64, InfiniteVariableRef})
        @test isa(pt * pt, GenericQuadExpr{Float64, PointVariableRef})
    end
    # test division
    @testset "Base.:/" begin
        # test some combos
        @test_throws ErrorException inf / par
        @test_throws ErrorException inf / inf
        @test_throws ErrorException pt / hold
    end
end

# Test operations between variables and GenericAffExpr
@testset "Variable--AffExpr" begin
    # initialize model and references and expressions
    m = InfiniteModel()
    @infinite_parameter(m, 0 <= par <= 1)
    @infinite_variable(m, inf(par))
    @point_variable(m, inf(0), pt)
    @hold_variable(m, hold)
    meas = MeasureRef(m, 1)
    m.meas_to_name[1] = "meas"
    aff1 = (inf + pt) - 2
    aff2 = hold - pt
    aff3 = (hold - meas) + 1
    # test addition
    @testset "Base.:+" begin
        # test infinite + aff1
        @test isa(inf + aff1, GenericAffExpr{Float64, GeneralVariableRef})
        @test (inf + aff1).constant == -2
        @test (inf + aff1).terms[pt] == 1
        @test (inf + aff1).terms[inf] == 2
        # test meas + aff2
        @test isa(meas + aff2, GenericAffExpr{Float64, MeasureFiniteVariableRef})
        @test (meas + aff2).terms[hold] == 1
        # test other combos
        @test isa(hold + aff3, GenericAffExpr{Float64, MeasureFiniteVariableRef})
        @test isa(par + aff2, GenericAffExpr{Float64, GeneralVariableRef})
        @test isa(par + aff1, GenericAffExpr{Float64, GeneralVariableRef})
        @test isa(par + aff3, GenericAffExpr{Float64, GeneralVariableRef})
        # test same
        @test isa(meas + (meas + meas), GenericAffExpr{Float64, MeasureRef})
        @test (meas + (meas + meas)).terms[meas] == 3
    end
    # test subtraction
    @testset "Base.:-" begin
        # test infinite - aff1
        @test isa(inf - aff1, GenericAffExpr{Float64, GeneralVariableRef})
        @test (inf - aff1).constant == 2
        @test (inf - aff1).terms[pt] == -1
        @test (inf - aff1).terms[inf] == 0
        # test meas + aff2
        @test isa(meas - aff2, GenericAffExpr{Float64, MeasureFiniteVariableRef})
        @test (meas - aff2).terms[hold] == -1
        # test other combos
        @test isa(hold - aff3, GenericAffExpr{Float64, MeasureFiniteVariableRef})
        @test isa(par - aff2, GenericAffExpr{Float64, GeneralVariableRef})
        @test isa(par - aff1, GenericAffExpr{Float64, GeneralVariableRef})
        @test isa(par - aff3, GenericAffExpr{Float64, GeneralVariableRef})
        # test same
        @test isa(meas - (meas + meas), GenericAffExpr{Float64, MeasureRef})
        @test (meas - (meas + meas)).terms[meas] == -1
    end
    # test multiplication
    @testset "Base.:*" begin
        # test infinite * aff1
        @test isa(inf * aff1, GenericQuadExpr{Float64, GeneralVariableRef})
        @test (inf * aff1).aff.constant == 0
        @test (inf * aff1).aff.terms[inf] == -2
        pair = UnorderedPair{GeneralVariableRef}(inf, pt)
        @test (inf * aff1).terms[pair] == 1
        pair = UnorderedPair{GeneralVariableRef}(inf, inf)
        @test (inf * aff1).terms[pair] == 1
        # test meas * aff2
        @test isa(meas * aff2, GenericQuadExpr{Float64, MeasureFiniteVariableRef})
        pair = UnorderedPair{MeasureFiniteVariableRef}(meas, hold)
        @test (meas * aff2).terms[pair] == 1
        @test length((meas * aff2).aff.terms) == 0
        # test other combos
        @test isa(hold * aff3, GenericQuadExpr{Float64, MeasureFiniteVariableRef})
        @test isa(par * aff2, GenericQuadExpr{Float64, GeneralVariableRef})
        @test isa(par * aff1, GenericQuadExpr{Float64, GeneralVariableRef})
        @test isa(par * aff3, GenericQuadExpr{Float64, GeneralVariableRef})
        # test same
        @test isa(meas * (meas + meas), GenericQuadExpr{Float64, MeasureRef})
        pair = UnorderedPair{MeasureRef}(meas, meas)
        @test (meas * (meas + meas)).terms[pair] == 2
    end
    # test division
    @testset "Base.:/" begin
        # test some combos
        @test_throws ErrorException inf / aff1
        @test_throws ErrorException inf / aff2
        @test_throws ErrorException pt / aff3
    end
end

# Test operations between GenericAffExpr and variable
@testset "AffExpr--Variable" begin
    # initialize model and references and expressions
    m = InfiniteModel()
    @infinite_parameter(m, 0 <= par <= 1)
    @infinite_variable(m, inf(par))
    @point_variable(m, inf(0), pt)
    @hold_variable(m, hold)
    meas = MeasureRef(m, 1)
    m.meas_to_name[1] = "meas"
    aff1 = (inf + pt) - 2
    aff2 = hold - pt
    aff3 = (hold - meas) + 1
    # test addition
    @testset "Base.:+" begin
        # test aff1 + inf
        @test isa(copy(aff1) + inf, GenericAffExpr{Float64, GeneralVariableRef})
        @test (copy(aff1) + inf).constant == -2
        @test (copy(aff1) + inf).terms[pt] == 1
        @test (copy(aff1) + inf).terms[inf] == 2
        # test aff2 + meas
        @test isa(copy(aff2) + meas, GenericAffExpr{Float64, MeasureFiniteVariableRef})
        @test (copy(aff2) + meas).terms[hold] == 1
        # test other combos
        @test isa(aff3 + hold, GenericAffExpr{Float64, MeasureFiniteVariableRef})
        @test isa(aff2 + par, GenericAffExpr{Float64, GeneralVariableRef})
        @test isa(aff1 + par, GenericAffExpr{Float64, GeneralVariableRef})
        @test isa(aff3 + par, GenericAffExpr{Float64, GeneralVariableRef})
        # test same
        @test isa((meas + meas) + meas, GenericAffExpr{Float64, MeasureRef})
        @test ((meas + meas) + meas).terms[meas] == 3
    end
    # test subtraction
    @testset "Base.:-" begin
        # test aff1 - inf
        @test isa(copy(aff1) - inf, GenericAffExpr{Float64, GeneralVariableRef})
        @test (copy(aff1) - inf).constant == -2
        @test (copy(aff1) - inf).terms[pt] == 1
        @test (copy(aff1) - inf).terms[inf] == 0
        # test aff2 - meas
        @test isa(copy(aff2) - meas, GenericAffExpr{Float64, MeasureFiniteVariableRef})
        @test (copy(aff2) - meas).terms[hold] == 1
        # test other combos
        @test isa(aff3 - hold, GenericAffExpr{Float64, MeasureFiniteVariableRef})
        @test isa(aff2 - par, GenericAffExpr{Float64, GeneralVariableRef})
        @test isa(aff1 - par, GenericAffExpr{Float64, GeneralVariableRef})
        @test isa(aff3 - par, GenericAffExpr{Float64, GeneralVariableRef})
        # test same
        @test isa((meas + meas) - meas, GenericAffExpr{Float64, MeasureRef})
        @test ((meas + meas) - meas).terms[meas] == 1
    end
    # test multiplication
    @testset "Base.:*" begin
        # test aff1 * inf
        @test isa(aff1 * inf, GenericQuadExpr{Float64, GeneralVariableRef})
        @test (aff1 * inf).aff.constant == 0
        @test (aff1 * inf).aff.terms[inf] == -2
        pair = UnorderedPair{GeneralVariableRef}(inf, pt)
        @test (aff1 * inf).terms[pair] == 1
        pair = UnorderedPair{GeneralVariableRef}(inf, inf)
        @test (aff1 * inf).terms[pair] == 1
        # test aff2 * meas
        @test isa(aff2 * meas, GenericQuadExpr{Float64, MeasureFiniteVariableRef})
        pair = UnorderedPair{MeasureFiniteVariableRef}(meas, hold)
        @test (aff2 * meas).terms[pair] == 1
        @test length((aff2 * meas).aff.terms) == 0
        # test other combos
        @test isa(aff3 * hold, GenericQuadExpr{Float64, MeasureFiniteVariableRef})
        @test isa(aff2 * par, GenericQuadExpr{Float64, GeneralVariableRef})
        @test isa(aff1 * par, GenericQuadExpr{Float64, GeneralVariableRef})
        @test isa(aff3 * par, GenericQuadExpr{Float64, GeneralVariableRef})
        # test same
        @test isa((meas + meas * meas), GenericQuadExpr{Float64, MeasureRef})
        pair = UnorderedPair{MeasureRef}(meas, meas)
        @test ((meas + meas) * meas).terms[pair] == 2
    end
    # test division
    @testset "Base.:/" begin
        # test some combos
        @test_throws ErrorException aff1 / inf
        @test_throws ErrorException aff2 / inf
        @test_throws ErrorException aff3 / pt
    end
end

# Test operators on GenericAffExpr and GenericAffExpr
@testset "AffExpr--AffExpr" begin
    # initialize model and references and expressions
    m = InfiniteModel()
    @infinite_parameter(m, 0 <= par <= 1)
    @infinite_variable(m, inf(par))
    @point_variable(m, inf(0), pt)
    @hold_variable(m, hold)
    @hold_variable(m, test[1:51])
    meas = MeasureRef(m, 1)
    m.meas_to_name[1] = "meas"
    aff1 = (inf + pt) - 2
    aff2 = hold - pt
    aff3 = (hold - meas) + 1
    aff4 = (inf + inf) + 3
    aff5 = 4pt - 42
    aff_big = sum(test[i] for i = 1:length(test))
    # test addition
    @testset "Base.:+" begin
        # test large input
        @test isa(aff_big + aff1, GenericAffExpr{Float64, GeneralVariableRef})
        @test isa(copy(aff1) + aff_big, GenericAffExpr{Float64, GeneralVariableRef})
        # test various mixed combination types
        @test isa(copy(aff1) + aff2, GenericAffExpr{Float64, GeneralVariableRef})
        @test isa(aff2 + aff3, GenericAffExpr{Float64, MeasureFiniteVariableRef})
        @test isa(copy(aff2) + aff5, GenericAffExpr{Float64, FiniteVariableRef})
        @test isa(aff4 + aff1, GenericAffExpr{Float64, GeneralVariableRef})
        # test results of general combinations
        @test (copy(aff1) + aff2).constant == -2
        @test (copy(aff1) + aff2).terms[inf] == 1
        @test (copy(aff1) + aff2).terms[pt] == 0
        @test (aff4 + aff1).constant == 1
        @test (aff4 + aff1).terms[inf] == 3
        @test (aff4 + aff1).terms[pt] == 1
        # test results of measurefinite combinations
        @test (aff2 + aff3).constant == 1
        @test (aff2 + aff3).terms[meas] == -1
        @test (aff2 + aff3).terms[hold] == 2
        @test (copy(aff3) + aff5).constant == -41
        @test (copy(aff3) + aff5).terms[meas] == -1
        @test (copy(aff3) + aff5).terms[pt] == 4
        # test results of fininite combinations
        @test (copy(aff2) + aff5).constant == -42
        @test (copy(aff2) + aff5).terms[pt] == 3
        @test (copy(aff2) + aff5).terms[hold] == 1
        # test pure additions
        @test isa(copy(aff4) + aff4, GenericAffExpr{Float64, InfiniteVariableRef})
        @test (copy(aff4) + aff4).constant == 6
        @test (copy(aff4) + aff4).terms[inf] == 4
        @test isa(copy(aff5) + aff5, GenericAffExpr{Float64, PointVariableRef})
        @test (copy(aff5) + aff5).constant == -84
        @test (copy(aff5) + aff5).terms[pt] == 8
    end
    # test subtraction
    @testset "Base.:-" begin
        # test various mixed combination types
        @test isa(copy(aff1) - aff2, GenericAffExpr{Float64, GeneralVariableRef})
        @test isa(aff2 - aff3, GenericAffExpr{Float64, MeasureFiniteVariableRef})
        @test isa(copy(aff2) - aff5, GenericAffExpr{Float64, FiniteVariableRef})
        @test isa(aff4 - aff1, GenericAffExpr{Float64, GeneralVariableRef})
        # test results of general combinations
        @test (copy(aff1) - aff2).constant == -2
        @test (copy(aff1) - aff2).terms[inf] == 1
        @test (copy(aff1) - aff2).terms[pt] == 2
        @test (aff4 - aff1).constant == 5
        @test (aff4 - aff1).terms[inf] == 1
        @test (aff4 - aff1).terms[pt] == -1
        # test results of measurefinite combinations
        @test (aff2 - aff3).constant == -1
        @test (aff2 - aff3).terms[meas] == 1
        @test (aff2 - aff3).terms[hold] == 0
        @test (copy(aff3) - aff5).constant == 43
        @test (copy(aff3) - aff5).terms[meas] == -1
        @test (copy(aff3) - aff5).terms[pt] == -4
        # test results of fininite combinations
        @test (copy(aff2) - aff5).constant == 42
        @test (copy(aff2) - aff5).terms[pt] == -5
        @test (copy(aff2) - aff5).terms[hold] == 1
        # test pure additions
        @test isa(copy(aff4) - aff4, GenericAffExpr{Float64, InfiniteVariableRef})
        @test (copy(aff4) - aff4).constant == 0
        @test (copy(aff4) - aff4).terms[inf] == 0
        @test isa(copy(aff5) - aff5, GenericAffExpr{Float64, PointVariableRef})
        @test (copy(aff5) - aff5).constant == 0
        @test (copy(aff5) - aff5).terms[pt] == 0
    end
    # test multiplication
    @testset "Base.:*" begin
        # test various mixed combination types
        @test isa(aff1 * aff2, GenericQuadExpr{Float64, GeneralVariableRef})
        @test isa(aff2 * aff3, GenericQuadExpr{Float64, MeasureFiniteVariableRef})
        @test isa(aff2 * aff5, GenericQuadExpr{Float64, FiniteVariableRef})
        @test isa(aff4 * aff1, GenericQuadExpr{Float64, GeneralVariableRef})
        # test results of general combinations
        @test (aff1 * aff2).aff.constant == 0
        @test (aff1 * aff2).aff.terms[hold] == -2
        pair = UnorderedPair{GeneralVariableRef}(inf, hold)
        @test (aff1 * aff2).terms[pair] == 1
        pair = UnorderedPair{GeneralVariableRef}(pt, pt)
        @test (aff1 * aff2).terms[pair] == -1
        @test (aff4 * aff1).aff.constant == -6
        @test (aff4 * aff1).aff.terms[inf] == -1
        pair = UnorderedPair{GeneralVariableRef}(inf, inf)
        @test (aff4 * aff1).terms[pair] == 2
        pair = UnorderedPair{GeneralVariableRef}(inf, pt)
        @test (aff4 * aff1).terms[pair] == 2
        # test results of measurefinite combinations
        @test (aff2 * aff3).aff.constant == 0
        @test (aff2 * aff3).aff.terms[pt] == -1
        pair = UnorderedPair{MeasureFiniteVariableRef}(pt, meas)
        @test (aff2 * aff3).terms[pair] == 1
        pair = UnorderedPair{MeasureFiniteVariableRef}(hold, hold)
        @test (aff2 * aff3).terms[pair] == 1
        @test (aff3 * aff5).aff.constant == -42
        @test (aff3 * aff5).aff.terms[meas] == 42
        pair = UnorderedPair{MeasureFiniteVariableRef}(hold, pt)
        @test (aff3 * aff5).terms[pair] == 4
        pair = UnorderedPair{MeasureFiniteVariableRef}(meas, pt)
        @test (aff3 * aff5).terms[pair] == -4
        # test results of fininite combinations
        @test (aff2 * aff5).aff.constant == 0
        @test (aff2 * aff5).aff.terms[hold] == -42
        pair = UnorderedPair{FiniteVariableRef}(hold, pt)
        @test (aff2 * aff5).terms[pair] == 4
        pair = UnorderedPair{FiniteVariableRef}(pt, pt)
        @test (aff2 * aff5).terms[pair] == -4
        # test pure additions
        @test isa(aff4 * aff4, GenericQuadExpr{Float64, InfiniteVariableRef})
        @test (aff4 * aff4).aff.constant == 9
        @test (aff4 * aff4).aff.terms[inf] == 12
        pair = UnorderedPair{InfiniteVariableRef}(inf, inf)
        @test (aff4 * aff4).terms[pair] == 4
        @test isa(aff5 * aff5, GenericQuadExpr{Float64, PointVariableRef})
        @test (aff5 * aff5).aff.constant == 42 * 42
        @test (aff5 * aff5).aff.terms[pt] == 2 * 4 * -42
        pair = UnorderedPair{PointVariableRef}(pt, pt)
        @test (aff5 * aff5).terms[pair] == 16
    end
    # test division
    @testset "Base.:/" begin
        # test some combos
        @test_throws ErrorException aff1 / aff2
        @test_throws ErrorException aff2 / aff4
        @test_throws ErrorException aff3 / aff5
        @test_throws ErrorException aff5 / aff5
    end
end

# Test operations for GenericQuadExpr--GeneralVariableRef
@testset "QuadExpr--Variable" begin
    # initialize model and references and expressions
    m = InfiniteModel()
    @infinite_parameter(m, 0 <= par <= 1)
    @infinite_variable(m, inf(par))
    @point_variable(m, inf(0), pt)
    @hold_variable(m, hold)
    @hold_variable(m, test[1:51])
    meas = MeasureRef(m, 1)
    m.meas_to_name[1] = "meas"
    aff1 = (inf + pt) - 2
    aff2 = hold - pt
    aff3 = (hold - meas) + 1
    aff5 = 4pt - 42
    quad1 = copy(aff1) * aff2 #(inf*hold - inf*pt + pt*hold - pt^2 - 2hold + 2pt)
    quad2 = (pt * aff3) - 3   #(pt*hold - pt*meas + pt - 3)
    quad3 = hold * aff2       #(hold^2 - hold*pt)
    quad4 = pt * aff5         #(4pt^2 - 42pt)
    # test addition
    @testset "Base.:+" begin
        # test mixed combo types
        @test isa(copy(quad1) + inf, GenericQuadExpr{Float64, GeneralVariableRef})
        @test isa(copy(quad2) + meas, GenericQuadExpr{Float64, MeasureFiniteVariableRef})
        @test isa(copy(quad3) + pt, GenericQuadExpr{Float64, FiniteVariableRef})
        @test isa(copy(quad4) + par, GenericQuadExpr{Float64, GeneralVariableRef})
        # test general combos
        @test (copy(quad1) + inf).aff.constant == 0
        @test (copy(quad1) + inf).aff.terms[pt] == 2
        @test (copy(quad1) + inf).aff.terms[inf] == 1
        pair = UnorderedPair{GeneralVariableRef}(inf, hold)
        @test (copy(quad1) + inf).terms[pair] == 1
        @test (copy(quad4) + par).aff.constant == 0
        @test (copy(quad4) + par).aff.terms[pt] == -42
        @test (copy(quad4) + par).aff.terms[par] == 1
        pair = UnorderedPair{GeneralVariableRef}(pt, pt)
        @test (copy(quad4) + par).terms[pair] == 4
        # test measurefinite combos
        @test (copy(quad2) + meas).aff.constant == -3
        @test (copy(quad2) + meas).aff.terms[pt] == 1
        @test (copy(quad2) + meas).aff.terms[meas] == 1
        pair = UnorderedPair{GeneralVariableRef}(pt, hold)
        @test (copy(quad2) + meas).terms[pair] == 1
        @test (copy(quad2) + pt).aff.constant == -3
        @test (copy(quad2) + pt).aff.terms[pt] == 2
        # test finite combos
        @test (copy(quad3) + pt).aff.constant == 0
        @test (copy(quad3) + pt).aff.terms[pt] == 1
        pair = UnorderedPair{GeneralVariableRef}(hold, pt)
        @test (copy(quad3) + pt).terms[pair] == -1
        # test same types
        @test isa(copy(quad4) + pt, GenericQuadExpr{Float64, PointVariableRef})
        @test (copy(quad4) + pt).aff.constant == 0
        @test (copy(quad4) + pt).aff.terms[pt] == -41
        pair = UnorderedPair{GeneralVariableRef}(pt, pt)
        @test (copy(quad4) + pt).terms[pair] == 4
    end
    # test subtraction
    @testset "Base.:-" begin
        # test mixed combo types
        @test isa(copy(quad1) - inf, GenericQuadExpr{Float64, GeneralVariableRef})
        @test isa(copy(quad2) - meas, GenericQuadExpr{Float64, MeasureFiniteVariableRef})
        @test isa(copy(quad3) - pt, GenericQuadExpr{Float64, FiniteVariableRef})
        @test isa(copy(quad4) - par, GenericQuadExpr{Float64, GeneralVariableRef})
        # test general combos
        @test (copy(quad1) - inf).aff.constant == 0
        @test (copy(quad1) - inf).aff.terms[pt] == 2
        @test (copy(quad1) - inf).aff.terms[inf] == -1
        pair = UnorderedPair{GeneralVariableRef}(inf, hold)
        @test (copy(quad1) - inf).terms[pair] == 1
        @test (copy(quad4) - par).aff.constant == 0
        @test (copy(quad4) - par).aff.terms[pt] == -42
        @test (copy(quad4) - par).aff.terms[par] == -1
        pair = UnorderedPair{GeneralVariableRef}(pt, pt)
        @test (copy(quad4) - par).terms[pair] == 4
        # test measurefinite combos
        @test (copy(quad2) - meas).aff.constant == -3
        @test (copy(quad2) - meas).aff.terms[pt] == 1
        @test (copy(quad2) - meas).aff.terms[meas] == -1
        pair = UnorderedPair{GeneralVariableRef}(pt, hold)
        @test (copy(quad2) - meas).terms[pair] == 1
        @test (copy(quad2) - pt).aff.constant == -3
        @test (copy(quad2) - pt).aff.terms[pt] == 0
        # test finite combos
        @test (copy(quad3) - pt).aff.constant == 0
        @test (copy(quad3) - pt).aff.terms[pt] == -1
        pair = UnorderedPair{GeneralVariableRef}(hold, pt)
        @test (copy(quad3) - pt).terms[pair] == -1
        # test same types
        @test isa(copy(quad4) - pt, GenericQuadExpr{Float64, PointVariableRef})
        @test (copy(quad4) - pt).aff.constant == 0
        @test (copy(quad4) - pt).aff.terms[pt] == -43
        pair = UnorderedPair{GeneralVariableRef}(pt, pt)
        @test (copy(quad4) - pt).terms[pair] == 4
    end
    # test multiplication
    @testset "Base.:*" begin
        # test some combos
        @test_throws ErrorException quad1 * inf
        @test_throws ErrorException quad2 * meas
        @test_throws ErrorException quad3 * par
        @test_throws ErrorException quad4 * pt
    end
    # test division
    @testset "Base.:/" begin
        # test some combos
        @test_throws ErrorException quad1 / inf
        @test_throws ErrorException quad2 / meas
        @test_throws ErrorException quad3 / par
        @test_throws ErrorException quad4 / pt
    end
end

# Test operations for GeneralVariableRef--GenericQuadExpr
@testset "Variable--QuadExpr" begin
    # initialize model and references and expressions
    m = InfiniteModel()
    @infinite_parameter(m, 0 <= par <= 1)
    @infinite_variable(m, inf(par))
    @point_variable(m, inf(0), pt)
    @hold_variable(m, hold)
    @hold_variable(m, test[1:51])
    meas = MeasureRef(m, 1)
    m.meas_to_name[1] = "meas"
    aff1 = (inf + pt) - 2
    aff2 = hold - pt
    aff3 = (hold - meas) + 1
    aff5 = 4pt - 42
    quad1 = copy(aff1) * aff2 #(inf*hold - inf*pt + pt*hold - pt^2 - 2hold + 2pt)
    quad2 = (pt * aff3) - 3   #(pt*hold - pt*meas + pt - 3)
    quad3 = hold * aff2       #(hold^2 - hold*pt)
    quad4 = pt * aff5         #(4pt^2 - 42pt)
    # test addition
    @testset "Base.:+" begin
        # test mixed combo types
        @test isa(inf + copy(quad1), GenericQuadExpr{Float64, GeneralVariableRef})
        @test isa(meas + copy(quad2), GenericQuadExpr{Float64, MeasureFiniteVariableRef})
        @test isa(pt + copy(quad3), GenericQuadExpr{Float64, FiniteVariableRef})
        @test isa(par + copy(quad4), GenericQuadExpr{Float64, GeneralVariableRef})
        # test general combos
        @test (inf + copy(quad1)).aff.constant == 0
        @test (inf + copy(quad1)).aff.terms[pt] == 2
        @test (inf + copy(quad1)).aff.terms[inf] == 1
        pair = UnorderedPair{GeneralVariableRef}(inf, hold)
        @test (inf + copy(quad1)).terms[pair] == 1
        @test (par + copy(quad4)).aff.constant == 0
        @test (par + copy(quad4)).aff.terms[pt] == -42
        @test (par + copy(quad4)).aff.terms[par] == 1
        pair = UnorderedPair{GeneralVariableRef}(pt, pt)
        @test (par + copy(quad4)).terms[pair] == 4
        # test measurefinite combos
        @test (meas + copy(quad2)).aff.constant == -3
        @test (meas + copy(quad2)).aff.terms[pt] == 1
        @test (meas + copy(quad2)).aff.terms[meas] == 1
        pair = UnorderedPair{GeneralVariableRef}(pt, hold)
        @test (meas + copy(quad2)).terms[pair] == 1
        @test (pt + copy(quad2)).aff.constant == -3
        @test (pt + copy(quad2)).aff.terms[pt] == 2
        # test finite combos
        @test (pt + copy(quad3)).aff.constant == 0
        @test (pt + copy(quad3)).aff.terms[pt] == 1
        pair = UnorderedPair{GeneralVariableRef}(hold, pt)
        @test (pt + copy(quad3)).terms[pair] == -1
        # test same types
        @test isa(pt + copy(quad4), GenericQuadExpr{Float64, PointVariableRef})
        @test (pt + copy(quad4)).aff.constant == 0
        @test (pt + copy(quad4)).aff.terms[pt] == -41
        pair = UnorderedPair{GeneralVariableRef}(pt, pt)
        @test (pt + copy(quad4)).terms[pair] == 4
    end
    # test subtraction
    @testset "Base.:-" begin
        # test mixed combo types
        @test isa(inf - copy(quad1), GenericQuadExpr{Float64, GeneralVariableRef})
        @test isa(meas - copy(quad2), GenericQuadExpr{Float64, MeasureFiniteVariableRef})
        @test isa(pt - copy(quad3), GenericQuadExpr{Float64, FiniteVariableRef})
        @test isa(par - copy(quad4), GenericQuadExpr{Float64, GeneralVariableRef})
        # test general combos
        @test (inf - copy(quad1)).aff.constant == 0
        @test (inf - copy(quad1)).aff.terms[pt] == -2
        @test (inf - copy(quad1)).aff.terms[inf] == 1
        pair = UnorderedPair{GeneralVariableRef}(inf, hold)
        @test (inf - copy(quad1)).terms[pair] == -1
        @test (par - copy(quad4)).aff.constant == 0
        @test (par - copy(quad4)).aff.terms[pt] == 42
        @test (par - copy(quad4)).aff.terms[par] == 1
        pair = UnorderedPair{GeneralVariableRef}(pt, pt)
        @test (par - copy(quad4)).terms[pair] == -4
        # test measurefinite combos
        @test (meas - copy(quad2)).aff.constant == 3
        @test (meas - copy(quad2)).aff.terms[pt] == -1
        @test (meas - copy(quad2)).aff.terms[meas] == 1
        pair = UnorderedPair{GeneralVariableRef}(pt, hold)
        @test (meas - copy(quad2)).terms[pair] == -1
        @test (pt - copy(quad2)).aff.constant == 3
        @test (pt - copy(quad2)).aff.terms[pt] == 0
        # test finite combos
        @test (pt - copy(quad3)).aff.constant == 0
        @test (pt - copy(quad3)).aff.terms[pt] == 1
        pair = UnorderedPair{GeneralVariableRef}(hold, pt)
        @test (pt - copy(quad3)).terms[pair] == 1
        # test same types
        @test isa(pt - copy(quad4), GenericQuadExpr{Float64, PointVariableRef})
        @test (pt - copy(quad4)).aff.constant == 0
        @test (pt - copy(quad4)).aff.terms[pt] == 43
        pair = UnorderedPair{GeneralVariableRef}(pt, pt)
        @test (pt - copy(quad4)).terms[pair] == -4
    end
    # test multiplication
    @testset "Base.:*" begin
        # test some combos
        @test_throws ErrorException inf * quad1
        @test_throws ErrorException meas * quad2
        @test_throws ErrorException par * quad3
        @test_throws ErrorException pt * quad4
    end
    # test division
    @testset "Base.:/" begin
        # test some combos
        @test_throws ErrorException inf / quad1
        @test_throws ErrorException meas / quad2
        @test_throws ErrorException par / quad3
        @test_throws ErrorException pt / quad4
    end
end

# Test operators on GenericAffExpr--GenericQuadExpr
@testset "AffExpr--QuadExpr" begin
    # initialize model and references and expressions
    m = InfiniteModel()
    @infinite_parameter(m, 0 <= par <= 1)
    @infinite_variable(m, inf(par))
    @point_variable(m, inf(0), pt)
    @hold_variable(m, hold)
    @hold_variable(m, test[1:51])
    meas = MeasureRef(m, 1)
    m.meas_to_name[1] = "meas"
    aff1 = (inf + pt) - 2
    aff2 = hold - pt
    aff3 = (hold - meas) + 1
    aff5 = 4pt - 42
    quad1 = copy(aff1) * aff2 #(inf*hold - inf*pt + pt*hold - pt^2 - 2hold + 2pt)
    quad2 = (pt * aff3) - 3   #(pt*hold - pt*meas + pt - 3)
    quad3 = hold * aff2       #(hold^2 - hold*pt)
    quad4 = pt * aff5         #(4pt^2 - 42pt)
    # test addition
    @testset "Base.:+" begin
        # test mixed combo types
        @test isa(copy(aff1) + quad1, GenericQuadExpr{Float64, GeneralVariableRef})
        @test isa(copy(aff2) + quad2, GenericQuadExpr{Float64, MeasureFiniteVariableRef})
        @test isa(copy(aff5) + quad3, GenericQuadExpr{Float64, FiniteVariableRef})
        @test isa(copy(aff3) + quad1, GenericQuadExpr{Float64, GeneralVariableRef})
        # test general combos
        @test (copy(aff1) + quad1).aff.constant == -2
        @test (copy(aff1) + quad1).aff.terms[pt] == 3
        @test (copy(aff1) + quad1).aff.terms[inf] == 1
        pair = UnorderedPair{GeneralVariableRef}(inf, hold)
        @test (copy(aff1) + quad1).terms[pair] == 1
        pair = UnorderedPair{GeneralVariableRef}(pt, pt)
        @test (copy(aff1) + quad1).terms[pair] == -1
        @test (copy(aff3) + quad1).aff.constant == 1
        @test (copy(aff3) + quad1).aff.terms[meas] == -1
        @test (copy(aff3) + quad1).aff.terms[hold] == -1
        pair = UnorderedPair{GeneralVariableRef}(inf, hold)
        @test (copy(aff3) + quad1).terms[pair] == 1
        # test measurefinite combos
        @test (copy(aff3) + quad2).aff.constant == -2
        @test (copy(aff3) + quad2).aff.terms[pt] == 1
        @test (copy(aff3) + quad2).aff.terms[hold] == 1
        pair = UnorderedPair{GeneralVariableRef}(pt, hold)
        @test (copy(aff3) + quad2).terms[pair] == 1
        @test (copy(aff3) + quad4).aff.constant == 1
        @test (copy(aff3) + quad4).aff.terms[pt] == -42
        @test (copy(aff3) + quad4).aff.terms[meas] == -1
        pair = UnorderedPair{GeneralVariableRef}(pt, pt)
        @test (copy(aff3) + quad4).terms[pair] == 4
        # test finite combos
        @test (copy(aff2) + quad3).aff.constant == 0
        @test (copy(aff2) + quad3).aff.terms[pt] == -1
        @test (copy(aff2) + quad3).aff.terms[hold] == 1
        pair = UnorderedPair{GeneralVariableRef}(hold, pt)
        @test (copy(aff2) + quad3).terms[pair] == -1
        # test same types
        @test isa(copy(aff5) + quad4, GenericQuadExpr{Float64, PointVariableRef})
        @test (copy(aff5) + quad4).aff.constant == -42
        @test (copy(aff5) + quad4).aff.terms[pt] == -38
        pair = UnorderedPair{GeneralVariableRef}(pt, pt)
        @test (copy(aff5) + quad4).terms[pair] == 4
    end
    # test subtraction
    @testset "Base.:-" begin
        # test mixed combo types
        @test isa(copy(aff1) - quad1, GenericQuadExpr{Float64, GeneralVariableRef})
        @test isa(copy(aff2) - quad2, GenericQuadExpr{Float64, MeasureFiniteVariableRef})
        @test isa(copy(aff5) - quad3, GenericQuadExpr{Float64, FiniteVariableRef})
        @test isa(copy(aff3) - quad1, GenericQuadExpr{Float64, GeneralVariableRef})
        # test general combos
        @test (copy(aff1) - quad1).aff.constant == -2
        @test (copy(aff1) - quad1).aff.terms[pt] == -1
        @test (copy(aff1) - quad1).aff.terms[inf] == 1
        pair = UnorderedPair{GeneralVariableRef}(inf, hold)
        @test (copy(aff1) - quad1).terms[pair] == -1
        pair = UnorderedPair{GeneralVariableRef}(pt, pt)
        @test (copy(aff1) - quad1).terms[pair] == 1
        @test (copy(aff3) - quad1).aff.constant == 1
        @test (copy(aff3) - quad1).aff.terms[meas] == -1
        @test (copy(aff3) - quad1).aff.terms[hold] == 3
        pair = UnorderedPair{GeneralVariableRef}(inf, hold)
        @test (copy(aff3) - quad1).terms[pair] == -1
        # test measurefinite combos
        @test (copy(aff3) - quad2).aff.constant == 4
        @test (copy(aff3) - quad2).aff.terms[pt] == -1
        @test (copy(aff3) - quad2).aff.terms[hold] == 1
        pair = UnorderedPair{GeneralVariableRef}(pt, hold)
        @test (copy(aff3) - quad2).terms[pair] == -1
        @test (copy(aff3) - quad4).aff.constant == 1
        @test (copy(aff3) - quad4).aff.terms[pt] == 42
        @test (copy(aff3) - quad4).aff.terms[meas] == -1
        pair = UnorderedPair{GeneralVariableRef}(pt, pt)
        @test (copy(aff3) - quad4).terms[pair] == -4
        # test finite combos
        @test (copy(aff2) - quad3).aff.constant == 0
        @test (copy(aff2) - quad3).aff.terms[pt] == -1
        @test (copy(aff2) - quad3).aff.terms[hold] == 1
        pair = UnorderedPair{GeneralVariableRef}(hold, pt)
        @test (copy(aff2) - quad3).terms[pair] == 1
        # test same types
        @test isa(copy(aff5) - quad4, GenericQuadExpr{Float64, PointVariableRef})
        @test (copy(aff5) - quad4).aff.constant == -42
        @test (copy(aff5) - quad4).aff.terms[pt] == 46
        pair = UnorderedPair{GeneralVariableRef}(pt, pt)
        @test (copy(aff5) - quad4).terms[pair] == -4
    end
    # test multiplication
    @testset "Base.:*" begin
        # test some combos
        @test_throws ErrorException aff1 * quad1
        @test_throws ErrorException aff2 * quad3
        @test_throws ErrorException aff3 * quad4
        @test_throws ErrorException aff5 * quad2
    end
    # test division
    @testset "Base.:/" begin
        # test some combos
        @test_throws ErrorException aff1 / quad1
        @test_throws ErrorException aff2 / quad3
        @test_throws ErrorException aff3 / quad4
        @test_throws ErrorException aff5 / quad2
    end
end

# Test operators on GenericQuadExpr--GenericAffExpr
@testset "QuadExpr--AffExpr" begin
    # initialize model and references and expressions
    m = InfiniteModel()
    @infinite_parameter(m, 0 <= par <= 1)
    @infinite_variable(m, inf(par))
    @point_variable(m, inf(0), pt)
    @hold_variable(m, hold)
    @hold_variable(m, test[1:51])
    meas = MeasureRef(m, 1)
    m.meas_to_name[1] = "meas"
    aff1 = (inf + pt) - 2
    aff2 = hold - pt
    aff3 = (hold - meas) + 1
    aff5 = 4pt - 42
    quad1 = copy(aff1) * aff2 #(inf*hold - inf*pt + pt*hold - pt^2 - 2hold + 2pt)
    quad2 = (pt * aff3) - 3   #(pt*hold - pt*meas + pt - 3)
    quad3 = hold * aff2       #(hold^2 - hold*pt)
    quad4 = pt * aff5         #(4pt^2 - 42pt)
    # test addition
    @testset "Base.:+" begin
        # test mixed combo types
        @test isa(copy(quad1) + copy(aff1), GenericQuadExpr{Float64, GeneralVariableRef})
        @test isa(copy(quad2) + copy(aff2) + quad2, GenericQuadExpr{Float64, MeasureFiniteVariableRef})
        @test isa(copy(quad3) + copy(aff5), GenericQuadExpr{Float64, FiniteVariableRef})
        @test isa(copy(quad1) + copy(aff3), GenericQuadExpr{Float64, GeneralVariableRef})
        # test general combos
        @test (copy(quad1)+ copy(aff1)).aff.constant == -2
        @test (copy(quad1) + copy(aff1)).aff.terms[pt] == 3
        @test (copy(quad1) + copy(aff1)).aff.terms[inf] == 1
        pair = UnorderedPair{GeneralVariableRef}(inf, hold)
        @test (copy(quad1) + copy(aff1)).terms[pair] == 1
        pair = UnorderedPair{GeneralVariableRef}(pt, pt)
        @test (copy(quad1) + copy(aff3)).terms[pair] == -1
        @test (copy(quad1) + copy(aff3)).aff.constant == 1
        @test (copy(quad1) + copy(aff3)).aff.terms[meas] == -1
        @test (copy(quad1) + copy(aff3)).aff.terms[hold] == -1
        pair = UnorderedPair{GeneralVariableRef}(inf, hold)
        @test (copy(quad1) + copy(aff3)).terms[pair] == 1
        # test measurefinite combos
        @test (copy(quad2) + copy(aff3)).aff.constant == -2
        @test (copy(quad2) + copy(aff3)).aff.terms[pt] == 1
        @test (copy(quad2) + copy(aff3)).aff.terms[hold] == 1
        pair = UnorderedPair{GeneralVariableRef}(pt, hold)
        @test (copy(quad2) + copy(aff3)).terms[pair] == 1
        @test (copy(quad4) + copy(aff3)).aff.constant == 1
        @test (copy(quad4) + copy(aff3)).aff.terms[pt] == -42
        @test (copy(quad4) + copy(aff3)).aff.terms[meas] == -1
        pair = UnorderedPair{GeneralVariableRef}(pt, pt)
        @test (copy(quad4) + copy(aff3)).terms[pair] == 4
        # test finite combos
        @test (copy(quad3) + copy(aff2)).aff.constant == 0
        @test (copy(quad3) + copy(aff2)).aff.terms[pt] == -1
        @test (copy(quad3) + copy(aff2)).aff.terms[hold] == 1
        pair = UnorderedPair{GeneralVariableRef}(hold, pt)
        @test (copy(quad3) + copy(aff2)).terms[pair] == -1
        # test same types
        @test isa(copy(quad4) + copy(aff5), GenericQuadExpr{Float64, PointVariableRef})
        @test (copy(quad4) + copy(aff5)).aff.constant == -42
        @test (copy(quad4) + copy(aff5)).aff.terms[pt] == -38
        pair = UnorderedPair{GeneralVariableRef}(pt, pt)
        @test (copy(quad4) + copy(aff5)).terms[pair] == 4
    end
    # test subtraction
    @testset "Base.:-" begin
        # test mixed combo types
        @test isa(copy(quad1) - copy(aff1), GenericQuadExpr{Float64, GeneralVariableRef})
        @test isa(copy(quad2) - copy(aff2), GenericQuadExpr{Float64, MeasureFiniteVariableRef})
        @test isa(copy(quad3) - copy(aff5), GenericQuadExpr{Float64, FiniteVariableRef})
        @test isa(copy(quad1) - copy(aff3), GenericQuadExpr{Float64, GeneralVariableRef})
        # test general combos
        @test (copy(quad1) - copy(aff1)).aff.constant == 2
        @test (copy(quad1) - copy(aff1)).aff.terms[pt] == 1
        @test (copy(quad1) - copy(aff1)).aff.terms[inf] == -1
        pair = UnorderedPair{GeneralVariableRef}(inf, hold)
        @test (copy(quad1) - copy(aff1)).terms[pair] == 1
        pair = UnorderedPair{GeneralVariableRef}(pt, pt)
        @test (copy(quad1) - copy(aff1)).terms[pair] == -1
        @test (copy(quad1) - copy(aff3)).aff.constant == -1
        @test (copy(quad1) - copy(aff3)).aff.terms[meas] == 1
        @test (copy(quad1) - copy(aff3)).aff.terms[hold] == -3
        pair = UnorderedPair{GeneralVariableRef}(inf, hold)
        @test (copy(quad1) - copy(aff3)).terms[pair] == 1
        # test measurefinite combos
        @test (copy(quad2) - copy(aff3)).aff.constant == -4
        @test (copy(quad2) - copy(aff3)).aff.terms[pt] == 1
        @test (copy(quad2) - copy(aff3)).aff.terms[hold] == -1
        pair = UnorderedPair{GeneralVariableRef}(pt, hold)
        @test (copy(quad2) - copy(aff3)).terms[pair] == 1
        @test (copy(quad4) - copy(aff3)).aff.constant == -1
        @test (copy(quad4) - copy(aff3)).aff.terms[pt] == -42
        @test (copy(quad4) - copy(aff3)).aff.terms[meas] == 1
        pair = UnorderedPair{GeneralVariableRef}(pt, pt)
        @test (copy(quad4)- copy(aff3)).terms[pair] == 4
        # test finite combos
        @test (copy(quad3) - copy(aff2)).aff.constant == 0
        @test (copy(quad3) - copy(aff2)).aff.terms[pt] == 1
        @test (copy(quad3) - copy(aff2)).aff.terms[hold] == -1
        pair = UnorderedPair{GeneralVariableRef}(hold, pt)
        @test (copy(quad3) - copy(aff2)).terms[pair] == -1
        # test same types
        @test isa(copy(quad4) - copy(aff5), GenericQuadExpr{Float64, PointVariableRef})
        @test (copy(quad4) - copy(aff5)).aff.constant == 42
        @test (copy(quad4) - copy(aff5)).aff.terms[pt] == -46
        pair = UnorderedPair{GeneralVariableRef}(pt, pt)
        @test (copy(quad4) - copy(aff5)).terms[pair] == 4
    end
    # test multiplication
    @testset "Base.:*" begin
        # test some combos
        @test_throws ErrorException quad1 * aff1
        @test_throws ErrorException quad3 * aff2
        @test_throws ErrorException quad4 * aff3
        @test_throws ErrorException quad2 * aff5
    end
    # test division
    @testset "Base.:/" begin
        # test some combos
        @test_throws ErrorException quad1 / aff1
        @test_throws ErrorException quad3 / aff2
        @test_throws ErrorException quad4 / aff3
        @test_throws ErrorException quad2 / aff5
    end
end

# Test operators on GenericQuadExpr--GenericQuadExpr
@testset "QuadExpr--QuadExpr" begin
    # initialize model and references and expressions
    m = InfiniteModel()
    @infinite_parameter(m, 0 <= par <= 1)
    @infinite_variable(m, inf(par))
    @point_variable(m, inf(0), pt)
    @hold_variable(m, hold)
    @hold_variable(m, test[1:51])
    meas = MeasureRef(m, 1)
    m.meas_to_name[1] = "meas"
    aff1 = (inf + pt) - 2
    aff2 = hold - pt
    aff3 = (hold - meas) + 1
    aff5 = 4pt - 42
    quad1 = copy(aff1) * aff2 #(inf*hold - inf*pt + pt*hold - pt^2 - 2hold + 2pt)
    quad2 = (pt * aff3) - 3   #(pt*hold - pt*meas + pt - 3)
    quad3 = hold * aff2       #(hold^2 - hold*pt)
    quad4 = pt * aff5         #(4pt^2 - 42pt)
    # test addition
    @testset "Base.:+" begin
        # test mixed combo types
        @test isa(copy(quad1) + quad2, GenericQuadExpr{Float64, GeneralVariableRef})
        @test isa(copy(quad2) + quad4, GenericQuadExpr{Float64, MeasureFiniteVariableRef})
        @test isa(copy(quad3) + quad4, GenericQuadExpr{Float64, FiniteVariableRef})
        @test isa(copy(quad4) + quad1, GenericQuadExpr{Float64, GeneralVariableRef})
        # test general combos
        @test (copy(quad1) + quad2).aff.constant == -3
        @test (copy(quad1) + quad2).aff.terms[pt] == 3
        @test (copy(quad1) + quad2).aff.terms[hold] == -2
        pair = UnorderedPair{GeneralVariableRef}(inf, hold)
        @test (copy(quad1) + quad2).terms[pair] == 1
        pair = UnorderedPair{GeneralVariableRef}(pt, hold)
        @test (copy(quad1) + quad2).terms[pair] == 2
        @test (copy(quad4) + quad1).aff.constant == 0
        @test (copy(quad4) + quad1).aff.terms[pt] == -40
        @test (copy(quad4) + quad1).aff.terms[hold] == -2
        pair = UnorderedPair{GeneralVariableRef}(pt, pt)
        @test (copy(quad4) + quad1).terms[pair] == 3
        # test measurefinite combos
        @test (copy(quad2) + quad4).aff.constant == -3
        @test (copy(quad2) + quad4).aff.terms[pt] == -41
        pair = UnorderedPair{GeneralVariableRef}(pt, hold)
        @test (copy(quad2) + quad4).terms[pair] == 1
        # test finite combos
        @test (copy(quad3) + quad4).aff.constant == 0
        @test (copy(quad3) + quad4).aff.terms[pt] == -42
        pair = UnorderedPair{GeneralVariableRef}(hold, pt)
        @test (copy(quad3) + quad4).terms[pair] == -1
        # test same types
        @test isa(copy(quad4) + quad4, GenericQuadExpr{Float64, PointVariableRef})
        @test (copy(quad4) + quad4).aff.constant == 0
        @test (copy(quad4) + quad4).aff.terms[pt] == -84
        pair = UnorderedPair{GeneralVariableRef}(pt, pt)
        @test (copy(quad4) + quad4).terms[pair] == 8
    end
    # test subtraction
    @testset "Base.:-" begin
        # test mixed combo types
        @test isa(copy(quad1) - quad2, GenericQuadExpr{Float64, GeneralVariableRef})
        @test isa(copy(quad2) - quad4, GenericQuadExpr{Float64, MeasureFiniteVariableRef})
        @test isa(copy(quad3) - quad4, GenericQuadExpr{Float64, FiniteVariableRef})
        @test isa(copy(quad4) - quad1, GenericQuadExpr{Float64, GeneralVariableRef})
        # test general combos
        @test (copy(quad1) - quad2).aff.constant == 3
        @test (copy(quad1) - quad2).aff.terms[pt] == 1
        @test (copy(quad1) - quad2).aff.terms[hold] == -2
        pair = UnorderedPair{GeneralVariableRef}(inf, hold)
        @test (copy(quad1) - quad2).terms[pair] == 1
        pair = UnorderedPair{GeneralVariableRef}(pt, hold)
        @test (copy(quad1) - quad2).terms[pair] == 0
        @test (copy(quad4) - quad1).aff.constant == 0
        @test (copy(quad4) - quad1).aff.terms[pt] == -44
        @test (copy(quad4) - quad1).aff.terms[hold] == 2
        pair = UnorderedPair{GeneralVariableRef}(pt, pt)
        @test (copy(quad4) - quad1).terms[pair] == 5
        # test measurefinite combos
        @test (copy(quad2) - quad4).aff.constant == -3
        @test (copy(quad2) - quad4).aff.terms[pt] == 43
        pair = UnorderedPair{GeneralVariableRef}(pt, hold)
        @test (copy(quad2) - quad4).terms[pair] == 1
        # test finite combos
        @test (copy(quad3) - quad4).aff.constant == 0
        @test (copy(quad3) - quad4).aff.terms[pt] == 42
        pair = UnorderedPair{GeneralVariableRef}(hold, pt)
        @test (copy(quad3) - quad4).terms[pair] == -1
        # test same types
        @test isa(copy(quad4) - quad4, GenericQuadExpr{Float64, PointVariableRef})
        @test (copy(quad4) - quad4).aff.constant == 0
        @test (copy(quad4) - quad4).aff.terms[pt] == 0
        pair = UnorderedPair{GeneralVariableRef}(pt, pt)
        @test (copy(quad4) - quad4).terms[pair] == 0
    end
    # test multiplication
    @testset "Base.:*" begin
        # test some combos
        @test_throws ErrorException quad1 * quad2
        @test_throws ErrorException quad3 * quad1
        @test_throws ErrorException quad4 * quad4
        @test_throws ErrorException quad2 * quad3
    end
    # test division
    @testset "Base.:/" begin
        # test some combos
        @test_throws ErrorException quad1 / quad2
        @test_throws ErrorException quad3 / quad1
        @test_throws ErrorException quad4 / quad4
        @test_throws ErrorException quad2 / quad3
    end
end


text:
abstract type AbstractFunctional end

mutable struct volumeAverage <: AbstractFunctional
end
function (obj::volumeAverage)(
              mesh::AbstractMesh{Tmsh},
              sbp::AbstractOperator,
              eqn::EllipticData{Tsol, Tres, Tdim},
              opts,
              val::Array{Tsol, 1}) where {Tmsh, Tsol, Tres, Tdim}
  @assert(length(val) == mesh.numDofPerNode)
  val[:] = 0.0

  w = sview(sbp.w, :)
  for elem = 1 : mesh.numEl
    jac = sview(mesh.jac, :, elem)
    q = sview(eqn.q, :, :, elem)
    for n = 1 : mesh.numNodesPerElement
      val[:] += w[n]/jac[n]*q[:, n]
    end
  end

  return nothing
end

mutable struct volumeEnergy <: AbstractFunctional
end
function (obj::volumeEnergy)(
              mesh::AbstractMesh{Tmsh},
              sbp::AbstractOperator,
              eqn::EllipticData{Tsol, Tres, Tdim},
              opts,
              val::Array{Tsol, 1}) where {Tmsh, Tsol, Tres, Tdim}
  @assert(length(val) == mesh.numDofPerNode)
  val[:] = 0.0

  w = sview(sbp.w, :)
  for elem = 1 : mesh.numEl
    jac = sview(mesh.jac, :, elem)
    q = sview(eqn.q, :, :, elem)
    for n = 1 : mesh.numNodesPerElement
      for dof = 1:mesh.numDofPerNode
        val[dof] += w[n]/jac[n]*q[dof, n]*q[dof,n]
      end
    end
  end

  return nothing
end
#

#
# IP Modification of target functional
#
function FunctionalModification(mesh::AbstractMesh{Tmsh},
                                sbp::AbstractOperator,
                                eqn::EllipticData{Tsol, Tres, Tdim},
                                otps::Array{Tsol, 1}) where {Tmsh, Tsol, Tres, Tdim}
  if haskey(opts, "FunctionalModification") && opts["FunctionalModification"] == "No"
    return nothing
  end

  for iBC = 1:mesh.numBC

    bc_func = mesh.bndry_funcs[iBC]
    if isNeumann(bc_func)
      continue
    end

    indx0 = mesh.bndry_offsets[iBC]
    indx1 = mesh.bndry_offsets[iBC+1] - 1

    for f = indx0:indx1

    end

  end
  return nothing
end

global const FunctionalDict = Dict{String, AbstractFunctional}(
  "volumeAverage" => volumeAverage(),
  "energy" => volumeEnergy()
)

function getFunctional(mesh::AbstractMesh,
                       sbp::AbstractOperator,
                       eqn::EllipticData,
                       opts)
  eqn.functional = FunctionalDict[opts["Functional"]]
end


text:
@testset "Paths" begin
  grid = RegularGrid{Float64}(100,100)

  for path in [LinearPath(), RandomPath(),
               ShiftedPath(LinearPath(), 0),
               SourcePath(1:3)]
    p = traverse(grid, path)
    @test length(p) == 100*100
  end

  @testset "LinearPath" begin
    p = traverse(grid, LinearPath())
    @test p == 1:100*100
  end

  @testset "RandomPath" begin
    p = traverse(grid, RandomPath())
    @test all(1 .≤ collect(p) .≤ 100*100)
  end

  @testset "SourcePath" begin
    grid = RegularGrid{Float64}(3,3)
    pset = PointSet(coordinates(grid))

    for sdomain in [grid, pset]
      p = traverse(sdomain, SourcePath([1,9]))
      @test collect(p) == [1,9,2,4,6,8,5,3,7]

      p = traverse(sdomain, SourcePath([1]))
      @test collect(p) == [1,2,4,5,3,7,6,8,9]
    end
  end

  @testset "ShiftedPath" begin
    grid = RegularGrid{Float64}(3,3)
    path = LinearPath()
    for offset in [0,1,-1]
      spath = ShiftedPath(path, offset)
      p  = traverse(grid, path)
      sp = traverse(grid, spath)
      @test length(sp) == 9
      @test collect(sp) == circshift(p, -offset)
    end
  end
end


text:
## -------------------------------------------------------------------
# bounds
let

    ps = Plots.Plot[]
    for ider in FLX_IDERS
        p = plot(title = ider, xlabel = "replica", ylabel = "flx")
        xticks =  (EXPS, string.(EXPS))

        for method in ALL_METHODS
            color = method_colors[method]

            bounds = DAT[method, :bounds, :flx, ider, EXPS]
            lb, ub = first.(bounds), last.(bounds)
            plot!(p, EXPS, [lb ub];
                fillrange = [ub lb],
                fillalpha = 0.2,
                label = "", color = :gray
            )

            ep_vals = DAT[method, :ep, :flx, ider, EXPS]
            ep_errs = DAT[method, :eperr, :flx, ider, EXPS]
            eplb, epub = ep_vals .- ep_errs, ep_vals .+ ep_errs
            plot!(p, EXPS, ep_vals;
                label = string(method), color, alpha = 0.5,
                fillrange = [epub eplb],
                lw = 1, xticks
            )
        end

        Hd_vals = DAT[:Hd, :flx, ider, EXPS]
        plot!(p, EXPS, Hd_vals;
            label = "exp", color = :black, alpha = 0.8, lw = 3, xticks)

        mysavefig(p, "bound_study"; ider)
        push!(ps, p)
    end
    mysavefig(ps, "bound_study")

end


text:
include("loadtrace.jl")

function modelmaker(attack_round, model_type, cipher_type, trace_data)
  num_unit = trace_data.num_unit
  plaintexts = trace_data.plaintexts
  ciphertexts = trace_data.ciphertexts
  unit_size = trace_data.unit_size
  unit_in_pt = div(length(plaintexts[1]),num_unit)
  keymax = 2^unit_size -1

  #the vector to store the model
  V = zeros(Int, num_unit, length(plaintexts), keymax+1)

  if attack_round == "first"
    if model_type == "HW"
      sbox = S_boxes[cipher_type] #S_boxes is a Dict() which returns the S_box of
      #the cipher_type which is to be provided in components
      #for almost all SPN ciphers the first round before Sbox is similar

      for i in 1:num_unit
        for keyGuess in 0:keymax
          for p in 1:length(plaintexts)
            pt_unit = parse(Int,plaintexts[p][(i-1)*unit_in_pt + 1 : i*unit_in_pt],16)
            xor_val = pt_unit$keyGuess
            V[i,p,keyGuess+1] = HW[sbox[xor_val+1]] #HW is Dict needed to be formed depending upon unit_size
          end
        end
      end
    elseif model_type == "HD"
      #the specification of model user needs to provide algo as to how to build model
    end

  elseif attack_round == "last"
    if model_type == "HD"
      #the last round may have different structure for different cipher_type
      #example last round of AES, PRESENT have different structure
      #hence the user may need to provide the algo for the model making of the last round
      #if needed for some cipher type the following is for AES-128
      if cipher_type == "AES"
        invsbox = invS_boxes["AES"]

        for i in 1:num_unit
          if i%4 == 0
              r = 4
          else
              r = i%4
          end
          c = trunc(Int8,ceil(i/4))

          for keyGuess in 0:keymax
            for ct in 1:length(ciphertexts)
              mat_c = transformToMatrix(ciphertexts[ct])
              cell = mat_c[r,c]
              #when I am trying to figure out the subkey corresponding to (r,c)
              #after taking xor with the subkey the value I get(xval) will actually
              #contribute in the leakage with the (r_,c_) in the cipher matrix such that
              #(r_,c_)---shifting-->(r,c) during encryption
              r_,c_ = invShiftRows(r,c)
              cell_ = mat_c[r_,c_]
              xval = cell$keyGuess
              b4_sbox = invsbox[xval+1]
              V[i, ct, keyGuess+1] = HW[b4_sbox$cell_]
            end
          end
        end
      #elseif cipher_type == "<cipher_name>"
        #the code for making model goes here
      end
    #elseif model_type == "<model_type(maybe HD or something)>"
      #the code goes here
    end
  end #the two standard places to perform attacks, can be extended to higher order DPA
V
end

max_stars_repo_name: udion/zypher

max_issues_repo_name: udion/zypher

max_forks_repo_name: udion/zypher


text:
"""
    Diagnostics

Accumulate mean fields and covariance statistics on the computational grid.

"""
module Diagnostics

using Dates
using FileIO
using JLD2
using MPI
using StaticArrays

using ..Atmos
using ..Atmos: thermo_state
using ..DGmethods: num_state, vars_state, num_aux, vars_aux, vars_diffusive, num_diffusive
using ..Mesh.Topologies
using ..Mesh.Grids
using ..MoistThermodynamics
using ..MPIStateArrays
using ..PlanetParameters
using ..VariableTemplates

export gather_diagnostics

include("diagnostic_vars.jl")

function extract_state(dg, localQ, ijk, e)
    bl = dg.balancelaw
    FT = eltype(localQ)
    nstate = num_state(bl, FT)
    l_Q = MArray{Tuple{nstate},FT}(undef)
    for s in 1:nstate
        l_Q[s] = localQ[ijk,s,e]
    end
    return Vars{vars_state(bl, FT)}(l_Q)
end

function extract_aux(dg, auxstate, ijk, e)
    bl = dg.balancelaw
    FT = eltype(auxstate)
    nauxstate = num_aux(bl, FT)
    l_aux = MArray{Tuple{nauxstate},FT}(undef)
    for s in 1:nauxstate
        l_aux[s] = auxstate[ijk,s,e]
    end
    return Vars{vars_aux(bl, FT)}(l_aux)
end

function extract_diffusion(dg, localdiff, ijk, e)
    bl = dg.balancelaw
    FT = eltype(localdiff)
    ndiff = num_diffusive(bl, FT)
    l_diff = MArray{Tuple{ndiff},FT}(undef)
    for s in 1:ndiff
        l_diff[s] = localdiff[ijk,s,e]
    end
    return Vars{vars_diffusive(bl, FT)}(l_diff)
end

# thermodynamic variables of interest
function vars_thermo(FT)
  @vars begin
    q_liq::FT
    q_ice::FT
    q_vap::FT
    T::FT
    θ_liq_ice::FT
    θ_dry::FT
    θ_v::FT
    e_int::FT
    h_m::FT
    h_t::FT
  end
end
num_thermo(FT) = varsize(vars_thermo(FT))
thermo_vars(array) = Vars{vars_thermo(eltype(array))}(array)

function compute_thermo!(FT, bl, state, k, ijk, ev, e, z, zvals, thermoQ, aux)
    e_tot = state.ρe / state.ρ
    ts = thermo_state(bl.moisture, bl.orientation, state, aux)
    e_int = internal_energy(ts)
    Phpart = PhasePartition(ts)

    th = thermo_vars(thermoQ[ijk,e])
    th.q_liq     = Phpart.liq
    th.q_ice     = Phpart.ice
    th.q_vap     = Phpart.tot-Phpart.liq-Phpart.ice
    th.T         = air_temperature(ts)
    th.θ_liq_ice = liquid_ice_pottemp(ts)
    th.θ_dry     = dry_pottemp(ts)
    th.θ_v       = virtual_pottemp(ts)
    th.e_int     = e_int

    # Moist and total henthalpy
    R_m          = gas_constant_air(ts)
    th.h_m       = e_int + R_m * th.T
    th.h_t       = e_tot + R_m * th.T

    return nothing
end

# horizontal averages
function vars_horzavg(FT)
  @vars begin
    ρ::FT
    ρu::FT
    ρv::FT
    ρw::FT
    e_tot::FT
    ρq_tot::FT
    q_liq::FT
    q_vap::FT
    θ_liq_ice::FT
    θ_dry::FT
    θ_v::FT
    e_int::FT
    h_m::FT
    h_t::FT
    qt_sgs::FT
    ht_sgs::FT
  end
end
num_horzavg(FT) = varsize(vars_horzavg(FT))
horzavg_vars(array) = Vars{vars_horzavg(eltype(array))}(array)

function compute_horzsums!(state, diffusive_flx, k, ijk, ev, e,
                           Nqk, nvertelem, MH, localaux, thermoQ,
                           horzsums, repdvsr, LWP)
    th = thermo_vars(thermoQ[ijk,e])
    hs = horzavg_vars(horzsums[k,ev])
    hs.ρ         += MH * state.ρ
    hs.ρu        += MH * state.ρu[1]
    hs.ρv        += MH * state.ρu[2]
    hs.ρw        += MH * state.ρu[3]
    hs.e_tot     += MH * state.ρe
    hs.ρq_tot    += MH * state.moisture.ρq_tot
    hs.q_liq     += MH * th.q_liq
    hs.q_vap     += MH * th.q_vap
    hs.θ_liq_ice += MH * th.θ_liq_ice
    hs.θ_dry     += MH * th.θ_dry
    hs.θ_v       += MH * th.θ_v
    hs.e_int     += MH * th.e_int
    hs.h_m       += MH * th.h_m
    hs.h_t       += MH * th.h_t
    hs.qt_sgs    += MH * diffusive_flx.moisture.ρd_q_tot[end]
    hs.ht_sgs    += MH * diffusive_flx.ρd_h_tot[end]

    repdvsr[Nqk*(ev-1)+k] += MH

    # liquid water path
    # this condition is also going to be used to get the number of points that
    # exist on a horizontal plane provided all planes have the same number of
    # points
    # TODO adjust for possibility of non equivalent horizontal slabs
    if ev == floor(nvertelem/2) && k == floor(Nqk/2)
        # TODO: uncomment the line below after rewriting the LWP assignment below using aux.∫dz...?
        # aux = extract_aux(dg, localaux, ijk, e)
        LWP[1] += MH * (localaux[ijk,1,e] + localaux[ijk,2,e])
    end

    return nothing
end

function compute_diagnosticsums!(state, k, ijk, ev, e, MH, zvals,
                                 thermoQ, horzavgs, dsums)
    th = thermo_vars(thermoQ[ijk,e])
    ha = horzavg_vars(horzavgs[k,ev])
    ds = diagnostic_vars(dsums[k,ev])

    u = state.ρu[1] / state.ρ
    v = state.ρu[2] / state.ρ
    w = state.ρu[3] / state.ρ
    q_tot = state.moisture.ρq_tot / state.ρ
    ũ = ha.ρu / ha.ρ
    ṽ = ha.ρv / ha.ρ
    w̃ = ha.ρw / ha.ρ
    ẽ = ha.e_tot / ha.ρ
    q̃_tot = ha.ρq_tot / ha.ρ

    # vertical coordinate
    ds.z        += MH * zvals[k,ev]

    # state and functions of state
    ds.u        += MH * ũ
    ds.v        += MH * ṽ
    ds.w        += MH * w̃
    ds.e_tot    += MH * ẽ
    ds.q_tot    += MH * ha.ρq_tot / ha.ρ
    ds.q_liq    += MH * ha.q_liq
    ds.thd      += MH * ha.θ_dry
    ds.thl      += MH * ha.θ_liq_ice
    ds.thv      += MH * ha.θ_v
    ds.e_int    += MH * ha.e_int
    ds.h_m      += MH * ha.h_m
    ds.h_t      += MH * ha.h_t
    ds.qt_sgs   += MH * ha.qt_sgs
    ds.ht_sgs   += MH * ha.ht_sgs

    # vertical fluxes
    ds.vert_eddy_mass_flx += MH * (w - w̃) * (state.ρ - ha.ρ)
    ds.vert_eddy_u_flx    += MH * (w - w̃) * (u - ha.ρu / ha.ρ)
    ds.vert_eddy_v_flx    += MH * (w - w̃) * (v - ha.ρv / ha.ρ)
    ds.vert_eddy_qt_flx   += MH * (w - w̃) * (q_tot - q̃_tot)
    ds.vert_qt_flx        += MH * w * q_tot
    ds.vert_eddy_ql_flx   += MH * (w - w̃) * (th.q_liq - ha.q_liq)
    ds.vert_eddy_qv_flx   += MH * (w - w̃) * (th.q_vap - ha.q_vap)
    ds.vert_eddy_thd_flx  += MH * (w - w̃) * (th.θ_dry - ha.θ_dry)
    ds.vert_eddy_thv_flx  += MH * (w - w̃) * (th.θ_v - ha.θ_v)
    ds.vert_eddy_thl_flx  += MH * (w - w̃) * (th.θ_liq_ice - ha.θ_liq_ice)

    # variances
    ds.uvariance += MH * (u - ũ)^2
    ds.vvariance += MH * (v - ṽ)^2
    ds.wvariance += MH * (w - w̃)^2

    # skewness
    ds.wskew     += MH * (w - w̃)^3

    # turbulent kinetic energy
    ds.TKE       += MH * 0.5 * (ds.uvariance + ds.vvariance + ds.wvariance)

    return nothing
end

macro visitQ(nhorzelem, nvertelem, Nqk, Nq, expr)
    return esc(
        quote
            for eh in 1:nhorzelem
                for ev in 1:nvertelem
                    e = ev + (eh - 1) * nvertelem
                    for k in 1:Nqk
                        for j in 1:Nq
                            for i in 1:Nq
                                ijk = i + Nq * ((j-1) + Nq * (k-1))
                                $expr
                            end
                        end
                    end
                end
        end
    end)
end

"""
    gather_diagnostics(mpicomm, dg, Q, current_time_string, out_dir)

Compute various diagnostic variables and write them to JLD2 files in `out_dir`,
indexed by `current_time_string`.
"""
function gather_diagnostics(mpicomm,
                            dg,
                            Q,
                            diagnostics_time_str,
                            sim_time_str,
                            out_dir)
    # make sure this time step is not already recorded
    try
        jldopen(joinpath(out_dir,
                "diagnostics-$(diagnostics_time_str).jld2"), "r") do file
            davgs = file[sim_time_str]
            @warn "diagnostics for time step $(sim_time_str) gathered already"
            return nothing
        end
    catch e
    end

    mpirank = MPI.Comm_rank(mpicomm)
    nranks = MPI.Comm_size(mpicomm)

    # extract grid information
    bl = dg.balancelaw
    grid = dg.grid
    topology = grid.topology
    N = polynomialorder(grid)
    Nq = N + 1
    Nqk = dimensionality(grid) == 2 ? 1 : Nq
    npoints = Nq * Nq * Nqk
    nrealelem = length(topology.realelems)
    nvertelem = topology.stacksize
    nhorzelem = div(nrealelem, nvertelem)

    # get the state, auxiliary and geo variables onto the host if needed
    if Array ∈ typeof(Q).parameters
        localQ    = Q.realdata
        localaux  = dg.auxstate.realdata
        localvgeo = grid.vgeo
        localdiff = dg.diffstate.realdata
    else
        localQ    = Array(Q.realdata)
        localaux  = Array(dg.auxstate.realdata)
        localvgeo = Array(grid.vgeo)
        localdiff = Array(dg.diffstate.realdata)
    end
    FT = eltype(localQ)

    nstate    = num_state(bl, FT)
    nauxstate = num_aux(bl, FT)
    ndiff     = num_diffusive(bl, FT)

    # vertical coordinates and thermo variables
    zvals = zeros(Nqk, nvertelem)
    thermoQ = [zeros(FT, num_thermo(FT)) for _ in 1:npoints, _ in 1:nrealelem]

    # divisor for horizontal averages
    l_repdvsr = zeros(FT, Nqk * nvertelem)

    # horizontal sums and the liquid water path
    l_LWP = zeros(FT, 1)
    horzsums = [zeros(FT, num_horzavg(FT)) for _ in 1:Nqk, _ in 1:nvertelem]

    # compute thermo variables and horizontal sums in a single pass
    @visitQ nhorzelem nvertelem Nqk Nq begin
        state = extract_state(dg, localQ, ijk, e)
        aux   = extract_aux(dg, localaux, ijk, e)

        z = localvgeo[ijk,grid.x3id,e]
        compute_thermo!(FT, bl, state, k, ijk, ev, e, z, zvals, thermoQ, aux)

        diffusive_flx = extract_diffusion(dg, localdiff, ijk, e)
        MH = localvgeo[ijk,grid.MHid,e]
        compute_horzsums!(state, diffusive_flx, k, ijk, ev, e, Nqk, nvertelem,
                          MH, localaux, thermoQ, horzsums, l_repdvsr, l_LWP)
    end

    # compute the full number of points on a slab
    repdvsr = MPI.Allreduce(l_repdvsr, +, mpicomm)

    # compute the horizontal and LWP averages
    horzavgs = [zeros(FT, num_horzavg(FT)) for _ in 1:Nqk, _ in 1:nvertelem]
    for ev in 1:nvertelem
        for k in 1:Nqk
            hsum = MPI.Allreduce(horzsums[k,ev], +, mpicomm)
            horzavgs[k,ev][1:end] = hsum ./ repdvsr[Nqk*(ev-1)+k]
        end
    end
    LWP = zero(FT)
    LWP = MPI.Reduce(l_LWP[1], +, 0, mpicomm)
    if mpirank == 0
        LWP /= repdvsr[1]
    end

    # compute the diagnostics with the previous computed variables
    dsums = [zeros(FT, num_diagnostic(FT)) for _ in 1:Nqk, _ in 1:nvertelem]
    @visitQ nhorzelem nvertelem Nqk Nq begin
        state = extract_state(dg, localQ, ijk, e)
        MH = localvgeo[ijk,grid.MHid,e]
        compute_diagnosticsums!(state, k, ijk, ev, e, MH, zvals,
                                thermoQ, horzavgs, dsums)
    end
    davgs = [zeros(FT, num_diagnostic(FT)) for _ in 1:Nqk, _ in 1:nvertelem]
    for ev in 1:nvertelem
        for k in 1:Nqk
            dsum = MPI.Reduce(dsums[k,ev], +, 0, mpicomm)
            if mpirank == 0
                davgs[k,ev][1:end] = dsum ./ repdvsr[Nqk*(ev-1)+k]
            end
        end
    end

    # write results
    if mpirank == 0
        jldopen(joinpath(out_dir,
                "diagnostics-$(diagnostics_time_str).jld2"), "a+") do file
            file[sim_time_str] = davgs
        end
        jldopen(joinpath(out_dir,
                "liquid_water_path-$(diagnostics_time_str).jld2"), "a+") do file
            file[sim_time_str] = LWP
        end
    end

    return nothing
end # function gather_diagnostics

end # module Diagnostics

max_stars_repo_name: emassoud/CLIMA

max_issues_repo_name: emassoud/CLIMA

max_forks_repo_name: emassoud/CLIMA