text:
#############################################################################
#Copyright (c) 2016-2019 Charles Le Losq
#
#The MIT License (MIT)
#
#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.
#
# This files contains functions to perform calculations of numeric integrale values on real dataset x y.
#
#
#############################################################################

"""

	trapz(x::Vector{Tx}, y::Vector{Ty}) where {Tx <: Number, Ty <: Number}

Trapezoidal integration. This function is particularly helpful to calculate the area under a portion of a spectrum, and can be used for various purposes (normalisation, area comparison, etc.).

Inputs
-----
	x: Vector{Float64}
		x values
	y: Vector{Float64}
		y values.

Outputs
-------
	area: Vector{Float64}
		trapezoidal integration value.

"""
function trapz(x::Vector{Tx}, y::Vector{Ty}) where {Tx <: Number, Ty <: Number}
    # Trapezoidal integration rule
    local n = length(x)
    if (length(y) != n)
        error("Vectors 'x', 'y' must be of same length")
    end
    r = zero(zero(Tx) + zero(Ty))
    if n == 1; return r; end
    for i in 2:n
        r += (x[i] - x[i-1]) * (y[i] + y[i-1])
    end
    trapz_int = r/2
    return trapz_int
end

"""

	bandarea(Amplitude::Array{Float64},HWHM::Array{Float64}; peak_shape = "Gaussian", error_switch = "no", eseAmplitude::Array{Float64} = [0.0], eseHWHM::Array{Float64} = [0.0])

This function allows calculating the area under a specific band, with different shapes. For now, only Gaussian bands are supported, but other band shapes will be added soon.

Inputs
------

	Amplitude: Array{Float64}
		peak amplitude
	HWHM: Array{Float64}
		peak half width at half maximum

Options
-------

	peak_shape: String
		shape of the peak. Only "Gaussian" is supported for now
	error_switch: String
		should be "yes" or "no". If "yes", the arrays containing the errors affecting the band amplitude and widhts should be provided in eseAmplitude and eseHWHM (see below);
	eseAmplitude: Array{Float64}
		array containing the errors affecting Amplitude
	eseHWHM: Array{Float64}
		array containing the errors affecting HWHM;

Outputs
-------
	area: Array{Float64}
		array containing peak areas

	if error_switch is set to "yes", then a second output is provided:

	esearea: Array{Float64}
		array that contains the propagated errors affecting the areas calculations.
"""
function bandarea(Amplitude::Array{Float64},HWHM::Array{Float64}; peak_shape = "Gaussian", error_switch = "no", eseAmplitude::Array{Float64} = [0.0], eseHWHM::Array{Float64} = [0.0])

	# first we check the desired peak shape, and apply the relevant calculation
	if peak_shape == "Gaussian"
    	area::Array{Float64} = sqrt(pi./log(2)).*Amplitude.*HWHM # Gaussian area, HWHM is the half-width at half-maximum
		if error_switch == "yes" # if errors are desired, perform the error calculation
			if size(eseAmplitude) != size(Amplitude) || size(eseHWHM) != size(HWHM) # error check
				error("Please check that you entered good arrays for the errors on the amplitude and widths of the bands")
			end
	        esearea::Array{Float64} = sqrt((pi./log(2).*HWHM).^2 .* eseAmplitude.^2 + (pi./log(2).*Amplitude).^2 .* eseHWHM.^2)
	    end
	else
		error("Not yet implemented.")
	end

	# Depending on the error switch, we output only the areas or also the errors
    if error_switch == "no"
		# a quick test to see if the output is an array or should be a Float64 number
		if size(area) == (1,)
        	return area[1] # we return a Float64 number
		else
			return area # we return an array
		end
	elseif error_switch =="yes"
		# a quick test to see if the output is an array or should be 2 Float64 numbers
		if size(area) == (1,)
        	return area[1], esearea[1] # we return two Float64 numbers
		else
			return area, esearea # we return an array
		end
	else
		error("error_switch should be equal to yes or no. Please select the appropriate value.")
	end

end

ext: jl


text:
"""
    `@df d x`

Convert every symbol in the expression `x` with the respective column in `d` if it exists.

If you want to avoid replacing the symbol, escape it with `^`.

`NA` values are replaced with `NaN` for columns of `Float64` and `""` or `Symbol()`
for strings and symbols respectively.

`x` can be either a plot command or a block of plot commands.
"""
macro df(d, x)
    esc(Expr(:call, df_helper(x), d))
end

"""
    `@df x`

Curried version of `@df d x`. Outputs an anonymous function `d -> @df d x`.
"""
macro df(x)
    esc(df_helper(x))
end

function df_helper(x)
    i = gensym()
    Expr(:(->), i, df_helper(i, x))
end

function df_helper(d, x)
    if isa(x, Expr) && x.head == :block # meaning that there were multiple plot commands
        commands = [df_helper(d, xx) for xx in x.args if !(isa(xx, Expr) && xx.head == :line || isa(xx, LineNumberNode))] # apply the helper recursively
        return Expr(:block, commands...)

    elseif isa(x, Expr) && x.head == :call # each function call is operated on alone
        syms = Any[]
        vars = Symbol[]
        plot_call = parse_table_call!(d, x, syms, vars)
        names = gensym()
        compute_vars = Expr(:(=), Expr(:tuple, Expr(:tuple, vars...), names),
            Expr(:call, :(StatsPlots.extract_columns_and_names), d, syms...))
        argnames = _argnames(names, x)
        if (length(plot_call.args) >= 2) && isa(plot_call.args[2], Expr) && (plot_call.args[2].head == :parameters)
            label_plot_call = Expr(:call, :(StatsPlots.add_label), plot_call.args[2], argnames,
                plot_call.args[1], plot_call.args[3:end]...)
        else
            label_plot_call = Expr(:call, :(StatsPlots.add_label), argnames, plot_call.args...)
        end
        return Expr(:block, compute_vars, label_plot_call)

    else
        error("Second argument ($x) can only be a block or function call")
    end
end

parse_table_call!(d, x, syms, vars) = x

function parse_table_call!(d, x::QuoteNode, syms, vars)
    new_var = gensym(x.value)
    push!(syms, x)
    push!(vars, new_var)
    return new_var
end

function parse_table_call!(d, x::Expr, syms, vars)
    if x.head == :. && length(x.args) == 2
        isa(x.args[2], QuoteNode) && return x
    elseif x.head == :call
        x.args[1] == :^ && length(x.args) == 2 && return x.args[2]
        if x.args[1] == :cols
            if length(x.args) == 1
                push!(x.args, :(StatsPlots.column_names(StatsPlots.getiterator($d))))
                return parse_table_call!(d, x, syms, vars)
            end
            range = x.args[2]
            new_vars = gensym("range")
            push!(syms, range)
            push!(vars, new_vars)
            return new_vars
        end
    elseif x.head==:braces # From Query: use curly brackets to simplify writing named tuples
        new_ex = Expr(:tuple, x.args...)

        for (j,field_in_NT) in enumerate(new_ex.args)
            if isa(field_in_NT, Expr) && field_in_NT.head==:(=)
                new_ex.args[j] = Expr(:(=), field_in_NT.args...)
            elseif field_in_NT isa QuoteNode
                new_ex.args[j] = Expr(:(=), field_in_NT.value, field_in_NT)
            elseif isa(field_in_NT, Expr)
                new_ex.args[j] = Expr(:(=), Symbol(filter(t -> t != ':', string(field_in_NT))), field_in_NT)
            elseif isa(field_in_NT, Symbol)
               new_ex.args[j] = Expr(:(=), field_in_NT, field_in_NT)
            end
        end
        return parse_table_call!(d, new_ex, syms, vars)
    end
    return Expr(x.head, (parse_table_call!(d, arg, syms, vars) for arg in x.args)...)
end

function column_names(t)
    s = schema(t)
    s === nothing ? propertynames(first(rows(t))) : s.names
end

not_kw(x) = true
not_kw(x::Expr) = !(x.head in [:kw, :parameters])

function insert_kw!(x::Expr, s::Symbol, v)
    index = isa(x.args[2], Expr) && x.args[2].head == :parameters ? 3 : 2
    x.args = vcat(x.args[1:index-1], Expr(:kw, s, v), x.args[index:end])
end

function _argnames(names, x::Expr)
    Expr(:vect, [_arg2string(names, s) for s in x.args[2:end] if not_kw(s)]...)
end

_arg2string(names, x) = stringify(x)
function _arg2string(names, x::Expr)
    if x.head == :call && x.args[1] == :cols
        return :(StatsPlots.compute_name($names, $(x.args[2])))
    elseif x.head == :call && x.args[1] == :hcat
        return hcat(stringify.(x.args[2:end])...)
    elseif x.head == :hcat
        return hcat(stringify.(x.args)...)
    else
        return stringify(x)
    end
end

stringify(x) = filter(t -> t != ':', string(x))

compute_name(names, i::Int) = names[i]
compute_name(names, i::Symbol) = i
compute_name(names, i) = reshape([compute_name(names, ii) for ii in i], 1, :)

# This function ensures that labels are passed to the plottig command.
function add_label(argnames, f, args...; kwargs...)
    i = findlast(t -> isa(t, Expr) || isa(t, AbstractArray), argnames)
    if (i === nothing)
        return f(args...; kwargs...)
    else
        return f(label = argnames[i], args...; kwargs...)
    end
end

get_col(s::Int, col_nt, names) = col_nt[names[s]]
get_col(s::Symbol, col_nt, names) = get(col_nt, s, s)
get_col(syms, col_nt, names) = hcat((get_col(s, col_nt, names) for s in syms)...)

"""
    extract_columns_and_names(df, syms...)

Extracts columns and their names (if the column number is an integer)
into a slightly complex `Tuple`.

The structure goes as `((columndata...), names)`.  This is unpacked by the [`@df`](@ref) macro into `gensym`'ed variables, which are passed to the plotting function.

!!! note
    If you want to extend the [`@df`](@ref) macro
    to work with your custom type, this is the
    function you should overload!
"""
function extract_columns_and_names(df, syms...)
    istable(df) || error("Only tables are supported")
    names = column_names(df)

    selected_cols = Symbol[]

    # get the appropriate name when passed an Integer
    add_sym!(s::Integer) = push!(selected_cols, names[s])
    # check for errors in Symbols
    add_sym!(s::Symbol) = s in names && push!(selected_cols, s)
    # recursively extract column names
    add_sym!(s) = foreach(add_sym!, s)

    foreach(add_sym!, syms)

    cols = columntable(select(df, unique(selected_cols)...))
    return Tuple(get_col(s, cols, names) for s in syms), names
end

convert_missing(el) = el
convert_missing(el::DataValue{T}) where {T} = get(el, missing)
convert_missing(el::DataValue{<:AbstractString}) = get(el, "")
convert_missing(el::DataValue{Symbol}) = get(el, Symbol())
convert_missing(el::DataValue{<:Real}) = get(convert(DataValue{Float64}, el), NaN)


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

@doc raw"""HorizontalPodAutoscalerCondition describes the state of a HorizontalPodAutoscaler at a certain point.

    IoK8sApiAutoscalingV2beta2HorizontalPodAutoscalerCondition(;
        lastTransitionTime=nothing,
        message=nothing,
        reason=nothing,
        status=nothing,
        type=nothing,
    )

    - lastTransitionTime::IoK8sApimachineryPkgApisMetaV1Time : lastTransitionTime is the last time the condition transitioned from one status to another
    - message::String : message is a human-readable explanation containing details about the transition
    - reason::String : reason is the reason for the condition's last transition.
    - status::String : status is the status of the condition (True, False, Unknown)
    - type::String : type describes the current condition
"""
mutable struct IoK8sApiAutoscalingV2beta2HorizontalPodAutoscalerCondition <: SwaggerModel
    lastTransitionTime::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1Time } # spec name: lastTransitionTime
    message::Any # spec type: Union{ Nothing, String } # spec name: message
    reason::Any # spec type: Union{ Nothing, String } # spec name: reason
    status::Any # spec type: Union{ Nothing, String } # spec name: status
    type::Any # spec type: Union{ Nothing, String } # spec name: type

    function IoK8sApiAutoscalingV2beta2HorizontalPodAutoscalerCondition(;lastTransitionTime=nothing, message=nothing, reason=nothing, status=nothing, type=nothing)
        o = new()
        validate_property(IoK8sApiAutoscalingV2beta2HorizontalPodAutoscalerCondition, Symbol("lastTransitionTime"), lastTransitionTime)
        setfield!(o, Symbol("lastTransitionTime"), lastTransitionTime)
        validate_property(IoK8sApiAutoscalingV2beta2HorizontalPodAutoscalerCondition, Symbol("message"), message)
        setfield!(o, Symbol("message"), message)
        validate_property(IoK8sApiAutoscalingV2beta2HorizontalPodAutoscalerCondition, Symbol("reason"), reason)
        setfield!(o, Symbol("reason"), reason)
        validate_property(IoK8sApiAutoscalingV2beta2HorizontalPodAutoscalerCondition, Symbol("status"), status)
        setfield!(o, Symbol("status"), status)
        validate_property(IoK8sApiAutoscalingV2beta2HorizontalPodAutoscalerCondition, Symbol("type"), type)
        setfield!(o, Symbol("type"), type)
        o
    end
end # type IoK8sApiAutoscalingV2beta2HorizontalPodAutoscalerCondition

const _property_map_IoK8sApiAutoscalingV2beta2HorizontalPodAutoscalerCondition = Dict{Symbol,Symbol}(Symbol("lastTransitionTime")=>Symbol("lastTransitionTime"), Symbol("message")=>Symbol("message"), Symbol("reason")=>Symbol("reason"), Symbol("status")=>Symbol("status"), Symbol("type")=>Symbol("type"))
const _property_types_IoK8sApiAutoscalingV2beta2HorizontalPodAutoscalerCondition = Dict{Symbol,String}(Symbol("lastTransitionTime")=>"IoK8sApimachineryPkgApisMetaV1Time", Symbol("message")=>"String", Symbol("reason")=>"String", Symbol("status")=>"String", Symbol("type")=>"String")
Base.propertynames(::Type{ IoK8sApiAutoscalingV2beta2HorizontalPodAutoscalerCondition }) = collect(keys(_property_map_IoK8sApiAutoscalingV2beta2HorizontalPodAutoscalerCondition))
Swagger.property_type(::Type{ IoK8sApiAutoscalingV2beta2HorizontalPodAutoscalerCondition }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiAutoscalingV2beta2HorizontalPodAutoscalerCondition[name]))}
Swagger.field_name(::Type{ IoK8sApiAutoscalingV2beta2HorizontalPodAutoscalerCondition }, property_name::Symbol) =  _property_map_IoK8sApiAutoscalingV2beta2HorizontalPodAutoscalerCondition[property_name]

function check_required(o::IoK8sApiAutoscalingV2beta2HorizontalPodAutoscalerCondition)
    (getproperty(o, Symbol("status")) === nothing) && (return false)
    (getproperty(o, Symbol("type")) === nothing) && (return false)
    true
end

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


text:
import Term: inspect, expressiontree, typestree, Dendogram, Tree
import Term.console: console_width

println("\nTesting logging, stdout temporarily disabled")
# @suppress_out begin
@testset "\e[34mINSPECT test" begin

    # define some types
    abstract type T1 end

    abstract type T2 <: T1 end

    """
        MyType

    It's just a useless type we've defined to provide an example of
    Term's awesome `inspect` functionality!
    """
    struct MyType <: T2
        x::Int
        y::String
        fn::Any
    end

    #  Now, what is MyType like?
    tofile(cleansprint(inspect, MyType), "./txtfiles/inspect_01.txt")
    if dotest
        @test cleansprint(inspect, MyType) == fromfile("./txtfiles/inspect_01.txt")
    else
        @test_nothrow inspect(MyType)
    end

    # Let's define some constructors and methods using MyType

    """
    constructors!
    """
    MyType(x::Int) = MyType(x, "no string", nothing)

    MyType(x::Int, y) = MyType(x, y, nothing)

    # methods
    useless_method(m::MyType) = m
    another_method(m1::MyType, m2::MyType) = print(m1, m2)

    tofile(cleansprint(inspect, MyType), "./txtfiles/inspect_02.txt")
    tofile(cleansprint(inspect, 1), "./txtfiles/inspect_03.txt")
    tofile(cleansprint(inspect, T2), "./txtfiles/inspect_04.txt")
    tofile(cleansprint(inspect, inspect), "./txtfiles/inspect_05.txt")

    if dotest
        @test cleansprint(inspect, MyType) == fromfile("./txtfiles/inspect_02.txt")
        @test cleansprint(inspect, 1) == fromfile("./txtfiles/inspect_03.txt")
        @test cleansprint(inspect, T2) == fromfile("./txtfiles/inspect_04.txt")
        @test cleansprint(inspect, inspect) == fromfile("./txtfiles/inspect_05.txt")
    else
        @test_nothrow inspect(MyType)
        @test_nothrow inspect(1)
        @test_nothrow inspect(T2)
        @test_nothrow inspect(inspect)
    end
end
# end

# define expressions
e1 = :(2x + 3y + 2)
e2 = :(2x + 3 + 2 + 2y)
e3 = :(2x^(3+y))
e4 = :(1 + 1 - 2x^2)
e5 = :(mod(22, 6))
e6 = :(2x^(3+y) + 2z)
e7 = :(2x + 3 * √(3x^2))
e8 = :(print(lstrip("test")))
expressions = (e1, e2, e3, e4, e5, e6, e7, e8)

# save expressions to file (for later comparisons)
# for (i, e) in enumerate(expressions)
#     tofile(string(Dendogram(e)), "./txtfiles/dendo_expr_$i.txt")
#     tofile(string(Tree(e)), "./txtfiles/tree_expr_$i.txt")
#     tofile(sprint(expressiontree, e), "./txtfiles/exptree_expr_$i.txt")
# end

@testset "Inspect: expressions" begin
    if !Sys.iswindows()
        # dendogram
        for (i, e) in enumerate(expressions)
            dendo = Dendogram(e)
            tree = Tree(e)

            @test dendo isa Dendogram
            @test fromfile("./txtfiles/dendo_expr_$i.txt") == cleanstring(dendo)

            @test tree isa Tree
            @test fromfile("./txtfiles/tree_expr_$i.txt") == cleanstring(tree)

            @test fromfile("./txtfiles/exptree_expr_$i.txt") == cleansprint(expressiontree, e)
        end
    end
end

@testset "Inspect: typestree" begin
    if console_width() >= 88
        @test cleansprint(typestree, Float64) == "\e[2m\e[34m╭───────────── \e[22m\e[38;2;255;167;38mTypes hierarchy\e[22m\e[39m\e[34m\e[2m\e[34m ───╮\e[22m\e[39m\e[0m\e[22m\e[39m\e[34m\n\e[2m\e[34m│\e[22m\e[39m                                 \e[0m\e[2m\e[34m│\e[22m\e[39m\e[0m\n\e[2m\e[34m│\e[22m\e[39m  \e[3m\e[38;5;10mAny\e[23m\e[39m                            \e[0m\e[2m\e[34m│\e[22m\e[39m\e[0m\n\e[2m\e[34m│\e[22m\e[39m \e[2m\e[3m\e[38;5;10m━━━━━\e[22m\e[23m\e[39m                           \e[0m\e[2m\e[34m│\e[22m\e[39m\e[0m\n\e[2m\e[34m│\e[22m\e[39m   \e[22m\e[2m\e[32m│   \e[22m\e[39m\e[22m                          \e[0m\e[2m\e[34m│\e[22m\e[39m\e[0m\n\e[2m\e[34m│\e[22m\e[39m   \e[22m\e[2m\e[32m├── \e[22m\e[39m\e[3m\e[38;2;255;238;88mReal\e[23m\e[39m\e[22m                      \e[0m\e[2m\e[34m│\e[22m\e[39m\e[0m\n\e[2m\e[34m│\e[22m\e[39m   \e[22m\e[2m\e[32m│   \e[22m\e[39m\e[2m\e[32m├── \e[22m\e[39m\e[3m\e[38;2;255;238;88mAbstractFloat\e[23m\e[39m\e[22m         \e[0m\e[2m\e[34m│\e[22m\e[39m\e[0m\n\e[2m\e[34m│\e[22m\e[39m   \e[22m\e[2m\e[32m│   \e[22m\e[39m\e[2m\e[32m│   \e[22m\e[39m\e[2m\e[32m├── \e[22m\e[39m\e[3m\e[38;2;255;238;88m\e[3m\e[4m\e[38;5;214mFloat64\e[23m\e[24m\e[39m\e[38;2;255;238;88m\e[23m\e[39m\e[22m           \e[0m\e[2m\e[34m│\e[22m\e[39m\e[0m\n\e[2m\e[34m│\e[22m\e[39m   \e[22m\e[2m\e[32m│   \e[22m\e[39m\e[2m\e[32m│   \e[22m\e[39m\e[2m\e[32m├── \e[22m\e[39mBigFloat\e[22m          \e[0m\e[2m\e[34m│\e[22m\e[39m\e[0m\n\e[2m\e[34m│\e[22m\e[39m   \e[22m\e[2m\e[32m│   \e[22m\e[39m\e[2m\e[32m│   \e[22m\e[39m\e[2m\e[32m├── \e[22m\e[39mFloat32\e[22m           \e[0m\e[2m\e[34m│\e[22m\e[39m\e[0m\n\e[2m\e[34m│\e[22m\e[39m   \e[22m\e[2m\e[32m│   \e[22m\e[39m\e[2m\e[32m│   \e[22m\e[39m\e[2m\e[32m└── \e[22m\e[39mFloat16\e[22m           \e[0m\e[2m\e[34m│\e[22m\e[39m\e[0m\n\e[2m\e[34m│\e[22m\e[39m   \e[22m\e[2m\e[32m│   \e[22m\e[39m\e[2m\e[32m│   \e[22m\e[39m\e[22m                      \e[0m\e[2m\e[34m│\e[22m\e[39m\e[0m\n\e[2m\e[34m│\e[22m\e[39m   \e[22m\e[2m\e[32m│   \e[22m\e[39m\e[2m\e[32m├── \e[22m\e[39mRational\e[22m              \e[0m\e[2m\e[34m│\e[22m\e[39m\e[0m\n\e[2m\e[34m│\e[22m\e[39m   \e[22m\e[2m\e[32m│   \e[22m\e[39m\e[2m\e[32m├── \e[22m\e[39mAbstractIrrational\e[22m    \e[0m\e[2m\e[34m│\e[22m\e[39m\e[0m\n\e[2m\e[34m│\e[22m\e[39m   \e[22m\e[2m\e[32m│   \e[22m\e[39m\e[2m\e[32m└── \e[22m\e[39mInteger\e[22m               \e[0m\e[2m\e[34m│\e[22m\e[39m\e[0m\n\e[2m\e[34m│\e[22m\e[39m   \e[22m\e[2m\e[32m│   \e[22m\e[39m\e[22m                          \e[0m\e[2m\e[34m│\e[22m\e[39m\e[0m\n\e[2m\e[34m│\e[22m\e[39m   \e[22m\e[2m\e[32m└── \e[22m\e[39mComplex\e[22m                   \e[0m\e[2m\e[34m│\e[22m\e[39m\e[0m\n\e[2m\e[34m│\e[22m\e[39m                                 \e[0m\e[2m\e[34m│\e[22m\e[39m\e[0m\n\e[2m\e[34m╰─────────────────────────────────╯\e[22m\e[39m\n"
    else
        @test_nothrow typestree(Float64)
    end
end


text:
using LearningFourierTransform
using Test

@testset "LearningFourierTransform.jl" begin

end


text:
# This file contains functions for the Radial Multiplier DEA model
"""
    AbstractRadialMultiplierDEAModel
An abstract type representing a radial DEA model.
"""
abstract type AbstractRadialMultiplierDEAModel <: AbstractTechnicalDEAModel end

"""
    RadialMultiplierDEAModel
An data structure representing a radial multiplier DEA model.
"""
struct RadialMultiplierDEAModel <: AbstractRadialMultiplierDEAModel
    n::Int64
    m::Int64
    s::Int64
    orient::Symbol
    rts::Symbol
    dmunames::Union{Vector{String},Nothing}
    eff::Vector
    v::Matrix
    u::Matrix
    omega::Vector
    Xtarget::Matrix
    Ytarget::Matrix
end

"""
    deam(X, Y)
Compute the radial multiplier model using data envelopment analysis for inputs X and outputs Y.

# Optional Arguments
- `orient=:Input`: chooses the radially oriented input mode. For the radially oriented output model choose `:Output`.
- `rts=:CRS`: chooses constant returns to scale. For variable returns to scale choose `:VRS`.
- `Xref=X`: Identifies the reference set of inputs against which the units are evaluated.
- `Yref=Y`: Identifies the reference set of outputs against which the units are evaluated.
- `names`: a vector of strings with the names of the decision making units.
"""
function deam(X::Union{Matrix,Vector}, Y::Union{Matrix,Vector};
    orient::Symbol = :Input, rts::Symbol = :CRS,
    Xref::Union{Matrix,Vector,Nothing} = nothing, Yref::Union{Matrix, Vector,Nothing} = nothing,
    names::Union{Vector{String},Nothing} = nothing,
    optimizer::Union{DEAOptimizer,Nothing} = nothing, progress::Bool = false)::RadialMultiplierDEAModel

    # Check parameters
    nx, m = size(X, 1), size(X, 2)
    ny, s = size(Y, 1), size(Y, 2)

    if Xref === nothing Xref = X end
    if Yref === nothing Yref = Y end

    nrefx, mref = size(Xref, 1), size(Xref, 2)
    nrefy, sref = size(Yref, 1), size(Yref, 2)

    if nx != ny
        throw(DimensionMismatch("number of rows in X and Y ($nx, $ny) are not equal"));
    end
    if nrefx != nrefy
        throw(DimensionMismatch("number of rows in Xref and Yref ($nrefx, $nrefy) are not equal"));
    end
    if m != mref
        throw(DimensionMismatch("number of columns in X and Xref ($m, $mref) are not equal"));
    end
    if s != sref
        throw(DimensionMismatch("number of columns in Y and Yref ($s, $sref) are not equal"));
    end

    # Default optimizer
    if optimizer === nothing
        optimizer = DEAOptimizer(:LP)
    end

    # Compute efficiency for each DMU
    n = nx
    nref = nrefx

    effi = zeros(n)
    vi = zeros(n, m)
    ui = zeros(n, s)
    omegai = zeros(n)

    p = progressbarmeter(n, desc = "Computing Radial DEA Model...", progress = progress)

    for i=1:n
        # Value of inputs and outputs to evaluate
        x0 = X[i,:]
        y0 = Y[i,:]

        # Create the optimization model
        deamodel = newdeamodel(optimizer)

        @variable(deamodel, v[1:m] >= 0)
        @variable(deamodel, u[1:s] >= 0)
        @variable(deamodel, omega)

        if orient == :Input
            @objective(deamodel, Max, sum(u[r] * y0[r] for r in 1:s) - omega)

            @constraint(deamodel, sum(v[i] * x0[i] for i in 1:m) == 1)
            @constraint(deamodel, [j in 1:nref], sum(u[r] * Yref[j,r] for r in 1:s) - sum(v[i] * Xref[j,i] for i in 1:m) - omega <= 0)

        elseif orient == :Output
            @objective(deamodel, Min, sum(v[i] * x0[i] for i in 1:m) + omega)

            @constraint(deamodel, sum(u[r] * y0[r] for r in 1:s) == 1)
            @constraint(deamodel, [j in 1:nref], sum(v[i] * Xref[j,i] for i in 1:m) - sum(u[r] * Yref[j,r] for r in 1:s) + omega >= 0)
        else
            throw(ArgumentError("`orient` must be :Input or :Output"));
        end

        # Add return to scale constraints
        if rts == :CRS
            @constraint(deamodel, omega == 0)
        elseif rts == :VRS
            # No contraint to add for variable returns to scale
        else
            throw(ArgumentError("`rts` must be :CRS or :VRS"));
        end

        # Optimize and return results
        JuMP.optimize!(deamodel)

        effi[i]  = JuMP.objective_value(deamodel)
        vi[i, :] = JuMP.value.(v)
        ui[i, :] = JuMP.value.(u)
        omegai[i] = JuMP.value(omega)

        # Check termination status
        if (termination_status(deamodel) != MOI.OPTIMAL) && (termination_status(deamodel) != MOI.LOCALLY_SOLVED)
            @warn ("DMU $i termination status: $(termination_status(deamodel)). Primal status: $(primal_status(deamodel)). Dual status: $(dual_status(deamodel))")
        end

        next!(p)
    end

    # Get first-stage X and Y targets
    if orient == :Input
        Xtarget = X .* effi
        Ytarget = Y
    elseif orient == :Output
        Xtarget = X
        Ytarget = Y .* effi
    end

    if typeof(Xtarget) <: AbstractVector    Xtarget = Xtarget[:,:]  end
    if typeof(Ytarget) <: AbstractVector    Ytarget = Ytarget[:,:]  end

    return RadialMultiplierDEAModel(n, m, s, orient, rts, names, effi, vi, ui, omegai, Xtarget, Ytarget)
end

function Base.show(io::IO, x::RadialMultiplierDEAModel)
    compact = get(io, :compact, false)

    n = nobs(x)
    m = ninputs(x)
    s = noutputs(x)
    dmunames = names(x)

    eff = efficiency(x)
    v = multipliers(x, :X)
    u = multipliers(x, :Y)

    if !compact
        print(io, "Radial DEA Model (Multiplier form)\n")
        print(io, "DMUs = ", n)
        print(io, "; Inputs = ", m)
        print(io, "; Outputs = ", s)
        print(io, "\n")
        print(io, "Orientation = ", string(x.orient))
        print(io, "; Returns to Scale = ", string(x.rts))
        print(io, "\n")

        show(io, CoefTable(hcat(eff, v, u), ["efficiency"; ["v$i" for i in 1:m ]; ["u$i" for i in 1:s ]], dmunames))
    end
end


text:
function layouttable(io::IO, printer::Function, table::Table, tchars::Dict, indent=0)
    celltext = map(r -> if r isa TableRow
                       [join(sprint.(printer, c.contents)) for c in r.cells]
                   else r end, table.rows)
    colwidths = map(r -> if !isa(r, TableHrule)
                        [textwidth(c) - ansi_escape_textwidth_offset(c) for c in r]
                    end, celltext) |>
                        rs -> filter(!isnothing, rs) |>
                        rs -> maximum(hcat(rs...); dims=2)
    for row in celltext
        print(io, ' '^indent)
        if row isa TableHrule
            print(io, tchars['>'],
                  join(repeat.(tchars['-'], colwidths .+ 2), tchars['+']),
                  tchars['<'], "\n")
        else
            print(io, tchars['|'])
            for (cell, fillwidth) in zip(row, colwidths)
                print(io, tchars['['],
                      rpad(cell, fillwidth + ansi_escape_textwidth_offset(cell)),
                      tchars[']'], tchars['|'])
            end
            row === last(celltext) || print(io, '\n')
        end
    end
end

function ansi_escape_textwidth_offset(s::AbstractString)
    escapeskip, i, slen = 0, 1, ncodeunits(s)
    while i < slen
        if s[i] == '\e' && i-1 < slen
            if s[i+1] == '[' # Control Sequence Introducer
                terminator = findfirst('m', @inbounds @view s[i:end])
                delta = if !isnothing(terminator) terminator-1 else 0 end
                escapeskip += delta
                i += delta
            elseif s[i+1] == ']' # Operating System Command
                terminator = findfirst("\e\\", @inbounds @view s[i:end])
                delta = if !isnothing(terminator) terminator.stop-1 else 0 end
                escapeskip += delta
                i += delta
            else
                @warn "Unrecognised ANSI code $(s[i+1]) in string."
            end
        end
        i = nextind(s, i)
    end
    escapeskip
end

"""
    wraplines(s::AbstractString, width::Integer, offset::Integer)
`s` is assumed to only contain significant whitespace, and no newlines
"""
function wraplines(s::AbstractString, width::Integer, offset::Integer=0)
    lines = SubString{String}[]
    i, lastwrap, slen = 1, 0, ncodeunits(s)
    mostrecentbreakoppotunity = 1
    while i < slen
        if s[i] == ' '
            mostrecentbreakoppotunity = i
        elseif s[i] == '\e'
            if s[i+1] == '[' # Control Sequence Introducer
                terminator = findfirst('m', @inbounds @view s[i:end])
                delta = if !isnothing(terminator) terminator else 0 end
                offset -= delta
                i += delta
            elseif s[i+1] == ']' # Operating System Command
                terminator = findfirst("\e\\", @inbounds @view s[i:end])
                delta = if !isnothing(terminator) terminator.stop else 0 end
                offset -= delta
                i += delta
            else
                @warn "Unrecognised ANSI code $(s[i+1]) in string."
                i = nextind(s, i)
            end
            continue
        elseif s[i] == '\n'
            push!(lines, @inbounds @view s[nextind(s, lastwrap):prevind(s, i)])
            lastwrap = i
            offset = 0
        elseif i - (lastwrap - offset) > width && mostrecentbreakoppotunity > 1
            if lastwrap == mostrecentbreakoppotunity
                nextbreak = findfirst(' ', @inbounds @view s[nextind(s, lastwrap):end])
                if isnothing(nextbreak)
                    mostrecentbreakoppotunity = slen
                else
                    mostrecentbreakoppotunity = lastwrap + nextbreak
                end
                i = mostrecentbreakoppotunity
            end
            push!(lines, @inbounds @view s[nextind(s, lastwrap):prevind(s, mostrecentbreakoppotunity)])
            lastwrap = mostrecentbreakoppotunity
            offset = 0
        end
        i = nextind(s, i)
    end
    if lastwrap < slen
        push!(lines, @inbounds @view s[nextind(s, lastwrap):end])
    end
    lines
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

"""
    @opaque ([type, ]args...) -> body

Marks a given closure as "opaque". Opaque closures capture the
world age of their creation (as opposed to their invocation).
This allows for more aggressive optimization of the capture
list, but trades off against the ability to inline opaque
closures at the call site, if their creation is not statically
visible.

An argument tuple type (`type`) may optionally be specified, to
specify allowed argument types in a more flexible way. In particular,
the argument type may be fixed length even if the function is variadic.

!!! warning
    This interface is experimental and subject to change or removal without notice.
"""
macro opaque(ex)
    esc(Expr(:opaque_closure, ex))
end

macro opaque(ty, ex)
    esc(Expr(:opaque_closure, ty, ex))
end

max_stars_repo_name: JuliaLabs/julia

max_issues_repo_name: JuliaLabs/julia

max_forks_repo_name: JuliaLabs/julia


text:
include("einsum3.jl")
include("./data_generation.jl")
include("./model.jl")
# include("../../GitHub/Einsum.jl/src/einsum2.jl")
import .ClusterTests
import .DataGeneration
# n = 500
# ang = pi/5
# noise = 0.05
# gaussian = false
# lambdamin = .1
# lambdamax = .5
# numlambda = 250
# lambdavec = collect(range(lambdamin, stop=lambdamax, length=numlambda))
# outfilename = "halfmoon_pi5_test5t.txt"
# maxit = 50000
# t0 = time()
# DataGeneration.run_test_lambda_halfmoon(n, ang, noise, gaussian,
#                                       lambdavec, 10,
#                                       outfilename, maxit)
# print(time() - t0)

"""
Apply sum-of-norms clustering to a mixture of Gaussians
"""
n = 50    # number of points
d = 2     # dimension
k = 2    # number of Gaussians
# sigma = 1.0    # standard deviation is sigma * sqrt(2)
theta = 0.7    # we only count Rand index for data points that are located at 0.5 * theta * min mean distance away from the respective means
lambdamulseq = 0.1:0.05:1.2    # the multiple of the lambda value which is supposed to recover Gaussians exactly from the theorem
meandistvec = 1:6    # the range of mean distance, mean distance is the ratio between the pairwise distance of Gaussian means and standard deviation
maxit = 50000    # maximum iteration
outfilename = "test_lambda3.txt"
t0 = time()
DataGeneration.run_test_lambda_rho_meandistance(n, d, k, theta, lambdamulseq, meandistvec, maxit, outfilename)
print(time() - t0)

max_stars_repo_name: taotolojiang/Certifying-clusters-SON

max_issues_repo_name: taotolojiang/Certifying-clusters-SON

max_forks_repo_name: taotolojiang/Certifying-clusters-SON


text:
#!/usr/local/bin/julia

include("cfrac.jl")

if length(ARGS) != 1
    print_with_color(:blue, "encode: \n")
    println("encode a file to continued fraction")
    println("usage: julia encode.jl ≪filename≫")
    exit(0)
end

#=
str = """
This is a sentence to be shadowed.
This sentence may be very long.
And it may contain unicode characters.
∫e^x ∂x = e^x+C
∀x(human(x)⇒mortal(x))
"""
=#

filename = ARGS[1]
inputfile = open(filename, "r")

str = join(readlines(inputfile), "\n")

arraytype{T,N}(::Array{T,N}) = T

acc = big(0)
arr = transcode(Cwchar_t, str)
u = typemax(arraytype(arr))
for x in arr
    acc *= big(u) + 1
    acc += x
end

k = Int(round(log2(acc))) + 21
g = Int(floor(log10(acc)))

setprecision(BigFloat, k) do
    v = BigFloat(acc)
    v /= big(10)^g

    m = k÷4
    infl = cfrac(v)
    l = take(m, infl)
    #=
    outputfile = open("[encoded]"*filename, "w")
    info("first iteration")
    println(outputfile, "continued-fraction encoded file with ", g, " decimal precision")
    println(l)
    close(outputfile)
    =#

    info("trying to decrease...")
    outputfile = open("[encoded]"*filename, "w")
    v′ = restore_cfrac(l)
    while v != v′
        m += g÷50
        l = take(m, infl)
        v′ = restore_cfrac(l)
        info("iter: trying ", m, " numbers")
    end
    @show str
    @show acc
    @show k, m, g
    @show v
    @show l
    @show v′
    println(outputfile, "continued-fraction encoded file with ", g, " decimal precision")
    println(l)
end

max_stars_repo_name: sclereid/collection

max_issues_repo_name: sclereid/collection

max_forks_repo_name: sclereid/collection


text:
"""
    NestedModel(loglike, prior_transform)
    NestedModel(loglike, priors::AbstractVector{<:Distribution})

`loglike` must be callable with a signature `loglike(::AbstractVector)` where the length of the vector must match the number of parameters in your model.

`prior_transform` must be a callable with a signature `prior_transform(::AbstractVector)` that returns the transformation from the unit-cube to parameter space. This is effectively the quantile or ppf of a statistical distribution. For convenience, if a vector of `Distribution` is provided (as a set of priors), a transformation function will automatically be constructed using `Distributions.quantile`.

**Note:**
`loglike` is the only function used for likelihood calculations. This means if you want your priors to be used for the likelihood calculations they must be manually included in the `loglike` function.
"""
struct NestedModel <: AbstractModel
    loglike::Function
    prior_transform::Function
end

function NestedModel(loglike, priors::AbstractVector{<:UnivariateDistribution})
    prior_transform(X) = quantile.(priors, X)
    return NestedModel(loglike, prior_transform)
end


text:
module MultiNest

import Base: run

export nested

# find the symbol to invoke MultiNest
# this function opens the MultiNest library, which is kept open to allow
# multiple runs
function nested_symbol()
    # common names for MultiNest library
    names = [ "libmultinest", "libnest3" ]

    # names for library when compiled with MPI support
    names_mpi = [ "libmultinest_mpi", "libnest3_mpi" ]

    # check whether `const MULTINEST_MPI = true` was defined in parent module
    mpi = isdefined(module_parent(current_module()), :MULTINEST_MPI) && module_parent(current_module()).MULTINEST_MPI

    # try find library on DL_LOAD_PATH
    lib = Libdl.find_library(mpi ? names_mpi : names, ASCIIString[])
    lib == "" && error("cannot find MultiNest library")

    # now open library
    dl = Libdl.dlopen_e(lib)
    dl == C_NULL && error("cannot open MultiNest library")

    # try to find symbol in library
    for sym in [ :__nested_MOD_nestrun, :nested_mp_nestrun_, :nested_mock ]
        (nestrun_ptr = Libdl.dlsym(dl, sym)) != C_NULL && return nestrun_ptr
    end

    # symbol not found
    error("cannot link MultiNest library, check symbol table")
end

# pointer to function that runs MultiNest
const nestrun_ptr = nested_symbol()

# convert to Fortran logical
macro logical(x)
    :(convert(Cint, $x != 0))
end

# this holds all the information to run MultiNest
immutable Nested
    ins::Cint
    mmodal::Cint
    ceff::Cint
    nlive::Cint
    tol::Cdouble
    efr::Cdouble
    ndim::Cint
    npar::Cint
    nclspar::Cint
    maxmodes::Cint
    updint::Cint
    ztol::Cdouble
    root::ASCIIString
    seed::Cint
    wrap::Vector{Cint}
    fb::Cint
    resume::Cint
    outfile::Cint
    initmpi::Cint
    logzero::Cdouble
    maxiter::Cint
    loglike::Function
    dumper::Function
    context::Any
end

# interface between MultiNest and Julia loglike function
function nested_loglike(
    cube_::Ptr{Cdouble},
    ndim_::Ptr{Cint},
    npar_::Ptr{Cint},
    lnew_::Ptr{Cdouble},
    nested_::Ptr{Void}
)
    ndim = unsafe_load(ndim_)
    npar = unsafe_load(npar_)
    cube = pointer_to_array(cube_, npar)
    nested = unsafe_pointer_to_objref(nested_)::Nested
    lnew = nested.loglike(cube, nested.context...)
    unsafe_store!(lnew_, lnew)
    return
end

# interface between MultiNest and Julia dumper function
function nested_dumper(
    nsamples_::Ptr{Cint},
    nlive_::Ptr{Cint},
    npar_::Ptr{Cint},
    physlive_::Ptr{Ptr{Cdouble}},
    posterior_::Ptr{Ptr{Cdouble}},
    paramconstr_::Ptr{Ptr{Cdouble}},
    maxloglike_::Ptr{Cdouble},
    logz_::Ptr{Cdouble},
    inslogz_::Ptr{Cdouble},
    logzerr_::Ptr{Cdouble},
    nested_::Ptr{Void}
)
    nsamples = unsafe_load(nsamples_)
    nlive = unsafe_load(nlive_)
    npar = unsafe_load(npar_)
    physlive = pointer_to_array(unsafe_load(physlive_), (nlive, npar+1))
    posterior = pointer_to_array(unsafe_load(posterior_), (nsamples, npar+2))
    paramconstr = pointer_to_array(unsafe_load(paramconstr_), (npar, 4))
    maxloglike = unsafe_load(maxloglike_)
    logz = unsafe_load(logz_)
    inslogz = unsafe_load(inslogz_)
    logzerr = unsafe_load(logzerr_)
    nested = unsafe_pointer_to_objref(nested_)::Nested
    nested.dumper(physlive, posterior, paramconstr, maxloglike, logz, inslogz, logzerr, nested.context...)
    return
end

# default dumper; does nothing
function default_dumper(
    physlive::Array{Cdouble, 2},
    posterior::Array{Cdouble, 2},
    paramconstr::Array{Cdouble, 2},
    maxloglike::Cdouble,
    logz::Cdouble,
    inslogz::Cdouble,
    logzerr::Cdouble,
    context...
)
    # noop
end

# wraparound can be given globally as Bool or as a vector for each dimension
const Wrap = Union{Bool, Vector{Bool}}

# create a Nested object from the various options given
function nested(
    loglike::Function,
    ndim::Integer,
    root::AbstractString;
    ins::Bool = true,
    mmodal::Bool = true,
    ceff::Bool = false,
    nlive::Integer = 1000,
    tol::AbstractFloat = 0.5,
    efr::AbstractFloat = 0.8,
    npar::Integer = ndim,
    nclspar::Integer = npar,
    maxmodes::Integer = 100,
    updint::Integer = 1000,
    ztol::AbstractFloat = -1E90,
    seed::Integer = -1,
    wrap::Wrap = false,
    fb::Bool = true,
    resume::Bool = false,
    outfile::Bool = true,
    initmpi::Bool = true,
    logzero::AbstractFloat = nextfloat(-Inf),
    maxiter::Integer = 0,
    dumper::Function = default_dumper,
    context = ()
    )

    # create root dir if necessary (avoid Fortran runtime error)
    root_dir = dirname(root)
    if !isdir(root_dir)
        mkdir(root_dir)
    end

    # convert arguments to the necessary representation
    ins_c = @logical ins
    mmodal_c = @logical mmodal
    ceff_c = @logical ceff
    root_c = rpad(root, 100, ' ')
    wrap_c = ndims(wrap) == 0 ? [ @logical wrap for i in 1:ndim ] : [ @logical w for w in wrap ]
    fb_c = @logical fb
    resume_c = @logical resume
    outfile_c = @logical outfile
    initmpi_c = @logical initmpi

    # sanity checks
    ndim > 0 || error("ndim must be positive")
    nlive > 0 || error("nlive must be positive")
    tol > 0 || error("tol must be positive")
    0 < efr <= 1 || error("efr must be in (0, 1]")
    npar >= ndim || error("npar must be greater than or equal to ndim")
    nclspar <= npar || error("nclspar must be less than or equal to npar")
    maxmodes > 0 || error("maxmodes must be positive")
    updint > 0 || error("updint must be positive")
    size(wrap_c) == (ndim,) || error("size of wrap does not match ndim")
    maxiter >= 0 || error("maxiter must be positive or zero for no maximum")

    # create Nested object
    Nested(ins_c, mmodal_c, ceff_c, nlive, tol, efr, ndim, npar, nclspar,
           maxmodes, updint, ztol, root_c, seed, wrap_c, fb_c, resume_c,
           outfile_c, initmpi_c, logzero, maxiter, loglike, dumper, context)
end

# run MultiNest
@eval begin
    function run(n::Nested)
        loglike_c = cfunction(nested_loglike, Void, (
            Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}, Ptr{Void}
        ))

        dumper_c = cfunction(nested_dumper, Void, (
            Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Ptr{Cdouble}},
            Ptr{Ptr{Cdouble}}, Ptr{Ptr{Cdouble}}, Ptr{Cdouble},
            Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Void}
        ))

        ccall($nestrun_ptr, Void, ( Ptr{Cint}, Ptr{Cint}, Ptr{Cint},
              Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint},
              Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}, Ptr{Uint8},
              Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint},
              Ptr{Cdouble}, Ptr{Cint}, Ptr{Void}, Ptr{Void}, Any ),
              &n.ins, &n.mmodal, &n.ceff, &n.nlive, &n.tol, &n.efr,
              &n.ndim, &n.npar, &n.nclspar, &n.maxmodes, &n.updint, &n.ztol,
              n.root, &n.seed, n.wrap, &n.fb, &n.resume, &n.outfile, &n.initmpi,
              &n.logzero, &n.maxiter, loglike_c, dumper_c, n)
    end
end

end


text:
"""This file contains basic abstractions for molecular spectroscopy
in Julia, based around the two types `Level` and `Transition`
corresponding to energy levels, and transitions between energy
levels respectively.
"""

"""Concrete type for an energy level. The three fields
correspond to the energy of the level, its degeneracy,
and a vector of integers encoding its quantum numbering.

For performance reasons, the encoding is stored as a
vector as opposed to a dictionary mapping with human
readable names. The `quantum_numbers` function provides
an interface to obtain this.
"""
@with_kw struct Level{T,U,S}
  e::T = 1f0
  g::U = 1
  encoding::Vector{S} = [0]
end

"""Concrete type for a transition object. The three fields
correspond to the frequency of the transition, the intensity
of the transition (such as the intrinsic linestrength, the
Einstein A coefficient), and the quantum number encoding
for the transition.

While not explicitly enforced, the encoding is expected to
contain an even number of elements corresponding to lower
and upper state quantum numbers.
"""
@with_kw struct Transition{T,U,S}
  ν::T = 1f0
  I::U = 1f0
  encoding::Vector{S} = [0,1]
  lower::Union{Nothing, Level} = nothing
  upper::Union{Nothing, Level} = nothing
end

# Aliases for collections of either types
Levels = Vector{Level}
Transitions = Vector{Transition}

e(l::Level) = l.e
g(l::Level) = l.g

ν(t::Transition) = t.ν
I(t::Transition) = t.I

# Generate a mapping between string quantum numbers and their value
function quantum_numbers(l::Level, names...)
  return Dict(Symbol(name) => value for (name, value) in zip(names, l.encoding))
end

function quantum_numbers(t::Transition, names...)
  encoding = Dict()
  for (name, l, u) in zip(names, lower_state_qnos(t), upper_state_qnos(t))
    encoding = merge(encoding, Dict(Symbol("$name", "_l") => l, Symbol("$name", "_u") => u))
  end
  return encoding
end

"""Function that converts the fields in a `Level` object into
a tuple of a single, promoted type.
"""
features(l::Level) = promote([l.e, l.g, l.encoding...]...)
features(t::Transition) = promote([t.ν, t.I, t.encoding...]...)

"""Function that converts a collection of spectroscopic types
into a contiguous matrix of promoted types, with the intention
of doing some kind of computation on them.
"""
features(ts) = mapreduce(x -> collect(features(x)), hcat, ts)

# function to compute the line intensity from Aij and partition function
intensity(t::Transition, z::Real) = t.I / z

"""Function to randomly generate linear molecule energy levels.
This is generally for testing purposes than for actual use.
"""
function make_linear_levels(n)
  e = [100f0 * i * (i + 1) for i ∈ 0:n-1]
  j = [[i, 0] for i ∈ 0:n-1]
  g = [2 * i + 1 for i ∈ 0:n-1]
  StructArray([Level(g=g[i], e=e[i], encoding=j[i]) for i ∈ eachindex(e,j,g)])
end

"""Find the index of an spectroscopic object from a vector
of objects that matches the encoding provided. In other words,
match the quantum numbers.
"""
function match_encoding(obj, enc)
  @inbounds for i ∈ eachindex(obj)
    if obj[i].encoding == enc
  	  return i
    end
  end
end

"""Function that processes an iterable of energy levels, and
creates a vector of transition objects from them, assuming that
transitions between each contiguous energy level.

Mostly for testing purposes.

A quick way to do this would be to pipe the output of
`make_linear_levels` into this function.
"""
function make_linear_transitions(levels)
  transitions = Transitions()
  sizehint!(transitions, length(levels))
  for i ∈ 1:length(levels)-1
    l, u = levels[i], levels[i+1]
    ν = abs(e(l) - e(u))
    encoding = vcat(l.encoding, u.encoding)
    push!(transitions, Transition(ν=ν, encoding=encoding, lower=l, upper=u))
  end
  return StructArray(transitions)
end

function upper_state_qnos(t::Transition)
  offset = length(t.encoding) ÷ 2
  return t.encoding[offset+1:length(t.encoding)]
end

function lower_state_qnos(t::Transition)
  offset = length(t.encoding) ÷ 2
  return t.encoding[1:offset]
end

"""Type for storing spectral data in the form of
frequency versus intensity arrays.
"""
@with_kw mutable struct Spectrum{T}
  ν::Vector{T}
  I::Vector{T}
  noise::Union{T, Vector{T}} = 1f0
end


text:
using DataStructures

mutable struct Compound
  word  :: String
  rules :: Dict{String,String}
end

function step!(c::Compound)
  pieces = map(i -> c.rules[c.word[i:i+1]], 1:length(c.word)-1)
  c.word = string(pieces...,c.word[end])
end

function parseLine(line::String) :: Pair{String,String}
  s = split(line," -> ")
  return Pair(string(s[1]), string(s[1][1],s[2]))
end

function readInput()
  lines = readlines("./inputs/day14.txt")

  word = lines[1]
  rules :: Dict{String,String} = Dict{String,String}()
  for line in lines[3:end]
    rule = parseLine(line)
    rules[rule.first] = rule.second
  end

  return Compound(word,rules)
end

myCompound = readInput()

nSteps = 40
for i in 1:nSteps
  println(i)
  step!(myCompound)
end

cnt = counter(myCompound.word)

display(maximum(x->x.second,cnt)-minimum(x->x.second,cnt))


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

mutable struct CreateTrustedIdProviderWithAccountParameters <: SwaggerModel
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    properties::Any # spec type: Union{ Nothing, CreateOrUpdateTrustedIdProviderProperties } # spec name: properties

    function CreateTrustedIdProviderWithAccountParameters(;name=nothing, properties=nothing)
        o = new()
        validate_property(CreateTrustedIdProviderWithAccountParameters, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(CreateTrustedIdProviderWithAccountParameters, Symbol("properties"), properties)
        setfield!(o, Symbol("properties"), properties)
        o
    end
end # type CreateTrustedIdProviderWithAccountParameters

const _property_map_CreateTrustedIdProviderWithAccountParameters = Dict{Symbol,Symbol}(Symbol("name")=>Symbol("name"), Symbol("properties")=>Symbol("properties"))
const _property_types_CreateTrustedIdProviderWithAccountParameters = Dict{Symbol,String}(Symbol("name")=>"String", Symbol("properties")=>"CreateOrUpdateTrustedIdProviderProperties")
Base.propertynames(::Type{ CreateTrustedIdProviderWithAccountParameters }) = collect(keys(_property_map_CreateTrustedIdProviderWithAccountParameters))
Swagger.property_type(::Type{ CreateTrustedIdProviderWithAccountParameters }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_CreateTrustedIdProviderWithAccountParameters[name]))}
Swagger.field_name(::Type{ CreateTrustedIdProviderWithAccountParameters }, property_name::Symbol) =  _property_map_CreateTrustedIdProviderWithAccountParameters[property_name]

function check_required(o::CreateTrustedIdProviderWithAccountParameters)
    (getproperty(o, Symbol("name")) === nothing) && (return false)
    (getproperty(o, Symbol("properties")) === nothing) && (return false)
    true
end

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


text:
function main_function13!(r::Robot)
    x = 0
    y = 0
    while !isborder(r, Sud)
      move!(r, Sud)
      y += 1
    end
    while !isborder(r, West)
      move!(r, West)
      x += 1
    end

    putmarkerss(r)

    for i in 1:x
      move!(r, Ost)
    end

    for i in 1:y
      move!(r, Nord)
    end

  end
  function putmarkerss(r::Robot)
    putmarker!(r)
    while !isborder(r, Nord) && !isborder(r, Ost)
      move!(r, Nord)
      move!(r, Ost)
      putmarker!(r)
    end
    while !isborder(r, West)
      move!(r, West)
    end
    putmarker!(r)
    while !isborder(r, Sud) && !isborder(r, Ost)
      move!(r, Sud)
      move!(r, Ost)
      putmarker!(r)
    end
    while !isborder(r, West)
      move!(r, West)
    end
  end


text:
import Base.*

export S¹
export U1

"""
    Represents a point in a one-sphere.
"""
abstract type S¹ end

"""
    Represents a point in U(1).

field: α.
"""
struct U1 <: S¹
    r::Float64
    U1(z::Complex) = begin
        @assert(isapprox(abs(z), 1), "The magnitude must be equal to 1, but it's $(abs(z)).")
        θ = angle(z) ≥ 0 ? angle(z) : π + angle(z)
        new(θ)
    end
    U1(α::Real) = begin
        #@assert(-pi ≤ α ≤ pi, "The phase angle must be in the interval [-π, π].")
        new(float(α))
    end
end

Base.angle(u::U1) = u.r
*(u1::U1, u2::U1) = U1(u1.r + u2.r)
*(u::U1, scale::Float64) = U1(u.r * scale)
*(scale::Float64, u::U1) = U1(u.r * scale)
Base.isapprox(u1::U1, u2::U1) = isapprox(u1.r, u2.r)

max_stars_repo_name: iamazadi/geometrizer

max_issues_repo_name: iamazadi/Porta

max_forks_repo_name: iamazadi/Porta


text:
import .ForwardDiff

_eps(::Type{<:ForwardDiff.Dual{<:Any, Real}}) = _eps(Real)
_eps(::Type{<:ForwardDiff.Dual{<:Any, <:Integer}}) = _eps(Real)

# AD implementations
function jacobian(
    b::Union{<:ADBijector{<:ForwardDiffAD}, Inverse{<:ADBijector{<:ForwardDiffAD}}},
    x::Real
)
    return ForwardDiff.derivative(b, x)
end
function jacobian(
    b::Union{<:ADBijector{<:ForwardDiffAD}, Inverse{<:ADBijector{<:ForwardDiffAD}}},
    x::AbstractVector{<:Real}
)
    return ForwardDiff.jacobian(b, x)
end


text:
@testset "Globals" begin

function kernel(X)
    ptr = AMDGPUnative.get_global_pointer(Val(:myglobal), Float32)
    Base.unsafe_store!(ptr, 3f0)
    nothing
end

hk = AMDGPUnative.rocfunction(kernel, Tuple{Int32})
gbl = HSARuntime.get_global(hk.mod.exe, :myglobal)
gbl_ptr = Base.unsafe_convert(Ptr{Float32}, gbl.ptr)

@test Base.unsafe_load(gbl_ptr) == 0f0
Base.unsafe_store!(gbl_ptr, 2f0)
@test Base.unsafe_load(gbl_ptr) == 2f0

wait(@roc groupsize=1 kernel(Int32(1)))
@test Base.unsafe_load(gbl_ptr) == 3f0

end


text:
#! format: off

abstract type AbstractHydroFormulation <: AbstractDeviceFormulation end
abstract type AbstractHydroDispatchFormulation <: AbstractHydroFormulation end
abstract type AbstractHydroUnitCommitment <: AbstractHydroFormulation end
abstract type AbstractHydroReservoirFormulation <: AbstractHydroDispatchFormulation end
struct HydroDispatchRunOfRiver <: AbstractHydroDispatchFormulation end
struct HydroDispatchReservoirBudget <: AbstractHydroReservoirFormulation end
struct HydroDispatchReservoirStorage <: AbstractHydroReservoirFormulation end
struct HydroDispatchPumpedStorage <: AbstractHydroReservoirFormulation end
struct HydroDispatchPumpedStoragewReservation <: AbstractHydroReservoirFormulation end
struct HydroCommitmentRunOfRiver <: AbstractHydroUnitCommitment end
struct HydroCommitmentReservoirBudget <: AbstractHydroUnitCommitment end
struct HydroCommitmentReservoirStorage <: AbstractHydroUnitCommitment end

get_variable_sign(_, ::Type{<:PSY.HydroGen}, ::AbstractHydroFormulation) = 1.0
########################### ActivePowerVariable, HydroGen #################################
get_variable_binary(::ActivePowerVariable, ::Type{<:PSY.HydroGen}, ::AbstractHydroFormulation) = false
get_variable_expression_name(::ActivePowerVariable, ::Type{<:PSY.HydroGen}) = :nodal_balance_active

get_variable_initial_value(::ActivePowerVariable, d::PSY.HydroGen, ::AbstractHydroFormulation) = PSY.get_active_power(d)

get_variable_lower_bound(::ActivePowerVariable, d::PSY.HydroGen, ::AbstractHydroFormulation) = PSY.get_active_power_limits(d).min
get_variable_upper_bound(::ActivePowerVariable, d::PSY.HydroGen, ::AbstractHydroFormulation) = PSY.get_active_power_limits(d).max

############## ActivePowerVariable, HydroDispatchRunOfRiver ####################
get_variable_lower_bound(::ActivePowerVariable, d::PSY.HydroGen, ::HydroDispatchRunOfRiver) = 0.0

############## ReactivePowerVariable, HydroGen ####################
get_variable_binary(::ReactivePowerVariable, ::Type{<:PSY.HydroGen}, ::AbstractHydroFormulation) = false
get_variable_expression_name(::ReactivePowerVariable, ::Type{<:PSY.HydroGen}) = :nodal_balance_reactive
get_variable_initial_value(::ReactivePowerVariable, d::PSY.HydroGen, ::AbstractHydroFormulation) = PSY.get_active_power(d)
get_variable_lower_bound(::ReactivePowerVariable, d::PSY.HydroGen, ::AbstractHydroFormulation) = PSY.get_active_power_limits(d).min
get_variable_upper_bound(::ReactivePowerVariable, d::PSY.HydroGen, ::AbstractHydroFormulation) = PSY.get_active_power_limits(d).max

############## EnergyVariable, HydroGen ####################
get_variable_binary(::EnergyVariable, ::Type{<:PSY.HydroGen}, ::AbstractHydroFormulation) = false
get_variable_initial_value(pv::EnergyVariable, d::PSY.HydroGen, ::AbstractHydroFormulation) = PSY.get_initial_storage(d)
get_variable_lower_bound(::EnergyVariable, d::PSY.HydroGen, ::AbstractHydroFormulation) = 0.0
get_variable_upper_bound(::EnergyVariable, d::PSY.HydroGen, ::AbstractHydroFormulation) = PSY.get_storage_capacity(d)

########################### EnergyVariableUp, HydroGen #################################

get_variable_binary(::EnergyVariableUp, ::Type{<:PSY.HydroGen}, ::AbstractHydroFormulation) = false

get_variable_initial_value(pv::EnergyVariableUp, d::PSY.HydroGen, ::AbstractHydroFormulation) = PSY.get_initial_storage(d).up

get_variable_lower_bound(::EnergyVariableUp, d::PSY.HydroGen, ::AbstractHydroFormulation) = 0.0
get_variable_upper_bound(::EnergyVariableUp, d::PSY.HydroGen, ::AbstractHydroFormulation) = PSY.get_storage_capacity(d).up

########################### EnergyVariableDown, HydroGen #################################

get_variable_binary(::EnergyVariableDown, ::Type{<:PSY.HydroGen}, ::AbstractHydroFormulation) = false

get_variable_initial_value(::EnergyVariableDown, d::PSY.HydroGen, ::AbstractHydroFormulation) = PSY.get_initial_storage(d).down

get_variable_lower_bound(::EnergyVariableDown, d::PSY.HydroGen, ::AbstractHydroFormulation) = 0.0
get_variable_upper_bound(::EnergyVariableDown, d::PSY.HydroGen, ::AbstractHydroFormulation) = PSY.get_storage_capacity(d).down

########################### ActivePowerInVariable, HydroGen #################################

get_variable_binary(::ActivePowerInVariable, ::Type{<:PSY.HydroGen}, ::AbstractHydroFormulation) = false
get_variable_expression_name(::ActivePowerInVariable, ::Type{<:PSY.HydroGen}) = :nodal_balance_active

get_variable_lower_bound(::ActivePowerInVariable, d::PSY.HydroGen, ::AbstractHydroFormulation) = 0.0
get_variable_upper_bound(::ActivePowerInVariable, d::PSY.HydroGen, ::AbstractHydroFormulation) = nothing
get_variable_sign(::ActivePowerInVariable, d::Type{<:PSY.HydroGen}, ::AbstractHydroFormulation) = -1.0

########################### ActivePowerOutVariable, HydroGen #################################

get_variable_binary(::ActivePowerOutVariable, ::Type{<:PSY.HydroGen}, ::AbstractHydroFormulation) = false
get_variable_expression_name(::ActivePowerOutVariable, ::Type{<:PSY.HydroGen}) = :nodal_balance_active

get_variable_lower_bound(::ActivePowerOutVariable, d::PSY.HydroGen, ::AbstractHydroFormulation) = 0.0
get_variable_upper_bound(::ActivePowerOutVariable, d::PSY.HydroGen, ::AbstractHydroFormulation) = nothing
get_variable_sign(::ActivePowerOutVariable, d::Type{<:PSY.HydroGen}, ::AbstractHydroFormulation) = 1.0

############## OnVariable, HydroGen ####################

get_variable_binary(::OnVariable, ::Type{<:PSY.HydroGen}, ::AbstractHydroFormulation) = true
get_variable_initial_value(::OnVariable, d::PSY.HydroGen, ::AbstractHydroFormulation) = PSY.get_active_power(d) > 0 ? 1.0 : 0.0

############## SpillageVariable, HydroGen ####################

get_variable_binary(::SpillageVariable, ::Type{<:PSY.HydroGen}, ::AbstractHydroFormulation) = false
get_variable_lower_bound(::SpillageVariable, d::PSY.HydroGen, ::AbstractHydroFormulation) = 0.0

############## ReserveVariable, HydroGen ####################

get_variable_binary(::ReserveVariable, ::Type{<:PSY.HydroGen}, ::AbstractHydroFormulation) = true
get_variable_binary(::ReserveVariable, ::Type{<:PSY.HydroPumpedStorage}, ::AbstractHydroFormulation) = true
############## EnergyShortageVariable, HydroGen ####################

get_variable_binary(::EnergyShortageVariable, ::Type{<:PSY.HydroGen}, ::AbstractHydroFormulation) = false
get_variable_lower_bound(::EnergyShortageVariable, d::PSY.HydroGen, ::AbstractHydroFormulation) = 0.0
get_variable_upper_bound(::EnergyShortageVariable, d::PSY.HydroGen, ::AbstractHydroFormulation) = PSY.get_storage_capacity(d)

############## EnergySurplusVariable, HydroGen ####################

get_variable_binary(::EnergySurplusVariable, ::Type{<:PSY.HydroGen}, ::AbstractHydroFormulation) = false
get_variable_upper_bound(::EnergySurplusVariable, d::PSY.HydroGen, ::AbstractHydroFormulation) = 0.0
get_variable_lower_bound(::EnergySurplusVariable, d::PSY.HydroGen, ::AbstractHydroFormulation) = - PSY.get_storage_capacity(d)

get_efficiency(v::T, var::Type{<:InitialConditionType}) where T <: PSY.HydroGen = (in = 1.0, out = 1.0)
get_efficiency(v::PSY.HydroPumpedStorage, var::Type{InitialEnergyLevelUp}) = (in = PSY.get_pump_efficiency(v), out = 1.0)
get_efficiency(v::PSY.HydroPumpedStorage, var::Type{InitialEnergyLevelDown}) = (in = 1.0, out = PSY.get_pump_efficiency(v))

#! format: on

"""
This function define the range constraint specs for the
reactive power for dispatch formulations.
"""
function DeviceRangeConstraintSpec(
    ::Type{<:RangeConstraint},
    ::Type{ReactivePowerVariable},
    ::Type{T},
    ::Type{<:AbstractHydroDispatchFormulation},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
    use_parameters::Bool,
    use_forecasts::Bool,
) where {T <: PSY.HydroGen}
    return DeviceRangeConstraintSpec(;
        range_constraint_spec = RangeConstraintSpec(;
            constraint_name = make_constraint_name(
                RangeConstraint,
                ReactivePowerVariable,
                T,
            ),
            variable_name = make_variable_name(ReactivePowerVariable, T),
            limits_func = x -> PSY.get_reactive_power_limits(x),
            constraint_func = device_range!,
            constraint_struct = DeviceRangeConstraintInfo,
        ),
    )
end

"""
This function define the range constraint specs for the
active power for dispatch Run of River formulations.
"""
function DeviceRangeConstraintSpec(
    ::Type{<:RangeConstraint},
    ::Type{ActivePowerVariable},
    ::Type{T},
    ::Type{<:AbstractHydroDispatchFormulation},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
    use_parameters::Bool,
    use_forecasts::Bool,
) where {T <: PSY.HydroGen}
    if !use_parameters && !use_forecasts
        return DeviceRangeConstraintSpec(;
            range_constraint_spec = RangeConstraintSpec(;
                constraint_name = make_constraint_name(
                    RangeConstraint,
                    ActivePowerVariable,
                    T,
                ),
                variable_name = make_variable_name(ActivePowerVariable, T),
                limits_func = x -> (min = 0.0, max = PSY.get_active_power(x)),
                constraint_func = device_range!,
                constraint_struct = DeviceRangeConstraintInfo,
            ),
        )
    end

    return DeviceRangeConstraintSpec(;
        timeseries_range_constraint_spec = TimeSeriesConstraintSpec(
            constraint_name = make_constraint_name(RangeConstraint, ActivePowerVariable, T),
            variable_name = make_variable_name(ActivePowerVariable, T),
            parameter_name = use_parameters ? ACTIVE_POWER : nothing,
            forecast_label = "max_active_power",
            multiplier_func = x -> PSY.get_max_active_power(x),
            constraint_func = use_parameters ? device_timeseries_param_ub! :
                              device_timeseries_ub!,
        ),
    )
end

"""
This function define the range constraint specs for the
active power for dispatch Reservoir formulations.
"""
function DeviceRangeConstraintSpec(
    ::Type{<:RangeConstraint},
    ::Type{ActivePowerVariable},
    ::Type{T},
    ::Type{<:AbstractHydroReservoirFormulation},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
    use_parameters::Bool,
    use_forecasts::Bool,
) where {T <: PSY.HydroGen}
    return DeviceRangeConstraintSpec(;
        range_constraint_spec = RangeConstraintSpec(;
            constraint_name = make_constraint_name(RangeConstraint, ActivePowerVariable, T),
            variable_name = make_variable_name(ActivePowerVariable, T),
            limits_func = x -> PSY.get_active_power_limits(x),
            constraint_func = device_range!,
            constraint_struct = DeviceRangeConstraintInfo,
        ),
    )
end

"""
This function define the range constraint specs for the
active power for commitment formulations (semi continuous).
"""
function DeviceRangeConstraintSpec(
    ::Type{<:RangeConstraint},
    ::Type{ActivePowerVariable},
    ::Type{T},
    ::Type{<:AbstractHydroUnitCommitment},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::Nothing,
    use_parameters::Bool,
    use_forecasts::Bool,
) where {T <: PSY.HydroGen}
    return DeviceRangeConstraintSpec(;
        range_constraint_spec = RangeConstraintSpec(;
            constraint_name = make_constraint_name(RangeConstraint, ActivePowerVariable, T),
            variable_name = make_variable_name(ActivePowerVariable, T),
            bin_variable_names = [make_variable_name(OnVariable, T)],
            limits_func = x -> PSY.get_active_power_limits(x),
            constraint_func = device_semicontinuousrange!,
            constraint_struct = DeviceRangeConstraintInfo,
        ),
    )
end

"""
This function define the range constraint specs for the
reactive power for commitment formulations (semi continuous).
"""
function DeviceRangeConstraintSpec(
    ::Type{<:RangeConstraint},
    ::Type{ReactivePowerVariable},
    ::Type{T},
    ::Type{<:AbstractHydroUnitCommitment},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::Nothing,
    use_parameters::Bool,
    use_forecasts::Bool,
) where {T <: PSY.HydroGen}
    return DeviceRangeConstraintSpec(;
        range_constraint_spec = RangeConstraintSpec(;
            constraint_name = make_constraint_name(
                RangeConstraint,
                ReactivePowerVariable,
                T,
            ),
            variable_name = make_variable_name(ReactivePowerVariable, T),
            bin_variable_names = [make_variable_name(OnVariable, T)],
            limits_func = x -> PSY.get_active_power_limits(x),
            constraint_func = device_semicontinuousrange!,
            constraint_struct = DeviceRangeConstraintInfo,
        ),
    )
end

function DeviceRangeConstraintSpec(
    ::Type{<:RangeConstraint},
    ::Type{ActivePowerOutVariable},
    ::Type{T},
    ::Type{<:HydroDispatchPumpedStorage},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
    use_parameters::Bool,
    use_forecasts::Bool,
) where {T <: PSY.HydroGen}
    return DeviceRangeConstraintSpec(;
        range_constraint_spec = RangeConstraintSpec(;
            constraint_name = make_constraint_name(
                RangeConstraint,
                ActivePowerOutVariable,
                T,
            ),
            variable_name = make_variable_name(ActivePowerOutVariable, T),
            limits_func = x -> PSY.get_active_power_limits(x),
            constraint_func = device_range!,
            constraint_struct = DeviceRangeConstraintInfo,
        ),
    )
end

function DeviceRangeConstraintSpec(
    ::Type{<:RangeConstraint},
    ::Type{ActivePowerInVariable},
    ::Type{T},
    ::Type{<:HydroDispatchPumpedStorage},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
    use_parameters::Bool,
    use_forecasts::Bool,
) where {T <: PSY.HydroGen}
    return DeviceRangeConstraintSpec(;
        range_constraint_spec = RangeConstraintSpec(;
            constraint_name = make_constraint_name(
                RangeConstraint,
                ActivePowerInVariable,
                T,
            ),
            variable_name = make_variable_name(ActivePowerInVariable, T),
            limits_func = x -> PSY.get_active_power_limits_pump(x),
            constraint_func = device_range!,
            constraint_struct = DeviceRangeConstraintInfo,
        ),
    )
end

function DeviceRangeConstraintSpec(
    ::Type{<:RangeConstraint},
    ::Type{ActivePowerOutVariable},
    ::Type{T},
    ::Type{<:HydroDispatchPumpedStoragewReservation},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
    use_parameters::Bool,
    use_forecasts::Bool,
) where {T <: PSY.HydroGen}
    return DeviceRangeConstraintSpec(;
        range_constraint_spec = RangeConstraintSpec(;
            constraint_name = make_constraint_name(
                RangeConstraint,
                ActivePowerOutVariable,
                T,
            ),
            variable_name = make_variable_name(ActivePowerOutVariable, T),
            bin_variable_names = [make_variable_name(ReserveVariable, T)],
            limits_func = x -> PSY.get_active_power_limits(x),
            constraint_func = reserve_device_semicontinuousrange!,
            constraint_struct = DeviceRangeConstraintInfo,
        ),
    )
end

function DeviceRangeConstraintSpec(
    ::Type{<:RangeConstraint},
    ::Type{ActivePowerInVariable},
    ::Type{T},
    ::Type{<:HydroDispatchPumpedStoragewReservation},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
    use_parameters::Bool,
    use_forecasts::Bool,
) where {T <: PSY.HydroGen}
    return DeviceRangeConstraintSpec(;
        range_constraint_spec = RangeConstraintSpec(;
            constraint_name = make_constraint_name(
                RangeConstraint,
                ActivePowerInVariable,
                T,
            ),
            variable_name = make_variable_name(ActivePowerInVariable, T),
            bin_variable_names = [make_variable_name(ReserveVariable, T)],
            limits_func = x -> PSY.get_active_power_limits_pump(x),
            constraint_func = reserve_device_semicontinuousrange!,
            constraint_struct = DeviceRangeConstraintInfo,
        ),
    )
end
######################## RoR constraints ############################

"""
This function define the range constraint specs for the
reactive power for Commitment Run of River formulation.
    `` P <= multiplier * P_max ``
"""
function commit_hydro_active_power_ub!(
    optimization_container::OptimizationContainer,
    devices,
    model::DeviceModel{V, W},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
) where {V <: PSY.HydroGen, W <: AbstractHydroUnitCommitment}
    use_parameters = model_has_parameters(optimization_container)
    use_forecasts = model_uses_forecasts(optimization_container)
    if use_parameters || use_forecasts
        spec = DeviceRangeConstraintSpec(;
            timeseries_range_constraint_spec = TimeSeriesConstraintSpec(
                constraint_name = make_constraint_name(COMMITMENT, V),
                variable_name = make_variable_name(ActivePowerVariable, V),
                parameter_name = use_parameters ? ACTIVE_POWER : nothing,
                forecast_label = "max_active_power",
                multiplier_func = x -> PSY.get_max_active_power(x),
                constraint_func = use_parameters ? device_timeseries_param_ub! :
                                  device_timeseries_ub!,
            ),
        )
        device_range_constraints!(optimization_container, devices, model, feedforward, spec)
    end
end

######################## Energy balance constraints ############################

"""
This function defines the constraints for the water level (or state of charge)
for the Hydro Reservoir.
"""
function DeviceEnergyBalanceConstraintSpec(
    ::Type{<:EnergyBalanceConstraint},
    ::Type{EnergyVariable},
    ::Type{H},
    ::Type{<:AbstractHydroFormulation},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
    use_parameters::Bool,
    use_forecasts::Bool,
) where {H <: PSY.HydroEnergyReservoir}
    return DeviceEnergyBalanceConstraintSpec(;
        constraint_name = make_constraint_name(ENERGY_CAPACITY, H),
        energy_variable = make_variable_name(ENERGY, H),
        initial_condition = InitialEnergyLevel,
        pout_variable_names = [
            make_variable_name(ACTIVE_POWER, H),
            make_variable_name(SPILLAGE, H),
        ],
        constraint_func = use_parameters ? energy_balance_param! : energy_balance!,
        parameter_name = INFLOW,
        forecast_label = "inflow",
        multiplier_func = x -> PSY.get_inflow(x) * PSY.get_conversion_factor(x),
    )
end

"""
This function defines the constraints for the water level (or state of charge)
for the HydroPumpedStorage.
"""
function DeviceEnergyBalanceConstraintSpec(
    ::Type{<:EnergyBalanceConstraint},
    ::Type{EnergyVariableUp},
    ::Type{H},
    ::Type{<:AbstractHydroFormulation},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
    use_parameters::Bool,
    use_forecasts::Bool,
) where {H <: PSY.HydroPumpedStorage}
    return DeviceEnergyBalanceConstraintSpec(;
        constraint_name = make_constraint_name(ENERGY_CAPACITY_UP, H),
        energy_variable = make_variable_name(ENERGY_UP, H),
        initial_condition = InitialEnergyLevelUp,
        pin_variable_names = [make_variable_name(ACTIVE_POWER_IN, H)],
        pout_variable_names = [
            make_variable_name(ACTIVE_POWER_OUT, H),
            make_variable_name(SPILLAGE, H),
        ],
        constraint_func = use_parameters ? energy_balance_param! : energy_balance!,
        parameter_name = INFLOW,
        forecast_label = "inflow",
        multiplier_func = x -> PSY.get_inflow(x) * PSY.get_conversion_factor(x),
    )
end

function DeviceEnergyBalanceConstraintSpec(
    ::Type{<:EnergyBalanceConstraint},
    ::Type{EnergyVariableDown},
    ::Type{H},
    ::Type{<:AbstractHydroFormulation},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
    use_parameters::Bool,
    use_forecasts::Bool,
) where {H <: PSY.HydroPumpedStorage}
    return DeviceEnergyBalanceConstraintSpec(;
        constraint_name = make_constraint_name(ENERGY_CAPACITY_DOWN, H),
        energy_variable = make_variable_name(ENERGY_DOWN, H),
        initial_condition = InitialEnergyLevelDown,
        pout_variable_names = [make_variable_name(ACTIVE_POWER_IN, H)],
        pin_variable_names = [
            make_variable_name(ACTIVE_POWER_OUT, H),
            make_variable_name(SPILLAGE, H),
        ],
        constraint_func = use_parameters ? energy_balance_param! : energy_balance!,
        parameter_name = OUTFLOW,
        forecast_label = "outflow",
        multiplier_func = x -> PSY.get_outflow(x) * PSY.get_conversion_factor(x),
    )
end

function energy_target_constraint!(
    optimization_container::OptimizationContainer,
    devices::IS.FlattenIteratorWrapper{T},
    model::DeviceModel{T, S},
    system_formulation::Type{<:PM.AbstractPowerModel},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
) where {T <: PSY.HydroGen, S <: AbstractHydroFormulation}
    key = ICKey(InitialEnergyLevel, T)
    parameters = model_has_parameters(optimization_container)
    use_forecast_data = model_uses_forecasts(optimization_container)
    time_steps = model_time_steps(optimization_container)
    target_forecast_label = "storage_target"
    constraint_infos_target = Vector{DeviceTimeSeriesConstraintInfo}(undef, length(devices))
    if use_forecast_data
        for (ix, d) in enumerate(devices)
            ts_vector_target =
                get_time_series(optimization_container, d, target_forecast_label)
            constraint_info_target = DeviceTimeSeriesConstraintInfo(
                d,
                x -> PSY.get_storage_capacity(x),
                ts_vector_target,
            )
            constraint_infos_target[ix] = constraint_info_target
        end
    else
        for (ix, d) in enumerate(devices)
            ts_vector_target =
                length(time_steps) == 1 ? [PSY.get_storage_target(d)] :
                vcat(zeros(time_steps[end - 1]), PSY.get_storage_target(d))
            constraint_info_target = DeviceTimeSeriesConstraintInfo(
                d,
                x -> PSY.get_storage_capacity(x),
                ts_vector_target,
            )
            constraint_infos_target[ix] = constraint_info_target
        end
    end

    if parameters
        energy_target_param!(
            optimization_container,
            constraint_infos_target,
            make_constraint_name(ENERGY_TARGET, T),
            (
                make_variable_name(ENERGY, T),
                make_variable_name(ENERGY_SHORTAGE, T),
                make_variable_name(ENERGY_SURPLUS, T),
            ),
            UpdateRef{T}(TARGET, target_forecast_label),
        )
    else
        energy_target!(
            optimization_container,
            constraint_infos_target,
            make_constraint_name(ENERGY_TARGET, T),
            (
                make_variable_name(ENERGY, T),
                make_variable_name(ENERGY_SHORTAGE, T),
                make_variable_name(ENERGY_SURPLUS, T),
            ),
        )
    end

    constraint_infos = Vector{DeviceRangeConstraintInfo}()
    for (ix, d) in enumerate(devices)
        op_cost = PSY.get_operation_cost(d)
        if PSY.get_energy_shortage_cost(op_cost) == 0.0
            dev_name = PSY.get_name(d)
            limits = (min = 0.0, max = 0.0)
            constraint_info = DeviceRangeConstraintInfo(dev_name, limits)
            push!(constraint_infos, constraint_info)
            @warn(
                "Device $dev_name has energy shortage cost set to 0.0, as a result the model will turnoff the EnergyShortageVariable to avoid infeasible/unbounded problem."
            )
        end
    end
    if !isempty(constraint_infos)
        device_range!(
            optimization_container,
            RangeConstraintSpecInternal(
                constraint_infos,
                make_constraint_name(RangeConstraint, EnergyShortageVariable, T),
                make_variable_name(EnergyShortageVariable, T),
                Vector{Symbol}(),
            ),
        )
    end

    return
end

########################## Make initial Conditions for a Model #############################
function initial_conditions!(
    optimization_container::OptimizationContainer,
    devices::IS.FlattenIteratorWrapper{H},
    device_formulation::AbstractHydroUnitCommitment,
) where {H <: PSY.HydroGen}
    status_initial_condition!(optimization_container, devices, device_formulation)
    output_initial_condition!(optimization_container, devices, device_formulation)
    duration_initial_condition!(optimization_container, devices, device_formulation)

    return
end

function initial_conditions!(
    optimization_container::OptimizationContainer,
    devices::IS.FlattenIteratorWrapper{H},
    device_formulation::AbstractHydroDispatchFormulation,
) where {H <: PSY.HydroGen}
    output_initial_condition!(optimization_container, devices, device_formulation)

    return
end

######################### Initialize Functions for Hydro #################################
function status_initial_condition!(
    optimization_container::OptimizationContainer,
    devices::IS.FlattenIteratorWrapper{T},
    ::D,
) where {T <: PSY.HydroGen, D <: AbstractHydroUnitCommitment}
    _make_initial_conditions!(
        optimization_container,
        devices,
        D(),
        OnVariable(),
        ICKey(DeviceStatus, T),
        _make_initial_condition_active_power,
        _get_variable_initial_value,
        # Doesn't require Cache
    )
end

function output_initial_condition!(
    optimization_container::OptimizationContainer,
    devices::IS.FlattenIteratorWrapper{T},
    ::D,
) where {T <: PSY.HydroGen, D <: AbstractHydroFormulation}
    _make_initial_conditions!(
        optimization_container,
        devices,
        D(),
        ActivePowerVariable(),
        ICKey(DevicePower, T),
        _make_initial_condition_active_power,
        _get_variable_initial_value,
        # Doesn't require Cache
    )

    return
end

function duration_initial_condition!(
    optimization_container::OptimizationContainer,
    devices::IS.FlattenIteratorWrapper{T},
    ::D,
) where {T <: PSY.HydroGen, D <: AbstractHydroUnitCommitment}
    for key in (ICKey(InitialTimeDurationOn, T), ICKey(InitialTimeDurationOff, T))
        _make_initial_conditions!(
            optimization_container,
            devices,
            D(),
            nothing,
            key,
            _make_initial_condition_active_power,
            _get_variable_initial_value,
            TimeStatusChange,
        )
    end

    return
end

function storage_energy_initial_condition!(
    optimization_container::OptimizationContainer,
    devices::IS.FlattenIteratorWrapper{T},
    ::D,
) where {T <: PSY.HydroGen, D <: AbstractHydroFormulation}
    key = ICKey(InitialEnergyLevel, T)
    _make_initial_conditions!(
        optimization_container,
        devices,
        D(),
        EnergyVariable(),
        key,
        _make_initial_condition_reservoir_energy,
        _get_variable_initial_value,
        StoredEnergy,
    )

    return
end

function storage_energy_initial_condition!(
    optimization_container::OptimizationContainer,
    devices::IS.FlattenIteratorWrapper{T},
    ::D,
) where {T <: PSY.HydroPumpedStorage, D <: AbstractHydroFormulation}
    key_up = ICKey(InitialEnergyLevelUp, T)
    _make_initial_conditions!(
        optimization_container,
        devices,
        D(),
        EnergyVariableUp(),
        key_up,
        _make_initial_condition_reservoir_energy_up,
        _get_variable_initial_value,
        StoredEnergy,
    )

    key_down = ICKey(InitialEnergyLevelDown, T)
    _make_initial_conditions!(
        optimization_container,
        devices,
        D(),
        EnergyVariableDown(),
        key_down,
        _make_initial_condition_reservoir_energy_down,
        _get_variable_initial_value,
        StoredEnergy,
    )

    return
end

########################## Addition to the nodal balances #################################

function NodalExpressionSpec(
    ::Type{T},
    ::Type{<:PM.AbstractPowerModel},
    use_forecasts::Bool,
) where {T <: PSY.HydroGen}
    return NodalExpressionSpec(
        "max_active_power",
        REACTIVE_POWER,
        use_forecasts ? x -> PSY.get_max_reactive_power(x) : x -> PSY.get_reactive_power(x),
        1.0,
        T,
    )
end

function NodalExpressionSpec(
    ::Type{T},
    ::Type{<:PM.AbstractActivePowerModel},
    use_forecasts::Bool,
) where {T <: PSY.HydroGen}
    return NodalExpressionSpec(
        "max_active_power",
        ACTIVE_POWER,
        use_forecasts ? x -> PSY.get_max_active_power(x) : x -> PSY.get_active_power(x),
        1.0,
        T,
    )
end

##################################### Water/Energy Budget Constraint ############################
energy_budget_constraints!(
    ::OptimizationContainer,
    ::IS.FlattenIteratorWrapper{<:PSY.HydroGen},
    ::DeviceModel{<:PSY.HydroGen, <:AbstractHydroFormulation},
    ::Type{<:PM.AbstractPowerModel},
    ::IntegralLimitFF,
) = nothing

"""
This function define the budget constraint for the
active power budget formulation.

`` sum(P[t]) <= Budget ``
"""
function energy_budget_constraints!(
    optimization_container::OptimizationContainer,
    devices::IS.FlattenIteratorWrapper{H},
    ::DeviceModel{H, <:AbstractHydroFormulation},
    ::Type{<:PM.AbstractPowerModel},
    ::Union{Nothing, AbstractAffectFeedForward},
) where {H <: PSY.HydroGen}
    forecast_label = "hydro_budget"
    constraint_data = Vector{DeviceTimeSeriesConstraintInfo}(undef, length(devices))
    for (ix, d) in enumerate(devices)
        ts_vector = get_time_series(optimization_container, d, forecast_label)
        @debug "time_series" ts_vector
        constraint_d =
            DeviceTimeSeriesConstraintInfo(d, x -> PSY.get_storage_capacity(x), ts_vector)
        constraint_data[ix] = constraint_d
    end

    if model_has_parameters(optimization_container)
        device_energy_budget_param_ub(
            optimization_container,
            constraint_data,
            make_constraint_name(ENERGY_BUDGET, H),
            UpdateRef{H}(ENERGY_BUDGET, forecast_label),
            make_variable_name(ACTIVE_POWER, H),
        )
    else
        device_energy_budget_ub(
            optimization_container,
            constraint_data,
            make_constraint_name(ENERGY_BUDGET),
            make_variable_name(ACTIVE_POWER, H),
        )
    end
end

"""
This function define the budget constraint (using params)
for the active power budget formulation.
"""
function device_energy_budget_param_ub(
    optimization_container::OptimizationContainer,
    energy_budget_data::Vector{DeviceTimeSeriesConstraintInfo},
    cons_name::Symbol,
    param_reference::UpdateRef,
    var_names::Symbol,
)
    time_steps = model_time_steps(optimization_container)
    resolution = model_resolution(optimization_container)
    inv_dt = 1.0 / (Dates.value(Dates.Second(resolution)) / SECONDS_IN_HOUR)
    variable_out = get_variable(optimization_container, var_names)
    set_name = [get_component_name(r) for r in energy_budget_data]
    constraint = add_cons_container!(optimization_container, cons_name, set_name)
    container =
        add_param_container!(optimization_container, param_reference, set_name, time_steps)
    multiplier = get_multiplier_array(container)
    param = get_parameter_array(container)
    for constraint_info in energy_budget_data
        name = get_component_name(constraint_info)
        for t in time_steps
            multiplier[name, t] = constraint_info.multiplier * inv_dt
            param[name, t] = add_parameter(
                optimization_container.JuMPmodel,
                constraint_info.timeseries[t],
            )
        end
        constraint[name] = JuMP.@constraint(
            optimization_container.JuMPmodel,
            sum([variable_out[name, t] for t in time_steps]) <= sum([multiplier[name, t] * param[name, t] for t in time_steps])
        )
    end

    return
end

"""
This function define the budget constraint
for the active power budget formulation.
"""
function device_energy_budget_ub(
    optimization_container::OptimizationContainer,
    energy_budget_constraints::Vector{DeviceTimeSeriesConstraintInfo},
    cons_name::Symbol,
    var_names::Symbol,
)
    time_steps = model_time_steps(optimization_container)
    variable_out = get_variable(optimization_container, var_names)
    names = [get_component_name(x) for x in energy_budget_constraints]
    constraint = add_cons_container!(optimization_container, cons_name, names)

    for constraint_info in energy_budget_constraints
        name = get_component_name(constraint_info)
        resolution = model_resolution(optimization_container)
        inv_dt = 1.0 / (Dates.value(Dates.Second(resolution)) / SECONDS_IN_HOUR)
        forecast = constraint_info.timeseries
        multiplier = constraint_info.multiplier * inv_dt
        constraint[name] = JuMP.@constraint(
            optimization_container.JuMPmodel,
            sum([variable_out[name, t] for t in time_steps]) <= multiplier * sum(forecast)
        )
    end

    return
end

##################################### Hydro generation cost ############################
function AddCostSpec(
    ::Type{T},
    ::Type{U},
    ::OptimizationContainer,
) where {T <: PSY.HydroGen, U <: AbstractHydroFormulation}
    # Hydro Generators currently have no OperationalCost
    cost_function = x -> (x === nothing ? 1.0 : PSY.get_variable(x))
    return AddCostSpec(;
        variable_type = ActivePowerVariable,
        component_type = T,
        fixed_cost = PSY.get_fixed,
        variable_cost = cost_function,
        multiplier = OBJECTIVE_FUNCTION_POSITIVE,
    )
end

############################
function AddCostSpec(
    ::Type{T},
    ::Type{U},
    ::OptimizationContainer,
) where {T <: PSY.HydroPumpedStorage, U <: AbstractHydroFormulation}
    # Hydro Generators currently have no OperationalCost
    cost_function = x -> (x === nothing ? 1.0 : PSY.get_variable(x))
    return AddCostSpec(;
        variable_type = ActivePowerOutVariable,
        component_type = T,
        fixed_cost = PSY.get_fixed,
        variable_cost = cost_function,
        multiplier = OBJECTIVE_FUNCTION_POSITIVE,
    )
end


text:
"""
    minimize_metabolic_adjustment_analysis(
        model::MetabolicModel,
        flux_ref::Union{Dict{String,Float64}, Vector{Float64}},
        optimizer;
        modifications = [],
        kwargs...
    )

Run minimization of metabolic adjustment (MOMA) on `model` with respect to
`flux_ref`, which is a vector of fluxes in the order of `reactions(model)`.
MOMA finds the shortest Euclidian distance between `flux_ref` and `model` with
`modifications`:
```
min Σᵢ (xᵢ - flux_refᵢ)²
s.t. S x = b
     xₗ ≤ x ≤ xᵤ
```
Because the problem has a quadratic objective, a QP solver is required. See
"Daniel, Vitkup & Church, Analysis of Optimality in Natural and Perturbed
Metabolic Networks, Proceedings of the National Academy of Sciences, 2002" for
more details.

Additional arguments are passed to [`flux_balance_analysis`](@ref).

Returns an optimized model that contains the resultant nearest flux.

# Example
```
model = load_model("e_coli_core.json")
flux_ref = flux_balance_analysis_vec(model, Gurobi.Optimizer)
optmodel = minimize_metabolic_adjustment(
    model,
    flux_ref,
    Gurobi.Optimizer;
    modifications = [change_constraint("PFL"; lb=0, ub=0)], # find flux of mutant that is closest to the wild type (reference) model
    )
value.(solution[:x])  # extract the flux from the optimizer
```
"""
minimize_metabolic_adjustment_analysis(
    model::MetabolicModel,
    flux_ref::Union{Dict{String,Float64},Vector{Float64}},
    optimizer;
    modifications = [],
    kwargs...,
) = flux_balance_analysis(
    model,
    optimizer;
    modifications = vcat([minimize_metabolic_adjustment(flux_ref)], modifications),
    kwargs...,
)

"""
    minimize_metabolic_adjustment(flux_ref::Vector{Float64})

An optimization model modification that implements the MOMA in
[`minimize_metabolic_adjustment_analysis`](@ref).
"""
minimize_metabolic_adjustment(flux_ref::Vector{Float64}) =
    (model, opt_model) -> begin
        length(opt_model[:x]) == length(flux_ref) || throw(
            DomainError(
                flux_ref,
                "length of the reference flux doesn't match the one in the optimization model",
            ),
        )
        @objective(opt_model, Min, sum((opt_model[:x] .- flux_ref) .^ 2))
    end

"""
    minimize_metabolic_adjustment(flux_ref_dict::Dict{String, Float64})

Overload of [`minimize_metabolic_adjustment`](@ref) that works with a
dictionary of fluxes.
"""
minimize_metabolic_adjustment(flux_ref_dict::Dict{String,Float64}) =
    (model, opt_model) ->
        minimize_metabolic_adjustment([flux_ref_dict[rid] for rid in reactions(model)])(
            model,
            opt_model,
        )

"""
    minimize_metabolic_adjustment_analysis_vec(args...; kwargs...)

Perform minimization of metabolic adjustment (MOMA) and return a vector of fluxes in the
same order as the reactions in `model`. Arguments are forwarded to
[`minimize_metabolic_adjustment`](@ref) internally.

This function is kept for backwards compatibility, use [`flux_vector`](@ref)
instead.
"""
minimize_metabolic_adjustment_analysis_vec(args...; kwargs...) =
    flux_vector(minimize_metabolic_adjustment_analysis(args...; kwargs...))

"""
    minimize_metabolic_adjustment_analysis_dict(model::MetabolicModel, args...; kwargs...)

Perform minimization of metabolic adjustment (MOMA) and return a dictionary mapping the
reaction IDs to fluxes. Arguments are forwarded to [`minimize_metabolic_adjustment`](@ref)
internally.

This function is kept for backwards compatibility, use [`flux_vector`](@ref)
instead.
"""
minimize_metabolic_adjustment_analysis_dict(model::MetabolicModel, args...; kwargs...) =
    flux_dict(model, minimize_metabolic_adjustment_analysis(model, args...; kwargs...))


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#    Coordinate transformations related with the geodetic and geocentric
#    coordinates.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
#
#   [1] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
#       Microcosm Press, Hawthorn, CA, USA.
#
#   [2] ESA Navipedia: http://www.navipedia.net/
#
#   [3] mu-blox ag (1999). Datum Transformations of GPS Positions. Application
#       Note.
#
#   [4] ISO TC 20/SC 14 N (2011). Geomagnetic Reference Models.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

export ECEFtoGeodetic, GeodetictoECEF

export GeodetictoGeocentric

"""
    ECEFtoGeodetic(r_e::AbstractVector)

Convert the vector `r_e` [m] represented in the Earth-Centered, Earth-Fixed
(ECEF) reference frame into Geodetic coordinates (WGS-84).

# Returns

* Latitude [rad].
* Longitude [rad].
* Altitude [m].

# Remarks

Based on algorithm in [3].

"""
function ECEFtoGeodetic(r_e::AbstractVector)
    # Auxiliary variables.
    X = r_e[1]
    Y = r_e[2]
    Z = r_e[3]

    # Auxiliary variables.
    p = sqrt(X^2 + Y^2)
    θ = atan(Z*a_wgs84, p*b_wgs84)
    e_wgs84² = e_wgs84^2

    # Compute Geodetic.
    lon = atan(Y, X)
    lat = atan(Z + el_wgs84^2*b_wgs84*sin(θ)^3,
               p -   e_wgs84²*a_wgs84*cos(θ)^3)

    sin_lat, cos_lat = sincos(lat)

    N = a_wgs84/sqrt(1 - e_wgs84²*sin_lat^2)

    # Avoid singularity if we are near the poles (~ 1 deg according to [1,
    # p.172]). Note that `cosd(1) = -0.01745240643728351`.
    if !(-0.01745240643728351 < cos_lat < 0.01745240643728351)
        h = p/cos_lat - N
    else
        h = Z/sin_lat - N*(1 - e_wgs84²)
    end

    (lat, lon, h)
end

"""
    GeodetictoECEF(lat::Number, lon::Number, h::Number)

Convert the latitude `lat` [rad], longitude `lon` [rad], and altitude `h` \\[m]
(WGS-84) into a vector represented on the Earth-Centered, Earth-Fixed (ECEF)
reference frame.

# Remarks

Based on algorithm in [3].

"""
function GeodetictoECEF(lat::Number, lon::Number, h::Number)
    # Auxiliary variables.
    sin_lat, cos_lat = sincos(lat)
    sin_lon, cos_lon = sincos(lon)

    # Radius of curvature [m].
    N = a_wgs84/sqrt(1 - e_wgs84^2*sin_lat^2 )

    # Compute the position in ECEF frame.
    SVector{3}((                     N + h)*cos_lat*cos_lon,
               (                     N + h)*cos_lat*sin_lon,
               ( (b_wgs84/a_wgs84)^2*N + h)*sin_lat)
end

"""
    GeodetictoGeocentric(ϕ_gd::Number, h::Number)

Compute the geocentric latitude and radius from the geodetic latitude `ϕ_gd`
(-π/2,π/2) \\[rad] and height above the reference ellipsoid `h` \\[m] (WGS-84).
Notice that the longitude is the same in both geocentric and geodetic
coordinates.

# Returns

* Geocentric latitude [rad].
* Radius from the center of the Earth [m].

# Remarks

Based on algorithm in [4, p. 3].

"""
function GeodetictoGeocentric(ϕ_gd::Number, h::Number)
    # Auxiliary variables to decrease computational burden.
    sin_ϕ_gd  = sin(ϕ_gd)
    sin2_ϕ_gd = sin_ϕ_gd^2
    e2_wgs84  = e_wgs84^2

    # Radius of curvature in the prime vertical [m].
    N    = a_wgs84/sqrt(1 - e2_wgs84*sin2_ϕ_gd )

    # Compute the geocentric latitude and radius from the Earth center.
    ρ    = (N + h)*cos(ϕ_gd)
    z    = (N * (1 - e2_wgs84) + h)*sin_ϕ_gd
    r    = sqrt(ρ^2 + z^2)
    ϕ_gc = asin(z/r)

    (ϕ_gc, r)
end


text:
using Test
using LinearAlgebra
using DataFrames
# Preprocess Tests:
include("preprocesstests.jl")
# Stats Tests:
include("stattests.jl")
# Models Tets:
include("modeltests.jl")