text:
# Mean-field variational update rules
@naiveVariationalRule(:node_type     => GeneralisedFilterX,
                      :outbound_type => Message{GaussianMeanPrecision},
                      :inbound_types => (Nothing, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution),
                      :name          => VariationalGFXOutNPPPPP)

@naiveVariationalRule(:node_type     => GeneralisedFilterX,
                      :outbound_type => Message{GaussianWeightedMeanPrecision},
                      :inbound_types => (ProbabilityDistribution, Nothing, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution),
                      :name          => VariationalGFXIn1PNPPPP)

@naiveVariationalRule(:node_type     => GeneralisedFilterX,
                      :outbound_type => Message{GaussianWeightedMeanPrecision},
                      :inbound_types => (ProbabilityDistribution, ProbabilityDistribution, Nothing, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution),
                      :name          => VariationalGFXIn2PPNPPP)

@naiveVariationalRule(:node_type     => GeneralisedFilterX,
                      :outbound_type => Message{GaussianWeightedMeanPrecision},
                      :inbound_types => (ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, Nothing, ProbabilityDistribution, ProbabilityDistribution),
                      :name          => VariationalGFXIn3PPPNPP)

@naiveVariationalRule(:node_type     => GeneralisedFilterX,
                      :outbound_type => Message{GaussianWeightedMeanPrecision},
                      :inbound_types => (ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, Nothing, ProbabilityDistribution),
                      :name          => VariationalGFXIn4PPPPNP)

@naiveVariationalRule(:node_type     => GeneralisedFilterX,
                      :outbound_type => Message{Gamma},
                      :inbound_types => (ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, Nothing),
                      :name          => VariationalGFXIn5PPPPPN)

# Structured updates
#TODO

max_stars_repo_name: wmkouw/gfx-node

max_issues_repo_name: wmkouw/GFX

max_forks_repo_name: wmkouw/GFX


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiNetworkingV1beta1IngressRule <: SwaggerModel
    host::Any # spec type: Union{ Nothing, String } # spec name: host
    http::Any # spec type: Union{ Nothing, IoK8sApiNetworkingV1beta1HTTPIngressRuleValue } # spec name: http

    function IoK8sApiNetworkingV1beta1IngressRule(;host=nothing, http=nothing)
        o = new()
        validate_property(IoK8sApiNetworkingV1beta1IngressRule, Symbol("host"), host)
        setfield!(o, Symbol("host"), host)
        validate_property(IoK8sApiNetworkingV1beta1IngressRule, Symbol("http"), http)
        setfield!(o, Symbol("http"), http)
        o
    end
end # type IoK8sApiNetworkingV1beta1IngressRule

const _property_map_IoK8sApiNetworkingV1beta1IngressRule = Dict{Symbol,Symbol}(Symbol("host")=>Symbol("host"), Symbol("http")=>Symbol("http"))
const _property_types_IoK8sApiNetworkingV1beta1IngressRule = Dict{Symbol,String}(Symbol("host")=>"String", Symbol("http")=>"IoK8sApiNetworkingV1beta1HTTPIngressRuleValue")
Base.propertynames(::Type{ IoK8sApiNetworkingV1beta1IngressRule }) = collect(keys(_property_map_IoK8sApiNetworkingV1beta1IngressRule))
Swagger.property_type(::Type{ IoK8sApiNetworkingV1beta1IngressRule }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiNetworkingV1beta1IngressRule[name]))}
Swagger.field_name(::Type{ IoK8sApiNetworkingV1beta1IngressRule }, property_name::Symbol) =  _property_map_IoK8sApiNetworkingV1beta1IngressRule[property_name]

function check_required(o::IoK8sApiNetworkingV1beta1IngressRule)
    true
end

function validate_property(::Type{ IoK8sApiNetworkingV1beta1IngressRule }, name::Symbol, val)
end

ext: jl


text:
@testset "logseries dist" begin
  @test logseriescdf(0.01)[1:3] ≈ [0.0, 0.994992, 0.999967] atol=1.0e-5
  @test logseriescdf(0.000000000001)[1:5] ≈ [0., 1., 1., 1., 1.] atol=1.0e-2
  #@test logseriesquantile(0.9, [0.25, 0.5, 0.75]) == [1, 2, 5]
  Random.seed!(43)
  w = logseriescdf(0.4)
  v = pmap(i -> findlast(w .< i), rand(1000000))
  @test mean(v) ≈ 1.304 atol=1.0e-2
  @test std(v) ≈ 0.687 atol=1.0e-2
  @test skewness(v) ≈ 3.1 atol=1.0e-2
  @test kurtosis(v) ≈ 13.5 atol=1.0
end
@testset "stable levy dist" begin
  Random.seed!(43)
  @test levyel(2., rand(), rand()) ≈ 0.18170339379413047
    Random.seed!(43)
  @test tiltedlevygen(0.2, 2.; rng = Random.GLOBAL_RNG) ≈ 0.00409133 atol=1.0e-5
  Random.seed!(43)
  @test tiltedlevygen(0.6, 2.; rng = Random.GLOBAL_RNG) ≈ 0.036822009 atol=1.0e-5
end
@testset "nested copulas data generators" begin
  @test Ginv(0.5, 0.5) ≈ 1.2732395447351625
  @test InvlaJ(4, 0.5) ≈ 0.7265625
  @test sampleInvlaJ(0.5, 0.5) == 1
  @test sampleInvlaJ(0.5, 0.8) == 8
  Random.seed!(43)
  w = logseriescdf(1-exp(-4.))
  @test elInvlaF(4., 2., w; rng = Random.GLOBAL_RNG) == 7
  Random.seed!(43)
  @test elInvlaF(4., .5, w; rng = Random.GLOBAL_RNG) == 16
  @test nestedfrankgen(4., 3., 1, w; rng = Random.GLOBAL_RNG) == 6
end


text:
# Simple tests to make sure Roots works with other
# packages. This is not part of runtests, as no expectation
# that CI should run these

using Roots
using Compat.Test
using Unitful
using SymEngine

@testset "Test composability with other packages" begin

    orders = [Order0(), Order1(), Order2(), Order5(),
              Order8(), Order16()]

    # Unitful
    @testset "find zero(s) with Unitful" begin

        s = u"s"; m = u"m"
        g = 9.8*m/s^2
        v0 = 10m/s
        y0 = 16m
        y(t) = -g*t^2 + v0*t + y0

        for order in orders
            @test find_zero(y, 1.8s, order) ≈  1.886053370668014s
        end

        for M in [Roots.Bisection(), Roots.A42(), Roots.AlefeldPotraShi()]
            @test find_zero(y, (1.8s, 1.9s), M) ≈ 1.886053370668014s
        end

        xrts = find_zeros(y, 0s, 10s)
        @test length(xrts) == 1
        @test xrts[1] ≈  1.886053370668014s

    end

    # SymEngine
    @testset "find zero(s) with SymEngine" begin

        m = s = 1.0
        g = 9.8*m/s^2
        v0 = 10m/s
        y0 = 16m
        y(t) = -g*t^2 + v0*t + y0

        @vars x
        for order in orders
            @test find_zero(y(x), 1.8, order) ≈  1.8860533706680143
        end

        for M in [Roots.Bisection(), Roots.A42(), Roots.AlefeldPotraShi()]
            @test find_zero(y(x), (1.8, 1.9), M) ≈ 1.8860533706680143
        end

        xrts = find_zeros(y(x), 0, 10)
        @test length(xrts) == 1
        @test xrts[1] ≈  1.8860533706680143
    end

end


text:
module TSSM

using Distributions
using HypothesisTests
import HypothesisTests.ksstats
using LinearAlgebra
using Random
using Statistics
using StatsBase

include("weight_matrix.jl")
include("summary_stats.jl")
include("distance_measure.jl")

export BlockBootstrapWeightMatrix, select_moments, weighted_moment_distance

end


text:
## A small driver file to run the bad LOO example and combine
#run like this
#julia -p 3 -L testHarness_Paper.jl tests_LOO.jl
#pass arguments for things via ARGS[1] is numRun per batch

numRuns = parse(Int, ARGS[1])

const param_path = "./Results/badLOOExample.csv"
n_grid = [2^i for i = 5:14]

tic()
a = @spawn test_ReadData(spath, numRuns, n_grid, 8675309000, param_path, Gamma_min=1, Gamma_max=10.)
b = @spawn test_ReadData(spath, numRuns, n_grid, 5164174290, param_path, Gamma_min=1, Gamma_max=10.)
c = @spawn test_ReadData(spath, numRuns, n_grid, 123456, param_path, Gamma_min=1, Gamma_max=10.)
d = @spawn test_ReadData(spath, numRuns, n_grid, 5167462266, param_path, Gamma_min=1, Gamma_max=10.)

######
file_a = fetch(a)
file_b = fetch(b)
file_c = fetch(c)
file_d = fetch(d)

time_stamp = toc()

##read everyone in, throw away a line
data, header = readdlm(file_a, ',', header=true)

data_t = readdlm(file_b, ',', skipstart=1)
data_t[:, 1] += numRuns
data = vcat(data, data_t)

data_t = readdlm(file_c, ',', skipstart=1)
data_t[:, 1] += 2numRuns
data = vcat(data, data_t)

data_t = readdlm(file_d, ',', skipstart=1)
data_t[:, 1] += 3numRuns
data = vcat(data, data_t)

#strip the name of file_a to make the numbers better
println("Filea \t", file_a)
println("spath \t", spath)

indx = search(file_a, spath)[end] + 1

f = open("$(spath)_$(file_a[indx:end])_full_$(4numRuns).csv", "w")
writedlm(f,  header)
writedlm(f,  data)
close(f)

println("Num Paths: \t $(numRuns) \t Time:", time_stamp )

max_stars_repo_path: Normal Model/testLOO.jl

max_issues_repo_path: Normal Model/testLOO.jl

max_forks_repo_path: Normal Model/testLOO.jl


text:
"""
    Layer(transformations::Tuple, data, positional::Tuple, named::NamedTuple)

Algebraic object encoding a single layer of a visualization. It is composed of a dataset,
positional and named arguments, as well as transformations to be applied to those.
`Layer` objects can be multiplied, yielding a novel `Layer` object, or added,
yielding a [`AlgebraOfGraphics.Layers`](@ref) object.
"""
struct Layer
    transformations::Tuple
    data::Any
    positional::Tuple
    named::NamedTuple
end

Layer(transformations::Tuple=()) = Layer(transformations, nothing, (), (;))

data(df) = Layer((), columns(df), (), (;))
mapping(args...; kwargs...) = Layer((), nothing, args, values(kwargs))

function Base.:*(l1::Layer, l2::Layer)
    t1, t2 = l1.transformations, l2.transformations
    d1, d2 = l1.data, l2.data
    p1, p2 = l1.positional, l2.positional
    n1, n2 = l1.named, l2.named
    transformations = (t1..., t2...)
    data = isnothing(d2) ? d1 : d2
    positional =(p1..., p2...)
    named = merge(n1, n2)
    return Layer(transformations, data, positional, named)
end


text:
module fock

import Base.==

using ..bases, ..states, ..operators, ..operators_sparse

export FockBasis, number, destroy, create, fockstate, coherentstate, qfunc

"""
Basis for a Fock space.

Arguments
---------

Nmin
    Minimal particle number included in this basis. By default it is zero, i.e.
    starting with the vacuum state.
Nmax
    Maximal particle number included in this basis.
"""
type FockBasis <: Basis
    shape::Vector{Int}
    Nmin::Int
    Nmax::Int
    function FockBasis(Nmin::Int, Nmax::Int)
        if Nmin < 0 || Nmax <= Nmin
            throw(DimensionMismatch())
        end
        new([Nmax-Nmin+1], Nmin, Nmax)
    end
end

FockBasis(Nmax::Int) = FockBasis(0, Nmax)

==(b1::FockBasis, b2::FockBasis) = b1.Nmin==b2.Nmin && b1.Nmax==b2.Nmax

"""
Number operator for the given Fock space.

Arguments
---------

b
    FockBasis of this operator.
"""
function number(b::FockBasis)
    N = b.Nmax-b.Nmin+1
    diag = Complex128[complex(x) for x=b.Nmin:b.Nmax]
    data = spdiagm(diag, 0, N, N)
    SparseOperator(b, data)
end

"""
Annihilation operator for the given Fock space.

Arguments
---------

b
    FockBasis of this operator.
"""
function destroy(b::FockBasis)
    N = b.Nmax-b.Nmin+1
    diag = Complex128[complex(sqrt(x)) for x=b.Nmin+1:b.Nmax]
    data = spdiagm(diag, 1, N, N)
    SparseOperator(b, data)
end

"""
Creation operator for the given Fock space.

Arguments
---------

b
    FockBasis of this operator.
"""
function create(b::FockBasis)
    N = b.Nmax-b.Nmin+1
    diag = Complex128[complex(sqrt(x)) for x=b.Nmin+1:b.Nmax]
    data = spdiagm(diag, -1, N, N)
    SparseOperator(b, data)
end

"""
Fock state for the given particle number.

Arguments
---------

b
    FockBasis for this state.
n
    Quantum number of the state.
"""
function fockstate(b::FockBasis, n::Int)
    @assert b.Nmin <= n <= b.Nmax
    basis_ket(b, n+1-b.Nmin)
end

"""
Coherent state :math:`| \\alpha \\rangle` for the given Fock basis.

Arguments
---------

b
    FockBasis for this state.
alpha
    Eigenvalue of annihilation operator.
"""
function coherentstate(b::FockBasis, alpha::Number)
    alpha = complex(alpha)
    x = zeros(Complex128, b.Nmax - b.Nmin + 1)
    if b.Nmin == 0
        x[1] = exp(-abs2(alpha)/2)
    else
        x[1] = exp(-abs2(alpha)/2)*alpha^b.Nmin/sqrt(factorial(b.Nmin))
    end
    for n=1:(b.Nmax-b.Nmin)
        x[n+1] = x[n]*alpha/sqrt(b.Nmin + n)
    end
    return Ket(b, x)
end

"""
Husimi Q representation :math:`\\frac{1}{\\pi} \\langle \\alpha | \\rho | \\alpha \\rangle`.
"""
function qfunc(rho::Operator, alpha::Complex128)
    psi = coherentstate(rho.basis_l, alpha)
    return real(dagger(psi)*rho*psi)/pi
end

function qfunc(rho::Operator, X::Vector{Float64}, Y::Vector{Float64})
    @assert rho.basis_l == rho.basis_r
    return Float64[qfunc(rho, complex(x,y)) for x=X, y=Y]
end

end # module


text:
# Test variable reference makers
@testset "Variable Reference Makers/Deleters" begin
    # initialize model and references
    m = InfiniteModel()
    @infinite_parameter(m, par1 in [0, 2])
    @infinite_parameter(m, par2 in [0, 2])
    @infinite_variable(m, inf1(par1))
    @infinite_variable(m, inf2(par1, par2))
    d1 = @deriv(inf1, par1)
    d2 = @deriv(inf2, par1)
    # test make_point_variable_ref (InfiniteModel)
    @testset "make_point_variable_ref (InfiniteModel)" begin
        # test with inf1
        pvref = GeneralVariableRef(m, 1, PointVariableIndex)
        @test make_point_variable_ref(m, inf1, Float64[0]) == pvref
        @test parameter_values(pvref) == (0.,)
        # test with inf2
        pvref = GeneralVariableRef(m, 2, PointVariableIndex)
        @test make_point_variable_ref(m, inf2, Float64[0, 0]) == pvref
        @test parameter_values(pvref) == (0., 0.)
        # test with d1
        pvref = GeneralVariableRef(m, 3, PointVariableIndex)
        @test make_point_variable_ref(m, d1, Float64[0]) == pvref
        @test parameter_values(pvref) == (0.,)
        # test with d2
        pvref = GeneralVariableRef(m, 4, PointVariableIndex)
        @test make_point_variable_ref(m, d2, Float64[0, 0]) == pvref
        @test parameter_values(pvref) == (0., 0.)
    end
    # test make_point_variable_ref with infinite parameter functions
    @testset "make_point_variable_ref (Parameter Function)" begin
        f = parameter_function(sin, par1)
        @test make_point_variable_ref(m, f, [0.]) == 0
        @test make_point_variable_ref(Model(), f, [0.]) == 0
    end
    # test add_measure_variable
    @testset "add_measure_variable" begin
        @test_throws ErrorException add_measure_variable(Model(), Bad(), Val(:some_key))
    end
    # test make_point_variable_ref (optmizer_model)
    @testset "make_point_variable_ref (optimizer_model)" begin
        opt_m = Model()
        opt_m.ext[:my_key] = 42
        @test_throws ErrorException make_point_variable_ref(opt_m, inf1, Float64[0])
    end
    # test make_reduced_variable_ref (InfiniteModel)
    @testset "make_reduced_variable_ref (InfiniteModel)" begin
        # test first addition
        rvref1 = GeneralVariableRef(m, 1, ReducedVariableIndex)
        @test make_reduced_variable_ref(m, inf2, [1], Float64[1]) == rvref1
        @test infinite_variable_ref(rvref1) == inf2
        @test eval_supports(rvref1) == Dict(1 => Float64(1))
        # test second addition
        rvref2 = GeneralVariableRef(m, 2, ReducedVariableIndex)
        @test make_reduced_variable_ref(m, inf2, [1], Float64[0]) == rvref2
        @test infinite_variable_ref(rvref2) == inf2
        @test eval_supports(rvref2) == Dict(1 => Float64(0))
        # test first addition with derivative
        rvref1 = GeneralVariableRef(m, 3, ReducedVariableIndex)
        @test make_reduced_variable_ref(m, d2, [1], Float64[1]) == rvref1
        @test infinite_variable_ref(rvref1) == d2
        @test eval_supports(rvref1) == Dict(1 => Float64(1))
        # test second addition with derivative
        rvref2 = GeneralVariableRef(m, 4, ReducedVariableIndex)
        @test make_reduced_variable_ref(m, d2, [1], Float64[0]) == rvref2
        @test infinite_variable_ref(rvref2) == d2
        @test eval_supports(rvref2) == Dict(1 => Float64(0))
    end
    # test make_reduced_variable_ref (optimizer_model)
    @testset "make_reduced_variable_ref (optimizer_model)" begin
        opt_m = Model()
        opt_m.ext[:my_key] = 42
        @test_throws ErrorException make_reduced_variable_ref(opt_m, inf2, [1], Float64[1])
    end
    # test delete_internal_reduced_variable (InfiniteModel)
    @testset "delete_internal_reduced_variable (InfiniteModel)" begin
        rvref1 = ReducedVariableRef(m, ReducedVariableIndex(1))
        rvref2 = ReducedVariableRef(m, ReducedVariableIndex(2))
        rvref3 = ReducedVariableRef(m, ReducedVariableIndex(3))
        rvref4 = ReducedVariableRef(m, ReducedVariableIndex(4))
        # test cannot delete
        push!(InfiniteOpt._constraint_dependencies(rvref1), ConstraintIndex(1))
        @test delete_internal_reduced_variable(m, rvref1) isa Nothing
        @test is_valid(m, rvref1)
        empty!(InfiniteOpt._constraint_dependencies(rvref1))
        # test can delete
        @test delete_internal_reduced_variable(m, rvref1) isa Nothing
        @test delete_internal_reduced_variable(m, rvref2) isa Nothing
        @test delete_internal_reduced_variable(m, rvref3) isa Nothing
        @test delete_internal_reduced_variable(m, rvref4) isa Nothing
        @test !is_valid(m, rvref1)
        @test !is_valid(m, rvref2)
        @test !is_valid(m, rvref3)
        @test !is_valid(m, rvref4)
    end
    # test delete_reduced_variable
    @testset "delete_reduced_variable" begin
        rvref = ReducedVariableRef(m, ReducedVariableIndex(1))
        warn = "'delete_reduced_variable' not extended for reduced variable type " *
              "`ReducedVariableRef` and optimizer model with key `bad`."
        @test_logs (:warn, warn) delete_reduced_variable(Model(), rvref, Val(:bad))
    end
    # test delete_internal_reduced_variable (optimizer_model)
    @testset "delete_internal_reduced_variable (optimizer_model)" begin
        rvref = ReducedVariableRef(m, ReducedVariableIndex(1))
        @test !is_valid(m, rvref)
        opt_m = Model()
        opt_m.ext[:my_key] = 42
        warn = "'delete_reduced_variable' not extended for reduced variable type " *
              "`ReducedVariableRef` and optimizer model with key `my_key`."
        @test_logs (:warn, warn) delete_internal_reduced_variable(opt_m, rvref)
    end
end

# Test measure expansion methods
@testset "expand_measure" begin
    # initialize model and references
    m = InfiniteModel()
    @infinite_parameter(m, 1 <= par1 <= 2)
    @infinite_parameter(m, 1 <= par2 <= 2)
    @infinite_parameter(m, 1 <= pars1[1:2] <= 2, independent = true)
    @infinite_parameter(m, 1 <= pars2[1:2] <= 2)
    @infinite_variable(m, inf1(par1))
    @infinite_variable(m, inf2(par1, par2))
    @infinite_variable(m, inf3(par2))
    @infinite_variable(m, inf4(pars1))
    @infinite_variable(m, inf5(pars1, pars2))
    @infinite_variable(m, inf6(pars2))
    @infinite_variable(m, inf7(par1, par2, pars1))
    @infinite_variable(m, inf8(pars1[1]))
    @infinite_variable(m, inf9(pars1[1], par1))
    @hold_variable(m, x)
    d1 = @deriv(inf1, par1)
    d2 = @deriv(inf2, par1)
    d3 = @deriv(inf3, par2)
    d4 = @deriv(inf4, pars1[1])
    d5 = @deriv(inf5, pars1[1])
    d6 = @deriv(inf6, pars2[1])
    d7 = @deriv(inf7, par1)
    d8 = @deriv(inf8, pars1[1])
    d9 = @deriv(inf9, pars1[1])
    f1 = parameter_function(sin, par1)
    f2 = parameter_function((a,b) -> 1, (par1, par2))
    f3 = parameter_function(cos, par2)
    f4 = parameter_function(a -> -1, (pars1,))
    f5 = parameter_function((a,b) -> 1, (pars1, pars2))
    f6 = parameter_function(a -> -1, (pars2,))
    f7 = parameter_function((a,b,c) -> 1, (par1, par2, pars1))
    f8 = parameter_function(cos, pars1[1])
    f9 = parameter_function((a,b) -> 1, (pars1[1], par1))
    # prepare measures
    w(t) = 3
    data1 = DiscreteMeasureData(par1, [0.5, 0.5], [1, 2])
    data2 = DiscreteMeasureData(par2, [0.5, 0.5], [1, 1], weight_function = w)
    data3 = DiscreteMeasureData(pars1, [1, 1], [[1, 1], [2, 2]])
    data4 = DiscreteMeasureData(pars2, [2, 2], [[1, 1], [1, 1]])
    data5 = DiscreteMeasureData([pars1[2], pars1[1]], [1, 1], [[1, 1.2], [2, 2]])
    data6 = DiscreteMeasureData([pars1[2]], [1, 1], [[1], [2]])
    data7 = DiscreteMeasureData(pars1[1], [0.5, 0.5], [1, 2])
    meas1 = measure(inf1, data1)
    meas2 = measure(2 * inf3 * x - par1, data2)
    # test expand_measure (infinite variable) with DiscreteMeasureData
    @testset "Infinite Variable (1D DiscreteMeasureData)" begin
        # test single param infinite var with measure param
        idx = length(InfiniteOpt._data_dictionary(m, PointVariable)) + 1
        pts = [GeneralVariableRef(m, idx, PointVariableIndex),
               GeneralVariableRef(m, idx + 1, PointVariableIndex)]
        @test InfiniteOpt.expand_measure(inf1, data1, m) == 0.5 * (pts[1] + pts[2])
        @test parameter_values.(pts) == [(1.,), (2.,)]
        # test single param infinite var without measure param
        @test InfiniteOpt.expand_measure(inf1, data2, m) == 3inf1
        # test single param infinite var with measure param and others
        idx = length(InfiniteOpt._data_dictionary(m, ReducedVariable)) + 1
        rv1 = GeneralVariableRef(m, idx, ReducedVariableIndex)
        rv2 = GeneralVariableRef(m, idx + 1, ReducedVariableIndex)
        @test InfiniteOpt.expand_measure(inf2, data1, m) == 0.5 * (rv1 + rv2)
        @test eval_supports(rv2) == Dict{Int, Float64}(1 => 2.)
        # test multi param infinite var with single element evaluation
        rv1 = GeneralVariableRef(m, idx + 2, ReducedVariableIndex)
        rv2 = GeneralVariableRef(m, idx + 3, ReducedVariableIndex)
        @test InfiniteOpt.expand_measure(inf7, data7, m) == 0.5 * (rv1 + rv2)
        @test eval_supports(rv2) == Dict{Int, Float64}(3 => 2.)
    end
    # test expand_measure (derivatives) with DiscreteMeasureData
    @testset "Derivative (1D DiscreteMeasureData)" begin
        # test single param infinite var with measure param
        idx = length(InfiniteOpt._data_dictionary(m, PointVariable)) + 1
        pts = [GeneralVariableRef(m, idx, PointVariableIndex),
               GeneralVariableRef(m, idx + 1, PointVariableIndex)]
        @test InfiniteOpt.expand_measure(d1, data1, m) == 0.5 * (pts[1] + pts[2])
        @test parameter_values.(pts) == [(1.,), (2.,)]
        # test single param infinite var without measure param
        @test InfiniteOpt.expand_measure(d1, data2, m) == 3d1
        # test single param infinite var with measure param and others
        idx = length(InfiniteOpt._data_dictionary(m, ReducedVariable)) + 1
        rv1 = GeneralVariableRef(m, idx, ReducedVariableIndex)
        rv2 = GeneralVariableRef(m, idx + 1, ReducedVariableIndex)
        @test InfiniteOpt.expand_measure(d2, data1, m) == 0.5 * (rv1 + rv2)
        @test eval_supports(rv2) == Dict{Int, Float64}(1 => 2.)
        # test multi param infinite var with single element evaluation
        rv1 = GeneralVariableRef(m, idx + 2, ReducedVariableIndex)
        rv2 = GeneralVariableRef(m, idx + 3, ReducedVariableIndex)
        @test InfiniteOpt.expand_measure(d7, data7, m) == 0.5 * (rv1 + rv2)
        @test eval_supports(rv2) == Dict{Int, Float64}(3 => 2.)
    end
    # test expand_measure (infinite parameter functions) with DiscreteMeasureData
    @testset "Parameter Function (1D DiscreteMeasureData)" begin
        # test single param infinite var with measure param
        @test InfiniteOpt.expand_measure(f1, data1, m) == 0.5 * (sin(1) + sin(2))
        # test single param infinite var without measure param
        @test InfiniteOpt.expand_measure(f1, data2, m) == 3 * f1
        # test single param infinite var with measure param and others
        idx = length(InfiniteOpt._data_dictionary(m, ReducedVariable)) + 1
        rv1 = GeneralVariableRef(m, idx, ReducedVariableIndex)
        rv2 = GeneralVariableRef(m, idx + 1, ReducedVariableIndex)
        @test InfiniteOpt.expand_measure(f2, data1, m) == 0.5 * (rv1 + rv2)
        @test eval_supports(rv2) == Dict{Int, Float64}(1 => 2.)
        # test multi param infinite var with single element evaluation
        rv1 = GeneralVariableRef(m, idx + 2, ReducedVariableIndex)
        rv2 = GeneralVariableRef(m, idx + 3, ReducedVariableIndex)
        @test InfiniteOpt.expand_measure(f7, data7, m) == 0.5 * (rv1 + rv2)
        @test eval_supports(rv2) == Dict{Int, Float64}(3 => 2.)
    end
    # test expand_measure (infinite variable) with Multi DiscreteMeasureData
    @testset "Infinite Variable (Multi DiscreteMeasureData)" begin
        # test infinite var with measure param
        idx = length(InfiniteOpt._data_dictionary(m, PointVariable)) + 1
        pts = [GeneralVariableRef(m, idx, PointVariableIndex),
               GeneralVariableRef(m, idx + 1, PointVariableIndex)]
        @test InfiniteOpt.expand_measure(inf4, data3, m) == pts[1] + pts[2]
        @test parameter_values(pts[1]) == (Float64[1., 1.],)
        # test infinite var without measure param
        @test InfiniteOpt.expand_measure(inf4, data4, m) == 4inf4
        # test infinite var with measure param and others
        idx = length(InfiniteOpt._data_dictionary(m, ReducedVariable)) + 1
        rv1 = GeneralVariableRef(m, idx, ReducedVariableIndex)
        rv2 = GeneralVariableRef(m, idx + 1, ReducedVariableIndex)
        @test InfiniteOpt.expand_measure(inf5, data3, m) == rv1 + rv2
        @test eval_supports(rv2) == Dict{Int, Float64}(1 => 2, 2 => 2)
        # test infinite var with measure param that is out of order
        idx = length(InfiniteOpt._data_dictionary(m, PointVariable)) + 1
        pts = [GeneralVariableRef(m, idx, PointVariableIndex),
               GeneralVariableRef(m, idx + 1, PointVariableIndex)]
        @test InfiniteOpt.expand_measure(inf4, data5, m) == pts[1] + pts[2]
        @test parameter_values(pts[1]) == (Float64[1.2, 1],)
        # test infinite var with measure param that doesn't overlap
        @test InfiniteOpt.expand_measure(inf8, data6, m) == 2inf8
        # test infinite var with measure param that doesn't overlap and there are others
        @test InfiniteOpt.expand_measure(inf9, data6, m) == 2inf9
        # test infinite variable has subset of multi params
        idx = length(InfiniteOpt._data_dictionary(m, PointVariable)) + 1
        pts = [GeneralVariableRef(m, idx, PointVariableIndex),
               GeneralVariableRef(m, idx + 1, PointVariableIndex)]
        @test InfiniteOpt.expand_measure(inf8, data3, m) == pts[1] + pts[2]
        @test parameter_values(pts[2]) == (2.,)
        # test making reduced vars with variable not containing all the measure prefs
        idx = length(InfiniteOpt._data_dictionary(m, ReducedVariable)) + 1
        rv1 = GeneralVariableRef(m, idx, ReducedVariableIndex)
        rv2 = GeneralVariableRef(m, idx + 1, ReducedVariableIndex)
        @test InfiniteOpt.expand_measure(inf9, data3, m) == rv1 + rv2
        @test eval_supports(rv2) == Dict{Int, Float64}(1 => 2)
    end
    # test expand_measure (derivatives) with Multi DiscreteMeasureData
    @testset "Derivative (Multi DiscreteMeasureData)" begin
        # test infinite var with measure param
        idx = length(InfiniteOpt._data_dictionary(m, PointVariable)) + 1
        pts = [GeneralVariableRef(m, idx, PointVariableIndex),
               GeneralVariableRef(m, idx + 1, PointVariableIndex)]
        @test InfiniteOpt.expand_measure(d4, data3, m) == pts[1] + pts[2]
        @test parameter_values(pts[1]) == (Float64[1., 1.],)
        # test infinite var without measure param
        @test InfiniteOpt.expand_measure(d4, data4, m) == 4d4
        # test infinite var with measure param and others
        idx = length(InfiniteOpt._data_dictionary(m, ReducedVariable)) + 1
        rv1 = GeneralVariableRef(m, idx, ReducedVariableIndex)
        rv2 = GeneralVariableRef(m, idx + 1, ReducedVariableIndex)
        @test InfiniteOpt.expand_measure(d5, data3, m) == rv1 + rv2
        @test eval_supports(rv2) == Dict{Int, Float64}(1 => 2, 2 => 2)
        # test infinite var with measure param that is out of order
        idx = length(InfiniteOpt._data_dictionary(m, PointVariable)) + 1
        pts = [GeneralVariableRef(m, idx, PointVariableIndex),
               GeneralVariableRef(m, idx + 1, PointVariableIndex)]
        @test InfiniteOpt.expand_measure(d4, data5, m) == pts[1] + pts[2]
        @test parameter_values(pts[1]) == (Float64[1.2, 1],)
        # test infinite var with measure param that doesn't overlap
        @test InfiniteOpt.expand_measure(d8, data6, m) == 2d8
        # test infinite var with measure param that doesn't overlap and there are others
        @test InfiniteOpt.expand_measure(d9, data6, m) == 2d9
        # test infinite variable has subset of multi params
        idx = length(InfiniteOpt._data_dictionary(m, PointVariable)) + 1
        pts = [GeneralVariableRef(m, idx, PointVariableIndex),
               GeneralVariableRef(m, idx + 1, PointVariableIndex)]
        @test InfiniteOpt.expand_measure(d8, data3, m) == pts[1] + pts[2]
        @test parameter_values(pts[2]) == (2.,)
        # test making reduced vars with variable not containing all the measure prefs
        idx = length(InfiniteOpt._data_dictionary(m, ReducedVariable)) + 1
        rv1 = GeneralVariableRef(m, idx, ReducedVariableIndex)
        rv2 = GeneralVariableRef(m, idx + 1, ReducedVariableIndex)
        @test InfiniteOpt.expand_measure(d9, data3, m) == rv1 + rv2
        @test eval_supports(rv2) == Dict{Int, Float64}(1 => 2)
    end
    # test expand_measure (infinite parameter functions) with Multi DiscreteMeasureData
    @testset "Parameter Function (Multi DiscreteMeasureData)" begin
        # test infinite var with measure param
        @test InfiniteOpt.expand_measure(f4, data3, m) == -2
        # test infinite var without measure param
        @test InfiniteOpt.expand_measure(f4, data4, m) == 4 * f4
        # test infinite var with measure param and others
        idx = length(InfiniteOpt._data_dictionary(m, ReducedVariable)) + 1
        rv1 = GeneralVariableRef(m, idx, ReducedVariableIndex)
        rv2 = GeneralVariableRef(m, idx + 1, ReducedVariableIndex)
        @test InfiniteOpt.expand_measure(f5, data3, m) == rv1 + rv2
        @test eval_supports(rv2) == Dict{Int, Float64}(1 => 2, 2 => 2)
        # test infinite var with measure param that is out of order
        @test InfiniteOpt.expand_measure(f4, data5, m) == -2
        # test infinite var with measure param that doesn't overlap
        @test InfiniteOpt.expand_measure(f8, data6, m) == 2 * f8
        # test infinite var with measure param that doesn't overlap and there are others
        @test InfiniteOpt.expand_measure(f9, data6, m) == 2 * f9
        # test infinite variable has subset of multi params
        @test InfiniteOpt.expand_measure(f8, data3, m) == cos(1) + cos(2)
        # test making reduced vars with variable not containing all the measure prefs
        idx = length(InfiniteOpt._data_dictionary(m, ReducedVariable)) + 1
        rv1 = GeneralVariableRef(m, idx, ReducedVariableIndex)
        rv2 = GeneralVariableRef(m, idx + 1, ReducedVariableIndex)
        @test InfiniteOpt.expand_measure(f9, data3, m) == rv1 + rv2
        @test eval_supports(rv2) == Dict{Int, Float64}(1 => 2)
    end
    @testset "_make_point_support (1D DiscreteMeasureData)" begin
        orig_prefs = parameter_list(inf2)
        support_dict = Dict{Int, Float64}(2 => 42)
        expected = Float64[23, 42]
        @test InfiniteOpt._make_point_support(orig_prefs, support_dict, 1, Float64(23)) == expected
    end
    # test expand_measure (reduced infinite variable)
    @testset "Reduced Variable (1D DiscreteMeasureData)" begin
        # test single param reduced var without measure param
        rv = make_reduced_variable_ref(m, inf2, [1], Float64[1])
        @test InfiniteOpt.expand_measure(rv, data1, m) == 0.5 * (rv + rv)
        # test single param reduced var with measure param
        idx = length(InfiniteOpt._data_dictionary(m, PointVariable)) + 1
        pts = [GeneralVariableRef(m, idx, PointVariableIndex),
               GeneralVariableRef(m, idx + 1, PointVariableIndex)]
        @test InfiniteOpt.expand_measure(rv, data2, m) == 1.5 * (pts[1] + pts[2])
        @test parameter_values(pts[1]) == (1., 1.)
        @test !is_valid(m, rv)
        # test single param reduced var with measure param and others
        rv = make_reduced_variable_ref(m, inf7, [1], Float64[1])
        idx = index(rv).value + 1
        rv1 = GeneralVariableRef(m, idx, ReducedVariableIndex)
        rv2 = GeneralVariableRef(m, idx + 1, ReducedVariableIndex)
        @test InfiniteOpt.expand_measure(rv, data2, m) == 1.5 * (rv1 + rv2)
        @test eval_supports(rv2) == Dict{Int, Float64}(1 => 1, 2 => 1)
        @test !is_valid(m, rv)
        # test single param reduced var partially from array element
        rv = make_reduced_variable_ref(m, inf7, [1], Float64[1])
        idx = index(rv).value + 1
        rv1 = GeneralVariableRef(m, idx, ReducedVariableIndex)
        rv2 = GeneralVariableRef(m, idx + 1, ReducedVariableIndex)
        @test InfiniteOpt.expand_measure(rv, data7, m) == 0.5 * (rv1 + rv2)
        @test eval_supports(rv2) == Dict{Int, Float64}(1 => 1, 3 => 2)
        @test !is_valid(m, rv)
    end
    @testset "_make_point_support (Multi DiscreteMeasureData)" begin
        orig_prefs = parameter_list(inf7)
        support_dict = Dict{Int, Float64}(1 => 23, 2 => 42)
        expected = Float64[23, 42, 3, 2]
        vals = Float64[3, 2]
        @test InfiniteOpt._make_point_support(orig_prefs, support_dict, [3, 4], vals) == expected
    end
    # test expand_measure (reduced infinite variable)
    @testset "Reduced Variable (Multi DiscreteMeasureData)" begin
        # test array param reduced var without measure param
        rv = make_reduced_variable_ref(m, inf5, [1, 2], Float64[1, 1])
        @test InfiniteOpt.expand_measure(rv, data3, m) == rv + rv
        # test array param reduced var with measure param
        idx = length(InfiniteOpt._data_dictionary(m, PointVariable)) + 1
        pts = [GeneralVariableRef(m, idx, PointVariableIndex),
               GeneralVariableRef(m, idx + 1, PointVariableIndex)]
        @test InfiniteOpt.expand_measure(rv, data4, m) == 2 * (pts[1] + pts[2])
        @test parameter_values(pts[2]) == (Float64[1, 1], Float64[1, 1])
        @test !is_valid(m, rv)
        # test array param reduced var with measure param and others
        rv = make_reduced_variable_ref(m, inf7, [1], Float64[1])
        idx = index(rv).value + 1
        rv1 = GeneralVariableRef(m, idx, ReducedVariableIndex)
        rv2 = GeneralVariableRef(m, idx + 1, ReducedVariableIndex)
        @test InfiniteOpt.expand_measure(rv, data3, m) == rv1 + rv2
        @test eval_supports(rv2) == Dict{Int, Float64}(1 => 1, 3 => 2, 4 => 2)
        @test !is_valid(m, rv)
        # test array param that is out of order and should make point variable
        rv = make_reduced_variable_ref(m, inf5, [3, 4], Float64[1, 1])
        idx = length(InfiniteOpt._data_dictionary(m, PointVariable)) + 1
        pts = [GeneralVariableRef(m, idx, PointVariableIndex),
               GeneralVariableRef(m, idx + 1, PointVariableIndex)]
        @test InfiniteOpt.expand_measure(rv, data5, m) == pts[1] + pts[2]
        @test parameter_values(pts[1]) == (Float64[1.2, 1], Float64[1, 1])
        @test !is_valid(m, rv)
        # test array param that with no others that partially overlap
        rv = make_reduced_variable_ref(m, inf5, [2, 3, 4], Float64[1, 1, 1])
        @test InfiniteOpt.expand_measure(rv, data6, m) == 2rv
        # test array param with others that partially overlap
        rv = make_reduced_variable_ref(m, inf5, [2], Float64[1])
        @test InfiniteOpt.expand_measure(rv, data6, m) == 2rv
        # test partial overlap
        rv = make_reduced_variable_ref(m, inf5, [2], Float64[1])
        idx = index(rv).value + 1
        rv1 = GeneralVariableRef(m, idx, ReducedVariableIndex)
        rv2 = GeneralVariableRef(m, idx + 1, ReducedVariableIndex)
        @test InfiniteOpt.expand_measure(rv, data3, m) == rv1 + rv2
        @test eval_supports(rv2) == Dict{Int, Float64}(2 => 1, 1 => 2)
        @test !is_valid(m, rv)
    end
    # test expand_measure (finite variable)
    @testset "Finite Variable (1D DiscreteMeasureData)" begin
        # test with single parameter measure
        @test InfiniteOpt.expand_measure(x, data1, m) == 0.5 * (x + x)
    end
    # test expand_measure (finite variable)
    @testset "Finite Variable (Multi DiscreteMeasureData)" begin
        # test with multi parameter measure
        @test InfiniteOpt.expand_measure(x, data3, m) == x + x
    end
    # test expand_measure (parameter)
    @testset "Infinite Parameter (1D DiscreteMeasureData)" begin
        # test with different parameter
        @test InfiniteOpt.expand_measure(par1, data2, m) == 1.5 * (par1 + par1)
        # test with same parameter
        @test InfiniteOpt.expand_measure(par1, data1, m) == 0.5 * (1. + 2.)
    end
    # test expand_measure (multi measure data with parameter)
    @testset "Infinite Parameter (Multi DiscreteMeasureData)" begin
        # test with different parameter
        @test InfiniteOpt.expand_measure(par1, data3, m) == par1 + par1
        # test with same parameter
        @test InfiniteOpt.expand_measure(pars1[2], data3, m) == 3.
    end
    # test expand_measure (Finite Expr) with DiscreteMeasureData
    @testset "Finite Expression (1D DiscreteMeasureData)"  begin
        # test AffExpr
        expr = 2x +3
        @test InfiniteOpt.expand_measure(expr, data2, m) == expr * 3
        # test QuadExpr
        expr = x^2 + 2x +3
        @test InfiniteOpt.expand_measure(expr, data2, m) == expr * 3
    end
    # test expand_measure (Finite Expr) with DiscreteMeasureData
    @testset "Finite Expression (Multi DiscreteMeasureData)"  begin
        # test AffExpr
        expr = 2x +3
        @test InfiniteOpt.expand_measure(expr, data3, m) == expr * 2
        # test QuadExpr
        expr = x^2 + 2x +3
        @test InfiniteOpt.expand_measure(expr, data4, m) == expr * 4
    end
    # test expand_measure (AffExpr)
    @testset "AffExpr (1D DiscreteMeasureData)" begin
        # test single param AffExpr, no measures
        expr = 2inf1 + par1 - x + 3par2 - 3
        idx = length(InfiniteOpt._data_dictionary(m, PointVariable)) + 1
        pts = [GeneralVariableRef(m, idx, PointVariableIndex),
               GeneralVariableRef(m, idx + 1, PointVariableIndex)]
        expected = 0.5 * (2pts[1] + 2pts[2] + 1 + 2 - x - x + 3par2 + 3par2 - 6)
        @test InfiniteOpt.expand_measure(expr, data1, m) == expected
        # test single param AffExpr, with measures
        expr = meas2 - x + par1
        idx = length(InfiniteOpt._data_dictionary(m, PointVariable)) + 1
        pts = [GeneralVariableRef(m, idx, PointVariableIndex),
               GeneralVariableRef(m, idx + 1, PointVariableIndex)]
        expected = 0.5 * (6 * (pts[1] * x + pts[2] * x) - 2x - 6)
        @test InfiniteOpt.expand_measure(expr, data1, m) == expected
    end
    # test expand_measure (AffExpr)
    @testset "AffExpr (Multi DiscreteMeasureData)" begin
        # test array param AffExpr, no measures
        expr = inf4 + par1 - x + 3pars1[2] - 1
        idx = length(InfiniteOpt._data_dictionary(m, PointVariable)) + 1
        pts = [GeneralVariableRef(m, idx, PointVariableIndex),
               GeneralVariableRef(m, idx + 1, PointVariableIndex)]
        expected = pts[1] + pts[2] + 2par1 - 2x + 9 - 2
        @test InfiniteOpt.expand_measure(expr, data3, m) == expected
        # test array param AffExpr, with measures
        expr = meas2 - x + par1 - inf4
        idx = length(InfiniteOpt._data_dictionary(m, PointVariable)) + 1
        pts = [GeneralVariableRef(m, idx, PointVariableIndex),
               GeneralVariableRef(m, idx + 1, PointVariableIndex),
               GeneralVariableRef(m, idx + 2, PointVariableIndex),
               GeneralVariableRef(m, idx + 3, PointVariableIndex)]
        expected = 6 * (pts[1] * x + pts[2] * x) - 4par1 - 2x - pts[3] - pts[4]
        @test InfiniteOpt.expand_measure(expr, data3, m) == expected
    end
    # test expand_measure (QuadExpr)
    @testset "QuadExpr (1D DiscreteMeasureData)" begin
        # test single param QuadExpr with both variables integrated or not
        expr = 2 * inf1 * inf2 - inf3 * inf4 + x + 2
        idx = length(InfiniteOpt._data_dictionary(m, PointVariable)) + 1
        pts = [GeneralVariableRef(m, idx, PointVariableIndex),
               GeneralVariableRef(m, idx + 1, PointVariableIndex)]
        idx = InfiniteOpt._data_dictionary(m, ReducedVariable).last_index + 1
        rv1 = GeneralVariableRef(m, idx, ReducedVariableIndex)
        rv2 = GeneralVariableRef(m, idx + 1, ReducedVariableIndex)
        expected = 0.5 * (2 * pts[1] * rv1 + - 2 * inf3 * inf4 + 2 * pts[2] *
                          rv2 + 2x + 4)
        @test InfiniteOpt.expand_measure(expr, data1, m) == expected
        # test single param QuadExpr with first variable not integrated
        expr = 3 * inf3 * inf1 + pars1[1] - 1
        idx = length(InfiniteOpt._data_dictionary(m, PointVariable)) + 1
        pts = [GeneralVariableRef(m, idx, PointVariableIndex),
               GeneralVariableRef(m, idx + 1, PointVariableIndex)]
        expected = 0.5 * (3 * inf3 * pts[1] + 3 * inf3 * pts[2] + 2pars1[1] - 2)
        @test InfiniteOpt.expand_measure(expr, data1, m) == expected
        # test single param QuadExpr with first variable integrated
        expr = 3 * inf2 * inf1 + pars1[1] + 1
        idx = InfiniteOpt._data_dictionary(m, ReducedVariable).last_index + 1
        rv1 = GeneralVariableRef(m, idx, ReducedVariableIndex)
        rv2 = GeneralVariableRef(m, idx + 1, ReducedVariableIndex)
        expected = 1.5 * (3 * rv1 * inf1 + 3 * rv2 * inf1 + 2pars1[1] + 2)
        @test InfiniteOpt.expand_measure(expr, data2, m) == expected
        # test single parameter with first quadratic term becomes number
        expr = par1 * inf1 + 2
        idx = length(InfiniteOpt._data_dictionary(m, PointVariable)) + 1
        pts = [GeneralVariableRef(m, idx, PointVariableIndex),
               GeneralVariableRef(m, idx + 1, PointVariableIndex)]
        expected = 0.5 * (pts[1] + 2 * pts[2] + 4)
        @test InfiniteOpt.expand_measure(expr, data1, m) == expected
        # test single parameter with first quadratic term becomes number, 2nd constant
        expr = par1 * x + 2
        expected = 0.5 * (x + 2 * x + 4)
        @test InfiniteOpt.expand_measure(expr, data1, m) == expected
        # test single parameter with both quadratic terms become numbers
        expr = par1 * par1 + 2
        expected = 0.5 * (1 + 2 * 2 + 4)
        @test InfiniteOpt.expand_measure(expr, data1, m) == Float64(expected)
        # test single parameter with 2nd quadratic term becomes number
        expr = inf1 * par1 + 2
        idx = length(InfiniteOpt._data_dictionary(m, PointVariable)) + 1
        pts = [GeneralVariableRef(m, idx, PointVariableIndex),
               GeneralVariableRef(m, idx + 1, PointVariableIndex)]
        expected = 0.5 * (pts[1] + 2 * pts[2] + 4)
        @test InfiniteOpt.expand_measure(expr, data1, m) == expected
        # test single parameter with 2nd quadratic term becomes number, 1st constant
        expr = x * par1 + 2
        expected = 0.5 * (x + 2 * x + 4)
        @test InfiniteOpt.expand_measure(expr, data1, m) == expected
    end
    # test expand_measure (QuadExpr)
    @testset "QuadExpr (Multi DiscreteMeasureData)" begin
        # test array param QuadExpr with both variables integrated
        expr = 2 * inf4 * inf5 - inf1 * inf2 + x + 2
        idx = length(InfiniteOpt._data_dictionary(m, PointVariable)) + 1
        pts = [GeneralVariableRef(m, idx, PointVariableIndex),
               GeneralVariableRef(m, idx + 1, PointVariableIndex)]
               idx = InfiniteOpt._data_dictionary(m, ReducedVariable).last_index + 1
               rv1 = GeneralVariableRef(m, idx, ReducedVariableIndex)
               rv2 = GeneralVariableRef(m, idx + 1, ReducedVariableIndex)
        expected = 2 * pts[1] * rv1 - 2 * inf1 * inf2 + 2 * pts[2] * rv2 + 2x + 4
        @test InfiniteOpt.expand_measure(expr, data3, m) == expected
        # test array param QuadExpr with first variable not integrated
        expr = 3 * inf1 * inf4 + par1 - 1
        idx = length(InfiniteOpt._data_dictionary(m, PointVariable)) + 1
        pts = [GeneralVariableRef(m, idx, PointVariableIndex),
               GeneralVariableRef(m, idx + 1, PointVariableIndex)]
        expected = 3 * inf1 * pts[1] + 3 * inf1 * pts[2] + 2par1 - 2
        @test InfiniteOpt.expand_measure(expr, data3, m) == expected
        # test array param QuadExpr with first variable integrated
        expr = 3 * inf5 * inf1 + pars1[1] + 1
        idx = InfiniteOpt._data_dictionary(m, ReducedVariable).last_index + 1
        rv1 = GeneralVariableRef(m, idx, ReducedVariableIndex)
        rv2 = GeneralVariableRef(m, idx + 1, ReducedVariableIndex)
        expected = 2 * (3 * rv1 * inf1 + 3 * rv2 * inf1 + 2pars1[1] + 2)
        @test InfiniteOpt.expand_measure(expr, data4, m) == expected
        # test array parameter with first quadratic term becomes number
        expr = pars1[1] * inf4 + 2
        idx = length(InfiniteOpt._data_dictionary(m, PointVariable)) + 1
        pts = [GeneralVariableRef(m, idx, PointVariableIndex),
               GeneralVariableRef(m, idx + 1, PointVariableIndex)]
        expected = pts[1] + 2 * pts[2] + 4
        @test InfiniteOpt.expand_measure(expr, data3, m) == expected
        # test array parameter with first quadratic term becomes number, 2nd constant
        expr = pars1[1] * x + 2
        expected = x + 2 * x + 4
        @test InfiniteOpt.expand_measure(expr, data3, m) == expected
        # test array parameter with both quadratic terms become numbers
        expr = pars1[1] * pars1[2] + 2
        expected = 1 + 2 * 2 + 4
        @test InfiniteOpt.expand_measure(expr, data3, m) == Float64(expected)
        # test array parameter with 2nd quadratic term becomes number
        expr = inf4 * pars1[1] + 2
        idx = length(InfiniteOpt._data_dictionary(m, PointVariable)) + 1
        pts = [GeneralVariableRef(m, idx, PointVariableIndex),
               GeneralVariableRef(m, idx + 1, PointVariableIndex)]
        expected = pts[1] + 2 * pts[2] + 4
        @test InfiniteOpt.expand_measure(expr, data3, m) == expected
        # test array parameter with 2nd quadratic term becomes number, 1st constant
        expr = x * pars1[1] + 2
        expected = x + 2 * x + 4
        @test InfiniteOpt.expand_measure(expr, data3, m) == expected
    end
    # test expand_measure (measure)
    @testset "Measure" begin
        # test single param measure
        idx = length(InfiniteOpt._data_dictionary(m, PointVariable)) + 1
        pts = [GeneralVariableRef(m, idx, PointVariableIndex),
               GeneralVariableRef(m, idx + 1, PointVariableIndex)]
        expected = 0.5 * (pts[1] + pts[2])
        @test InfiniteOpt.expand_measure(meas1, data1, m) == expected
        # test another single param
        idx = length(InfiniteOpt._data_dictionary(m, PointVariable)) + 1
        pts = [GeneralVariableRef(m, idx, PointVariableIndex),
               GeneralVariableRef(m, idx + 1, PointVariableIndex)]
        expected = 0.5 * (6 * pts[1] * x + 6 * pts[2] * x - 3 * (1 + 2))
        @test InfiniteOpt.expand_measure(meas2, data1, m) == expected
        # test array param measure
        idx = length(InfiniteOpt._data_dictionary(m, PointVariable)) + 1
        pts = [GeneralVariableRef(m, idx, PointVariableIndex),
               GeneralVariableRef(m, idx + 1, PointVariableIndex)]
        expected = pts[1] + pts[2]
        @test InfiniteOpt.expand_measure(meas1, data3, m) == expected
    end
    # test expand_measure (FunctionalDiscreteMeasureData)
    @testset "FunctionalDiscreteMeasureData" begin
        # test without generative supports
        coef(a) = ones(length(a))
        data = FunctionalDiscreteMeasureData(par1, coef, 2, All, is_expect = true)
        @test InfiniteOpt.expand_measure(x, data, m) == 2x
        # test with generative supports
        coef2(a) = ones(length(a) + 1)
        info = UniformGenerativeInfo([0.5], InternalLabel)
        data = FunctionalDiscreteMeasureData(par1, coef2, 2, All, generative_support_info = info)
        InfiniteOpt._set_generative_support_info(dispatch_variable_ref(par1), info)
        @test InfiniteOpt.expand_measure(x, data, m) == 3x
    end
    # test expand_measure (other)
    @testset "Other" begin
        # prepare test
        @variable(Model(), y)
        # test it
        @test_throws ErrorException InfiniteOpt.expand_measure(y, BadData(), m)
    end
end

# Test analytic_expansion
@testset "analytic_expansion" begin
    # initialize model and references
    m = InfiniteModel()
    @infinite_parameter(m, 1 <= par1 <= 2)
    @infinite_parameter(m, 1 <= par2 <= 2)
    @infinite_parameter(m, 1 <= pars1[1:2] <= 2, num_supports = 2)
    @infinite_variable(m, inf3(par2))
    @hold_variable(m, x)
    # test with 1D data
    @testset "1D Discrete Data" begin
        # test with bounds
        data = DiscreteMeasureData(par1, [1], [1], lower_bound = 1,
                                   upper_bound = 1.5)
        @test InfiniteOpt.analytic_expansion(2x + inf3, data, m) == 0.5 * (2x + inf3)
        # test as expectation
        coef(a) = ones(length(a))
        data = FunctionalDiscreteMeasureData(par2, coef, 0, WeightedSample,
                                             is_expect = true)
        @test InfiniteOpt.analytic_expansion(inf3, data, m) == inf3
        # test other
        data = DiscreteMeasureData(par1, [1, 0.5], [1, 1.5])
        @test InfiniteOpt.analytic_expansion(x, data, m) == 1.5 * x
    end
    # test with multi dimensional data
    @testset "Multi Discrete Data" begin
        # test with bounds
        coef(a) = ones(length(a))
        data = FunctionalDiscreteMeasureData(pars1, coef, 0, All,
                    lower_bounds = [1, 1], upper_bounds = [1.5, 2])
        @test InfiniteOpt.analytic_expansion(x, data, m) == 0.5x
        # test as expectation
        data = DiscreteMeasureData(pars1, [1], [[1, 1]], is_expect = true)
        @test InfiniteOpt.analytic_expansion(inf3, data, m) == inf3
        # test other
        coef(a) = ones(size(a, 2))
        data = FunctionalDiscreteMeasureData(pars1, coef, 2, All)
        @test InfiniteOpt.analytic_expansion(x, data, m) == 2x
    end
    # test with fallback
    @testset "Fallback" begin
        data = TestData(par1, 1, 2)
        @test_throws ErrorException InfiniteOpt.analytic_expansion(x, data, m)
    end
end

# Test expand_measures
@testset "expand_measures" begin
    # initialize model and references
    m = InfiniteModel()
    @infinite_parameter(m, 1 <= par1 <= 2)
    @infinite_parameter(m, 1 <= par2 <= 2)
    @infinite_parameter(m, 1 <= pars1[1:2] <= 2)
    @infinite_variable(m, inf1(par1) >= 1)
    @infinite_variable(m, inf2(par1, par2))
    @infinite_variable(m, inf3(par2))
    @infinite_variable(m, inf4(pars1))
    @hold_variable(m, x)
    # prepare measures
    data1 = DiscreteMeasureData(par1, [0.5, 0.5], [1, 2], is_expect = true)
    data2 = DiscreteMeasureData(pars1, [1, 1], [[1, 1], [2, 2]])
    meas1 = measure(inf1 + 3x - inf2 + inf3 - 2, data1)
    meas2 = measure(2 * inf4 * x - pars1[2] + inf2, data2)
    # test expand_measures (measure)
    @testset "MeasureRef" begin
        # test the first measure
        idx = length(InfiniteOpt._data_dictionary(m, PointVariable)) + 1
        pts = [GeneralVariableRef(m, idx, PointVariableIndex),
               GeneralVariableRef(m, idx + 1, PointVariableIndex)]
        idx = InfiniteOpt._data_dictionary(m, ReducedVariable).last_index + 1
        rv1 = GeneralVariableRef(m, idx, ReducedVariableIndex)
        rv2 = GeneralVariableRef(m, idx + 1, ReducedVariableIndex)
        expected = 0.5 * (pts[1] + pts[2] + 6x - rv1 - rv2 + 2inf3 - 4)
        @test InfiniteOpt.expand_measures(meas1, m) == expected
        # test the second measure
        idx = length(InfiniteOpt._data_dictionary(m, PointVariable)) + 1
        pts = [GeneralVariableRef(m, idx, PointVariableIndex),
               GeneralVariableRef(m, idx + 1, PointVariableIndex)]
        expected = 2 * pts[1] * x + 2 * pts[2] * x - 3 + 2inf2
        @test InfiniteOpt.expand_measures(meas2, m) == expected
        # test analytic expansion
        meas3 = measure(x, data1)
        @test InfiniteOpt.expand_measures(meas3, m) == x
    end
    # test expand_measures (GeneralVariableRef)
    @testset "GeneralVariableRef" begin
        @test InfiniteOpt.expand_measures(x, m) == x
        @test InfiniteOpt.expand_measures(par1, m) == par1
        @test InfiniteOpt.expand_measures(inf1, m) == inf1
    end
    # test expand_measures (Finite Expression)
    @testset "Expression (Finite)" begin
        expr = 2x + 2
        @test InfiniteOpt.expand_measures(expr, m) == expr
        expr = x^2 + 2x + 2
        @test InfiniteOpt.expand_measures(expr, m) == expr
    end
    # test expand_measures (GenericAffExpr)
    @testset "AffExpr (General)" begin
        # test returning AffExpr
        expr = 2inf4 + x + 2meas1
        idx = length(InfiniteOpt._data_dictionary(m, PointVariable)) + 1
        pts = [GeneralVariableRef(m, idx, PointVariableIndex),
               GeneralVariableRef(m, idx + 1, PointVariableIndex)]
        idx = InfiniteOpt._data_dictionary(m, ReducedVariable).last_index + 1
        rv1 = GeneralVariableRef(m, idx, ReducedVariableIndex)
        rv2 = GeneralVariableRef(m, idx + 1, ReducedVariableIndex)
        expected = 2inf4 + 7x + pts[1] + pts[2] - rv1 - rv2 + 2inf3 - 4
        @test InfiniteOpt.expand_measures(expr, m) == expected
        # test returning QuadExpr
        expr = x + 3par1 + meas2
        idx = length(InfiniteOpt._data_dictionary(m, PointVariable)) + 1
        pts = [GeneralVariableRef(m, idx, PointVariableIndex),
               GeneralVariableRef(m, idx + 1, PointVariableIndex)]
        expected = x + 3par1 + 2 * pts[1] * x + 2 * pts[2] * x - 3 + 2inf2
        @test InfiniteOpt.expand_measures(expr, m) == expected
    end
    # test expand_measures (GenericQuadExpr)
    @testset "QuadExpr (General)" begin
        # prepare first measure expansion
        expr = 2 * meas1 * x - x + 1
        idx = length(InfiniteOpt._data_dictionary(m, PointVariable)) + 1
        pts = [GeneralVariableRef(m, idx, PointVariableIndex),
               GeneralVariableRef(m, idx + 1, PointVariableIndex)]
        idx = InfiniteOpt._data_dictionary(m, ReducedVariable).last_index + 1
        rv1 = GeneralVariableRef(m, idx, ReducedVariableIndex)
        rv2 = GeneralVariableRef(m, idx + 1, ReducedVariableIndex)
        m1 = 0.5 * (pts[1] + pts[2] + 6x - rv1 - rv2 + 2inf3 - 4)
        # test with this thorough expression
        @test InfiniteOpt.expand_measures(expr, m) == 2 * m1 * x - x + 1
    end
    # test expand_measures (Fallback)
    @testset "Fallback" begin
        @test_throws ErrorException InfiniteOpt.expand_measures(1, m)
    end
end

# Test user expansion methods
@testset "User Methods" begin
    # initialize model and references
    m = InfiniteModel()
    @infinite_parameter(m, 1 <= par1 <= 2)
    @infinite_parameter(m, 1 <= par2 <= 2)
    @infinite_parameter(m, 1 <= pars1[1:2] <= 2)
    @infinite_variable(m, inf1(par1) >= 1)
    @infinite_variable(m, inf2(par1, par2) >= 0)
    @infinite_variable(m, inf3(par2))
    @infinite_variable(m, inf4(pars1))
    @hold_variable(m, x >= 0)
    # prepare measures
    data1 = DiscreteMeasureData(par1, [0.5, 0.5], [1, 2], lower_bound = 1, upper_bound = 2)
    data2 = DiscreteMeasureData(pars1, [1, 1], [[1, 1], [2, 2]])
    data3 = DiscreteMeasureData(par2, [2, 2], [1, 2])
    meas1 = measure(inf1 + 3x - inf2 + inf3 - 2, data1)
    meas2 = measure(2 * inf4 * x - pars1[2] + inf2, data2)
    # test expand
    @testset "expand" begin
        # test the first measure
        idx = length(InfiniteOpt._data_dictionary(m, PointVariable)) + 1
        pts = [GeneralVariableRef(m, idx, PointVariableIndex),
               GeneralVariableRef(m, idx + 1, PointVariableIndex)]
        idx = InfiniteOpt._data_dictionary(m, ReducedVariable).last_index + 1
        rv1 = GeneralVariableRef(m, idx, ReducedVariableIndex)
        rv2 = GeneralVariableRef(m, idx + 1, ReducedVariableIndex)
        expected = 0.5 * (pts[1] + pts[2] + 6x - rv1 - rv2 + 2inf3 - 4)
        @test expand(meas1) == expected
        @test parameter_values(pts[2]) == (2.,)
        @test eval_supports(rv1) == Dict{Int, Float64}(1 => 1)
        # test the second measure
        idx = length(InfiniteOpt._data_dictionary(m, PointVariable)) + 1
        pts = [GeneralVariableRef(m, idx, PointVariableIndex),
               GeneralVariableRef(m, idx + 1, PointVariableIndex)]
        expected = 2 * pts[1] * x + 2 * pts[2] * x - 3 + 2inf2
        @test expand(meas2) == expected
        @test parameter_values(pts[2]) == (Float64[2, 2],)
        # test analytic
        meas = measure(x, data1)
        @test is_analytic(meas)
        @test expand(meas) == 1*x
        # test integral with only 1 support
        delete(m, meas1)
        delete(m, meas)
        delete_supports(par1)
        add_supports(par1, 1)
        @test_throws ErrorException expand(integral(inf1, par1))
    end

    # test expand_all_measures!
    @testset "expand_all_measures!" begin
        # prepare the model for testing
        @hold_variable(m, y, parameter_bounds = (par2 == 1))
        set_objective_function(m, x + measure(inf1 * par1 + 3, data1))
        @constraint(m, c1, inf1 + x >= 42.)
        @constraint(m, c2, 2x - measure(measure(inf1 * x + par1 + inf2, data1),
                                        data3) == 0.)
        @BDconstraint(m, c3(par2 in [1, 1.5]), measure(inf2, data1) + inf1 >= 0.)
        @constraint(m, c4, measure(inf4 + y, data2) <= 3.)
        # prepare comparison expressions
        idx = length(InfiniteOpt._data_dictionary(m, PointVariable)) + 1
        pts = [GeneralVariableRef(m, idx, PointVariableIndex),
               GeneralVariableRef(m, idx + 1, PointVariableIndex)]
        obj = x + 0.5pts[1] + pts[2] + 3
        pts = [GeneralVariableRef(m, idx + 2, PointVariableIndex),
               GeneralVariableRef(m, idx + 3, PointVariableIndex),
               GeneralVariableRef(m, idx + 4, PointVariableIndex),
               GeneralVariableRef(m, idx + 5, PointVariableIndex),
               GeneralVariableRef(m, idx + 6, PointVariableIndex),
               GeneralVariableRef(m, idx + 7, PointVariableIndex)]
        c2_expected = 2x - (2 * pts[1] * x + 2 * pts[2] * x + pts[3] + pts[4] +
                      pts[5] + pts[6])
        idx2 = InfiniteOpt._data_dictionary(m, ReducedVariable).last_index + 3
        rv1 = GeneralVariableRef(m, idx2, ReducedVariableIndex)
        rv2 = GeneralVariableRef(m, idx2 + 1, ReducedVariableIndex)
        c3_expected = 0.5rv1 + 0.5rv2 + inf1
        pts = [GeneralVariableRef(m, idx + 8, PointVariableIndex),
               GeneralVariableRef(m, idx + 9, PointVariableIndex)]
        c4_expected = pts[1] + pts[2] + 2y
        # test the expansion
        @test isa(expand_all_measures!(m), Nothing)
        @test objective_function(m) == obj
        @test is_valid(m, c2)
        @test is_valid(m, c3)
        @test is_valid(m, c4)
        @test name(c1) == "c1"
        @test name(c2) == "c2"
        @test name(c3) == "c3"
        @test name(c4) == "c4"
        @test constraint_object(c1).func == inf1 + x
        @test constraint_object(c2).func == c2_expected
        @test constraint_object(c3).func == c3_expected
        @test constraint_object(c4).func == c4_expected
        @test constraint_object(c2).func isa JuMP.GenericQuadExpr{Float64, GeneralVariableRef}
        @test constraint_object(c3).func isa JuMP.GenericAffExpr{Float64, GeneralVariableRef}
        @test constraint_object(c4).func isa JuMP.GenericAffExpr{Float64, GeneralVariableRef}
        @test constraint_object(c1).set == MOI.GreaterThan(42.)
        @test constraint_object(c2).set == MOI.EqualTo(6.)
        @test constraint_object(c3).set == MOI.GreaterThan(0.)
        @test constraint_object(c4).set == MOI.LessThan(3.)
        @test isa(constraint_object(c1), ScalarConstraint)
        @test isa(constraint_object(c2), ScalarConstraint)
        @test isa(constraint_object(c3), BoundedScalarConstraint)
        @test isa(constraint_object(c4), BoundedScalarConstraint)
        @test constraint_object(c3).bounds.intervals == Dict(par2 => IntervalSet(1, 1.5))
        @test constraint_object(c4).bounds.intervals == Dict(par2 => IntervalSet(1, 1))
    end
end


text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

struct Relation <: AbstractRelation
    iid::String
    type::RelationType
end

function Relation(t::Proto.Thing)
    iid = bytes2hex(t.iid)
    isempty(iid) && throw(TypeDBClientException(CONCEPT_MISSING_IID))
    return Relation(iid, instantiate(t._type))
end

# TODO seems unnecessary?
as_relation(r::Relation) = r

function add_player(trans::AbstractCoreTransaction,
                    relation::AbstractRelation,
                    role::AbstractRoleType,
                    player_thing::AbstractThing)

    proto_thing = get_proto_thing(ConceptManager(trans), player_thing.iid)
    player_req = RelationRequestBuilder.add_player_req(relation.iid, proto(role), proto_thing)

    execute(trans, player_req)

    return nothing
end

function remove_player(transaction::AbstractCoreTransaction,
                            relation::AbstractRelation,
                            role_type::AbstractRoleType,
                            player_thing::AbstractThing)

    proto_player = get_proto_thing(ConceptManager(transaction), player_thing.iid)
    rem_player_req = RelationRequestBuilder.remove_player_req(relation.iid, proto(role_type), proto_player)
    execute(transaction, rem_player_req)

    return nothing
end

function get_players(transaction::AbstractCoreTransaction,
                        relation::AbstractRelation,
                        roles::Optional{Vector{<:AbstractRoleType}}= nothing)

    role_types = roles !== nothing ? proto.(roles) : nothing

    play_req = RelationRequestBuilder.get_players_req(relation.iid, role_types)
    player_res = stream(transaction, play_req)

    return instantiate.(collect(Iterators.flatten(
        r.thing_res_part.relation_get_players_res_part.things for r in player_res)))
end

function get_relations(
    transaction::AbstractCoreTransaction,
    thing::AbstractThing,
    role_types::Optional{AbstractVector{<:AbstractRoleType}}=nothing)

    proto_role_types = role_types !== nothing ? proto.(role_types) : nothing

    rel_req = ThingRequestBuilder.get_relations_req(thing.iid, proto_role_types)
    res_rel = stream(transaction, rel_req)
    return return instantiate.(collect(Iterators.flatten(
        r.thing_res_part.thing_get_relations_res_part.relations for r in res_rel)))
end


text:
dxt1 .- [-0.40429154961292091
-0.00028833951030255518
3.1210170895092109
-0.00035116191491866704
-2.3890914916989998
-0.00041003907445329613]

0.0
5.421010862427522e-20
0.0
0.0
0.0
0.000121699564150741

j20::Float64 = PZ90Ellipsoid.j20,
mu::Float64 = PZ90Ellipsoid.gm_km,
ae::Float64 = PZ90Ellipsoid.a_km

j20 = -1.08262575e-3
mu = 398600.4418
ae = 6378136.0 / 1000.0

initialState = [12023.740191250332
    -0.4042915496129209
 14643.53682860572
     3.121017089509211
 17095.20605469
    -2.389091491699]

accel = [ -1.8369219039096033e-9
  3.084880401099731e-10
 -9.313225746155e-10]


text:
"""
    countBinaryOnes(ind::Individual)

documentation
"""
function countBinaryOnes(genotype::BinaryGenotype)
    genotype = genotype._representation
    tam = length(genotype)
    count=0
    for i=1:tam
        if genotype[i]==1
            count+=1
        end
    end
    return count
end # function

"""
    countBinaryOnes(ind::Individual)

documentation
"""
function countBinaryZeros(genotype::BinaryGenotype)
    genotype = genotype._representation
    genLen = length(genotype)
    count=0
    for i=1:genLen
        if genotype[i]==0
            count+=1
        end
    end
    return count
end # function

"""
    llegar50(genotype::IntegerGenotype)

documentation
"""
function llegar50(genotype::IntegerGenotype)
    genotype = genotype._representation
    genLen = length(genotype)
    acc = 0;
    for i=1:genLen
        acc+=genotype[i]
    end
    return abs(acc-50)
end # function

function checkVelocidad(genotype::IntegerGenotype)
    genotype = genotype._representation

    function probar1(iteraciones)

        for i=1:iteraciones
            randomIndexSelection2(genotype[1],genotype[2],rng)
        end

    end

    function probar2(iteraciones)

        for i=1:iteraciones
            randomIndexSelection(genotype[1],genotype[2],rng)
            #randomIndexSelection2(100,75)
        end

    end
    val, a, bytes, gctime, memallocs = @timed probar1(10000)
    val, b, bytes, gctime, memallocs = @timed probar2(10000)

    abs(a-b)
end

function nNodes(canonicalgp::GPGenotype, gpInfo)
    length(canonicalgp._representation)
end

max_stars_repo_name: SergioGmezM/GeneticJulia

max_issues_repo_name: SergioGmezM/GeneticJulia

max_forks_repo_name: SergioGmezM/GeneticJulia


text:
"""
GenX: An Configurable Capacity Expansion Model
Copyright (C) 2021,  Massachusetts Institute of Technology
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
A complete copy of the GNU General Public License v2 (GPLv2) is available
in LICENSE.txt.  Users uncompressing this from an archive may not have
received this license file.  If not, see <http://www.gnu.org/licenses/>.
"""

@doc raw"""
	storage_all(EP::Model, inputs::Dict, Reserves::Int, OperationWrapping::Int)

Sets up variables and constraints common to all storage resources. See ```storage()``` in ```storage.jl``` for description of constraints.
"""
function storage_all(EP::Model, inputs::Dict, Reserves::Int, OperationWrapping::Int)
	# Setup variables, constraints, and expressions common to all storage resources
	println("Storage Core Resources Module")

	dfGen = inputs["dfGen"]

	G = inputs["G"]     # Number of resources (generators, storage, DR, and DERs)
	T = inputs["T"]     # Number of time steps (hours)
	Z = inputs["Z"]     # Number of zones

	STOR_ALL = inputs["STOR_ALL"]
	STOR_SHORT_DURATION = inputs["STOR_SHORT_DURATION"]

	START_SUBPERIODS = inputs["START_SUBPERIODS"]
	INTERIOR_SUBPERIODS = inputs["INTERIOR_SUBPERIODS"]

	hours_per_subperiod = inputs["hours_per_subperiod"] #total number of hours per subperiod

	### Variables ###

	# Storage level of resource "y" at hour "t" [MWh] on zone "z" - unbounded
	@variable(EP, vS[y in STOR_ALL, t=1:T] >= 0);

	# Energy withdrawn from grid by resource "y" at hour "t" [MWh] on zone "z"
	@variable(EP, vCHARGE[y in STOR_ALL, t=1:T] >= 0);

	### Expressions ###

	# Energy losses related to technologies (increase in effective demand)
	@expression(EP, eELOSS[y in STOR_ALL], sum(inputs["omega"][t]*EP[:vCHARGE][y,t] for t in 1:T) - sum(inputs["omega"][t]*EP[:vP][y,t] for t in 1:T))

	## Objective Function Expressions ##

	#Variable costs of "charging" for technologies "y" during hour "t" in zone "z"
	@expression(EP, eCVar_in[y in STOR_ALL,t=1:T], inputs["omega"][t]*dfGen[!,:Var_OM_Cost_per_MWh_In][y]*vCHARGE[y,t])

	# Sum individual resource contributions to variable charging costs to get total variable charging costs
	@expression(EP, eTotalCVarInT[t=1:T], sum(eCVar_in[y,t] for y in STOR_ALL))
	@expression(EP, eTotalCVarIn, sum(eTotalCVarInT[t] for t in 1:T))
	EP[:eObj] += eTotalCVarIn

	## Power Balance Expressions ##

	# Term to represent net dispatch from storage in any period
	@expression(EP, ePowerBalanceStor[t=1:T, z=1:Z],
		sum(EP[:vP][y,t]-EP[:vCHARGE][y,t] for y in intersect(dfGen[dfGen.Zone.==z,:R_ID],STOR_ALL)))

	EP[:ePowerBalance] += ePowerBalanceStor

	### Constraints ###

	## Storage energy capacity and state of charge related constraints:

	# Links state of charge in first time step with decisions in last time step of each subperiod
	# We use a modified formulation of this constraint (cSoCBalLongDurationStorageStart) when operations wrapping and long duration storage are being modeled

	if OperationWrapping ==1 && !isempty(inputs["STOR_LONG_DURATION"])
		@constraint(EP, cSoCBalStart[t in START_SUBPERIODS, y in STOR_SHORT_DURATION], EP[:vS][y,t] ==
			EP[:vS][y,t+hours_per_subperiod-1]-(1/dfGen[!,:Eff_Down][y]*EP[:vP][y,t])
			+(dfGen[!,:Eff_Up][y]*EP[:vCHARGE][y,t])-(dfGen[!,:Self_Disch][y]*EP[:vS][y,t+hours_per_subperiod-1]))
	else
		@constraint(EP, cSoCBalStart[t in START_SUBPERIODS, y in STOR_ALL], EP[:vS][y,t] ==
			EP[:vS][y,t+hours_per_subperiod-1]-(1/dfGen[!,:Eff_Down][y]*EP[:vP][y,t])
			+(dfGen[!,:Eff_Up][y]*EP[:vCHARGE][y,t])-(dfGen[!,:Self_Disch][y]*EP[:vS][y,t+hours_per_subperiod-1]))
	end

	@constraints(EP, begin

		# Max and min constraints on energy storage capacity built (as proportion to discharge power capacity)
		[y in STOR_ALL], EP[:eTotalCapEnergy][y] >= dfGen[!,:Min_Duration][y] * EP[:eTotalCap][y]
		[y in STOR_ALL], EP[:eTotalCapEnergy][y] <= dfGen[!,:Max_Duration][y] * EP[:eTotalCap][y]

		# Maximum energy stored must be less than energy capacity
		[y in STOR_ALL, t in 1:T], EP[:vS][y,t] <= EP[:eTotalCapEnergy][y]

		# energy stored for the next hour
		cSoCBalInterior[t in INTERIOR_SUBPERIODS, y in STOR_ALL], EP[:vS][y,t] ==
			EP[:vS][y,t-1]-(1/dfGen[!,:Eff_Down][y]*EP[:vP][y,t])+(dfGen[!,:Eff_Up][y]*EP[:vCHARGE][y,t])-(dfGen[!,:Self_Disch][y]*EP[:vS][y,t-1])
	end)

	# Storage discharge and charge power (and reserve contribution) related constraints:
	if Reserves == 1
		EP = storage_all_reserves(EP, inputs)
	else
		# Note: maximum charge rate is also constrained by maximum charge power capacity, but as this differs by storage type,
		# this constraint is set in functions below for each storage type

		# Maximum discharging rate must be less than power rating OR available stored energy in the prior period, whichever is less
		# wrapping from end of sample period to start of sample period for energy capacity constraint
		@constraints(EP, begin
			[y in STOR_ALL, t=1:T], EP[:vP][y,t] <= EP[:eTotalCap][y]
			[y in STOR_ALL, t in INTERIOR_SUBPERIODS], EP[:vP][y,t] <= EP[:vS][y,t-1]*dfGen[!,:Eff_Down][y]
			[y in STOR_ALL, t in START_SUBPERIODS], EP[:vP][y,t] <= EP[:vS][y,t+hours_per_subperiod-1]*dfGen[!,:Eff_Down][y]
		end)
	end
	# From co2 Policy module
	@expression(EP, eELOSSByZone[z=1:Z],
		sum(EP[:eELOSS][y] for y in intersect(inputs["STOR_ALL"], dfGen[dfGen[!,:Zone].==z,:R_ID]))
	)

	return EP
end

function storage_all_reserves(EP::Model, inputs::Dict)

	dfGen = inputs["dfGen"]
	T = inputs["T"]

	START_SUBPERIODS = inputs["START_SUBPERIODS"]
	INTERIOR_SUBPERIODS = inputs["INTERIOR_SUBPERIODS"]
	hours_per_subperiod = inputs["hours_per_subperiod"]

	STOR_ALL = inputs["STOR_ALL"]

	STOR_REG_RSV = intersect(STOR_ALL, inputs["REG"], inputs["RSV"]) # Set of storage resources with both REG and RSV reserves

	STOR_REG = intersect(STOR_ALL, inputs["REG"]) # Set of storage resources with REG reserves
	STOR_RSV = intersect(STOR_ALL, inputs["RSV"]) # Set of storage resources with RSV reserves

	STOR_NO_RES = setdiff(STOR_ALL, STOR_REG, STOR_RSV) # Set of storage resources with no reserves

	STOR_REG_ONLY = setdiff(STOR_REG, STOR_RSV) # Set of storage resources only with REG reserves
	STOR_RSV_ONLY = setdiff(STOR_RSV, STOR_REG) # Set of storage resources only with RSV reserves

	if !isempty(STOR_REG_RSV)
		# Storage units charging can charge faster to provide reserves down and charge slower to provide reserves up
		@constraints(EP, begin
			# Maximum storage contribution to reserves is a specified fraction of installed discharge power capacity
			[y in STOR_REG_RSV, t=1:T], EP[:vREG][y,t] <= dfGen[!,:Reg_Max][y]*EP[:eTotalCap][y]
			[y in STOR_REG_RSV, t=1:T], EP[:vRSV][y,t] <= dfGen[!,:Rsv_Max][y]*EP[:eTotalCap][y]

			# Actual contribution to regulation and reserves is sum of auxilary variables for portions contributed during charging and discharging
			[y in STOR_REG_RSV, t=1:T], EP[:vREG][y,t] == EP[:vREG_charge][y,t]+EP[:vREG_discharge][y,t]
			[y in STOR_REG_RSV, t=1:T], EP[:vRSV][y,t] == EP[:vRSV_charge][y,t]+EP[:vRSV_discharge][y,t]

			# Maximum charging rate plus contribution to reserves up must be greater than zero
			# Note: when charging, reducing charge rate is contributing to upwards reserve & regulation as it drops net demand
			[y in STOR_REG_RSV, t=1:T], EP[:vCHARGE][y,t]-EP[:vREG_charge][y,t]-EP[:vRSV_charge][y,t] >= 0

			# Maximum discharging rate and contribution to reserves down must be greater than zero
			# Note: when discharging, reducing discharge rate is contributing to downwards regulation as it drops net supply
			[y in STOR_REG_RSV, t=1:T], EP[:vP][y,t]-EP[:vREG_discharge][y,t] >= 0

			# Maximum charging rate plus contribution to regulation down must be less than available storage capacity
			## Made change to let the model run and not have key error issue for time -Sam (04/20/2021)
			[y in STOR_REG_RSV, t in START_SUBPERIODS], EP[:vCHARGE][y,t]+EP[:vREG_charge][y,t] <= EP[:eTotalCapEnergy][y]-EP[:vS][y,t+hours_per_subperiod-1]
			[y in STOR_REG_RSV, t in INTERIOR_SUBPERIODS], EP[:vCHARGE][y,t]+EP[:vREG_charge][y,t] <= EP[:eTotalCapEnergy][y]-EP[:vS][y,t-1]
			# Note: maximum charge rate is also constrained by maximum charge power capacity, but as this differs by storage type,
			# this constraint is set in functions below for each storage type

			# Maximum discharging rate and contribution to reserves up must be less than power rating OR available stored energy in prior period, whichever is less
			# wrapping from end of sample period to start of sample period for energy capacity constraint
			[y in STOR_REG_RSV, t=1:T], EP[:vP][y,t]+EP[:vREG_discharge][y,t]+EP[:vRSV_discharge][y,t] <= EP[:eTotalCap][y]
			[y in STOR_REG_RSV, t in INTERIOR_SUBPERIODS], EP[:vP][y,t]+EP[:vREG_discharge][y,t]+EP[:vRSV_discharge][y,t] <= EP[:vS][y,t-1]
			[y in STOR_REG_RSV, t in START_SUBPERIODS], EP[:vP][y,t]+EP[:vREG_discharge][y,t]+EP[:vRSV_discharge][y,t] <= EP[:vS][y,t+hours_per_subperiod-1]
		end)
	end
	if !isempty(STOR_REG_ONLY)
		# Storage units charging can charge faster to provide reserves down and charge slower to provide reserves up
		@constraints(EP, begin
			# Maximum storage contribution to reserves is a specified fraction of installed capacity
			[y in STOR_REG_ONLY, t=1:T], EP[:vREG][y,t] <= dfGen[!,:Reg_Max][y]*EP[:eTotalCap][y]

			# Actual contribution to regulation and reserves is sum of auxilary variables for portions contributed during charging and discharging
			[y in STOR_REG_ONLY, t=1:T], EP[:vREG][y,t] == EP[:vREG_charge][y,t]+EP[:vREG_discharge][y,t]

			# Maximum charging rate plus contribution to reserves up must be greater than zero
			# Note: when charging, reducing charge rate is contributing to upwards reserve & regulation as it drops net demand
			[y in STOR_REG_ONLY, t=1:T], EP[:vCHARGE][y,t]-EP[:vREG_charge][y,t] >= 0

			# Maximum discharging rate and contribution to reserves down must be greater than zero
			# Note: when discharging, reducing discharge rate is contributing to downwards regulation as it drops net supply
			[y in STOR_REG_ONLY, t=1:T], EP[:vP][y,t] - EP[:vREG_discharge][y,t] >= 0

			# Maximum charging rate plus contribution to regulation down must be less than available storage capacity
			[y in STOR_REG_ONLY, t in START_SUBPERIODS], EP[:vCHARGE][y,t]+EP[:vREG_charge][y,t] <= EP[:eTotalCapEnergy][y]-EP[:vS][y,t+hours_per_subperiod-1]
			[y in STOR_REG_ONLY, t in INTERIOR_SUBPERIODS], EP[:vCHARGE][y,t]+EP[:vREG_charge][y,t] <= EP[:eTotalCapEnergy][y]-EP[:vS][y,t-1]
			# Note: maximum charge rate is also constrained by maximum charge power capacity, but as this differs by storage type,
			# this constraint is set in functions below for each storage type

			# Maximum discharging rate and contribution to reserves up must be less than power rating OR available stored energy in prior period, whichever is less
			# wrapping from end of sample period to start of sample period for energy capacity constraint
			[y in STOR_REG_ONLY, t=1:T], EP[:vP][y,t] + EP[:vREG_discharge][y,t] <= EP[:eTotalCap][y]
			[y in STOR_REG_ONLY, t in INTERIOR_SUBPERIODS], EP[:vP][y,t]+EP[:vREG_discharge][y,t] <= EP[:vS][y,t-1]
			[y in STOR_REG_ONLY, t in START_SUBPERIODS], EP[:vP][y,t]+EP[:vREG_discharge][y,t]<= EP[:vS][y,t+hours_per_subperiod-1]
		end)
	end
	if !isempty(STOR_RSV_ONLY)
		# Storage units charging can charge faster to provide reserves down and charge slower to provide reserves up
		@constraints(EP, begin
			# Maximum storage contribution to reserves is a specified fraction of installed capacity
			[y in STOR_RSV_ONLY, t=1:T], EP[:vRSV][y,t] <= dfGen[!,:Rsv_Max][y]*EP[:eTotalCap][y]

			# Actual contribution to regulation and reserves is sum of auxilary variables for portions contributed during charging and discharging
			[y in STOR_RSV_ONLY, t=1:T], EP[:vRSV][y,t] == EP[:vRSV_charge][y,t]+EP[:vRSV_discharge][y,t]

			# Maximum charging rate plus contribution to reserves up must be greater than zero
			# Note: when charging, reducing charge rate is contributing to upwards reserve & regulation as it drops net demand
			[y in STOR_RSV_ONLY, t=1:T], EP[:vCHARGE][y,t]-EP[:vRSV_charge][y,t] >= 0

			# Note: maximum charge rate is also constrained by maximum charge power capacity, but as this differs by storage type,
			# this constraint is set in functions below for each storage type

			# Maximum discharging rate and contribution to reserves up must be less than power rating OR available stored energy in prior period, whichever is less
			# wrapping from end of sample period to start of sample period for energy capacity constraint
			[y in STOR_RSV_ONLY, t=1:T], EP[:vP][y,t]+EP[:vRSV_discharge][y,t] <= EP[:eTotalCap][y]
			[y in STOR_RSV_ONLY, t in INTERIOR_SUBPERIODS], EP[:vP][y,t]+EP[:vRSV_discharge][y,t] <= EP[:vS][y,t-1]
			[y in STOR_RSV_ONLY, t in START_SUBPERIODS], EP[:vP][y,t]+EP[:vRSV_discharge][y,t] <= EP[:vS][y,t+hours_per_subperiod-1]
		end)
	end
	if !isempty(STOR_NO_RES)
		# Maximum discharging rate must be less than power rating OR available stored energy in prior period, whichever is less
		# wrapping from end of sample period to start of sample period for energy capacity constraint
		@constraints(EP, begin
			[y in STOR_NO_RES, t=1:T], EP[:vP][y,t] <= EP[:eTotalCap][y]
			[y in STOR_NO_RES, t in INTERIOR_SUBPERIODS], EP[:vP][y,t] <= EP[:vS][y,t-1]
			[y in STOR_NO_RES, t in START_SUBPERIODS], EP[:vP][y,t] <= EP[:vS][y,t+hours_per_subperiod-1]
		end)
	end

	return EP
end

max_stars_repo_name: Betristor/LTESOM

max_issues_repo_name: Betristor/LTESOM

max_forks_repo_name: Betristor/LTESOM


text:
using HDF5

"""
Need to have method `get_data`
"""
abstract type TimeSeries{N,T} end

struct BoundaryTimeSeries{T} <: TimeSeries{2,T}
    ts    :: T
    field :: Symbol

    function BoundaryTimeSeries(foldername::String, field::Symbol)
        ts = OpenPMDTimeSeries(foldername, "boundary_")
        new{typeof(ts)}(ts, field)
    end
end

struct XiTimeSeries{T} <: TimeSeries{2,T}
    ts :: T

    function XiTimeSeries(foldername::String)
        ts = OpenPMDTimeSeries(foldername, "gauge_")
        new{typeof(ts)}(ts)
    end
end

struct BulkTimeSeries{T} <: TimeSeries{3,T}
    ts        :: T
    field     :: Symbol
    component :: Int

    function BulkTimeSeries(foldername::String, field::Symbol, c::Int)
        ts = OpenPMDTimeSeries(foldername, "bulk_")
        new{typeof(ts)}(ts, field, c)
    end
end

struct ConstrainedTimeSeries{T} <: TimeSeries{3,T}
    ts        :: T
    field     :: Symbol
    component :: Int

    function ConstrainedTimeSeries(foldername::String, field::Symbol, c::Int)
        ts = OpenPMDTimeSeries(foldername, "constrained_")
        new{typeof(ts)}(ts, field, c)
    end
end

struct VEVTimeSeries{T} <: TimeSeries{2,T}
    ts  :: T
    vev :: Symbol

    function VEVTimeSeries(foldername::String, vev::Symbol)
        ts = OpenPMDTimeSeries(foldername, "boundary_")
        new{typeof(ts)}(ts, vev)
    end
end

struct AHTimeSeries{T} <: TimeSeries{2,T}
    ts            :: T
    ts_diag       :: T
    property :: Symbol

    function AHTimeSeries(foldername::String, property::Symbol)
        ts       = OpenPMDTimeSeries(foldername, "constrained_")
        ts_diag  = OpenPMDTimeSeries(foldername, "diagnostics_")
        new{OpenPMDTimeSeries}(ts, ts_diag, property)
    end
end

function get_data(ff::BoundaryTimeSeries, it::Int)
    f, chart = get_field(ff.ts, it=it, field=String(ff.field))
    _, x, y  = chart[:]
    f[1,:,:], [x, y]
end

function get_data(xi::XiTimeSeries, it::Int)
    f, chart = get_field(xi.ts, it=it, field="xi")
    _, x, y  = chart[:]
    f[1,:,:], [x, y]
end

function get_data(ff::BulkTimeSeries, it::Int)
    field = "$(ff.field) c=$(ff.component)"
    f, chart = get_field(ff.ts, it=it, field=field)
    u, x, y = chart[:]
    f, [u, x, y]
end

function get_data(ff::ConstrainedTimeSeries, it::Int)
    field = "$(ff.field) c=$(ff.component)"
    f, chart = get_field(ff.ts, it=it, field=field)
    u, x, y = chart[:]
    f, [u, x, y]
end

function get_data(ff::VEVTimeSeries, it::Int)
    if ff.vev == :energy
        f, chart = get_energy(ff.ts, it=it)
    elseif ff.vev == :Jx
        f, chart = get_Jx(ff.ts, it=it)
    elseif ff.vev == :Jy
        f, chart = get_Jy(ff.ts, it=it)
    elseif ff.vev == :px
        f, chart = get_px(ff.ts, it=it)
    elseif ff.vev == :py
        f, chart = get_py(ff.ts, it=it)
    elseif ff.vev == :pz
        f, chart = get_pz(ff.ts, it=it)
    elseif ff.vev == :Jx
        f, chart = get_Jx(ff.ts, it=it)
    elseif ff.vev == :pxy
        f, chart = get_pxy(ff.ts, it=it)
    elseif ff.vev == :Ophi
        f, chart = get_Ophi(ff.ts, it=it)
    else
        error("Unknown VEV")
    end
    _, x, y  = chart[:]
    f[1,:,:], [x, y]
end

function get_data(ff::AHTimeSeries, it::Int)
    if ff.property == :temperature
        f, chart = get_temperature(ff.ts, ff.ts_diag, it=it)
    elseif ff.property == :entropy
        f, chart = get_entropy(ff.ts, ff.ts_diag, it=it)
    else
        error("Unknown property")
    end
    _, x, y = chart[:]
    f[1,:,:], [x,y]
end

function Base.size(ff::TimeSeries)
    Nt  = length(ff.ts.iterations)
    it0 = ff.ts.iterations[1]
    f0, = get_data(ff, it0)
    size_ = size(f0)
    Nt, size_...
end

@inline Base.firstindex(ff::TimeSeries, d::Int) = 1
@inline Base.lastindex(ff::TimeSeries, d::Int) = size(ff)[d]

function Base.getindex(ff::TimeSeries, a::Int, idx::Vararg)
    it = ff.ts.iterations[a]
    f, = get_data(ff, it)
    f[idx...]
end

function Base.getindex(ff::TimeSeries, aa, idx::Vararg)
    it  = ff.ts.iterations[aa[1]]
    f0, = get_data(ff, it)

    Na    = length(aa)
    size_ = size(f0[idx...])
    f     = zeros(Na, size_...)

    slicer = [Colon() for _ in 1:ndims(f0)]
    for (i,a) in enumerate(aa)
        it  = ff.ts.iterations[a]
        f0, = get_data(ff, it)
        f[i,slicer...] .= f0[idx...]
    end
    f
end

function Base.getindex(ff::TimeSeries, ::Colon, idx::Vararg)
    Na  = length(ff.ts.iterations)
    getindex(ff, 1:Na, idx...)
end

function Jecco.get_coords(ff::TimeSeries{N}, a::Int, idx::Vararg) where{N}
    it = ff.ts.iterations[a]
    f, coords = get_data(ff, it)
    t = ff.ts.current_t
    coord = [coords[a][idx[a]] for a in 1:N]
    t, coord...
end

function Jecco.get_coords(ff::TimeSeries{N}, aa, idx::Vararg) where{N}
    it = ff.ts.iterations[aa[1]]
    f0, coords = get_data(ff, it)
    t0 = ff.ts.current_t
    coord = [coords[a][idx[a]] for a in 1:N]

    Na = length(aa)
    t  = Vector{typeof(t0)}(undef, Na)
    for (i,a) in enumerate(aa)
        it   = ff.ts.iterations[a]
        get_data(ff, it)
        t0   = ff.ts.current_t
        t[i] = t0
    end
    t, coord...
end

function Jecco.get_coords(ff::TimeSeries{N}, ::Colon, idx::Vararg) where{N}
    Na = length(ff.ts.iterations)
    get_coords(ff, 1:Na, idx...)
end

function convert_to_mathematica(dirname::String; outfile::String="data_mathematica.h5")
    ts = OpenPMDTimeSeries(dirname, prefix="boundary_")

    iterations = ts.iterations

    Nt = length(iterations)
    t  = zeros(Nt)

    it = 0
    en0, chart = get_energy(ts, it=it)

    _, x, y = chart[:]
    Nx = length(x)
    Ny = length(y)

    en   = zeros(Nt,Nx,Ny)
    Jx   = zeros(Nt,Nx,Ny)
    Jy   = zeros(Nt,Nx,Ny)
    px   = zeros(Nt,Nx,Ny)
    py   = zeros(Nt,Nx,Ny)
    pz   = zeros(Nt,Nx,Ny)
    pxy  = zeros(Nt,Nx,Ny)
    Ophi = zeros(Nt,Nx,Ny)

    for (idx,it) in enumerate(iterations)
        en[idx,:,:]   .= get_energy(ts, it=it)[1][1,:,:]
        Jx[idx,:,:]   .= get_Jx(ts, it=it)[1][1,:,:]
        Jy[idx,:,:]   .= get_Jy(ts, it=it)[1][1,:,:]
        px[idx,:,:]   .= get_px(ts, it=it)[1][1,:,:]
        py[idx,:,:]   .= get_py(ts, it=it)[1][1,:,:]
        pz[idx,:,:]   .= get_pz(ts, it=it)[1][1,:,:]
        pxy[idx,:,:]  .= get_pxy(ts, it=it)[1][1,:,:]
        Ophi[idx,:,:] .= get_Ophi(ts, it=it)[1][1,:,:]

        t[idx] = ts.current_t
    end

    # store in an array suitable for Mathematica
    T_m = zeros(11, Nt*Nx*Ny)

    n = 1
    for i in 1:Nt
        for j in 1:Nx
            for k in 1:Ny
                T_m[:,n] = [t[i] x[j] y[k] en[i,j,k] Jx[i,j,k] Jy[i,j,k] px[i,j,k] pxy[i,j,k] py[i,j,k] pz[i,j,k] Ophi[i,j,k]]
                n += 1
            end
        end
    end

    output   = abspath(dirname, outfile)
    fid      = h5open(output, "w")
    group_st = g_create(fid, "data")
    Jecco.write_dataset(group_st, "VEVs", T_m)
    close(fid)
end

function convert_to_mathematica_local(dirname::String; outfile::String="local_data_mathematica.h5")#,
                                #phi0, oophiM2)

    ts = OpenPMDTimeSeries(dirname, prefix="boundary_")

    iterations = ts.iterations

    Nt = length(iterations)
    t  = zeros(Nt)

    it = 0
    en0, chart = get_energy(ts, it=it)

    _, x, y = chart[:]
    Nx = length(x)
    Ny = length(y)

    ut        = zeros(Nt,Nx,Ny)
    ux        = zeros(Nt,Nx,Ny)
    uy        = zeros(Nt,Nx,Ny)
    en_local  = zeros(Nt,Nx,Ny)
    p1_local  = zeros(Nt,Nx,Ny)
    p2_local  = zeros(Nt,Nx,Ny)

    for (idx,it) in enumerate(iterations)
        en   = get_energy(ts, it=it)[1][1,:,:]
        Jx   = get_Jx(ts, it=it)[1][1,:,:]
        Jy   = get_Jy(ts, it=it)[1][1,:,:]
        px   = get_px(ts, it=it)[1][1,:,:]
        py   = get_py(ts, it=it)[1][1,:,:]
        pxy  = get_pxy(ts, it=it)[1][1,:,:]

        t[idx] = ts.current_t

        for j in 1:Ny
            for i in 1:Nx
                ut[idx,i,j],ux[idx,i,j],uy[idx,i,j],en_local[idx,i,j],p1_local[idx,i,j],p2_local[idx,i,j] = AdS5_3_1.compute_local_VEVs(en[i,j],Jx[i,j],Jy[i,j],px[i,j],py[i,j],pxy[i,j])
            end
        end
    end

    # store in an array suitable for Mathematica
    T_m = zeros(9, Nt*Nx*Ny)

    n = 1
    for i in 1:Nt
        for j in 1:Nx
            for k in 1:Ny
                T_m[:,n] = [t[i] x[j] y[k] ut[i,j,k] ux[i,j,k] uy[i,j,k] en_local[i,j,k] p1_local[i,j,k] p2_local[i,j,k]]
                n += 1
            end
        end
    end

    output   = abspath(dirname, outfile)
    fid      = h5open(output, "w")
    group_st = g_create(fid, "data")
    Jecco.write_dataset(group_st, "Local VEVs", T_m)
    close(fid)

end


text:
using BinaryBuilder

name = "utf8proc"
version = v"2.6.1"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/JuliaStrings/utf8proc.git",
              "3203baa7374d67132384e2830b2183c92351bffc"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/utf8proc*

if [[ "${target}" == *-mingw* ]]; then
    make -j${nproc} libutf8proc.a
    mkdir -p ${libdir} ${prefix}/lib ${includedir}
    cp utf8proc.h ${includedir}
    ar x libutf8proc.a
    cc -shared -o "${libdir}/libutf8proc.${dlext}" *.o
    cp libutf8proc.a ${prefix}/lib
else
    make -j${nproc}
    make install prefix=${prefix} libdir=${libdir}
fi
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(; experimental=true)

# The products that we will ensure are always built
products = [
    LibraryProduct("libutf8proc", :libutf8proc),
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil