text:
# setup
# using Pkg
# Pkg.pkg"add Makie#master AbstractPlotting#master GLMakie#master VideoIO#master Observables DataStructures"

using Makie, VideoIO, Dates, Observables

# include(joinpath(@__DIR__, "testvideo.jl"))

# testvideo = joinpath(tempdir(), "test.mp4")
#
# # create the test video if it doesn't exist
# if !isfile(testvideo)
#     createtestvideo(testvideo)
# end

testvideo = joinpath(tempdir(), "testvideo.mp4")
# f = VideoIO.openvideo(avf)
_correctimg(img) = rotr90(img)

# open the video
f = VideoIO.openvideo(testvideo)

duration = VideoIO.get_duration(f.avin.io)

_img = _correctimg(read(f))
seekstart(f)

img = Observable(_img)
current = Observable(0.0)
on(current) do t
    try
        seek(f, t)
    catch ex
        if !isa(ex, EOFError)
            throw(ex)
        end
    end
end
correctcurrent = lift(img) do _
    gettime(f)
end
timestamp = lift(correctcurrent) do cc
    Time(0,0,0,1) + Millisecond(round(Int, 1000cc))
end
# all the stuff needed for display
rsc() = Scene(;camera = campixel!, raw = true, backgroundcolor = :black)
pixelaspectratio = VideoIO.aspect_ratio(f)
h = f.height
w = round(typeof(h), f.width * pixelaspectratio)
AbstractPlotting.set_theme!(backgroundcolor = :black)
scene = Scene(resolution = (w, h))
Makie.image!(scene, 1:w, 1:h, img, show_axis = false, scale_plot = false)
# Makie.rotate!(scene, -0.5π)
update_limits!(scene)
update_cam!(scene)
Makie.update!(scene)

# the slider, updating currentframe and getting its string from timestamps
slidersteps = range(0, duration, length = 100)
slider_h = slider!(rsc(), slidersteps, start = 0.0, valueprinter = _ -> "lolno", textcolor = :white)

lift(slider_h[end][:value]) do t
    current[] = t
    # if !eof(f)
    img[] = _correctimg(read(f))
    # end
    nothing
end

# this is for a separate display of the time stamp, due to the fact that we can not calculate with the time stamp is from the frame number (which is what the slider has)
timestamp_h = text!(rsc(), lift(string, timestamp), color = :white)

# I have to have a forward button becuase apparently the stuff in the lift gets evaluated, so to have the movie start in frame #1 I need to go forward once and the backwards...
fwdbutton = button!(rsc(), ">", textcolor = :white)
lift(fwdbutton[end][:clicks]) do _
    # slider_h[end].value[] = slider_h[end].value[] + steplen
    if !eof(f)
        img[] = _correctimg(read(f))
    end
    #=if correctcurrent[] - slider_h[end][:value][] > step(slidersteps)/2
        @info "do"
        slider_h[end][:value][] += step(slidersteps) # this doesn't work
    end=#
    nothing
end
#=bckbutton = button!(rsc(), "<", textcolor = :white)
lift(bckbutton[end][:clicks]) do _
    t2 = correctcurrent[]
    seek(f, max(t2 - 1, 0.0))
    # read(f)
    # t0 = 0.0
    # t1 = 0.0
    # while t1 < t2
    #     t0 = gettime(f)
    #     read(f) # TODO FIXME inefficient
    #     t1 = gettime(f)
    # end
    current[] = max(t2 - 1, 0.0)
    img[] = read(f) |> _correctimg
end=#

# done!

sc = Scene(resolution = (Int64(w), h+Int64(100)))
sc.center=false

seekstart(f)
img[] = _correctimg(read(f))

hbox(vbox(slider_h, fwdbutton, timestamp_h), scene; parent=sc)

# setup keyboard controls
kb = on(sc.events.keyboardbuttons) do kb
    if ispressed(sc, Keyboard.right)
        fwdbutton[end][:clicks][] += 1
    elseif ispressed(sc,Keyboard.left)
        bckbutton[end][:clicks][] += 1
    end
end

lastmpos = Node(Point2f0(0e0, 0e0))

to_screen(scene, mpos) = Point2f0(mpos) .- Point2f0(minimum(pixelarea(scene)[]))
mouseposition(scene) = to_world(scene, to_screen(scene, events(scene).mouseposition[]))

mb = on(scene.events.mousebuttons) do mb
    if ispressed(scene, Mouse.left)
        lastmpos[] = mouseposition(scene)
    end
end

return sc, lastmpos

Observables.off(scene.events.mousebuttons, mb)


text:
"""
    make_point_variable_ref(ivref::InfiniteVariableRef)::PointVariableRef

Internal method for measure expansion that returns a `PointVariableRef` with
a unique index. This is useful for writing extensions to [`expand_measure`](@ref).
"""
function make_point_variable_ref(ivref::InfiniteVariableRef)::PointVariableRef
    inf_model = JuMP.owner_model(ivref)
    index = inf_model.next_var_index += 1
    return PointVariableRef(inf_model, index)
end

## Helper function for making place holder infinite variables
"""
    make_reduced_variable_ref(ivref::InfiniteVariableRef, removed_index::Int,
                              support::Union{Number,
                              JuMPC.SparseAxisArray{<:Number}}
                              )::ReducedInfiniteVariableRef

Internal method for measure expansion that returns a `ReducedInfiniteVariableRef`
and that is connected with `ivref.model`. This method is intended for intial
reduction of an infinite variable. Here `removed_index` is the index of the
reduced parameter in `parameter_refs(ivref)` and `support` is the value that will
support that particular infinite parameter. This is useful for writing
extensions to [`expand_measure`](@ref).
"""
function make_reduced_variable_ref(ivref::InfiniteVariableRef, removed_index::Int,
                                support::Union{Number,
                                JuMPC.SparseAxisArray{<:Number}}
                                )::ReducedInfiniteVariableRef
    inf_model = JuMP.owner_model(ivref)
    index = inf_model.next_var_index += 1
    inf_model.reduced_info[index] = ReducedInfiniteInfo(ivref,
                                                 Dict(removed_index => support))
    if haskey(inf_model.infinite_to_reduced, JuMP.index(ivref))
        push!(inf_model.infinite_to_reduced[JuMP.index(ivref)], index)
    else
        inf_model.infinite_to_reduced[JuMP.index(ivref)] = [index]
    end
    return ReducedInfiniteVariableRef(inf_model, index)
end

"""
    make_reduced_variable_ref(ivref::InfiniteVariableRef,
                              supports::Dict)::ReducedInfiniteVariableRef

Internal method for measure expansion that returns a `ReducedInfiniteVariableRef`
and that is connected with `ivref.model`. This method is intended for further
reduction of an infinite variable based on an existing reduced infinite variable.
Here `supports` is a dictionary mapping the indices of parameters in
`parameter_refs(ivref)` to their supports values. This is useful for writing
extensions to [`expand_measure`](@ref).
"""
function make_reduced_variable_ref(ivref::InfiniteVariableRef,
                                    supports::Dict)::ReducedInfiniteVariableRef
    inf_model = JuMP.owner_model(ivref)
    index = inf_model.next_var_index += 1
    inf_model.reduced_info[index] = ReducedInfiniteInfo(ivref, copy(supports))
    if haskey(inf_model.infinite_to_reduced, JuMP.index(ivref))
        push!(inf_model.infinite_to_reduced[JuMP.index(ivref)], index)
    else
        inf_model.infinite_to_reduced[JuMP.index(ivref)] = [index]
    end
    return ReducedInfiniteVariableRef(inf_model, index)
end

## Make helper functions for extracting evaluated parameter values
# DiscreteMeasureData
function _get_param_value_list(pref::ParameterRef, data::DiscreteMeasureData)
    return supports(data)
end

# MultiDiscreteMeasureData
function _get_param_value_list(pref::ParameterRef,
                               data::MultiDiscreteMeasureData)
    key = first(filter(p -> p[2] == pref, parameter_refs(data).data))[1]
    return [supports(data)[i][key] for i in eachindex(supports(data))]
end

"""
    expand_measure(expr, data::AbstractMeasureData,
                   write_model::JuMP.AbstractModel,
                   point_mapper::Function)::JuMP.AbstractJuMPScalar

Return the finite reformulation of a measure containing a variable/parameter
expression `expr` with measure data `data`. Here `write_model` is the target
model where this expanded expression will be used. Thus, any variables that need
to be created will be added to `write_model`. The function `point_mapper` should
be called to make and map point variables via the syntax:
`point_mapper(write_model, pvref, ivref, support)`. Here `ivref` denotes the
`InfiniteVariableRef` that the point variable will be based on at `support` where
`support` is a `Tuple` of values that matches the format of `parameter_refs(ivref)`.
`pvref` is the pre-made `ParameterRef` that should be created via
[`make_point_variable_ref`](@ref). [`make_reduced_variable_ref`](@ref) may also
be helpful. Note this is intended as an internal function,
but will need to be extended for unsupported `expr` types and for user-defined
measure data types.
"""
function expand_measure end

# InfiniteVariableRef
function expand_measure(ivref::InfiniteVariableRef,
                         data::Union{DiscreteMeasureData,
                                     MultiDiscreteMeasureData},
                         trans_model::JuMP.AbstractModel,
                         point_mapper::Function)::JuMP.AbstractJuMPScalar
    # figure out the parameter groups
    group = group_id(first(data.parameter_ref))
    groups = _group.(parameter_refs(ivref))
    # prepare return AffExpr and get necessary information
    aff = zero(JuMP.GenericAffExpr{Float64, GeneralVariableRef})
    # treat variable as constant if doesn't have measure parameter
    if !(group in groups)
        for i in eachindex(data.supports)
            JuMP.add_to_expression!(aff, data.coefficients[i] *
                                    data.weight_function(data.supports[i]),
                                    ivref)
        end
    # convert variable into point variables if its only parameter is the
    # measure parameter
    elseif length(parameter_refs(ivref)) == 1
        for i in eachindex(data.supports)
            pvref = make_point_variable_ref(ivref)
            point_mapper(trans_model, pvref, ivref, (data.supports[i],))
            JuMP.add_to_expression!(aff, data.coefficients[i] *
                                    data.weight_function(data.supports[i]),
                                    pvref)
        end
    # make reduced variables if the variable contains other parameters
    else
        tuple_loc = findfirst(isequal(group), groups)
        for i in eachindex(data.supports)
            rvref = make_reduced_variable_ref(ivref, tuple_loc, data.supports[i])
            JuMP.add_to_expression!(aff, data.coefficients[i] *
                                    data.weight_function(data.supports[i]),
                                    rvref)
        end
    end
    return aff
end

# ReducedInfiniteVariableRef
function expand_measure(rvref::ReducedInfiniteVariableRef,
                         data::Union{DiscreteMeasureData,
                                     MultiDiscreteMeasureData},
                         trans_model::JuMP.AbstractModel,
                         point_mapper::Function)::JuMP.AbstractJuMPScalar
    # figure out the parameters used by the reduced infinite variable
    orig_prefs = parameter_refs(infinite_variable_ref(rvref))
    # figure out the parameter groups
    group = group_id(first(data.parameter_ref))
    groups = _group.(parameter_refs(rvref))
    # prepare return AffExpr and get necessary information
    aff = zero(JuMP.GenericAffExpr{Float64, GeneralVariableRef})
    # treat variable as constant if doesn't have measure parameter
    if !(group in groups)
        for i in eachindex(data.supports)
            JuMP.add_to_expression!(aff, data.coefficients[i] *
                                    data.weight_function(data.supports[i]),
                                    rvref)
        end
    # convert variable into point variables if its only parameter is the
    # measure parameter
    elseif length(parameter_refs(rvref)) == 1
        tuple_loc = findfirst(isequal(group), _group.(orig_prefs))
        for i in eachindex(data.supports)
            pvref = make_point_variable_ref(infinite_variable_ref(rvref))
            _reduced_info(rvref).eval_supports[tuple_loc] = data.supports[i]
            support = Tuple(eval_supports(rvref)[j] for j in eachindex(eval_supports(rvref)))
            point_mapper(trans_model, pvref, infinite_variable_ref(rvref), support)
            JuMP.add_to_expression!(aff, data.coefficients[i] *
                                    data.weight_function(data.supports[i]),
                                    pvref)
        end
    # make reduced variables if the variable contains other parameters
    else
        tuple_loc = findfirst(isequal(group), _group.(orig_prefs))
        for i in eachindex(data.supports)
            new_rvref = make_reduced_variable_ref(infinite_variable_ref(rvref),
                                               eval_supports(rvref))
            _reduced_info(new_rvref).eval_supports[tuple_loc] = data.supports[i]
            JuMP.add_to_expression!(aff, data.coefficients[i] *
                                    data.weight_function(data.supports[i]),
                                    new_rvref)
        end
    end
    return aff
end

# FiniteVariableRef
function expand_measure(vref::FiniteVariableRef,
                         data::Union{DiscreteMeasureData,
                                     MultiDiscreteMeasureData},
                         trans_model::JuMP.AbstractModel,
                         point_mapper::Function)::JuMP.AbstractJuMPScalar
    aff = zero(JuMP.GenericAffExpr{Float64, GeneralVariableRef})
    # treat the variable as a constant
    for i in eachindex(data.supports)
        JuMP.add_to_expression!(aff, data.coefficients[i] *
                                data.weight_function(data.supports[i]),
                                vref)
    end
    return aff
end

# ParameterRef with scalar data
function expand_measure(pref::ParameterRef,
                         data::DiscreteMeasureData,
                         trans_model::JuMP.AbstractModel,
                         point_mapper::Function)::JuMP.AbstractJuMPScalar
    aff = zero(JuMP.GenericAffExpr{Float64, GeneralVariableRef})
    # replace the parameter with its value if it is the measure parameter
    if data.parameter_ref == pref
        for i in eachindex(data.supports)
            JuMP.add_to_expression!(aff, data.coefficients[i] *
                                    data.weight_function(data.supports[i]) *
                                    data.supports[i])
        end
    # treat the parameter as a constant otherwise
    else
        for i in eachindex(data.supports)
            JuMP.add_to_expression!(aff, data.coefficients[i] *
                                    data.weight_function(data.supports[i]),
                                    pref)
        end
    end
    return aff
end

# ParameterRef with vector data
function expand_measure(pref::ParameterRef,
                         data::MultiDiscreteMeasureData,
                         trans_model::JuMP.AbstractModel,
                         point_mapper::Function)::JuMP.AbstractJuMPScalar
    aff = zero(JuMP.GenericAffExpr{Float64, GeneralVariableRef})
    # determine if pref is part of the measure parameters
    pref_dict = filter(p -> p[2] == pref, data.parameter_ref.data)
    # replace the parameter with its value if it is the measure parameter
    if length(pref_dict) != 0
        for i in eachindex(data.supports)
            JuMP.add_to_expression!(aff, data.coefficients[i] *
                                    data.weight_function(data.supports[i]) *
                                    data.supports[i][collect(keys(pref_dict))[1]])
        end
    # treat the parameter as a constant otherwise
    else
        for i in eachindex(data.supports)
            JuMP.add_to_expression!(aff, data.coefficients[i] *
                                    data.weight_function(data.supports[i]),
                                    pref)
        end
    end
    return aff
end

# GenericAffExpr
function expand_measure(expr::JuMP.GenericAffExpr,
                         data::Union{DiscreteMeasureData,
                                     MultiDiscreteMeasureData},
                         trans_model::JuMP.AbstractModel,
                         point_mapper::Function)::JuMP.AbstractJuMPScalar
    # need to use a quadratic expression in case contains measures with
    # quadratic expressions
    quad = zero(JuMP.GenericQuadExpr{Float64, GeneralVariableRef})
    # expand each variable independently and add all together
    for (var, coef) in expr.terms
        JuMP.add_to_expression!(quad, coef, expand_measure(var, data,
                                                            trans_model,
                                                            point_mapper))
    end
    # expand over the constant
    if expr.constant != 0
        for i in eachindex(data.supports)
            JuMP.add_to_expression!(quad, data.coefficients[i] *
                                    data.weight_function(data.supports[i]) *
                                    expr.constant)
        end
    end
    # return affexpr if appropriate
    if length(quad.terms) == 0
        return quad.aff
    else
        return quad
    end
end

# GenericQuadExpr
function expand_measure(expr::JuMP.GenericQuadExpr,
                         data::Union{DiscreteMeasureData,
                                     MultiDiscreteMeasureData},
                         trans_model::JuMP.AbstractModel,
                         point_mapper::Function)::JuMP.AbstractJuMPScalar
    quad = zero(JuMP.GenericQuadExpr{Float64, GeneralVariableRef})
    # convert the GenericAffExpr
    quad.aff = expand_measure(expr.aff, data, trans_model, point_mapper)
    for (pair, coef) in expr.terms
        # expand on both variables
        expr_a = expand_measure(pair.a, data, trans_model, point_mapper)
        expr_b = expand_measure(pair.b, data, trans_model, point_mapper)
        vars_a = collect(keys(expr_a.terms))
        vars_b = collect(keys(expr_b.terms))
        alt_terms = false
        # check for case that a variable was a parameter converted to a number
        if length(vars_a) == 0
            vars_a = _get_param_value_list(pair.a, data)
            terms = [data.coefficients[i] * data.weight_function(data.supports[i]) for i in eachindex(data.supports)]
            alt_terms = true
        end
        if length(vars_b) == 0
            vars_b = _get_param_value_list(pair.b, data)
            terms = [data.coefficients[i] * data.weight_function(data.supports[i]) for i in eachindex(data.supports)]
            alt_terms = true
        end
        # combine both variable expressions using the coefficients from one
        # of them
        if length(vars_a) == length(vars_b)
            # are same length therefore have same coefficients
            for i in eachindex(vars_a)
                if alt_terms
                    JuMP.add_to_expression!(quad, coef * terms[i], vars_a[i],
                                            vars_b[i])
                else
                    JuMP.add_to_expression!(quad, coef * expr_a.terms[vars_a[i]],
                                            vars_a[i], vars_b[i])
                end
            end
        elseif length(vars_a) == 1
            # var_a was effectively a constant and var_b was't
            for i in eachindex(vars_b)
                if alt_terms
                    JuMP.add_to_expression!(quad, coef * terms[i], vars_a[1],
                                            vars_b[i])
                else
                    JuMP.add_to_expression!(quad, coef * expr_b.terms[vars_b[i]],
                                            vars_a[1], vars_b[i])
                end
            end
        else
            # var_b was effectively a constant and var_a was't
            for i in eachindex(vars_a)
                if alt_terms
                    JuMP.add_to_expression!(quad, coef * terms[i], vars_a[i],
                                            vars_b[1])
                else
                    JuMP.add_to_expression!(quad, coef * expr_a.terms[vars_a[i]],
                                            vars_a[i], vars_b[1])
                end
            end
        end
    end
    return quad
end

# MeasureRef
function expand_measure(mref::MeasureRef,
                         data::Union{DiscreteMeasureData,
                                     MultiDiscreteMeasureData},
                         trans_model::JuMP.AbstractModel,
                         point_mapper::Function)::JuMP.AbstractJuMPScalar
    # determine function and data of the inner measure
    deeper_func = measure_function(mref)
    deeper_data = measure_data(mref)
    # expand the inner measure (note this is recursive for nested measures)
    new_func = expand_measure(deeper_func, deeper_data, trans_model,
                               point_mapper)
    # expand current level with the inner measure now expanded
    return expand_measure(new_func, data, trans_model, point_mapper)
end

# Catch all method for undefined behavior
function expand_measure(expr, data::AbstractMeasureData,
                         trans_model::JuMP.AbstractModel,
                         point_mapper::Function)
    expr_type = typeof(expr)
    data_type = typeof(data)
    error("Undefined behavior to expand expression of type `$expr_type` with " *
          "measure data `$data_type`. If this functionality is needed consider " *
          "extending `expand_measure`.")
end

# Map temp point variable references to actual variables by adding a variable to
# the infinite model. This is used by expand(mref) to expand measures within an
# infinite model
function _add_mapped_point_variable(model::InfiniteModel,
                                    pvref::PointVariableRef,
                                    ivref::InfiniteVariableRef, support::Tuple)
    # build variable
    var = JuMP.build_variable(error, _variable_info(ivref), Point,
                              infinite_variable_ref = ivref,
                              parameter_values = support)
    # add the variable completely (note the reference is already made)
    desired_index = JuMP.index(pvref)
    curr_index = model.next_var_index
    model.next_var_index = desired_index - 1
    JuMP.add_variable(model, var)
    model.next_var_index = curr_index
    return
end

"""
    expand(mref::MeasureRef)::JuMP.AbstractJuMPScalar

Return a JuMP scalar function containing the explicit expansion of the measure
`mref`. This expansion is done according to the measure data. Note that
variables are added to the model as necessary to accomodate the expansion (i.e.,
point variables and reduced infinite variables are made as needed). Errors if
expansion is undefined for the measure data and/or the measure expression. If
desired this can be used in combination with [`measure`](@ref) to expand measures
on the fly.

This is useful for extensions that employ a custom optimizer model since it
can be used evaluate measures before expressions are translated to the new model.
This method can also be extended to handle custom measure data types by extending
```expand_measure``` which should be of the form:
```julia
expand_measure(expr::AbstractJuMPScalar, data::AbstractMeasureData,
                model::JuMP.AbstractModel, point_mapper::Function)
```
See [`expand_measure`](@ref) for more information.

**Example**
```julia-repl
julia> tdata = DiscreteMeasureData(t, [0.5, 0.5], [0, 1])

julia> expr = expand(measure(g + z + T - h - 2, tdata))
0.5 g(0) + 0.5 g(1) + z + 0.5 T(0, x) + 0.5 T(1, x) - h(x) - 2
```
"""
function expand(mref::MeasureRef)::JuMP.AbstractJuMPScalar
    return expand_measure(measure_function(mref), measure_data(mref),
                           JuMP.owner_model(mref), _add_mapped_point_variable)
end

## Helper functions for expanding the measure references in expressions
# MeasureRef
function _expand_measures(mref::MeasureRef,
                          expand_model::JuMP.AbstractModel,
                          point_mapper::Function)::JuMP.AbstractJuMPScalar
    return expand_measure(measure_function(mref), measure_data(mref),
                           expand_model, point_mapper)
end

# GenericAffExpr
function _expand_measures(expr::JuMP.GenericAffExpr{C, <:GeneralVariableRef},
                          expand_model::JuMP.AbstractModel,
                          point_mapper::Function)::JuMP.AbstractJuMPScalar where {C}
    # use a QuadExpr in case measures contain quadratic espressions
    quad = zero(JuMP.GenericQuadExpr{C, GeneralVariableRef})
    quad.aff.constant = expr.constant
    # add the variables to the expr, converting measures into expanded exprs
    for (var, coef) in expr.terms
        if isa(var, MeasureRef)
            JuMP.add_to_expression!(quad, coef, expand_measure(measure_function(var),
                                                                measure_data(var),
                                                                expand_model,
                                                                point_mapper))
        else
            JuMP.add_to_expression!(quad, coef, var)
        end
    end
    # return a AffExpr if there are no quadratic terms
    if length(quad.terms) == 0
        return quad.aff
    else
        return quad
    end
end

# GenericQuadExpr
function _expand_measures(expr::JuMP.GenericQuadExpr{C, <:GeneralVariableRef},
                          expand_model::JuMP.AbstractModel,
                          point_mapper::Function)::JuMP.GenericQuadExpr where {C}
    quad = zero(JuMP.GenericQuadExpr{C, GeneralVariableRef})
    quad.aff = _expand_measures(expr.aff, expand_model, point_mapper)
    # add the quadratic terms to the expr, converting measures into expanded exprs
    # note that this will error if the expanded forms are not linear or quadratic
    for (pair, coef) in expr.terms
        var_a = pair.a
        var_b = pair.b
        if isa(var_a, MeasureRef)
            var_a = expand_measure(measure_function(var_a), measure_data(var_a),
                                    expand_model, point_mapper)
        end
        if isa(var_b, MeasureRef)
            var_b = expand_measure(measure_function(var_b), measure_data(var_b),
                                    expand_model, point_mapper)
        end
        JuMP.add_to_expression!(quad, convert(JuMP.GenericQuadExpr{C,
                                GeneralVariableRef}, coef * var_a * var_b))
    end
    return quad
end

"""
    expand_all_measures!(model::InfiniteModel)

Expand all of the measures used in the objective and/or constraints of `model`.
The objective and constraints are updated accordingly. Note that
variables are added to the model as necessary to accomodate the expansion (i.e.,
point variables and reduced infinite variables are made as needed). Errors if
expansion is undefined for the measure data and/or the measure expression. Also
errors if the expanded objective function is not finite.

This is useful for extensions that employ a custom optimizer model since it
can be used evaluate measures before `model` is translated into the new model.
This method can also be extended to handle custom measure data types by extending
```expand_measure``` which should be of the form:
```julia
expand_measure(expr::AbstractJuMPScalar, data::AbstractMeasureData,
                model::JuMP.AbstractModel, point_mapper::Function)
```
See [`expand_measure`](@ref) for more information.

**Example**
```julia-repl
julia> print(model)
Min measure(g(t)*t) + z
Subject to
 T(t, xi) >= 0.0
 z >= 0.0
 g(t) + z >= 42.0
 measure(T(t, xi)) >= 0.0, for all xi in [-1, 1]
 t in [0, 6]
 xi in Normal(μ=0.0, σ=1.0)

julia> expand_all_measures!(model)

julia> print(model)
Min 3 g(6) + z
Subject to
 T(t, xi) >= 0.0
 z >= 0.0
 g(t) + z >= 42.0
 0.5 T(0, xi) + 0.5 T(6, xi) >= 0.0, for all xi in [-1, 1]
 t in [0, 6]
 xi in Normal(μ=0.0, σ=1.0)
```
"""
function expand_all_measures!(model::InfiniteModel)
    # expand the objective if it contains measures
    if JuMP.objective_function_type(model) <: MeasureExpr
        new_obj = _possible_convert(FiniteVariableRef,
                         _expand_measures(JuMP.objective_function(model), model,
                                          _add_mapped_point_variable))
        isa(new_obj, InfiniteExpr) && error("Objective is not finite, ensure " *
                                            "all infinite variables/parameters " *
                                            "in measures are evaluated " *
                                            "completely.")
        JuMP.set_objective_function(model, new_obj)
    end
    # expand all of the constraints that contain measures
    for cindex in sort(unique(vcat(values(model.meas_to_constrs)...)))
        # expand the expression
        new_func = _possible_convert(FiniteVariableRef,
                                 _expand_measures(model.constrs[cindex].func,
                                                  model,
                                                  _add_mapped_point_variable))
        # get the necessary info
        cref = _make_constraint_ref(model, cindex)
        name = JuMP.name(cref)
        set = model.constrs[cindex].set
        curr_index = model.next_constr_index
        # delete the old cosntraint and replace it with the expanded version
        model.next_constr_index = cindex - 1
        if isa(model.constrs[cindex], BoundedScalarConstraint) && isa(new_func,
                                                                   InfiniteExpr)
            orig_bounds = model.constrs[cindex].orig_bounds
            JuMP.delete(model, cref)
            JuMP.add_constraint(model, JuMP.build_constraint(error, new_func,
                                set; parameter_bounds = orig_bounds), name)
        else
            JuMP.delete(model, cref)
            JuMP.add_constraint(model, JuMP.build_constraint(error, new_func,
                                set), name)
        end
        model.next_constr_index = curr_index
    end
    return
end

ext: jl


text:
# Configuration for nearly neutral simuation
function dfe_mixd( x::Int64 )
  dfe_mixed( x::Int64; adv_probability=0.01,
    alpha_disadv=dfe_disadv_alpha, beta_disadv=dfe_disadv_beta,
    alpha_adv=dfe_adv_alpha, beta_adv=dfe_adv_beta )
end

const nn_simtype = 1    # nearly neutral infinite alleles
type_str = "M2"
const mu_list_flag = false
const N_list     = [25,100,400]     # popsize
const N_mu_list= [1.0,2.0]   # Mutation rate as a multiple of 1.0/N
const L = 1           # number of loci
const ngens = 40000      # Generations after burn-in
const burn_in= 2.0    # generations of burn_in as a multiple of 1/mu
const use_poplist=true   # Using poplist does additional correctness checking, but increases memory usage
dfe = dfe_mixd
dfe_str = "dfe_mixd advprob=0.01 alpha_adv=1.0 alpha_disadv=0.2 beta_adv=0.01 beta_disadv=0.5"
dfe_adv_prob=0.01
dfe_adv_alpha=1.0
dfe_disadv_alpha=0.2
dfe_adv_beta=0.01
dfe_disadv_beta=0.5

max_stars_repo_name: um-tech-evolution/nearly-neutral

max_issues_repo_name: um-tech-evolution/nearly-neutral

max_forks_repo_name: um-tech-evolution/nearly-neutral


text:
# Observer pattern (can be an observer or something that is observable)
mutable struct ObserverMixin
  observers::Vector
  observables::Vector
end

ObserverMixin() = ObserverMixin([], [])

# observable patterns
function add_observer!(ob::Observer, obsv::T) where {T}
  if ~in(obsv, get_observers(ob))
    push!(get_observers(ob), obsv)
  end

  return ob
end

function notify_observers!(ob::Observer)
  for obsv in get_observers(ob)
    update!(obsv)
  end

  return ob
end


text:
type IntArray{w,T<:Unsigned,n} <: AbstractArray{T,n}
    buffer::Buffer{w,T}
    size::NTuple{n,Int}
    function IntArray(buffer::Buffer{w,T}, size::NTuple{n,Int})
        if w > bitsof(T)
            error("w = $w cannot be encoded with $T")
        end
        new(buffer, size)
    end
end

# call this function when creating an array
function call{w,T,n}(::Type{IntArray{w,T}}, dims::NTuple{n,Int}, mmap::Bool=false)
    return IntArray{w,T,n}(Buffer{w,T}(prod(dims), mmap), dims)
end

function call{w,T}(::Type{IntArray{w,T}}, len::Integer, mmap::Bool=false)
    return IntArray{w,T}((len,), mmap)
end

function call{w,T}(::Type{IntArray{w,T}}, I::Integer...)
    return IntArray{w,T}(I)
end

function call{w,T,n}(::Type{IntArray{w,T,n}}, dims::NTuple{n,Int}, mmap::Bool=false)
    return IntArray{w,T}(dims, mmap)
end

function convert{w,T,n}(::Type{IntArray{w,T,n}}, array::AbstractArray{T,n})
    iarray = IntArray{w,T}(size(array))
    @inbounds for i in eachindex(array)
        iarray[i] = array[i]
    end
    return iarray
end

function convert{w,T,n}(::Type{IntArray{w}}, array::AbstractArray{T,n})
    return convert(IntArray{w,T,n}, array)
end

linearindexing{T<:IntArray}(::Type{T}) = Base.LinearFast()

size(array::IntArray) = array.size
length(array::IntArray) = prod(array.size)
sizeof(array::IntArray) = sizeof(array.buffer.data)

@inline function getindex{w,T}(array::IntArray{w,T}, i::Integer)
    checkbounds(array, i)
    return unsafe_getindex(array, i)
end

@inline function unsafe_getindex{w,T}(array::IntArray{w,T}, i::Integer)
    return array.buffer[i] % T
end

# when I removed type parameters, array[i] fell into an infinite recursive call...
function getindex{w,T}(array::IntArray{w,T}, i::Integer, j::Integer...)
    checkbounds(array, i, j...)
    return unsafe_getindex(array, sub2ind(array.size, i, j...))
end

@inline function setindex!(array::IntArray, x::Unsigned, i::Integer)
    checkbounds(array, i)
    return unsafe_setindex!(array, x, i)
end

@inline function unsafe_setindex!{w,T}(array::IntArray{w,T}, x::Integer, i::Integer)
    array.buffer[i] = x % T
    return array
end

function setindex!{w,T}(array::IntArray{w,T}, x::Integer, i::Integer, j::Integer...)
    checkbounds(array, i, j...)
    return unsafe_setindex!(array, x, sub2ind(array.size, i, j...))
end

function similar{w,T<:Unsigned}(array::IntArray{w}, ::Type{T}, dims::Dims)
    n = length(dims)
    IntArray{w,T,n}(dims)
end

function fill!{w,T}(array::IntArray{w,T}, x::Integer)
    if x == 0
        fill0!(array.buffer)
    elseif x == (1 << w) - 1
        fill1!(array.buffer)
    else
        fill!(array.buffer, x % T)
    end
    return array
end

function copy!{w}(a::IntArray{w}, b::IntArray{w})
    len_a = length(a)
    len_b = length(b)
    if len_a < len_b
        throw(BoundsError())
    elseif len_a == len_b
        copy!(a.buffer.data, b.buffer.data)
    else
        for i in 1:len_b
            a[i] = b[i]
        end
    end
    return a
end


text:
"Julia in reverse: Parametric Inversion in pure Julia"
module ParametricInversion

const PI = ParametricInversion

using Spec
using IRTools: Statement, varargs!, insertafter!, xcall, var, IR, Block, Branch, blocks, arguments
using IRTools.Inner: Variable
import IRTools
using IRTools: varargs!, inlineable!, pis!, slots!, IR, var, xcall
using IRTools.Inner: argnames!, update!, argument!
using Mjolnir

export invertinvoke

include("util.jl")
using .Util

include("irutil.jl")
include("address.jl")

include("cfg.jl")
include("duplify.jl")

include("places.jl")
include("choose.jl")

include("reorient.jl")
include("invert.jl")
include("cycle.jl")

include("param.jl")
include("var.jl")
include("primitives.jl")
include("passes.jl")
include("traits.jl")
include("pgf.jl")

end


text:
#==
Servables!
==#
abstract type Component end
#==
 text/html
    Components
==#
"""
### html(::String) -> ::Function
------------------
Creates a servable from the provided string, which should be HTML.
#### example
"""
function html(hypertxt::String)
    return(http -> hypertxt)
end

"""
### html_file(URI::String) -> ::Function
------------------
Creates a servable which will read and return the file denoted by its path in
URI.
#### example
"""
function html_file(URI::String)
    return(http -> HTTP.Response(200, read(URI)))
end

"""
### file(URI::String) -> ::Function
------------------
Creates a servable which will read and return the file denoted by its path in
URI.
#### example
"""
function file(URI::String)
    return(http -> HTTP.Response(200, read(URI)))
end

"""
### html(::String) -> ::Function
------------------
Creates a servable from the provided string, which should be CSS.
#### example
"""
function css(css::String)
    return(http -> "")
end

"""
### css_file(URI::String) -> ::Function
------------------
Creates a servable which will read and return the file denoted by its path in
URI.
#### example
"""
function css_file(URI::String)
    http -> """<link rel="stylesheet" href="$URI">"""
end

"""
### html(::String) -> ::Function
------------------
Creates a servable from the provided string, which should be JavaScript.
#### example
"""
function js(js::String)
    return(http -> "")
end

"""
### js_file(URI::String) -> ::Function
------------------
Creates a servable which will read and return the file denoted by its path in
URI.
#### example
"""
function js_file(URI::String)
    http -> """"""
end
#==
Functions
==#
"""
### fn(::Function) -> ::Function
------------------
Turns any function into a servable. Functions can optionally take the single
    positional argument "http."
#### example

"""
function fn(f::Function)
    m = first(methods(f))
    if m.nargs > 2 | m.nargs < 1
        throw(ArgumentError("Expected either 1 or 2 arguments."))
    elseif m.nargs == 2
        http -> f(http)
    else
        http -> f()
    end
end

    #==
    Interactive
        Components
        ==#
abstract type FormComponent <: Component end
"""
### Button
name::String
action::String
label::String
f::Function
onAction::Function
html::String
------------------
A button is a form component that allows Toolips.jl to communicate with button
clicks on a web-page. The **onAction** function denotes what the button is
to do when clicked. Name will also become the id inside of the Document. \
------------------

"""
mutable struct Button <: FormComponent
    name::String
    action::String
    label::String
    f::Function
    onAction::Function
    html::String
    function Button(name::String; onAction::Function = http -> "",
         label = "Button", value = "none")
        action = "/connect/$name"
        html = """<button name="$name" value="$value">$label</button>"""
        f(http) = """<form action="$action" method="post">""" * html * "</form>"
        new(name, action, label, f, onAction, html)
    end
end

mutable struct TextArea <: FormComponent
    name::String
    text::String
    maxlength::Int64
    rows::Int64
    cols::Int64
    f::Function
    html::String
    onAction::Function
    function TextArea(name::String; maxlength::Int64 = 25, rows::Int64 = 25,
        cols::Int64 = 50, text = "~", onAction = http -> "")
        html = """
        <textarea id="$name" name="$name" rows="$rows" cols="$cols" maxlength = "$maxlength">
        $text
        </textarea>"""
        f(http) = """<form action="$action">
        $html
        </form>"""
        set_txt(to::String) = begin
            upperlower = split(html, "~")
        end
        new(name, text, maxlength, rows, cols, f, html, onAction)
    end
end

mutable struct TextBox <: FormComponent
    name::String
    text::String
    maxlength::Int64
    f::Function
    html::String
    onAction::Function
    function Text(name::String; maxlength::Int64 = 25, text::String = "",
        onAction = http -> "")
        html = """
        <input type = "text" id = "$name" name = "$name" maxlength = "$maxlength">
        $text
        </input>
        """
        f(http) = """<form action="$action">
        $html
        </form>"""
        new(name, text, mexlength, f, html, onAction)
    end
end

mutable struct RadioSet <: FormComponent
    name::String
    setdict::Dict
    f::Function
    html::String
    onAction::Function
    function RadioSet(name::String, setdict::Dict; onAction = http -> "",
        multiple = false)
        if multiple == true
            multiple = "multiple"
        else
            multiple = ""
        end
        htmlopen = """<select id="$name" name="$name" $multiple>"""
        for option in setdict
            label = setdic[option]
            htmlopen = htmlopen * """<option value="$option">$label</option>"""
        end
        html = htmlopen * "</select>"
        f(http) = """<form action="$action">
        $html
        </form>"""
        new(name, setdict, f, html, onAction)
    end
end

mutable struct Slider <: FormComponent
    name::String
    range::UnitRange
    f::Function
    html::String
    onAction::Function
    function Slider(name::String; range = 0:100, onAction = http -> "")
        min, max = range[1], range[2]
        html = """<input type="range" id="$name" name="$name"
               min="$min" max="$max">"""
        f(http) = """<form action="$action">
               $html
               </form>"""
        new(name, range, f, html, onAction)
    end
end

mutable struct Form <: FormComponent
    action::String
    f::Function
    html::String
    components::AbstractArray
    onAction::Function
    function Form(components...; onAction::Any = http -> "", action::String = "")
        html = """<form action="$action" method = "GET">"""
        components = [c for c in components]
        for comp in components
            html = html * comp.html
        end
        html = html * "</form>"
        f(http) = html
        new(action, f, html, components, onAction)
    end
end
#==
Page
==#
mutable struct Page <: Component
    f::Function
    components::AbstractVector
    add::Function
    # Document
    title::String
    icon::String
    function Page(components = [], title = "Toolips Webapp"; icon = "/")
        f(http) = generate_page(http, title, components, icon)
        add(x::Function) = push!(components, x)
        add(x::Component) = push!(components, x)
        new(f, components, add, title, icon)
    end
end

function generate_page(http, title, components, icon = "/")
    html = "" * "<title>$title</title>" * ""
    write(http, html)
    body = ""
    for comp in components
        try
            body = body * comp.f(http)
        catch
            body = body * comp(http)
        end
    end
    body
end
#==
Other
    Components
==#
abstract type ComponentPart <: Component end
mutable struct Context <: ComponentPart
    codestrings::AbstractArray
    update::Function
    fillRect::Function
    strokeRect::Function
    beginPath::Function
    closePath::Function
    moveTo::Function
    lineTo::Function
    fill::Function
    stroke::Function
    arc::Function
    rect::Function
    function Context(ctx::String = "2d", name::String = "canvas")
        codestrings = []
        push!(codestrings, """var canvas = document.getElementById("$name");
            var $namectx = canvas.getContext($ctx);""")
        namectx = "$name" * ctx
        update() = script = join(codestrings)
        rect(x::Integer, y::Integer, w::Integer, h::Integer) =
            push!(codestrings, "$namectx.rect($x, $y, $w, $h);")
        fillRect(x::Integer, y::Integer, w::Integer, h::Integer) = begin
            push!(codestrings, "$namectx.fillRect($x, $y, $w, $h);")
        end
        strokeRect(x::Integer, y::Integer, w::Integer, h::Integer) =
            push!(codestrings, "$namectx.strokeRect($x, $y, $w, $h);")
        arc(x, y, r, sa, ea, ccw) = begin
            push!(codestrings, "$namectx.arc($x, $y, $r, $sa, $ea, $ccw);")
        end
        beginPath() = begin
            push!(codestrings, "$namectx.beginPath();")
        end
        closePath() = begin
            push!(codestrings, "$namectx.closePath();")
        end
        moveTo(x::Integer, y::Integer) = begin
            push!(codestrings, "$namectx.moveTo($x, $y);")
        end
        lineTo(x::Integer, y::Integer) = begin
            push!(codestrings, "$namectx.lineTo($x, $y);")
        end
        fill() = begin
            push!(codestrings, "$namectx.fill();")
        end
        stroke() = begin
            push!(codestrings, "$namectx.stroke();")
        end
        new(codestrings, update, fillRect, strokeRect, beginPath, closePath, moveTo,
        lineTo, fill, stroke, arc, rect)
    end
end
mutable struct Canvas <: Component
    name::String
    width::Int64
    height::Int64
    html::String
    ctx::Context
    context::Function
    f::Function
    function Canvas(name = "canvas"; width = 200, height = 200, mode = "2d")
        ctx = Context(mode, name)
        html = """<canvas id="$name" width="$width" height="$height"></canvas>"""
        context(f::Function) = f(ctx)
        f(http) = """<canvas id="$name" width="$width" height="$height"script = '""" * join(ctx.codestrings) * "'</script></canvas>"
        new(name, width, height, html, ctx, context, f)
    end
end

abstract type ListComponent <: Component end

mutable struct List <: ListComponent
    name::String
    label::String
    html::String
    f::Function
    href::String
    function List(name::String = "list"; label = "hello world!", href = "")
        if href != ""
            href = "href='$href'"
        end
        html = """
$label
"""
        f(http) = html
        new(name, label, html, f, href)
    end
end

mutable struct UnorderedList
    name::String
    html::String
    f::Function
    lists::Vector{List}
    function UnorderedList(name::String = "ul", comps::Array{List} = [])
        html = "
"

        f(http) = "
" * join([l.f(http) for l in comps]) * "
"
        new(name, html, f, comps)
    end
     function UnorderedList(name::String = "ul")
        UnorderedList(name, comps)
    end
         function UnorderedList(comps::Vector{List})
        name = ""
        UnorderedList(name, comps)
    end
end

mutable struct A <: ComponentPart
    name::String
    href::String
    f::Function
    html::String
    function A(name = "a"; href = "#", label = "a")
        f(http) = "$label"
        html = "$label"
        new(name, href, f, html)
    end
end

mutable struct DropDown <: Component
    name::String
    html::String
    href::String
    As::AbstractArray{Component}
    f::Function
    function DropDrown(name::String = "dropdown", As::A ...;
        label::String = "dropdown", href::String = "#")
        html = "

"
        f(http) = begin
            """

              <button class="dropbtn">$label
                </button>

                """ * join([a.f(http) for a in As]) * "

"
            end
            new(name, html, href, As, f)
        end
end
include("frontend.jl")
include("methods.jl")
include("../server/serve.jl")


text:
# ---
# title: 831. Masking Personal Information
# id: problem831
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: String
# link: <https://leetcode.com/problems/masking-personal-information/description/>
# hidden: true
# ---
#
# We are given a personal information string `S`, which may represent either
# **an email address** or **a phone number.**
#
# We would like to mask this personal information according to the following
# rules:
#
#
# _**1\. Email address:**_
#
# We define a  **name** to be a string of `length ≥ 2` consisting of only
# lowercase letters `a-z` or uppercase letters `A-Z`.
#
# An email address starts with a name, followed by the symbol `'@'`, followed by
# a name, followed by the dot `'.'` and followed by a name.
#
# All email addresses are guaranteed to be valid and in the format of
# `"name1@name2.name3".`
#
# To mask an email, **all names must be converted to lowercase** and **all
# letters between the first and last letter of the first name** must be replaced
# by 5 asterisks `'*'`.
#
#
# _**2\. Phone number:**_
#
# A phone number is a string consisting of only the digits `0-9` or the
# characters from the set `{'+', '-', '(', ')', ' '}.` You may assume a phone
# number contains 10 to 13 digits.
#
# The last 10 digits make up the local number, while the digits before those
# make up the country code. Note that the country code is optional. We want to
# expose only the last 4 digits and mask all other digits.
#
# The local number should be formatted and masked as `"***-***-1111", `where `1`
# represents the exposed digits.
#
# To mask a phone number with country code like `"+111 111 111 1111"`, we write
# it in the form `"+***-***-***-1111".`  The `'+'` sign and the first `'-'` sign
# before the local number should only exist if there is a country code.  For
# example, a 12 digit phone number mask should start with `"+**-"`.
#
# Note that extraneous characters like `"(", ")", " "`, as well as extra dashes
# or plus signs not part of the above formatting scheme should be removed.
#
#
#
# Return the correct "mask" of the information provided.
#
#
#
# **Example 1:**
#
#
#
#     Input: "LeetCode@LeetCode.com"
#     Output: "l*****e@leetcode.com"
#     **Explanation:  **All names are converted to lowercase, and the letters between the
#                  first and last letter of the first name is replaced by 5 asterisks.
#                  Therefore, "leetcode" -> "l*****e".
#
#
# **Example 2:**
#
#
#
#     Input: "AB@qq.com"
#     Output: "a*****b@qq.com"
#     **Explanation:  **There must be 5 asterisks between the first and last letter
#                  of the first name "ab". Therefore, "ab" -> "a*****b".
#
#
# **Example 3:**
#
#
#
#     Input: "1(234)567-890"
#     Output: "***-***-7890"
#     Explanation:  10 digits in the phone number, which means all digits make up the local number.
#
#
# **Example 4:**
#
#
#
#     Input: "86-(10)12345678"
#     Output: "+**-***-***-5678"
#     Explanation:  12 digits, 2 digits for country code and 10 digits for local number.
#
#
# **Notes:**
#
#   1. `S.length <= 40`.
#   2. Emails have length at least 8.
#   3. Phone numbers have length at least 10.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
module DiscontinuousFESpacesTests

using Test
using Gridap.Arrays
using Gridap.ReferenceFEs
using Gridap.Geometry
using Gridap.FESpaces

domain =(0,1,0,1)
partition = (3,3)
model = CartesianDiscreteModel(domain,partition)

order = 3

reffe = ReferenceFE(lagrangian,Float64,order)
V = FESpace(model,reffe,conformity=:L2)
@test isa(V,UnconstrainedFESpace)
@test isa(V.cell_dofs_ids,Table{Int32})
test_single_field_fe_space(V)
@test num_free_dofs(V) == num_cells(model)*(order+1)^2

U = V
f(x) = sin(pi*x[1])*cos(2*pi*x[2])
fh = interpolate(f,U)
uh = FEFunction(V,rand(num_free_dofs(V)))

reffe = ReferenceFE(lagrangian,Float64,order;space=:S)
V = FESpace(model,reffe,conformity=:L2)
@test isa(V,UnconstrainedFESpace)
@test isa(V.cell_dofs_ids,Table{Int32})
test_single_field_fe_space(V)

U = V
fh = interpolate(f,U)
uh = FEFunction(V,rand(num_free_dofs(V)))

end # module


text:
# singleton for indicating if no default arguments are present
struct NoDefault end
const NO_DEFAULT = NoDefault()

"""
    @addlogprob!(ex)

Add the result of the evaluation of `ex` to the joint log probability.
"""
macro addlogprob!(ex)
    return quote
        acclogp!($(esc(:(_varinfo))), $(esc(ex)))
    end
end

"""
    getargs_dottilde(x)

Return the arguments `L` and `R`, if `x` is an expression of the form `L .~ R` or
`(~).(L, R)`, or `nothing` otherwise.
"""
getargs_dottilde(x) = nothing
function getargs_dottilde(expr::Expr)
    return MacroTools.@match expr begin
        (.~)(L_, R_) => (L, R)
        (~).(L_, R_) => (L, R)
        # Julia 1.6: see https://github.com/TuringLang/Turing.jl/issues/1525
        (L_ .~ R_) => (L, R)
        x_ => nothing
    end
end

"""
    getargs_tilde(x)

Return the arguments `L` and `R`, if `x` is an expression of the form `L ~ R`, or `nothing`
otherwise.
"""
getargs_tilde(x) = nothing
function getargs_tilde(expr::Expr)
    return MacroTools.@match expr begin
        (~)(L_, R_) => (L, R)
        x_ => nothing
    end
end

############################################
# Julia 1.2 temporary fix - Julia PR 33303 #
############################################
function to_namedtuple_expr(syms, vals=syms)
    if length(syms) == 0
        nt = :(NamedTuple())
    else
        nt_type = Expr(:curly, :NamedTuple,
            Expr(:tuple, QuoteNode.(syms)...),
            Expr(:curly, :Tuple, [:(Core.Typeof($x)) for x in vals]...)
        )
        nt = Expr(:call, :($(DynamicPPL.namedtuple)), nt_type, Expr(:tuple, vals...))
    end
    return nt
end

if VERSION == v"1.2"
    @eval function namedtuple(::Type{NamedTuple{names, T}}, args::Tuple) where {names, T <: Tuple}
        if length(args) != length(names)
            throw(ArgumentError("Wrong number of arguments to named tuple constructor."))
        end
        # Note T(args) might not return something of type T; e.g.
        # Tuple{Type{Float64}}((Float64,)) returns a Tuple{DataType}
        $(Expr(:splatnew, :(NamedTuple{names,T}), :(T(args))))
    end
else
    function namedtuple(::Type{NamedTuple{names, T}}, args::Tuple) where {names, T <: Tuple}
        return NamedTuple{names, T}(args)
    end
end

#####################################################
# Helper functions for vectorize/reconstruct values #
#####################################################

vectorize(d::UnivariateDistribution, r::Real) = [r]
vectorize(d::MultivariateDistribution, r::AbstractVector{<:Real}) = copy(r)
vectorize(d::MatrixDistribution, r::AbstractMatrix{<:Real}) = copy(vec(r))

# NOTE:
# We cannot use reconstruct{T} because val is always Vector{Real} then T will be Real.
# However here we would like the result to be specifric type, e.g. Array{Dual{4,Float64}, 2},
# otherwise we will have error for MatrixDistribution.
# Note this is not the case for MultivariateDistribution so I guess this might be lack of
# support for some types related to matrices (like PDMat).
reconstruct(d::UnivariateDistribution, val::AbstractVector) = val[1]
reconstruct(d::MultivariateDistribution, val::AbstractVector) = copy(val)
function reconstruct(d::MatrixDistribution, val::AbstractVector)
    return reshape(copy(val), size(d))
end
function reconstruct!(r, d::Distribution, val::AbstractVector)
    return reconstruct!(r, d, val)
end
function reconstruct!(r, d::MultivariateDistribution, val::AbstractVector)
    r .= val
    return r
end
function reconstruct(d::Distribution, val::AbstractVector, n::Int)
    return reconstruct(d, val, n)
end
function reconstruct(d::UnivariateDistribution, val::AbstractVector, n::Int)
    return copy(val)
end
function reconstruct(d::MultivariateDistribution, val::AbstractVector, n::Int)
    return copy(reshape(val, size(d)[1], n))
end
function reconstruct(d::MatrixDistribution, val::AbstractVector, n::Int)
    tmp = reshape(val, size(d)[1], size(d)[2], n)
    orig = [tmp[:, :, i] for i in 1:size(tmp, 3)]
    return orig
end
function reconstruct!(r, d::Distribution, val::AbstractVector, n::Int)
    return reconstruct!(r, d, val, n)
end
function reconstruct!(r, d::MultivariateDistribution, val::AbstractVector, n::Int)
    r .= val
    return r
end

# Uniform random numbers with range 4 for robust initializations
# Reference: https://mc-stan.org/docs/2_19/reference-manual/initialization.html
randrealuni(rng::Random.AbstractRNG) = 4 * rand(rng) - 2
randrealuni(rng::Random.AbstractRNG, args...) = 4 .* rand(rng, args...) .- 2

const Transformable = Union{PositiveDistribution,UnitDistribution,TransformDistribution,
                            SimplexDistribution,PDMatDistribution}
istransformable(dist) = false
istransformable(::Transformable) = true

#################################
# Single-sample initialisations #
#################################

inittrans(rng, dist::UnivariateDistribution) = invlink(dist, randrealuni(rng))
inittrans(rng, dist::MultivariateDistribution) = invlink(dist, randrealuni(rng, size(dist)[1]))
inittrans(rng, dist::MatrixDistribution) = invlink(dist, randrealuni(rng, size(dist)...))

################################
# Multi-sample initialisations #
################################

inittrans(rng, dist::UnivariateDistribution, n::Int) = invlink(dist, randrealuni(rng, n))
function inittrans(rng, dist::MultivariateDistribution, n::Int)
    return invlink(dist, randrealuni(rng, size(dist)[1], n))
end
function inittrans(rng, dist::MatrixDistribution, n::Int)
    return invlink(dist, [randrealuni(rng, size(dist)...) for _ in 1:n])
end


text:
mutable struct UniformBlock
	ubo::GLuint
	index::GLuint
	binding::GLuint
	size::UInt32
	name::Symbol

	UniformBlock() = new(0, 0, 0, 0)
end

isvalid(block::UniformBlock) = block.ubo != 0

function init(block::UniformBlock, program::GLuint, blockIndex::GLuint)
	block.index = blockIndex
	val = GLint[1]

	block.binding = block.index
	glUniformBlockBinding(program, block.index, block.binding)

	glGetActiveUniformBlockiv(program, block.index, GL_UNIFORM_BLOCK_DATA_SIZE, val)
	@assert glGetError() == GL_NO_ERROR
	block.size = val[1]

	glGetActiveUniformBlockiv(program, block.index, GL_UNIFORM_BLOCK_NAME_LENGTH, val)
	@assert glGetError() == GL_NO_ERROR
	nameBuf = Array{UInt8}(undef, val[1])
	glGetActiveUniformBlockName(program, block.index, length(nameBuf), C_NULL, nameBuf)
	@assert glGetError() == GL_NO_ERROR
	block.name = Symbol(nameBuf[1:end-1])

	ubo = GLuint[0]
	glGenBuffers(1, ubo)
	block.ubo = ubo[1]

	data = zeros(UInt8, block.size)
	setbufferdata(block, data)
end

function done(block::UniformBlock)
	if isvalid(block)
		ubo = GLuint[block.ubo]
		glDeleteBuffers(1, ubo)
		block.ubo = 0
	end
end

function setbufferdata(block::UniformBlock, data)
	@assert sizeof(data) == block.size
	glBindBuffer(GL_UNIFORM_BUFFER, block.ubo)
	glBufferData(GL_UNIFORM_BUFFER, sizeof(data), data, GL_DYNAMIC_DRAW)
end

function apply(block::UniformBlock)
	glBindBufferBase(GL_UNIFORM_BUFFER, block.binding, block.ubo)
end

mutable struct UniformVar
	name::Symbol
	varType::DataType
	location::GLint
	offset::Int32
	arraySize::Int32
	arrayStride::Int32
	matrixStride::Int32
	blockId::Int32
end

inblock(var::UniformVar) = var.blockId >= 0

function getptr(var::UniformVar, buffer::Vector{UInt8}, index::Int)
	@assert 1 <= index <= var.arraySize
	ptr = pointer(buffer) + var.offset + var.arrayStride * (index - 1)
	@assert convert(UInt64, pointer(buffer)) <= convert(UInt64, ptr + sizeof(var.varType)) <= convert(UInt64, pointer(buffer) + length(buffer))
	return convert(Ptr{var.varType}, ptr)
end

function setvalue(var::UniformVar, buffer::Vector{UInt8}, value, index::Int = 1)
	ptr = getptr(var, buffer, index)
	unsafe_store!(ptr, value)
end

function setvalue(var::UniformVar, buffer::Vector{UInt8}, array::Array, index::Int = 1)
	ptr = getptr(var, buffer, index)
	# @assert ismatrix(var.varType)
	if sizeof(var.varType) != sizeof(array)
		@info("$(var.name) has size $(sizeof(var.varType)), trying to set value of size $(sizeof(array))")
		@assert sizeof(var.varType) == sizeof(array)
	end

	if !ismatrix(var.varType) || size(var.varType, 1) == 4
		unsafe_copyto!(convert(Ptr{UInt8}, ptr), convert(Ptr{UInt8}, pointer(array)), sizeof(var.varType))
	else
		elType = eltype(var.varType)
		elSize = sizeof(elType)
		rowSize = size(var.varType, 1)
		@assert rowSize < 4
		pSrc = convert(Ptr{elType}, pointer(array))
		pDst = convert(Ptr{elType}, ptr)
		for col = 1:length(array)/rowSize
			unsafe_copyto!(pDst, pSrc, rowSize)
			pDst += 4elSize
			pSrc += rowSize * elSize
		end
		@assert pSrc == pointer(array) + length(array) * sizeof(eltype(array))
	end
end

load_ptr(p::Ptr{T}) where T = unsafe_load(p)

function getvalue(var::UniformVar, buffer::Vector{UInt8}, index::Int = 1)
	ptr = getptr(var, buffer, index)

	if ismatrix(var.varType)
		result = Array{eltype(var.varType)}(undef, size(var.varType))
		@assert sizeof(result) == sizeof(var.varType)
		unsafe_copyto!(convert(Ptr{UInt8}, pointer(result)), convert(Ptr{UInt8}, ptr), sizeof(var.varType))
		return result
	else
		# for some reason unsafe_load can't be called directly here, wrap it in a function
		# return unsafe_load(ptr)
		return load_ptr(ptr)
	end
end


text:
function connect!(c, j, i, σ = 1e-6)
    W = sparse(c.I, c.J, c.W, length(c.rowptr) - 1, length(c.colptr) - 1)
    W[i, j] = σ * randn(Float32)
    c.rowptr, c.colptr, c.I, c.J, c.index, c.W = dsparse(W)
    c.tpre, c.tpost, c.Apre, c.Apost = zero(c.W), zero(c.W), zero(c.W), zero(c.W)
    return nothing
end

function dsparse(A)
    At = sparse(A')
    colptr = A.colptr
    rowptr = At.colptr
    I = rowvals(A)
    V = nonzeros(A)
    J = zero(I)
    # FIXME: Breaks when A is empty
    for j in 1:(length(colptr) - 1)
        J[colptr[j]:(colptr[j+1] - 1)] .= j
    end
    index = zeros(size(I))
    coldown = zeros(eltype(index), length(colptr) - 1)
    for i in 1:(length(rowptr) - 1)
        for st in rowptr[i]:(rowptr[i+1] - 1)
            j = At.rowval[st]
            index[st] = colptr[j] + coldown[j]
            coldown[j] += 1
        end
    end
    # Test.@test At.nzval == A.nzval[index]
    rowptr, colptr, I, J, index, V
end

function record!(obj)
    for (key, val) in obj.records
        if isa(key, Tuple)
            sym, ind = key
            push!(val, getindex(getfield(obj, sym),ind))
        else
            push!(val, copy(getfield(obj, key)))
        end
    end
end

function monitor(obj, keys)
    for key in keys
        if isa(key, Tuple)
            sym, ind = key
        else
            sym = key
        end
        typ = typeof(getfield(obj, sym))
        obj.records[key] = Vector{typ}()
    end
end

function monitor(objs::Array, keys)
    for obj in objs
        monitor(obj, keys)
    end
end

function getrecord(p, sym)
    key = sym
    for (k,val) in p.records
        isa(k, Tuple) && k[1] == sym && (key = k)
    end
    p.records[key]
end

function clear_records(obj)
    for (key, val) in obj.records
        empty!(val)
    end
end

@inline function exp32(x::Float32)
    x = ifelse(x < -10f0, -32f0, x)
    x = 1f0 + x / 32f0
    x *= x; x *= x; x *= x; x *= x; x *= x
    return x
end

@inline function exp256(x::Float32)
    x = ifelse(x < -10f0, -256f0, x)
    x = 1.0f0 + x / 256.0f0
    x *= x; x *= x; x *= x; x *= x
    x *= x; x *= x; x *= x; x *= x
    return x
end

macro symdict(x...)
    ex = Expr(:block)
    push!(ex.args, :(d = Dict{Symbol,Any}()))
    for p in x
        push!(ex.args, :(d[$(QuoteNode(p))] = $(esc(p))))
    end
    push!(ex.args, :(d))
    return ex
end

snn_kw_str_param(x::Symbol) = (x,)
function snn_kw_str_param(x::Expr)
    if x.head == :(<:)
        return (x.args...,)
    elseif x.head == :(=)
        if x.args[1] isa Expr && x.args[1].head == :(<:)
            return (x.args[1].args..., x.args[2])
        elseif x.args[1] isa Symbol
            return (x.args[1], Any, x.args[2])
        end
    end
    error("Can't handle param Expr: $x")
end
snn_kw_str_field(x::Symbol) = (x,)
function snn_kw_str_field(x::Expr)
    if x.head == :(::)
        return (x.args...,)
    elseif x.head == :(=)
        return (x.args[1].args[1:2]..., x.args[2])
    end
    error("Can't handle field Expr: $x")
end
function snn_kw_str_kws(x::Tuple)
    if 1 <= length(x) <= 2
        return x[1]
    elseif length(x) == 3
        return Expr(:kw, x[1], x[3])
    end
end
function snn_kw_str_kws_types(x::Tuple)
    if 1 <= length(x) <= 2
        return Expr(:(::), x[1], x[2])
    elseif length(x) == 3
        return Expr(:kw, x[1], x[3])
    end
end
struct KwStrSentinel end
function snn_kw_str_sentinels(x)
    if length(x) == 1
        return (x[1], Any, :(KwStrSentinel()))
    elseif length(x) == 2
        return (x[1], Any, :(KwStrSentinel()))
    else
        return x
    end
end
snn_kw_str_sentinel_check(x) =
    :(if $(x[1]) isa KwStrSentinel
        $(x[1]) = $(length(x) > 1 ? x[2] : Any)
    end)
"A minimal implementation of `Base.@kwdef` with default type parameter support"
macro snn_kw(str)
    str_abs = nothing
    if str.args[2] isa Expr && str.args[2].head == :(<:)
        # Lower abstract type
        str_abs = str.args[2].args[2]
        str.args[2] = str.args[2].args[1]
    end
    if str.args[2] isa Symbol
        # No type params
        str_name = str.args[2]
        str_params = []
    else
        # Has type params
        str_name = str.args[2].args[1]
        str_params = map(snn_kw_str_param, str.args[2].args[2:end])
    end
    @assert str_name isa Symbol
    @assert str_abs isa Union{Symbol,Nothing}
    str_fields = map(snn_kw_str_field, filter(x->!(x isa LineNumberNode),
                                          str.args[3].args))

    # Remove default type params
    if length(str_params) > 0
        idx = 1
        for idx in 2:length(str.args[2].args)
            param = str_params[idx-1]
            if length(param) == 1
                str.args[2].args[idx] = param[1]
            else
                str.args[2].args[idx] = Expr(:(<:), param[1:2]...)
            end
        end
    end

    # Remove default field values
    idx = 1
    subidx = 1
    for idx in 1:length(str.args[3].args)
        if !(str.args[3].args[idx] isa LineNumberNode)
            field = str_fields[subidx]
            if length(field) == 1
                str.args[3].args[idx] = field[1]
            else
                str.args[3].args[idx] = Expr(:(::), field[1:2]...)
            end
            subidx += 1
        end
    end

    # Replace abstract type
    if str_abs !== nothing
        str.args[2] = Expr(:(<:), str.args[2], str_abs)
    end

    # Use sentinels to track if type param kwargs are assigned
    ctor_params = snn_kw_str_sentinels.(str_params)
    ctor_params_bodies = snn_kw_str_sentinel_check.(str_params)

    # Constructor accepts field values and type params as kwargs
    ctor_kws = Expr(:parameters,
                    map(snn_kw_str_kws, str_fields)...,
                    map(snn_kw_str_kws_types, ctor_params)...)
    ctor_sig = Expr(:call, str_name, ctor_kws)
    ctor_call = if length(str_params) > 0
        Expr(:curly, str_name, first.(str_params)...)
    else
        str_name
    end
    ctor_body = Expr(:block,
                     ctor_params_bodies...,
                     Expr(:call, ctor_call, first.(str_fields)...))
    ctor = Expr(:function, ctor_sig, ctor_body)

    return quote
        $(esc(str))
        $(esc(ctor))
    end
end


text:
"""
Abstract type for time series storage implementations.

All subtypes must implement:
- add_time_series_reference!
- check_read_only
- clear_time_series!
- deserialize_time_series
- get_compression_settings
- get_num_time_series
- is_read_only
- remove_time_series!
- serialize_time_series!
"""
abstract type TimeSeriesStorage end

const DEFAULT_COMPRESSION = false

@scoped_enum(CompressionTypes, BLOSC = 0, DEFLATE = 1,)

"""
Provides customization of HDF5 compression settings.
Refer to the HDF5.jl and HDF5 documention for more information.
"""
struct CompressionSettings
    "Controls whether compression is enabled."
    enabled::Bool
    "Specifies the type of compression to use."
    type::CompressionTypes
    "Supported values are 0-9. Higher values deliver better compression ratios but take longer."
    level::Int
    "Controls whether to enable the shuffle filter. Used with DEFLATE."
    shuffle::Bool
end

function CompressionSettings(;
    enabled = DEFAULT_COMPRESSION,
    type = CompressionTypes.DEFLATE,
    level = 3,
    shuffle = true,
)
    return CompressionSettings(enabled, type, level, shuffle)
end

function make_time_series_storage(;
    in_memory = false,
    filename = nothing,
    directory = nothing,
    compression = CompressionSettings(),
)
    if in_memory
        storage = InMemoryTimeSeriesStorage()
    elseif !isnothing(filename)
        storage = Hdf5TimeSeriesStorage(; filename = filename, compression = compression)
    else
        storage =
            Hdf5TimeSeriesStorage(true; directory = directory, compression = compression)
    end

    return storage
end

const COMPONENT_name_DELIMITER = "__"

function make_component_name(component_uuid::UUIDs.UUID, name::AbstractString)
    return string(component_uuid) * COMPONENT_name_DELIMITER * name
end

function deserialize_component_name(component_name::AbstractString)
    data = split(component_name, COMPONENT_name_DELIMITER)
    component = UUIDs.UUID(data[1])
    name = data[2]
    return component, name
end

function serialize(storage::TimeSeriesStorage, file_path::AbstractString)
    if storage isa Hdf5TimeSeriesStorage
        if abspath(get_file_path(storage)) == abspath(file_path)
            if !is_read_only(storage)
                error("Attempting to overwrite identical time series file")
            end

            @debug "Skip time series serialization because the paths are identical" _group =
                LOG_GROUP_TIME_SERIES
            return
        end

        # The data is currently in a temp file, so we can just make a copy.
        copy_file(get_file_path(storage), file_path)
    elseif storage isa InMemoryTimeSeriesStorage
        convert_to_hdf5(storage, file_path)
    else
        error("unsupported type $(typeof(storage))")
    end

    @info "Serialized time series data to $file_path."
end


text:
# Adapted from ParameterHandling.jl with the following license.
#=
Copyright (c) 2020 Invenia Technical Computing Corporation

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
=#

"""
    flatten(x)

Returns a "flattened" representation of `x` as a vector of real numbers, and a function
`unflatten` that takes a vector of reals of the same length and returns an object of the
same type as `x`.

`unflatten` is the inverse of `flatten`, so
```julia
julia> x = (randn(5), 5.0, (a=5.0, b=randn(2, 3)));

julia> v, unflatten = flatten(x);

julia> x == unflatten(v)
true
```
"""
function flatten end

maybeflatten(x::Real) = x
maybeflatten(x) = flatten(x)

function flatten(x::Real)
    v = [x]
    unflatten_to_Real(v) = only(v)
    return v, unflatten_to_Real
end

flatten(x::Vector{<:Real}) = (identity.(x), identity)

function flatten(x::AbstractVector)
    x_vecs_and_backs = map(val -> flatten(val), identity.(x))
    x_vecs, backs = first.(x_vecs_and_backs), last.(x_vecs_and_backs)
    function Vector_from_vec(x_vec)
        sz = _cumsum(map(_length, x_vecs))
        x_Vec = [backs[n](x_vec[sz[n] - _length(x_vecs[n]) + 1:sz[n]]) for n in eachindex(x)]
        return x_Vec
    end
    return reduce(vcat, x_vecs), Vector_from_vec
end

function flatten(x::AbstractArray)
    x_vec, from_vec = flatten(vec(identity.(x)))
    Array_from_vec(x_vec) = reshape(from_vec(x_vec), size(x))
    return identity.(x_vec), Array_from_vec
end

# Zygote can return a sparse vector co-tangent
# even if the input is a vector. This is causing
# issues in the rrule definition of Unflatten
flatten(x::SparseVector) = flatten(Array(x))

function flatten(x::JuMP.Containers.DenseAxisArray)
    x_vec, from_vec = flatten(vec(identity.(x.data)))
    Array_from_vec(x_vec) = JuMP.Containers.DenseAxisArray(reshape(from_vec(x_vec), size(x)), axes(x)...)
    return identity.(x_vec), Array_from_vec
end

function flatten(x::SparseMatrixCSC)
    x_vec, from_vec = flatten(x.nzval)
    Array_from_vec(x_vec) = SparseMatrixCSC(x.m, x.n, x.colptr, x.rowval, from_vec(x_vec))
    return identity.(x_vec), Array_from_vec
end

function flatten(x::Tuple)
    x_vecs_and_backs = map(val -> flatten(val), x)
    x_vecs, x_backs = first.(x_vecs_and_backs), last.(x_vecs_and_backs)
    lengths = map(_length, x_vecs)
    sz = _cumsum(lengths)
    function unflatten_to_Tuple(v)
        map(x_backs, lengths, sz) do x_back, l, s
            return x_back(v[s - l + 1:s])
        end
    end
    return reduce(vcat, x_vecs), unflatten_to_Tuple
end

function flatten(x::NamedTuple)
    x_vec, unflatten = flatten(values(x))
    function unflatten_to_NamedTuple(v)
        v_vec_vec = unflatten(v)
        return typeof(x)(v_vec_vec)
    end
    return identity.(x_vec), unflatten_to_NamedTuple
end

function flatten(d::AbstractDict, ks = collect(keys(d)))
    _d = OrderedDict(k => d[k] for k in ks)
    d_vec, unflatten = flatten(identity.(collect(values(_d))))
    function unflatten_to_Dict(v)
        v_vec_vec = unflatten(v)
        return _build_ordered_dict(v_vec_vec, keys(_d))
    end
    return identity.(d_vec), unflatten_to_Dict
end
function _build_ordered_dict(vals, keys)
    OrderedDict(key => vals[n] for (n, key) in enumerate(keys))
end
function ChainRulesCore.rrule(::typeof(_build_ordered_dict), vals, keys)
    _build_ordered_dict(vals, keys), Δ -> begin
        NoTangent(), values(Δ), NoTangent()
    end
end

function flatten(x)
    v, un = flatten(ntfromstruct(x))
    return identity.(v), Unflatten(x, y -> structfromnt(typeof(x), un(y)))
end

function zygote_flatten(::Real, x::Real)
    v = [x]
    unflatten_to_Real(v) = only(v)
    return v, unflatten_to_Real
end

zygote_flatten(::Vector{<:Real}, x::Vector{<:Real}) = (identity.(x), identity)

function zygote_flatten(x1::AbstractVector, x2::AbstractVector)
    x_vecs_and_backs = map((val) -> zygote_flatten(val[1], val[2]), zip(identity.(x1), identity.(x2)))
    x_vecs, backs = first.(x_vecs_and_backs), last.(x_vecs_and_backs)
    function Vector_from_vec(x_vec)
        sz = _cumsum(map(_length, x_vecs))
        x_Vec = [backs[n](x_vec[sz[n] - _length(x_vecs[n]) + 1:sz[n]]) for n in eachindex(x2)]
        return x_Vec
    end
    return reduce(vcat, x_vecs), Vector_from_vec
end

function zygote_flatten(x1::AbstractArray, x2::AbstractArray)
    x_vec, from_vec = zygote_flatten(vec(identity.(x1)), vec(identity.(x2)))
    Array_from_vec(x_vec) = reshape(from_vec(x_vec), size(x2))
    return identity.(x_vec), Array_from_vec
end

# Zygote can return a sparse vector co-tangent
# even if the input is a vector. This is causing
# issues in the rrule definition of Unflatten
zygote_flatten(x1::SparseVector, x2::SparseVector) = zygote_flatten(Array(x1), Array(x2))

function zygote_flatten(x1::JuMP.Containers.DenseAxisArray, x2::NamedTuple)
    x_vec, from_vec = zygote_flatten(vec(identity.(x1.data)), vec(identity.(x2.data)))
    Array_from_vec(x_vec) = JuMP.Containers.DenseAxisArray(reshape(from_vec(x_vec), size(x2)), axes(x2)...)
    return identity.(x_vec), Array_from_vec
end

function zygote_flatten(x1::JuMP.Containers.DenseAxisArray, x2::JuMP.Containers.DenseAxisArray)
    x_vec, from_vec = zygote_flatten(vec(identity.(x1.data)), vec(identity.(x2.data)))
    Array_from_vec(x_vec) = JuMP.Containers.DenseAxisArray(reshape(from_vec(x_vec), size(x2)), axes(x2)...)
    return identity.(x_vec), Array_from_vec
end

function zygote_flatten(x1::SparseMatrixCSC, x2::SparseMatrixCSC)
    x_vec, from_vec = zygote_flatten(x1.nzval, x2.nzval)
    Array_from_vec(x_vec) = SparseMatrixCSC(x.m, x.n, x.colptr, x.rowval, from_vec(x_vec))
    return identity.(x_vec), Array_from_vec
end

function zygote_flatten(x1::Tuple, x2::Tuple)
    x_vecs_and_backs = map(val -> zygote_flatten(val[1], val[2]), zip(x1, x2))
    x_vecs, x_backs = first.(x_vecs_and_backs), last.(x_vecs_and_backs)
    lengths = map(_length, x_vecs)
    sz = _cumsum(lengths)
    function unflatten_to_Tuple(v)
        map(x_backs, lengths, sz) do x_back, l, s
            return x_back(v[s - l + 1:s])
        end
    end
    return reduce(vcat, x_vecs), unflatten_to_Tuple
end

function zygote_flatten(x1, x2::Tangent)
    zygote_flatten(x1, ntfromstruct(x2).backing)
end
function zygote_flatten(x1, x2::NamedTuple)
    zygote_flatten(ntfromstruct(x1), x2)
end

function zygote_flatten(x1::NamedTuple, x2::NamedTuple)
    x_vec, unflatten = zygote_flatten(values(x1), values(x2))
    function unflatten_to_NamedTuple(v)
        v_vec_vec = unflatten(v)
        return typeof(x)(v_vec_vec)
    end
    return identity.(x_vec), unflatten_to_NamedTuple
end

function zygote_flatten(d1::AbstractDict, d2::AbstractDict, ks = collect(keys(d2)))
    _d1 = OrderedDict(k => d1[k] for k in ks)
    _d2 = OrderedDict(k => d2[k] for k in ks)
    d_vec, unflatten = zygote_flatten(identity.(collect(values(_d1))), identity.(collect(values(_d2))))
    function unflatten_to_Dict(v)
        v_vec_vec = unflatten(v)
        return OrderedDict(key => v_vec_vec[n] for (n, key) in enumerate(ks))
    end
    return identity.(d_vec), unflatten_to_Dict
end

function zygote_flatten(x1, x2)
    v, un = zygote_flatten(ntfromstruct(x1), ntfromstruct(x2))
    return identity.(v), Unflatten(x1, y -> structfromnt(typeof(x2), un(y)))
end

_length(x) = length(x)
_length(::Nothing) = 0

function ChainRulesCore.rrule(::typeof(flatten), x)
    d_vec, un = flatten(x)
    return (d_vec, un), Δ -> begin
        (NoTangent(), un(Δ[1]), NoTangent())
    end
end
function ChainRulesCore.rrule(::typeof(flatten), d::AbstractDict, ks)
    _d = OrderedDict(k => d[k] for k in ks)
    d_vec, un = flatten(_d, ks)
    return (d_vec, un), Δ -> begin
        (NoTangent(), un(Δ[1]), NoTangent())
    end
end

struct Unflatten{X, F} <: Function
    x::X
    unflatten::F
end
(f::Unflatten)(x) = f.unflatten(x)

_zero(x::Real) = zero(x)
_zero(x::AbstractArray) = _zero.(x)
_zero(x::AbstractDict) = Dict(keys(x) .=> map(_zero, values(x)))
_zero(x::NamedTuple) = map(_zero, x)
_zero(x::Tuple) = map(_zero, x)
_zero(x) = structfromnt(typeof(x), _zero(ntfromstruct(x)))

function _merge(d1, d2::AbstractDict)
    _d = OrderedDict(k => _zero(v) for (k, v) in d1)
    return sort!(merge(_d, OrderedDict(d2)))
end
_merge(::Any, d2) = d2

function ChainRulesCore.rrule(un::Unflatten, v)
    x = un(v)
    return x, Δ -> begin
        _Δ = _merge(x, Δ)
        return (NoTangent(), zygote_flatten(un.x, _Δ)[1])
    end
end

function flatten(::Nothing)
    return Float64[], _ -> nothing
end
function flatten(::NoTangent)
    return Float64[], _ -> NoTangent()
end
function flatten(::ZeroTangent)
    return Float64[], _ -> ZeroTangent()
end
function flatten(::Tuple{})
    return Float64[], _ -> ()
end

function zygote_flatten(x, ::Nothing)
    t = flatten(x)
    return zero(t[1]), Base.tail(t)
end
function zygote_flatten(x, ::NoTangent)
    t = flatten(x)
    return zero(t[1]), Base.tail(t)
end
function zygote_flatten(x, ::ZeroTangent)
    t = flatten(x)
    return zero(t[1]), Base.tail(t)
end
function zygote_flatten(::Any, ::Tuple{})
    return Float64[], _ -> ()
end

macro constructor(T)
    return flatten_expr(T, T)
end
macro constructor(T, C)
    return flatten_expr(T, C)
end
flatten_expr(T, C) = quote
    function flatten(x::$(esc(T)))
        v, un = flatten(ntfromstruct(x))
        return identity.(v), Unflatten(x, y -> structfromnt($(esc(C)), un(y)))
    end
    function zygote_flatten(x1::$(esc(T)), x2::$(esc(T)))
        v, un = zygote_flatten(ntfromstruct(x1), ntfromstruct(x2))
        return identity.(v), Unflatten(x2, y -> structfromnt($(esc(C)), un(y)))
    end
    _zero(x::$(esc(T))) = structfromnt($(esc(C)), _zero(ntfromstruct(x)))
end

_cumsum(x) = cumsum(x)
if VERSION < v"1.5"
    _cumsum(x::Tuple) = (_cumsum(collect(x))..., )
end


text:
const ROOT_GET = """
/:
  get:
    description: Testing swagger markdown!
    responses:
      '200':
        description: Returns a mysterious string.
"""

@testset "@swagger writes markdowns to files" begin
    tmp = SwaggerMarkdown.TMP
    @swagger """
    /:
      get:
        description: Testing swagger markdown!
        responses:
          '200':
            description: Returns a mysterious string.
    """

    @assert isdir(SwaggerMarkdown.TMP)
    content = ""
    for filename in readdir(tmp)
        ext = split(filename, ".")[end]
        if (ext == "yml")
            io = open(joinpath(tmp, filename), "r")
            content *= read(io, String)
            close(io)
        end
    end
    cleanup()

    @test ROOT_GET == content

end


text:
module Dhall

using Pkg.Artifacts: @artifact_str
using JSON: JSON
using FileIO: @format_str, File, Stream, FileIO

const PATHS = Dict(
        "dhall"         => joinpath(artifact"dhall", "bin", "dhall"),
        "dhall-to-json" => joinpath(artifact"dhall-json", "bin", "dhall-to-json"),
        "json-to-dhall" => joinpath(artifact"dhall-json", "bin", "json-to-dhall"),
        "dhall-to-yaml" => joinpath(artifact"dhall-json", "bin", "dhall-to-yaml"),
       )

function load(f::File{format"dhall"}; kw...)
    open(f) do stream
        load(stream; kw...)
    end
end
function load(st::Stream{format"dhall"}; kw...)
    io = FileIO.stream(st)
    load_via_json(io; kw...)
end

function path_from_io(io::IO)
    re = r"<file (.*)>"
    match(re, io.name)[1]
end

function load_via_json(io::IO; dir=nothing)
    if dir == nothing
        dir = dirname(path_from_io(io))
    end
    exe = Dhall.PATHS["dhall-to-json"]
    io_out = Pipe()
    p = pipeline(Cmd(`$exe`, dir=dir), stdin=io, stdout=io_out)
    run(p)
    JSON.parse(io_out)
end

function load_via_json(path::AbstractString; kw...)
    open(path) do io_in
        load_via_json(io_in; kw...)
    end
end

end#module


text:
"""
Summary statistics.
"""

using CSV
using DataFrames

function sumstats(fh; alleles = false, dropna = true)
    # Parses .sumstats files. See docs/file_formats_sumstats.txt
    #dtype_dict = {"SNP": str,   "Z": float, "N": float, "A1": str, "A2": str}
    compression = get_compression(fh)
    usecols = ["SNP", "Z", "N"]
    if alleles usecols += ["A1", "A2"] end

    try
        x = read_csv(fh, delim_whitespace = True, na_values = ".")
    catch e
        println("Improperly formatted sumstats file:", str(e.args))
    end

    if dropna
        x = x.dropna(how = "any")
    end

    return x
end

""" Estimate h² and partitioned h² """
function estimate_h2()
    args = deepcopy(LDScore.args)

    # TODO: more checking

    M_annot, w_ld_cname, ref_ld_cnames, sumstats, novar_cols = _read_ld_sumstats()

    ref_ld = sumstats[!, "LD"]
    ref_ld = reshape(ref_ld, (size(ref_ld)[1], 1))
    n_snp = size(sumstats)[1]
    n_blocks = min(n_snp, args["n-blocks"])
    n_annot = size(ref_ld_cnames)[1]
    χ²_max = args["χ²-max"]
    old_weights = false

    if n_annot == 1
        if args["two-step"] == nothing && args["intercept-h²"] == nothing
            args["two-step"] = 30
        end
    else
        old_weights = true
        if args["χ²-max"] == nothing
            χ²_max = max(0.001 * max(sumstats[!, "N"]), 80)
        end
    end

    s(x) = reshape(x, (n_snp, 1))
    χ² = s(sumstats[!, "Z"] .^ 2)

    if χ²_max != nothing
        ii = vec(χ² < χ²_max)
        sumstats = sumstats.ix[ii, :]

        n_snp = sum(ii) # TODO double check that this works
        ref_ld = [sumstats[ref_ld_cnames]]
        χ² = reshape(χ²[ii], (n_snp, 1))
    end

    ĥ² = Hsq(
        χ², ref_ld, s(sumstats[!, w_ld_cname]), s(sumstats[!, "N"]), M_annot,
        n_blocks = n_blocks, intercept = args["intercept-h²"],
        slow = false, step1_ii = args["two-step"],
        old_weights = old_weights,
    )

    if args["overlap-annot"]
        overlap_matrix, M_tot = _read_annot(args)

        df_results = ĥ²._overlap_output(
            ref_ld_cnames, overlap_matrix, M_annot, M_tot,
            args["print-coefficients"],
        )
        df_results.to_csv(
            args["out"] * ".results",
            sep = "\t", index = false,
        )
    end

    return ĥ²
end # function estimate_h2

# TODO
function _read_M(args, n_annot) end
function _check_variance(M_annot, ref_ld) end
function _read_w_ld(args) end
function _merge_and_log(ld, sumstats, noun) end

function _read_sumstats(args, fh; alleles=false, dropna=false)
    # Parser summary statistics
    @info "Reading summary statistics from" fh
    sumstats = parse_sumstats(fh, alleles=alleles, dropna=dropna)
    @info "Read summary statistics for" size(sumstats) "SNPs."
    m = size(sumstats)
    sumstats = sumstats.drop_duplicates(subset="SNP")
    if m > size(sumstats)
        @info "Dropped SNPs with duplicated rs numbers:" m-size(sumstats)
    end
    return sumstats
end

function _read_ref_ld()
    # Read reference LD Scores
    ref_ld = _read_chr_split_files(
        LDScore.args["ref-ld-chr"], LDScore.args["ref-ld"],
        "reference panel LD Score", ldscore_fromlist,
    )
    @info "Num. of SNPs read for reference panel LD Scores:" size(ref_ld)
    return ref_ld
end

function _load_testset_1()
    df_ref_id = DataFrames.DataFrame(CSV.File(LDScore.args["ref-ld"] * ".l2.ldscore"))
    df_h2 = DataFrames.DataFrame(CSV.File(LDScore.args["h²"]))

    M_annot = reshape([[155881.2526]], (1, 1))
    w_ld_cname = "CHR"
    ref_ld_cnames = ["LD_0"]
    sumstats = DataFrames.innerjoin(df_ref_id, df_h2, on = "SNP")
    novar_cols = Dict([("LD_0", false)])
    return (M_annot, w_ld_cname, ref_ld_cnames, sumstats, novar_cols)
end

#= TODO: Just loading testing set now,
    need to actually implement to load generic datasets =#
function _read_ld_sumstats(alleles = false, dropna = true)
    return _load_testset_1()
end


text:
# *********************************************************************************
# REopt, Copyright (c) 2019-2020, Alliance for Sustainable Energy, LLC.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification,
# are permitted provided that the following conditions are met:
#
# Redistributions of source code must retain the above copyright notice, this list
# of conditions and the following disclaimer.
#
# Redistributions in binary form must reproduce the above copyright notice, this
# list of conditions and the following disclaimer in the documentation and/or other
# materials provided with the distribution.
#
# Neither the name of the copyright holder nor the names of its contributors may be
# used to endorse or promote products derived from this software without specific
# prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
# IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
# OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
# OF THE POSSIBILITY OF SUCH DAMAGE.
# *********************************************************************************
using Xpress

@testset "Thermal loads" begin
    m = Model(optimizer_with_attributes(Xpress.Optimizer, "OUTPUTLOG" => 0))
    results = run_reopt(m, "./scenarios/thermal_load.json")

    @test round(results["ExistingBoiler"]["year_one_boiler_fuel_consumption_mmbtu"], digits=0) ≈ 2905

    data = JSON.parsefile("./scenarios/thermal_load.json")
    data["DomesticHotWaterLoad"]["fuel_loads_mmbtu_per_hour"] = repeat([0.5], 8760)
    data["SpaceHeatingLoad"]["fuel_loads_mmbtu_per_hour"] = repeat([0.5], 8760)
    s = Scenario(data)
    inputs = REoptInputs(s)
    m = Model(optimizer_with_attributes(Xpress.Optimizer, "OUTPUTLOG" => 0))
    results = run_reopt(m, inputs)

    @test round(results["ExistingBoiler"]["year_one_boiler_fuel_consumption_mmbtu"], digits=0) ≈ 8760
end
#=
add a time-of-export rate that is greater than retail rate for the month of January,
check to make sure that PV does NOT export unless the site load is met first for the month of January.
=#
@testset "Do not allow_simultaneous_export_import" begin
    model = Model(optimizer_with_attributes(Xpress.Optimizer, "OUTPUTLOG" => 0))
    data = JSON.parsefile("./scenarios/monthly_rate.json")

    # create wholesale_rate with compensation in January > retail rate
    jan_rate = data["ElectricTariff"]["monthly_energy_rates"][1]
    data["ElectricTariff"]["wholesale_rate"] =
        append!(repeat([jan_rate + 0.1], 31 * 24), repeat([0.0], 8760 - 31*24))
    data["ElectricTariff"]["monthly_demand_rates"] = repeat([0], 12)
    data["ElectricUtility"] = Dict("allow_simultaneous_export_import" => false)

    s = Scenario(data)
    inputs = REoptInputs(s)
    results = run_reopt(model, inputs)

    @test all(x == 0.0 for (i,x) in enumerate(results["ElectricUtility"]["year_one_to_load_series_kw"][1:744])
              if results["PV"]["year_one_to_grid_series_kw"][i] > 0)
end

@testset "Solar and Storage w/BAU" begin
    m1 = Model(optimizer_with_attributes(Xpress.Optimizer, "OUTPUTLOG" => 0))
    m2 = Model(optimizer_with_attributes(Xpress.Optimizer, "OUTPUTLOG" => 0))
    results = run_reopt([m1,m2], "./scenarios/pv_storage.json")

    @test results["PV"]["size_kw"] ≈ 216.6667 atol=0.01
    @test results["PV"]["lcoe_per_kwh"] ≈ 0.0483 atol = 0.001
    @test results["Financial"]["lcc"] ≈ 1.240037e7 rtol=1e-5
    @test results["Financial"]["lcc_bau"] ≈ 12766397 rtol=1e-5
    @test results["Storage"]["size_kw"] ≈ 55.9 atol=0.1
    @test results["Storage"]["size_kwh"] ≈ 78.9 atol=0.1
    proforma_npv = REoptLite.npv(results["Financial"]["offtaker_annual_free_cashflows"] -
        results["Financial"]["offtaker_annual_free_cashflows_bau"], 0.081)
    @test results["Financial"]["npv"] ≈ proforma_npv rtol=0.0001
end

@testset "Outage with Generator, outate simulator, BAU critical load outputs" begin
    m1 = Model(optimizer_with_attributes(Xpress.Optimizer, "OUTPUTLOG" => 0))
    m2 = Model(optimizer_with_attributes(Xpress.Optimizer, "OUTPUTLOG" => 0))
    p = REoptInputs("./scenarios/generator.json")
    results = run_reopt([m1,m2], p)
    @test results["Generator"]["size_kw"] ≈ 8.13 atol=0.01
    @test (sum(results["Generator"]["year_one_to_load_series_kw"][i] for i in 1:9) +
           sum(results["Generator"]["year_one_to_load_series_kw"][i] for i in 13:8760)) == 0
    @test results["ElectricLoad"]["bau_critical_load_met"] == false
    @test results["ElectricLoad"]["bau_critical_load_met_time_steps"] == 0

    simresults = simulate_outages(results, p)
    @test simresults["resilience_hours_max"] == 11
end

@testset "Minimize Unserved Load" begin
    m = Model(optimizer_with_attributes(Xpress.Optimizer, "OUTPUTLOG" => 0))
    results = run_reopt(m, "./scenarios/outage.json")

    @test results["expected_outage_cost"] ≈ 0
    @test sum(results["unserved_load_per_outage"]) ≈ 0
    @test value(m[:binMGTechUsed]["Generator"]) == 1
    @test value(m[:binMGTechUsed]["PV"]) == 0
    @test value(m[:binMGStorageUsed]) == 1
    @test results["Financial"]["lcc"] ≈ 7.3879557e7 atol=5e4

    #=
    Scenario with $0/kWh value_of_lost_load_per_kwh, 12x169 hour outages, 1kW load/hour, and min_resil_timesteps = 168
    - should meet 168 kWh in each outage such that the total unserved load is 12 kWh
    =#
    m = Model(optimizer_with_attributes(Xpress.Optimizer, "OUTPUTLOG" => 0))
    results = run_reopt(m, "./scenarios/nogridcost_minresilhours.json")
    @test sum(results["unserved_load_per_outage"]) ≈ 12

    # testing dvUnserved load, which would output 100 kWh for this scenario before output fix
    m = Model(optimizer_with_attributes(Xpress.Optimizer, "OUTPUTLOG" => 0))
    results = run_reopt(m, "./scenarios/nogridcost_multiscenario.json")
    @test sum(results["unserved_load_per_outage"]) ≈ 60

end

@testset "Multiple Sites" begin
    m = Model(optimizer_with_attributes(Xpress.Optimizer, "OUTPUTLOG" => 0))
    ps = [
        REoptInputs("./scenarios/pv_storage.json"),
        REoptInputs("./scenarios/monthly_rate.json"),
    ];
    results = run_reopt(m, ps)
    @test results[3]["Financial"]["lcc"] + results[10]["Financial"]["lcc"] ≈ 1.240037e7 + 437169.0 rtol=1e-5
end

@testset "MPC" begin
    model = Model(optimizer_with_attributes(Xpress.Optimizer, "OUTPUTLOG" => 0))
    r = run_mpc(model, "./scenarios/mpc.json")
    @test maximum(r["ElectricUtility"]["to_load_series_kw"][1:15]) <= 98.0
    @test maximum(r["ElectricUtility"]["to_load_series_kw"][16:24]) <= 97.0
    @test sum(r["PV"]["to_grid_series_kw"]) ≈ 0
end

@testset "Complex Incentives" begin
    """
    This test was compared against the API test:
        reo.tests.test_reopt_url.EntryResourceTest.test_complex_incentives
    when using the hardcoded levelization_factor in this package's REoptInputs function.
    The two LCC's matched within 0.00005%. (The Julia pkg LCC is  1.0971991e7)
    """
    m = Model(optimizer_with_attributes(Xpress.Optimizer, "OUTPUTLOG" => 0))
    results = run_reopt(m, "./scenarios/incentives.json")
    @test results["Financial"]["lcc"] ≈ 1.0968526e7 atol=5e4
end

@testset verbose = true "Rate Structures" begin

    @testset "Tiered Energy" begin
        m = Model(optimizer_with_attributes(Xpress.Optimizer, "OUTPUTLOG" => 0))
        results = run_reopt(m, "./scenarios/tiered_rate.json")
        @test results["ElectricTariff"]["year_one_energy_cost"] ≈ 2342.88
    end

    @testset "Lookback Demand Charges" begin
        m = Model(optimizer_with_attributes(Xpress.Optimizer, "OUTPUTLOG" => 0))
        results = run_reopt(m, "./scenarios/lookback_rate.json")
        @test results["ElectricTariff"]["year_one_demand_cost"] ≈ 721.99
    end

    @testset "Blended tariff" begin
        model = Model(optimizer_with_attributes(Xpress.Optimizer, "OUTPUTLOG" => 0))
        results = run_reopt(model, "./scenarios/no_techs.json")
        @test results["ElectricTariff"]["year_one_energy_cost"] ≈ 1000.0
        @test results["ElectricTariff"]["year_one_demand_cost"] ≈ 136.99
    end

    @testset "Coincident Peak Charges" begin
        model = Model(optimizer_with_attributes(Xpress.Optimizer, "OUTPUTLOG" => 0))
        results = run_reopt(model, "./scenarios/coincident_peak.json")
        @test results["ElectricTariff"]["year_one_coincident_peak_cost"] ≈ 15.0
        @test results["ElectricTariff"]["lifecycle_coincident_peak_cost"] ≈ 15.0 * 12.94887 atol=0.1
    end

    # # tiered monthly demand rate  TODO: expected results?
    # m = Model(optimizer_with_attributes(Xpress.Optimizer, "OUTPUTLOG" => 0))
    # data = JSON.parsefile("./scenarios/tiered_rate.json")
    # data["ElectricTariff"]["urdb_label"] = "59bc22705457a3372642da67"
    # s = Scenario(data)
    # inputs = REoptInputs(s)
    # results = run_reopt(m, inputs)

    # TODO test for tiered TOU demand rates
end

@testset "Wind" begin
    m = Model(optimizer_with_attributes(Xpress.Optimizer, "OUTPUTLOG" => 0))
    results = run_reopt(m, "./scenarios/wind.json")
    @test results["Wind"]["size_kw"] ≈ 3752 atol=0.1
    @test results["Financial"]["lcc"] ≈ 8.591017e6 rtol=1e-5
    #=
    0.5% higher LCC in this package as compared to API ? 8,591,017 vs 8,551,172
    - both have zero curtailment
    - same energy to grid: 5,839,317 vs 5,839,322
    - same energy to load: 4,160,683 vs 4,160,677
    - same city: Boulder
    - same total wind prod factor

    REoptLite.jl has:
    - bigger turbine: 3752 vs 3735
    - net_capital_costs_plus_om: 8,576,590 vs. 8,537,480

    TODO: will these discrepancies be addressed once NMIL binaries are added?
    =#
end

@testset "Multiple PVs" begin
    m1 = Model(optimizer_with_attributes(Xpress.Optimizer, "OUTPUTLOG" => 0))
    m2 = Model(optimizer_with_attributes(Xpress.Optimizer, "OUTPUTLOG" => 0))
    results = run_reopt([m1,m2], "./scenarios/multiple_pvs.json")

    ground_pv = results["PV"][findfirst(pv -> pv["name"] == "ground", results["PV"])]
    roof_west = results["PV"][findfirst(pv -> pv["name"] == "roof_west", results["PV"])]
    roof_east = results["PV"][findfirst(pv -> pv["name"] == "roof_east", results["PV"])]

    @test ground_pv["size_kw"] ≈ 15 atol=0.1
    @test roof_west["size_kw"] ≈ 7 atol=0.1
    @test roof_east["size_kw"] ≈ 4 atol=0.1
    @test ground_pv["lifecycle_om_cost_bau"] ≈ 782.0 atol=0.1
    @test roof_west["lifecycle_om_cost_bau"] ≈ 782.0 atol=0.1
    @test ground_pv["average_annual_energy_produced_kwh_bau"] ≈ 8844.19 atol=0.1
    @test roof_west["average_annual_energy_produced_kwh_bau"] ≈ 7440.1 atol=0.1
    @test ground_pv["average_annual_energy_produced_kwh"] ≈ 26533.54 atol=0.1
    @test roof_west["average_annual_energy_produced_kwh"] ≈ 10416.52 atol=0.1
    @test roof_east["average_annual_energy_produced_kwh"] ≈ 6482.37 atol=0.1
end

## equivalent REopt Lite API Post for test 2:
#   NOTE have to hack in API levelization_factor to get LCC within 5e-5 (Mosel tol)
# {"Scenario": {
#     "Site": {
#         "longitude": -118.1164613,
#         "latitude": 34.5794343,
#         "roof_squarefeet": 5000.0,
#         "land_acres": 1.0,
#     "PV": {
#         "macrs_bonus_pct": 0.4,
#         "installed_cost_per_kw": 2000.0,
#         "tilt": 34.579,
#         "degradation_pct": 0.005,
#         "macrs_option_years": 5,
#         "federal_itc_pct": 0.3,
#         "module_type": 0,
#         "array_type": 1,
#         "om_cost_per_kw": 16.0,
#         "macrs_itc_reduction": 0.5,
#         "azimuth": 180.0,
#         "federal_rebate_per_kw": 350.0,
#         "dc_ac_ratio": 1.1
#     },
#     "LoadProfile": {
#         "doe_reference_name": "RetailStore",
#         "annual_kwh": 10000000.0,
#         "city": "LosAngeles"
#     },
#     "Storage": {
#         "total_rebate_per_kw": 100.0,
#         "macrs_option_years": 5,
#         "can_grid_charge": true,
#         "macrs_bonus_pct": 0.4,
#         "macrs_itc_reduction": 0.5,
#         "total_itc_pct": 0,
#         "installed_cost_per_kw": 1000.0,
#         "installed_cost_per_kwh": 500.0,
#         "replace_cost_per_kw": 460.0,
#         "replace_cost_per_kwh": 230.0
#     },
#     "ElectricTariff": {
#         "urdb_label": "5ed6c1a15457a3367add15ae"
#     },
#     "Financial": {
#         "escalation_pct": 0.026,
#         "offtaker_discount_pct": 0.081,
#         "owner_discount_pct": 0.081,
#         "analysis_years": 20,
#         "offtaker_tax_pct": 0.4,
#         "owner_tax_pct": 0.4,
#         "om_cost_escalation_pct": 0.025
#     }
# }}}

max_stars_repo_name: NREL/REoptLite

max_issues_repo_name: NREL/REoptLite

max_forks_repo_name: NREL/REoptLite


text:
import Pkg

function get_version(filename::AbstractString)::VersionNumber
    project = Pkg.TOML.parsefile(filename)
    return VersionNumber(project["version"])
end

function get_version()::VersionNumber
    filename = joinpath(pwd(), "Project.toml")
    return get_version(filename)
end

function get_upstream_version(repo_url::AbstractString;
                              master_branch::AbstractString)::VersionNumber
    original_directory = pwd()
    result = with_cloned_repo(repo_url) do
        run(`git checkout $(master_branch)`)
        result = get_version()
    end
    cd(original_directory)
    return result
end

function get_head_version(repo_url::AbstractString;
                          build_directory::AbstractString)::VersionNumber
    original_directory = pwd()
    cd(build_directory)
    result = get_version()
    cd(original_directory)
    return result
end