text:
# select: select particular columns
# df |> select(:Chrom, :Pos, :Val)

# could be simplified by using
# NamedTuple{(columns to select)}(intuple)

export select

using Tables

# I is type of rows iterator
# O is output eltype
#
# having columns to select in type signature allows us to use
# @generated iterate function

struct Select{I,O}
    rows::I
end

"""
    select(rows, columns::Symbol...)
    rows |> select(columns::Symbol...)

Select `columns` from genome ordered stream `rows`.
"""
function select(rows, columns::Symbol...)
    sch = Tables.schema(rows)
    col2type = Dict(zip(sch.names, sch.types))
    outtypes = [col2type[a] for a in columns]
    O = NamedTuple{Tuple(columns), Tuple{outtypes...}}

    Select{typeof(rows), O}(rows)
end

select(columns::Symbol...) = rows -> select(rows, columns...)

Base.IteratorEltype(::Type{Select{I,O}}) where {I,O} = Base.HasEltype()
Base.eltype(::Type{Select{I,O}}) where {I,O} = O
Base.IteratorSize(::Type{Select{I,O}}) where {I,O} = Base.IteratorSize(I)
Base.length(r::Select{I,O}) where {I,O} = Base.length(r.rows)
Base.size(r::Select{I,O}) where {I,O} = Base.size(r.rows)

Tables.istable(r::Select{I,O}) where {I,O} = true
Tables.rowaccess(r::Select{I,O}) where {I,O} = true
Tables.schema(r::Select{I,O}) where {I,O} = Tables.Schema(O)

function Base.iterate(s::Select{I,O}) where {I,O}
    elt_s = iterate(s.rows)
    iterate(s, elt_s)
end

@generated function Base.iterate(s::Select{I,O}, state) where {I,O}
    # want to have
    # O(Tuple([elt[:col1], elt[:col2],...]))
    #
    # splatting only works as argument
    # :($(vals...)) creates error
    # but :($(vals...),)  or :(Tuple($(vals...))) works

    vals = []
    for n in O.parameters[1]
        push!(vals, :(elt[$(Meta.quot(n))]))
    end

    quote
        elt_s = state
        elt_s == nothing && return nothing
        elt, is = elt_s
        O(($(vals...),)), iterate(s.rows, is)
    end
end

# should we just do a more powerful select?
#
# @select :x, :y, c = 2*:y, _
#
# where _ is a fill-in for everything
# every symbol mentioned explicitly is excluded from _

#
# use macro first to create function for all specified cols
# then pass names (:x, :y, :c, :_) and  (x, y, c) = func(row)
#

# general strategy to process command
#
# @select :a, b = 23 - :b, _, d = 12
#
# 1. process AST using macro stage
#
#    (:a, :b, :_, :d), row -> (row[:a], 23 - row[:b], 12)
#    names           , func
#
#    macro returns function
#    rows -> pass2(rows, names, func)
#
# 2. type inference
#
#    rsch = Tables.schema(rows)
#    rtypes = rsch.types
#    rnames = rsch.names
#    ftypes = Base.return_types(func, (eltype(rows),))[1]
#    fnames = [x for x in names if x != :(_)]
#    ridx = [i for i = 1:length(rnames) if ~(rnames[i] in names)]
#    _names = rnames[ridx]
#    _types = rtypes[idx]
#
#    onames = []
#    otypes = []
#    idx = 1
#    for n in names
#      if n == :(_)
#        # insert _names, _types
#      else
#        push!(names, n)
#        push!(types, ftypes[i])
#        i += 1
#      end
#    end
#
#    pass2 returns a structure that is parameterized on
#    input iter type, fnames, onames, otypes
#
# 3. generated function
#
#    fvals = func(row)
#    name2val = onames => row[onames]
#    name2val[fnames] = fvals[fnameidx]
#

# NOTE: step 2 and 3 are basically mutate in mutate.jl


text:
using Test
using SPGBox
using Unitful
using ReverseDiff

@testset "simple quadratic" begin

    f(x) = x[1]^2 + (x[2] - 1)^2

    function g!(g, x)
        g[1] = 2 * x[1]
        g[2] = 2 * (x[2] - 1)
    end

    x = [10.0, 18.0]
    R = spgbox!(f, g!, x)
    @test R.f ≈ 0.0
    @test R.x ≈ [0.0, 1.0]

    x = [10.0, 18.0]
    R = spgbox!(f, g!, x, lower = [2, -Inf])
    @test R.f ≈ 4.0
    @test R.x ≈ [2.0, 1.0]

    x = [10.0, 18.0]
    R = spgbox!(f, g!, x, lower = [-Inf, 2.0])
    @test R.f ≈ 1.0
    @test R.x ≈ [0.0, 2.0]

    x = [10.0, 18.0]
    R = spgbox!(f, g!, x, upper = [+Inf, -2])
    @test R.f ≈ 9.0
    @test R.x ≈ [0.0, -2.0]

    x = [10.0, 18.0]
    R = spgbox!(f, g!, x, upper = [-2, +Inf])
    @test R.f ≈ 4.0
    @test R.x ≈ [-2.0, 1.0]

    # just testing the interface
    x = [10.0, 18.0]
    R = spgbox!(f, g!, x, lower = [-Inf, 2.0])
    @test R.f ≈ 1.0
    @test R.x ≈ [0.0, 2.0]

    x = [10.0, 18.0]
    R = spgbox!(f, g!, x, upper = [-2, +Inf])
    @test R.f ≈ 4.0
    @test R.x ≈ [-2.0, 1.0]

    x = [10.0, 18.0]
    R = spgbox!(f, g!, x, lower = [-Inf, 2.0], upper = [-2.0, +Inf])
    @test R.f ≈ 5.0
    @test R.x ≈ [-2.0, 2.0]

    x = [10.0, 18.0]
    R = spgbox!(f, g!, [-Inf, -2.0], [+Inf, 2], x)
    @test R.f ≈ 0.0
    @test R.x ≈ [0.0, 1.0]

    # Test the mutating call
    x = [10.0, 18.0]
    R = spgbox(f, g!, [-Inf, -2.0], [+Inf, 2], x)
    @test x == [10.0, 18.0]
    @test R.f ≈ 0.0
    @test R.x ≈ [0.0, 1.0]

    x = [10.0, 18.0]
    R = spgbox(f, g!, x, lower = [-Inf, -2.0], upper = [+Inf, 2])
    @test x == [10.0, 18.0]
    @test R.f ≈ 0.0
    @test R.x ≈ [0.0, 1.0]

    #BigFloat support
    x = BigFloat[10.0, 18.0]
    R = spgbox!(f, g!, x)
    @test R.f ≈ 0.0
    @test R.x ≈ [0.0, 1.0]
    @test eltype(R.f) == eltype(R.x) == BigFloat

    #Float32 support
    x = Float32[10.0, 18.0]
    R = spgbox!(f, g!, x)
    @test R.f ≈ 0.0
    @test R.x ≈ [0.0, 1.0]
    @test eltype(R.f) == eltype(R.x) == Float32

    # Matrix input support
    x = [ 10.0 18.0 ]
    R = spgbox!(f, g!, x)
    @test R.f ≈ 0.0
    @test R.x ≈ [0.0 1.0]

    # Unitful support
    f_units(x) = x[1]^2 + (x[2] - oneunit(eltype(x)))^2
    function g_units!(g, x)
        g[1] = 2 * x[1]
        g[2] = 2 * (x[2] - oneunit(eltype(x)))
    end
    x = [10.0, 18.0]u"nm"
    R = spgbox!(f_units, g_units!, x)
    @test R.f ≈ 0.0u"nm^2"
    @test R.x ≈ [0.0, 1.0]u"nm"
    x = [10.0, 18.0]u"nm"
    R = spgbox!(f_units, g_units!, x, lower = [-Inf, 2.0]u"nm")
    @test R.f ≈ 1.0u"nm^2"
    @test R.x ≈ [0.0, 2.0]u"nm"
    R = spgbox!(f_units, g_units!, x, lower = [-Inf, 2.0]u"nm")
    @test R.f ≈ 1.0u"nm^2"
    @test R.x ≈ [0.0, 2.0]u"nm"
    R = spgbox!(f_units, g_units!, x, upper = [-5.0, +Inf]u"nm")
    @test R.f ≈ 25.0u"nm^2"
    @test R.x ≈ [-5.0, 1.0]u"nm"
    R = spgbox!(f_units, g_units!, x, lower = [-Inf, 2.0]u"nm", upper = [-5.0, +Inf]u"nm")
    @test R.f ≈ 26.0u"nm^2"
    @test R.x ≈ [-5.0, 2.0]u"nm"

    #
    # With a single function to compute function and gradient
    #

    function fg!(g,x)
        g[1] = 2 * x[1]
        g[2] = 2 * (x[2] - 1)
        fx = x[1]^2 + (x[2] - 1)^2
        return fx
    end

    x = [10.0, 18.0]
    R = spgbox!(fg!, x)
    @test R.f ≈ 0.0
    @test R.x ≈ [0.0, 1.0]

    x = [10.0, 18.0]
    R = spgbox!(fg!, x, lower = [-Inf, 2.0], upper = [-2.0, +Inf])
    @test R.f ≈ 5.0
    @test R.x ≈ [-2.0, 2.0]

    x = [10.0, 18.0]
    R = spgbox(fg!, x)
    @test R.f ≈ 0.0
    @test R.x ≈ [0.0, 1.0]

    x = [10.0, 18.0]
    R = spgbox!(fg!, x, lower = [-Inf, 2.0], upper = [-2.0, +Inf])
    @test R.f ≈ 5.0
    @test R.x ≈ [-2.0, 2.0]

    #
    # An example that verify that the line search was fixed
    #
    x0 = [i / 9 for i = 1:8]
    fH368(x) = -1.0*sum(x.^2)*sum(x.^4) + sum(x.^3)^2;
    function gH368!(g, x)
        ReverseDiff.gradient!(g, fH368, x)
    end
    lower = zeros(8)
    upper = ones(8)
    R = spgbox!(fH368, gH368!, x0, lower = lower, upper = upper)
    g = similar(x0)
    gH368!(g, R.x)
    @test SPGBox.pr_gradnorm(g, R.x, lower, upper) <= 1.0e-5

end

ext: jl


text:
# MIT License

# Copyright (c) Microsoft Corporation.

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

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

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

module GlassCat

using Polynomials
using Plots
using StringEncodings
using Unitful
using StaticArrays
using Base: @.
import Unitful: Length, Temperature, Quantity, Units

include("constants.jl")

include("GlassTypes.jl")
export GlassID, info, glassid, glassname, glassforid
include("Air.jl")
export Air, isair

# include built glass cat source files
include("data/AGFGlassCat.jl")
include("data/CARGILLE.jl")

# include functionality for managing runtime (dynamic) glass cats: MIL_GLASSES and MODEL_GLASSES
include("runtime.jl")
export glassfromMIL, modelglass

# include functions for searching the glass cats
include("search.jl")
export glasscatalogs, glassnames, findglass

include("utilities.jl")
export plot_indices, index, polyfit_indices, absairindex, absorption

end # module
export GlassCat


text:
"""
A statically sized `Int`.
Use `StaticInt(N)` instead of `Val(N)` when you want it to behave like a number.
"""
struct StaticInt{N} <: Integer
    StaticInt{N}() where {N} = new{N::Int}()
end

const Zero = StaticInt{0}
const One = StaticInt{1}

Base.@pure StaticInt(N::Int) = StaticInt{N}()
StaticInt(N::Integer) = StaticInt(convert(Int, N))
StaticInt(::StaticInt{N}) where {N} = StaticInt{N}()
StaticInt(::Val{N}) where {N} = StaticInt{N}()
# Base.Val(::StaticInt{N}) where {N} = Val{N}()
Base.convert(::Type{T}, ::StaticInt{N}) where {T<:Number,N} = convert(T, N)
Base.Bool(x::StaticInt{N}) where {N} = Bool(N)
Base.BigInt(x::StaticInt{N}) where {N} = BigInt(N)
Base.Integer(x::StaticInt{N}) where {N} = x
(::Type{T})(x::StaticInt{N}) where {T<:Integer,N} = T(N)
(::Type{T})(x::Int) where {T<:StaticInt} = StaticInt(x)
Base.convert(::Type{StaticInt{N}}, ::StaticInt{N}) where {N} = StaticInt{N}()

Base.promote_rule(::Type{<:StaticInt}, ::Type{T}) where {T <: AbstractIrrational} = promote_rule(Int, T)
Base.promote_rule(::Type{T}, ::Type{<:StaticInt}) where {T <: AbstractIrrational} = promote_rule(T, Int)
for (S,T) ∈ [(:Complex,:Real), (:Rational, :Integer), (:(Base.TwicePrecision),:Any)]
    @eval Base.promote_rule(::Type{$S{T}}, ::Type{<:StaticInt}) where {T <: $T} = promote_rule($S{T}, Int)
end
Base.promote_rule(::Type{Union{Nothing,Missing}}, ::Type{<:StaticInt}) = Union{Nothing, Missing, Int}
Base.promote_rule(::Type{T}, ::Type{<:StaticInt}) where {T >: Union{Missing,Nothing}} = promote_rule(T, Int)
Base.promote_rule(::Type{T}, ::Type{<:StaticInt}) where {T >: Nothing} = promote_rule(T, Int)
Base.promote_rule(::Type{T}, ::Type{<:StaticInt}) where {T >: Missing} = promote_rule(T, Int)
for T ∈ [:Bool, :Missing, :BigFloat, :BigInt, :Nothing, :Any]
# let S = :Any
    @eval begin
        Base.promote_rule(::Type{S}, ::Type{$T}) where {S <: StaticInt} = promote_rule(Int, $T)
        Base.promote_rule(::Type{$T}, ::Type{S}) where {S <: StaticInt} = promote_rule($T, Int)
    end
end
Base.promote_rule(::Type{<:StaticInt}, ::Type{<:StaticInt}) = Int
Base.:(%)(::StaticInt{N}, ::Type{Integer}) where {N} = N

Base.eltype(::Type{T}) where {T<:StaticInt} = Int
Base.iszero(::Zero) = true
Base.iszero(::StaticInt) = false
Base.isone(::One) = true
Base.isone(::StaticInt) = false
Base.zero(::Type{T}) where {T<:StaticInt} = Zero()
Base.one(::Type{T}) where {T<:StaticInt} = One()

for T = [:Real, :Rational, :Integer]
    @eval begin
        @inline Base.:(+)(i::$T, ::Zero) = i
        @inline Base.:(+)(i::$T, ::StaticInt{M}) where {M} = i + M
        @inline Base.:(+)(::Zero, i::$T) = i
        @inline Base.:(+)(::StaticInt{M}, i::$T) where {M} = M + i
        @inline Base.:(-)(i::$T, ::Zero) = i
        @inline Base.:(-)(i::$T, ::StaticInt{M}) where {M} = i - M
        @inline Base.:(*)(i::$T, ::Zero) = Zero()
        @inline Base.:(*)(i::$T, ::One) = i
        @inline Base.:(*)(i::$T, ::StaticInt{M}) where {M} = i * M
        @inline Base.:(*)(::Zero, i::$T) = Zero()
        @inline Base.:(*)(::One, i::$T) = i
        @inline Base.:(*)(::StaticInt{M}, i::$T) where {M} = M * i
    end
end
@inline Base.:(+)(::Zero, ::Zero) = Zero()
@inline Base.:(+)(::Zero, ::StaticInt{M}) where {M} = StaticInt{M}()
@inline Base.:(+)(::StaticInt{M}, ::Zero) where {M} = StaticInt{M}()

@inline Base.:(-)(::StaticInt{M}) where {M} = StaticInt{-M}()
@inline Base.:(-)(::StaticInt{M}, ::Zero) where {M} = StaticInt{M}()

@inline Base.:(*)(::Zero, ::Zero) = Zero()
@inline Base.:(*)(::One, ::Zero) = Zero()
@inline Base.:(*)(::Zero, ::One) = Zero()
@inline Base.:(*)(::One, ::One) = One()
@inline Base.:(*)(::StaticInt{M}, ::Zero) where {M} = Zero()
@inline Base.:(*)(::Zero, ::StaticInt{M}) where {M} = Zero()
@inline Base.:(*)(::StaticInt{M}, ::One) where {M} = StaticInt{M}()
@inline Base.:(*)(::One, ::StaticInt{M}) where {M} = StaticInt{M}()
for f ∈ [:(+), :(-), :(*), :(/), :(÷), :(%), :(<<), :(>>), :(>>>), :(&), :(|), :(⊻)]
    @eval @generated Base.$f(::StaticInt{M}, ::StaticInt{N}) where {M,N} = Expr(:call, Expr(:curly, :StaticInt, $f(M, N)))
end
for f ∈ [:(==), :(!=), :(<), :(≤), :(>), :(≥)]
    @eval begin
        @inline Base.$f(::StaticInt{M}, ::StaticInt{N}) where {M,N} = $f(M, N)
        @inline Base.$f(::StaticInt{M}, x::Int) where {M} = $f(M, x)
        @inline Base.$f(x::Int, ::StaticInt{M}) where {M} = $f(x, M)
    end
end

@inline function maybe_static(f::F, g::G, x) where {F, G}
    L = f(x)
    isnothing(L) ? g(x) : StaticInt(L)
end
@inline static_length(x) = maybe_static(known_length, length, x)
@inline static_first(x) = maybe_static(known_first, first, x)
@inline static_last(x) = maybe_static(known_last, last, x)
@inline static_step(x) = maybe_static(known_step, step, x)

Base.UnitRange{T}(start::StaticInt, stop) where {T<:Real} = UnitRange{T}(T(start), stop)
Base.UnitRange{T}(start, stop::StaticInt) where {T<:Real} = UnitRange{T}(start, T(stop))
function Base.UnitRange{T}(start::StaticInt, stop::StaticInt) where {T<:Real}
    return UnitRange{T}(T(start), T(stop))
end

Base.UnitRange(start::StaticInt, stop) = UnitRange(Int(start), stop)
Base.UnitRange(start, stop::StaticInt) = UnitRange(start, Int(stop))
function Base.UnitRange(start::StaticInt, stop::StaticInt)
    return UnitRange(Int(start), Int(stop))
end


text:
# This file was generated, do not modify it. # hide
using Gnuplot
x = y = -15:0.33:15
fz(x,y) = sin.(sqrt.(x.^2 + y.^2))./sqrt.(x.^2+y.^2)
fxy = [fz(x,y) for x in x, y in y]
@gsp x y fxy "w pm3d" "set ylabel 'y' textcolor rgb 'red'"
@gsp :- "set xlabel 'x' textcolor rgb 'blue'"
save(term="pngcairo size 800,800", output="./code/plt3d_ex3_1.png") # hide

max_stars_repo_name: lazarusA/gnuplot-examples

max_issues_repo_name: lazarusA/gnuplot-examples

max_forks_repo_name: lazarusA/gnuplot-examples


text:
@testset "M" begin
    @testset "m2eul" begin
        ticam = [0.49127379678135830 0.50872620321864170 0.70699908539882417;
                 -0.50872620321864193 -0.49127379678135802 0.70699908539882428;
                 0.70699908539882406 -0.70699908539882439 0.01745240643728360]
        kappa, ang2, ang1 = m2eul(ticam, 3, 1, 3)
        alpha = ang1 + 1.5 * π
        delta = π/2 - ang2
        expected = [315.000000, 1.000000, 45.000000]
        actual = rad2deg.([alpha, delta, kappa])
        @testset for i in eachindex(actual, expected)
            @test actual[i] ≈ expected[i]
        end
        @test_throws ArgumentError m2eul(randn(3,4), 1, 2, 3)
        @test_throws SpiceError m2eul(randn(3,3), 1, 1, 1)
    end
    @testset "matchi" begin
        string = "  ABCDEFGHIJKLMNOPQRSTUVWXYZ  "
        wstr = "*"
        wchr = "%"
        @test matchi(string, "*A*", wstr, wchr)
        @test !matchi(string, "A%D*", wstr, wchr)
        @test matchi(string, "A%C*", wstr, wchr)
        @test !matchi(string, "%A*", wstr, wchr)
        @test matchi(string, "%%CD*Z", wstr, wchr)
        @test !matchi(string, "%%CD", wstr, wchr)
        @test matchi(string, "A*MN*Y*Z", wstr, wchr)
        @test !matchi(string, "A*MN*Y*%Z", wstr, wchr)
        @test matchi(string, "*BCD*Z*", wstr, wchr)
        @test !matchi(string, "*bdc*z*", wstr, wchr)
        @test matchi(string, " *bcD*Z*", wstr, wchr)
        @test_throws SpiceError matchi("", "", wstr, wchr)
    end
    @testset "matchw" begin
        string = "  ABCDEFGHIJKLMNOPQRSTUVWXYZ  "
        wstr = '*'
        wchr = '%'
        @test matchw(string, "*A*", wstr, wchr)
        @test !matchw(string, "A%D*", wstr, wchr)
        @test matchw(string, "A%C*", wstr, wchr)
        @test !matchw(string, "%A*", wstr, wchr)
        @test matchw(string, "%%CD*Z", wstr, wchr)
        @test !matchw(string, "%%CD", wstr, wchr)
        @test matchw(string, "A*MN*Y*Z", wstr, wchr)
        @test !matchw(string, "A*MN*Y*%Z", wstr, wchr)
        @test matchw(string, "*BCD*Z*", wstr, wchr)
        @test !matchw(string, "*bdc*z*", wstr, wchr)
        @test matchw(string, " *BCD*Z*", wstr, wchr)
        @test_throws SpiceError matchw("", "", wstr, wchr)
    end
    @testset "maxd" begin
        values = [4.0, 3.0, 5.0]
        @test SPICE._maxd(values...) == max(values...)
    end
    @testset "maxi" begin
        values = [4, 3, 5]
        @test SPICE._maxi(values...) == max(values...)
    end
    @testset "mequ" begin
        a1 = randn(3, 3)
        a2 = copy(a1)
        b1 = randn(3, 3)
        b2 = copy(b1)
        SPICE._mequ(a1, b1)
        b2 .= a2
        @test b1 == b2
    end
    @testset "mequg" begin
        a1 = randn(12, 15)
        a2 = copy(a1)
        b1 = randn(12, 15)
        b2 = copy(b1)
        SPICE._mequg(a1, b1)
        b2 .= a2
        @test b1 == b2
    end
    @testset "mind" begin
        values = [4.0, 3.0, 5.0]
        @test SPICE._mind(values...) == min(values...)
    end
    @testset "mini" begin
        values = [4, 3, 5]
        @test SPICE._mini(values...) == min(values...)
    end
    @testset "mtxm" begin
        a = randn(3, 3)
        b = randn(3, 3)
        @test SPICE._mtxm(a, b) ≈ a' * b
    end
    @testset "mtxmg" begin
        a = randn(2, 4)
        b = randn(2, 3)
        @test SPICE._mtxmg(a, b) ≈ a' * b
    end
    @testset "mtxv" begin
        a = randn(3, 3)
        v = randn(3)
        @test SPICE._mtxv(a, v) ≈ a' * v
    end
    @testset "mtxvg" begin
        a = randn(3, 2)
        v = randn(3)
        @test SPICE._mtxvg(a, v) ≈ a' * v
    end
    @testset "mxm" begin
        a = randn(3, 3)
        b = randn(3, 3)
        @test SPICE._mxm(a, b) ≈ a * b
    end
    @testset "mxmg" begin
        a = randn(4, 2)
        b = randn(2, 3)
        @test SPICE._mxmg(a, b) ≈ a * b
    end
    @testset "mxmt" begin
        a = randn(3, 3)
        b = randn(3, 3)
        @test SPICE._mxmt(a, b) ≈ a * b'
    end
    @testset "mxmtg" begin
        a = randn(4, 2)
        b = randn(3, 2)
        @test SPICE._mxmtg(a, b) ≈ a * b'
    end
    @testset "mxv" begin
        a = randn(3, 3)
        v = randn(3)
        @test SPICE._mxv(a, v) ≈ a * v
    end
    @testset "mxvg" begin
        a = randn(2, 3)
        v = randn(3)
        @test SPICE._mxvg(a, v) ≈ a * v
    end
    @testset "m2q" begin
        r = rotate(π/2, 3)
        act = m2q(r)
        exp = [sqrt(2)/2.0, 0.0, 0.0, -sqrt(2)/2.0]
        @testset for i in eachindex(act, exp)
            @test act[i] ≈ exp[i]
        end
        @test_throws ArgumentError m2q(randn(3,4))
        @test_throws SpiceError m2q(fill(1.0, 3, 3))
    end
end


text:
# Used by BuildExecutable.jl to build a self-contained executable
import Playground
main() = Playground.main(ARGS)


text:
# # Implement this example
# https://theoreticalecology.wordpress.com/2010/09/17/metropolis-hastings-mcmc-in-r/
using Distributions
using StatsBase
using Parameters
using DynamicHMC
using ContinuousTransformations
using MCMCDiagnostics
using DiffWrappers
using PyPlot

# # Creating test data
true_a = 5.0
true_b = 0.0
true_sd = 10.0
sample_size = 31

# create independent x-values
x = -(sample_size - 1) / 2 : (sample_size - 1) / 2
# create dependent values according to ax + b + N(0, sd)
y = true_a .* x + true_b + rand(Normal(0, true_sd), sample_size)
plot(x, y, "o")

struct LinRegProblem{T}
    x::Vector{T}
end

function likelihood(prob::LinRegProblem, θ)
    a, b, sd = θ
    return sum(logpdf(Normal(a * prob.x[i] + b, sd), y[i]) for i in 1:length(y))
end

function priors(θ)
    a, b, sd = θ
    return logpdf(Uniform(0, 10), a) + logpdf(Normal(0, 5), b) + logpdf(Uniform(0, 30), sd)
end

prob = LinRegProblem(collect(x))
likelihood(prob, [true_a, true_b, true_sd])

slope_values(x) = likelihood(prob, [x, true_b, true_sd])
slope_likelihoods = slope_values.(3:0.05:7)
begin
    plot(3:0.05:7, slope_likelihoods)
    xlabel("values of slope parameter a")
    ylabel("Log Likelihood")
end

"posterier of `LinRegProblem`"
function (prob::LinRegProblem)(θ)
    return likelihood(prob, θ) + priors(θ)
end

prob([true_a, true_b, true_sd])

# for the priors I will want:
# in the domain for the Uniform priors
# All of ℝ for the Normal
θ_transform = TransformationTuple((bridge(ℝ, ℝ), bridge(ℝ, ℝ), bridge(ℝ, ℝ⁺)))
tprob = TransformLogLikelihood(prob, θ_transform)
tprob([true_a, true_b, true_sd])

tprob∇ = ForwardGradientWrapper(tprob, fill(0.0, length(tprob)))
tprob∇(ones(length(tprob∇)))

@time samp, NUTS_tuned = NUTS_init_tune_mcmc(tprob∇, fill(0.0, length(tprob)), 20000)

NUTS_statistics(samp)
ESS = squeeze(mapslices(effective_sample_size, ungrouping_map(Array, get_position, samp), 1), 1)

# Note the use of ungrouping_map, a utility function that maps the collects posterior results,
# and groups then in tuples of vectors or arrays. You could do this manually, but you have to
# use get_position to extract the posterior position for each point. Also, in order to do inference,
# you would want to transform the "raw" values in $\mathbb{R}^n$ using the parameter transformation.
(a, b, sd) = ungrouping_map(Array, θ_transform ∘ get_position, samp)

# Posterior means are fairly close to the parameters (note that they will not be the same,
# because of sampling variation):
mean(a), true_a
mean(b), true_b
mean(sd), true_sd

median(a), true_a
median(b), true_b
median(sd), true_sd

linreg(x, y)

# # Plot some results
# ## Posterior
function plot_posterior(known_value, posterior_draws, varname; nbins = 20)
    #postr = fit(Histogram, posterior_draws; closed = :right, args...)
    plt[:hist](posterior_draws, bins = nbins, label = "posterior")
    xlabel(varname)
    axvline([known_value], label = "known value", c = "red")
    legend()
end

begin
    figure()
    subplot(3, 1, 1)
    plot_posterior(true_a, a, L"a"; nbins = 20)
    subplot(3, 1, 2)
    plot_posterior(true_b, b, L"b"; nbins = 20)
    subplot(3, 1, 3)
    plot_posterior(true_sd, sd, L"\sigma"; nbins = 20)
    tight_layout()
end

# ## MCMC sampling
begin
    figure()
    subplot(3, 1, 1)
    plot(a, label = L"a")
    axhline([true_a], c = "red")
    legend()
    subplot(3, 1, 2)
    plot(b, label = L"b")
    axhline([true_b], c = "red")
    legend()
    subplot(3, 1, 3)
    plot(sd, label = L"\sigma")
    axhline([true_sd], c = "red")
    legend()
    tight_layout()
end

max_stars_repo_name: gabrielgellner/LearnHMC

max_issues_repo_name: gabrielgellner/LearnHMC

max_forks_repo_name: gabrielgellner/LearnHMC


text:
using Oceananigans
using Oceananigans.Units

using Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities:
    CATKEVerticalDiffusivity,
    SurfaceTKEFlux,
    MixingLength

using OceanTurbulenceParameterEstimation.Observations: FieldTimeSeriesCollector
using JLD2
using GLMakie
using Printf
using Statistics

#dir = "two_day_suite/highres"
dir = "six_day_suite"
filename = "weak_wind_strong_cooling"
xy_filepath = joinpath(dir, filename * "_xy_slice.jld2")

cases = [
         "strong_wind_no_rotation",
         "free_convection",
         "weak_wind_strong_cooling",
         "strong_wind_weak_cooling",
         "strong_wind",
        ]

k_case = findfirst(c -> c == filename, cases)

#####
##### Metadata
#####

xy_file = jldopen(xy_filepath)

Nx = xy_file["grid/Nx"]
Ny = xy_file["grid/Ny"]
Nz = xy_file["grid/Nz"]
Hx = xy_file["grid/Hx"]
Hy = xy_file["grid/Hy"]
Hz = xy_file["grid/Hz"]
Lx = xy_file["grid/Lx"]
Ly = xy_file["grid/Ly"]
Lz = xy_file["grid/Lz"]

iterations = parse.(Int, keys(xy_file["timeseries/t"]))
times = [xy_file["timeseries/t/$i"] for i in iterations]
@show Nt = length(times)

Qᵀ = xy_file["parameters/temperature_flux"]
Qᵘ = xy_file["parameters/momentum_flux"]
f  = xy_file["parameters/coriolis_parameter"]
α  = xy_file["parameters/thermal_expansion_coefficient"]
f  = xy_file["parameters/coriolis_parameter"]
dTdz = xy_file["parameters/dθdz_deep"]

close(xy_file)

chop(a::AbstractArray{<:Any, 1}, H=Hz) = a[1+H:end-H]
chop(a::AbstractArray{<:Any, 2}, H=Hz) = a[1+H:end-H, 1+H:end-H]

function extract_slices(dir, filename; name="w")
    xz_filepath = joinpath(dir, filename * "_xz_slice.jld2")
    yz_filepath = joinpath(dir, filename * "_yz_slice.jld2")
    xy_filepath = joinpath(dir, filename * "_xy_slice.jld2")
    statistics_filepath = joinpath(dir, filename * "_instantaneous_statistics.jld2")
    catke_filepath = "catke_simulation$(dir).jld2"

    all_iterations = []
    slices = []
    for (filepath, dims) in zip((yz_filepath, xz_filepath, xy_filepath), (1, 2, 3))
        file = jldopen(filepath)
        iterations = parse.(Int, keys(file["timeseries/t"]))
        push!(slices, [chop(dropdims(file["timeseries/$name/$i"]; dims)) for i in iterations])
        push!(all_iterations, iterations)
        close(file)
    end

    file = jldopen(statistics_filepath)
    iterations = parse.(Int, keys(file["timeseries/t"]))
    push!(all_iterations, iterations)
    statistics = [chop(file["timeseries/$name/$i"][1, 1, :]) for i in iterations]
    close(file)

    file = jldopen(catke_filepath)
    iterations = parse.(Int, keys(file["timeseries/t"]))
    push!(all_iterations, iterations)
    catke_conv_adj    = [chop(file["timeseries/$name/$i"][1, k_case, :], 0) for i in iterations]
    catke_constant_Pr = [chop(file["timeseries/$name/$i"][2, k_case, :], 0) for i in iterations]
    catke_variable_Pr = [chop(file["timeseries/$name/$i"][3, k_case, :], 0) for i in iterations]
    close(file)

    return (; yz=slices[1], xz=slices[2], xy=slices[3], statistics, all_iterations,
            catke_conv_adj, catke_constant_Pr, catke_variable_Pr)
end

#####
##### 3D visualization
#####

grid = RectilinearGrid(size=(Nx, Ny, Nz), halo=(Hx, Hy, Hz), x=(0, Lx), y=(0, Ly), z=(-Lz, 0))

e = extract_slices(dir, filename, name="e")
u = extract_slices(dir, filename, name="u")
v = extract_slices(dir, filename, name="v")
T = extract_slices(dir, filename, name="T")

x, y, z = nodes((Center, Center, Center), grid)

x_xz = repeat(x, 1, Nz)
y_xz = 0.995 * Ly * ones(Nx, Nz)
z_xz = repeat(reshape(z, 1, Nz), Nx, 1)

x_yz = 0.995 * Lx * ones(Ny, Nz)
y_yz = repeat(y, 1, Nz)
z_yz = repeat(reshape(z, 1, Nz), grid.Ny, 1)

# Slight displacements to "stitch" the cube together
x_xy = x
y_xy = y
z_xy = - 0.001 * Lz * ones(grid.Nx, grid.Ny)

#####
##### Animate!
#####

fig = Figure(resolution=(1800, 800))

azimuth = 0.86
elevation = 0.38
perspectiveness = 0.6
aspect = :data
xlabel = "x (m)"
ylabel = "y (m)"
zlabel = "z (m)"

row = 1:6
bottom = row[end] + 2
ax_e = Axis3(fig[row, 1:4]; aspect, xlabel, ylabel, zlabel, azimuth, elevation, perspectiveness)
ax_T = Axis3(fig[row, 5:8]; aspect, xlabel, ylabel, zlabel, azimuth, elevation, perspectiveness)

ax_T_avg = Axis(fig[3:6, 9:10], xlabel="Temperature (ᵒC)", ylabel="z (m)")
ax_u_avg = Axis(fig[3:6, 11:12], xlabel="Velocity components (m s⁻¹)", ylabel="z (m)")

slider = Slider(fig[bottom+2, :], range=1:Nt, startvalue=1)
n = slider.value #Observable(1)

title = @lift @sprintf("Ocean surface boundary layer turbulence forced by cooling and light winds for %.1f hours",
                       times[$n] / hour)
lab = Label(fig[0, :], title, textsize=24)

###
### Vertical velocity
###

e_yz = @lift e.yz[$n]
e_xz = @lift e.xz[$n]
e_xy = @lift e.xy[$n]

T_yz = @lift T.yz[$n]
T_xz = @lift T.xz[$n]
T_xy = @lift T.xy[$n]

T_avg = @lift T.statistics[$n]
e_avg = @lift e.statistics[$n]
u_avg = @lift u.statistics[$n]
v_avg = @lift v.statistics[$n]

T_catke_conv_adj = @lift T.catke_conv_adj[$n]
e_catke_conv_adj = @lift e.catke_conv_adj[$n]
u_catke_conv_adj = @lift u.catke_conv_adj[$n]
v_catke_conv_adj = @lift v.catke_conv_adj[$n]

T_catke_constant_Pr = @lift T.catke_constant_Pr[$n]
e_catke_constant_Pr = @lift e.catke_constant_Pr[$n]
u_catke_constant_Pr = @lift u.catke_constant_Pr[$n]
v_catke_constant_Pr = @lift v.catke_constant_Pr[$n]

T_catke_variable_Pr = @lift T.catke_variable_Pr[$n]
e_catke_variable_Pr = @lift e.catke_variable_Pr[$n]
u_catke_variable_Pr = @lift u.catke_variable_Pr[$n]
v_catke_variable_Pr = @lift v.catke_variable_Pr[$n]

colormap_T = :oslo
colormap_e = :solar

e_max = maximum(maximum.(abs, e.xz))
colorrange_e = (0, e_max/20)

# Compute maximum temperature for colorrange with a moving average
ΔT = 0.04 # Absolute range
Δn = 20 # Half-width of moving average
Tf_max = maximum(T.xy[end])
T0_max = maximum(T.xy[1])

colorrange_T = @lift begin
    if $n < Δn + 1
        T_max = T0_max
    elseif $n < Nt - Δn
        T_max = mean([maximum(T.xy[nn]) for nn in $n-Δn:$n+Δn])
    else
        T_max = Tf_max
    end

    T_min = T_max - ΔT

    (T_min, T_max)
end

function box!(ax, yz, xz, xy; colormap, colorrange)
    sfc = surface!(ax, x_xz, y_xz, z_xz; color=xz, colormap, colorrange)
          surface!(ax, x_yz, y_yz, z_yz; color=yz, colormap, colorrange)
          surface!(ax, x_xy, y_xy, z_xy; color=xy, colormap, colorrange)
    return sfc
end

sfc_e = box!(ax_e, e_yz, e_xz, e_xy, colormap=colormap_e, colorrange=colorrange_e)
sfc_T = box!(ax_T, T_yz, T_xz, T_xy, colormap=colormap_T, colorrange=colorrange_T)

cp_e = fig[bottom, 2:3] = Colorbar(fig, sfc_e, flipaxis=false, vertical=false, label="Turbulent kinetic energy (m² s⁻²)")
cp_T = fig[bottom, 6:7] = Colorbar(fig, sfc_T, flipaxis=false, vertical=false, label="Temperature (ᵒC)")

lines!(ax_T_avg, T_catke_constant_Pr, z, linewidth=2, color=(:navy, 0.8), label="CATKE constant Pr")
lines!(ax_T_avg, T_catke_variable_Pr, z, linewidth=2, color=(:slateblue2, 0.8), label="CATKE variable Pr")
lines!(ax_T_avg, T_catke_conv_adj, z, linewidth=2, color=(:skyblue1, 0.8), label="CATKE w conv adj")

lines!(ax_u_avg, u_catke_constant_Pr, z, linewidth=2, color=(:navy, 0.8), label="u, CATKE constant Pr")
lines!(ax_u_avg, u_catke_variable_Pr, z, linewidth=2, color=(:slateblue2, 0.8), label="u, CATKE variable Pr")
lines!(ax_u_avg, u_catke_conv_adj, z, linewidth=2, color=(:skyblue1, 0.8), label="u, CATKE w conv adj")

lines!(ax_u_avg, v_catke_constant_Pr, z, linewidth=1, color=(:red4, 0.8), label="v, CATKE constant Pr")
lines!(ax_u_avg, v_catke_variable_Pr, z, linewidth=1, color=(:darkorange3, 0.8), label="v, CATKE variable Pr")
lines!(ax_u_avg, v_catke_conv_adj, z, linewidth=1, color=(:coral1, 0.8), label="v, CATKE w conv adj")

lines!(ax_T_avg, T_avg, z, linewidth=5, color=(:gray21, 0.6), label="LES")
lines!(ax_u_avg, u_avg, z, linewidth=5, color=(:gray21, 0.6), label="u, LES")
lines!(ax_u_avg, v_avg, z, linewidth=3, color=(:gray62, 0.6), label="v, LES")

xlims!(ax_u_avg, -0.2, 0.3)
fig[5, 13] = Legend(fig, ax_u_avg)
fig[7, 13] = Legend(fig, ax_T_avg)

hideydecorations!(ax_u_avg, grid=false)
hidespines!(ax_u_avg, :r, :t, :l)
hidespines!(ax_T_avg, :r, :t)

colgap!(fig.layout, Relative(0))
colgap!(fig.layout, 10, Relative(0.03))
rowgap!(fig.layout, Relative(0))

display(fig)

record(fig, joinpath(dir, filename * "_one_row_intro.mp4"), 1:Nt; framerate=16) do nn
    @info "Drawing frame $nn of $Nt..."
    n[] = nn
end

max_stars_repo_name: glwagner/LocalOceanClosureCalibration


text:
using Reactive, SignalView, GeometryTypes, Colors, GLWindow, GLVisualize

window = glscreen()

# create a nice looking signal graph
t = Signal(1.0)
color_a = Signal(RGBA{Float32}(1,0,0,1))
color_b = Signal(RGBA{Float32}(0,0,1,1))

colormap = map(color_a, color_b) do a, b
    _a, _b = RGB{Float32}(a), RGB{Float32}(b)
    cmap = map(RGBA{Float32}, linspace(_a, _b, 20))
    reshape(cmap, (1, length(cmap)))
end

# generate some pretty data
function xy_data(x,y,i)
    Float32(sin(1.3*x*i)*cos(0.9*y)+cos(.8*x)*sin(1.9*y)+cos(y*.2*x))
end
const N = 128
const range = linspace(-5f0, 5f0, N)
const heightmap = Float32[xy_data(x, y, value(t)) for x=range, y=range]

val = map(t) do v
    v/10f0
end
"""
some nice looking surface function
"""
function contourdata(t)
    for i=1:size(heightmap, 1)
        @simd for j=1:size(heightmap, 2)
            @inbounds heightmap[i,j] = xy_data(t, range[i], range[j])
        end
    end
    heightmap
end

hmap_color = map(colormap, val) do c, t
    hmap = contourdata(t)
    cmap = vec(c)
    mini, maxi = first(range), last(range)
    # you can pass along a tuple of pairs, which will map to the keyword args of visualize
    ((hmap, :surface), (
        :color_map=>cmap,
        :color_norm=>Vec2f0(-1,2),
        :ranges => (range, range),
        :boundingbox=>AABB{Float32}(Vec3f0(mini, mini, -1f0), Vec3f0(maxi-mini, maxi-mini, 3f0))
    ))
end

# view it as a contourplot as well!
map(hmap_color) do hmap_color
    reinterpret(Intensity{1, Float32}, hmap_color[1][1])
end

view_signal(t, window)
renderloop(window)


text:
"""
For each pixel, NOAA produces monthly estimates of radiance using the mean of radiance measured on days considered free of clouds. Radiance on months with no cloud-free observations are marked as zero, these should be marked as missing or or not available. The mark_missing function uses the radiance and the cloud-free observations image to marks missing wherever there no 0 cloud-free observations.
```julia
radiance = rand(1:10.0, 10, 10)
cloud = rand(0:5, 10, 10)
mark_missing(radiance, cloud)
```
"""
function mark_missing(radiance, clouds::Array{T, 2}) where T <: Any
    radiance = Array{Union{Float16, Missing}}(radiance)
    for i in 1:size(clouds)[1]
        for j in 1:size(clouds)[2]
            if clouds[i,j] == 0
                radiance[i,j]= missing
            end
        end
    end
    return radiance
end

"""
The mark_missing function also works on datacubes.
```julia
radiance = rand(1:10.0, 10, 10, 10)
cloud = rand(0:5, 10, 10, 10)
mark_missing(radiance, cloud)
```
Wherever the number of cloud-free observations is 0, radiance will be marked as missing.
"""
function mark_missing(radiance_datacube, clouds_datacube::Array{T, 3}) where T <: Any
    radiance_datacube = Array{Union{Float16, Missing}}(radiance_datacube)
    for i in 1:size(clouds_datacube)[3]
        radiance_datacube[:, :, i] = mark_missing(radiance_datacube[:, :, i], clouds_datacube[:, :, i])
    end
    return radiance_datacube
end


text:
module DynamicallyLoadedEmbedding

# Well, this is pretty much a dummy, at least for now. All the action
# happens in the `c` and `test` directories.

end


text:
#!/usr/bin/env julia

export pe1

pe1(n::Integer=1000)=union(3:3:(n-1),5:5:(n-1))|>sum

if !haskey(ENV,"PROJECT_EULER_WITHOUT_RESULT")
    pe1()|>print
end

max_stars_repo_path: 1. Multiples of 3 and 5/1.jl

max_stars_repo_name: miRoox/ProjectEuler

max_issues_repo_path: 1. Multiples of 3 and 5/1.jl

max_issues_repo_name: miRoox/ProjectEuler

max_forks_repo_path: 1. Multiples of 3 and 5/1.jl

max_forks_repo_name: miRoox/ProjectEuler


text:
using FluidInducedAseismicSlip
using Test

@testset "FluidInducedAseismicSlip.jl" begin
    # Test the calculation of λ for different values of T
    @test lambda_analytical_gs(0.1, 500)[1] ≈ 3.6403634658646724
    @test lambda_analytical_gs(0.2, 500)[1] ≈ 1.9124952153633206
    @test lambda_analytical_gs(0.3, 500)[1] ≈ 1.3474176336951285
    @test lambda_analytical_gs(0.4, 500)[1] ≈ 1.0252408283009076
    @test lambda_analytical_gs(0.5, 500)[1] ≈ 0.7916049887046924
    @test lambda_analytical_gs(0.6, 500)[1] ≈ 0.601228108399377
    @test lambda_analytical_gs(0.7, 500)[1] ≈ 0.4351313800638411
    @test lambda_analytical_gs(0.8, 500)[1] ≈ 0.2833769094180259
    @test lambda_analytical_gs(0.9, 500)[1] ≈ 0.1398124625560765
end

max_stars_repo_name: ajacquey/FluidInducedAseismicSlip

max_issues_repo_name: ajacquey/FluidInducedAseismicSlip

max_forks_repo_name: ajacquey/FluidInducedAseismicSlip


text:
module Generators

using Random

using Games: Player, NormalFormGame

include("bimatrix_generators.jl")
export blotto_game, ranking_game, sgc_game, unit_vector_game, tournament_game

end # module


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule libfacedetection_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("libfacedetection")
JLLWrappers.@generate_main_file("libfacedetection", UUID("d63ee80a-a26b-5956-a0ed-fe8a48e7c1c8"))
end  # module libfacedetection_jll


text:
const POWER_SYSTEM_DESCRIPTOR_FILE =
    joinpath(dirname(pathof(PowerSystems)), "descriptors", "power_system_inputs.json")

struct PowerSystemTableData
    basepower::Float64
    category_to_df::Dict{InputCategory, DataFrames.DataFrame}
    timeseries_metadata_file::Union{String, Nothing}
    directory::String
    user_descriptors::Dict
    descriptors::Dict
    generator_mapping::Dict{NamedTuple, DataType}
end

function PowerSystemTableData(
    data::Dict{String, Any},
    directory::String,
    user_descriptors::Union{String, Dict},
    descriptors::Union{String, Dict},
    generator_mapping::Union{String, Dict};
    timeseries_metadata_file = joinpath(directory, "timeseries_pointers"),
)
    category_to_df = Dict{InputCategory, DataFrames.DataFrame}()
    categories = [
        ("branch", BRANCH::InputCategory),
        ("bus", BUS::InputCategory),
        ("dc_branch", DC_BRANCH::InputCategory),
        ("gen", GENERATOR::InputCategory),
        ("load", LOAD::InputCategory),
        ("reserves", RESERVE::InputCategory),
        ("storage", STORAGE::InputCategory),
    ]

    if !haskey(data, "bus")
        throw(DataFormatError("key 'bus' not found in input data"))
    end

    if !haskey(data, "basepower")
        @warn "key 'basepower' not found in input data; using default=$(DEFAULT_BASE_MVA)"
    end
    basepower = get(data, "basepower", DEFAULT_BASE_MVA)

    for (label, category) in categories
        val = get(data, label, nothing)
        if isnothing(val)
            @warn "key '$label' not found in input data, set to nothing"
        else
            category_to_df[category] = val
        end
    end

    if !isfile(timeseries_metadata_file)
        if isfile(string(timeseries_metadata_file, ".json"))
            timeseries_metadata_file = string(timeseries_metadata_file, ".json")
        elseif isfile(string(timeseries_metadata_file, ".csv"))
            timeseries_metadata_file = string(timeseries_metadata_file, ".csv")
        else
            timeseries_metadata_file = nothing
        end
    end

    if user_descriptors isa AbstractString
        user_descriptors = _read_config_file(user_descriptors)
    end

    if descriptors isa AbstractString
        descriptors = _read_config_file(descriptors)
    end

    if generator_mapping isa AbstractString
        generator_mapping = get_generator_mapping(generator_mapping)
    end

    return PowerSystemTableData(
        basepower,
        category_to_df,
        timeseries_metadata_file,
        directory,
        user_descriptors,
        descriptors,
        generator_mapping,
    )
end

"""
     PowerSystemTableData(directory::AbstractString,
                          basepower::Float64,
                          user_descriptor_file::AbstractString;
                          descriptor_file=POWER_SYSTEM_DESCRIPTOR_FILE)

Reads in all the data stored in csv files
The general format for data is
    folder:
        gen.csv
        branch.csv
        bus.csv
        ..
        load.csv

# Arguments
- `directory::AbstractString`: directory containing CSV files
- `basepower::Float64`: base power for System
- `user_descriptor_file::AbstractString`: customized input descriptor file
- `descriptor_file=POWER_SYSTEM_DESCRIPTOR_FILE`: PowerSystems descriptor file
- `generator_mapping_file=GENERATOR_MAPPING_FILE`: generator mapping configuration file
"""
function PowerSystemTableData(
    directory::AbstractString,
    basepower::Float64,
    user_descriptor_file::AbstractString;
    descriptor_file = POWER_SYSTEM_DESCRIPTOR_FILE,
    generator_mapping_file = GENERATOR_MAPPING_FILE,
    timeseries_metadata_file = joinpath(directory, "timeseries_pointers"),
)
    files = readdir(directory)
    REGEX_DEVICE_TYPE = r"(.*?)\.csv"
    REGEX_IS_FOLDER = r"^[A-Za-z]+$"
    data = Dict{String, Any}()

    if length(files) == 0
        error("No files in the folder")
    else
        data["basepower"] = basepower
    end

    encountered_files = 0
    for d_file in files
        try
            if match(REGEX_IS_FOLDER, d_file) != nothing
                @info "Parsing csv files in $d_file ..."
                d_file_data = Dict{String, Any}()
                for file in readdir(joinpath(directory, d_file))
                    if match(REGEX_DEVICE_TYPE, file) != nothing
                        @info "Parsing csv data in $file ..."
                        encountered_files += 1
                        fpath = joinpath(directory, d_file, file)
                        raw_data = DataFrames.DataFrame(CSV.File(fpath))
                        d_file_data[split(file, r"[.]")[1]] = raw_data
                    end
                end

                if length(d_file_data) > 0
                    data[d_file] = d_file_data
                    @info "Successfully parsed $d_file"
                end

            elseif match(REGEX_DEVICE_TYPE, d_file) != nothing
                @info "Parsing csv data in $d_file ..."
                encountered_files += 1
                fpath = joinpath(directory, d_file)
                raw_data = DataFrames.DataFrame(CSV.File(fpath))
                data[split(d_file, r"[.]")[1]] = raw_data
                @info "Successfully parsed $d_file"
            end
        catch ex
            @error "Error occurred while parsing $d_file" exception = ex
            throw(ex)
        end
    end
    if encountered_files == 0
        error("No csv files or folders in $directory")
    end

    return PowerSystemTableData(
        data,
        directory,
        user_descriptor_file,
        descriptor_file,
        generator_mapping_file,
        timeseries_metadata_file = timeseries_metadata_file,
    )
end

"""
Return the custom name stored in the user descriptor file.

Throws DataFormatError if a required value is not found in the file.
"""
function get_user_field(
    data::PowerSystemTableData,
    category::InputCategory,
    field::AbstractString,
)
    if !haskey(data.user_descriptors, category)
        throw(DataFormatError("Invalid category=$category"))
    end

    try
        for item in data.user_descriptors[category]
            if item["name"] == field
                return Symbol(item["custom_name"])
            end
        end
    catch
        (err)
        if err == KeyError
            msg = "Failed to find category=$category field=$field in input descriptors $err"
            throw(DataFormatError(msg))
        else
            throw(err)
        end
    end

    msg = "Failed to find category=$category field=$field in input descriptors"
    throw(DataFormatError(msg))
end

"""Return a vector of user-defined fields for the category."""
function get_user_fields(data::PowerSystemTableData, category::InputCategory)
    if !haskey(data.user_descriptors, category)
        throw(DataFormatError("Invalid category=$category"))
    end

    return [x["name"] for x in data.user_descriptors[category]]
end

"""Return the dataframe for the category."""
function get_dataframe(data::PowerSystemTableData, category::InputCategory)
    @assert haskey(data.category_to_df, category)
    return data.category_to_df[category]
end

"""
    iterate_rows(data::PowerSystemTableData, category; na_to_nothing=true)

Return a NamedTuple of parameters from the descriptor file for each row of a dataframe,
making type conversions as necessary.

Refer to the PowerSystems descriptor file for field names that will be created.
"""
function iterate_rows(data::PowerSystemTableData, category; na_to_nothing = true)
    df = get_dataframe(data, category)
    field_infos = _get_field_infos(data, category, names(df))

    Channel() do channel
        for row in eachrow(df)
            obj = _read_data_row(data, row, field_infos; na_to_nothing = na_to_nothing)
            put!(channel, obj)
        end
    end
end

"""
    System(data::PowerSystemTableData)

Construct a System from PowerSystemTableData data.

# Arguments
- `forecast_resolution::Union{DateTime, Nothing}=nothing`: only store forecasts that match
  this resolution.
- `time_series_in_memory::Bool=false`: Store time series data in memory instead of HDF5 file
- `runchecks::Bool=true`: Validate struct fields.

Throws DataFormatError if forecasts with multiple resolutions are detected.
- A forecast has a different resolution than others.
- A forecast has a different horizon than others.

"""
function System(
    data::PowerSystemTableData;
    forecast_resolution = nothing,
    time_series_in_memory = false,
    runchecks = true,
)
    sys = System(
        data.basepower;
        time_series_in_memory = time_series_in_memory,
        runchecks = runchecks,
    )

    loadzone_csv_parser!(sys, data)
    bus_csv_parser!(sys, data)
    load_csv_parser!(sys, data)

    # Services and forecasts must be last.
    parsers = (
        (get_dataframe(data, BRANCH::InputCategory), branch_csv_parser!),
        (get_dataframe(data, DC_BRANCH::InputCategory), dc_branch_csv_parser!),
        (get_dataframe(data, GENERATOR::InputCategory), gen_csv_parser!),
        (get_dataframe(data, RESERVE::InputCategory), services_csv_parser!),
    )

    for (val, parser) in parsers
        if !isnothing(val)
            parser(sys, data)
        end
    end

    if !isnothing(data.timeseries_metadata_file)
        add_forecasts!(sys, data.timeseries_metadata_file; resolution = forecast_resolution)
    end

    check!(sys)
    return sys
end

"""
    bus_csv_parser!(sys::System, bus_raw::DataFrames.DataFrame)

Add buses and areas to the System from the raw data.

"""
function bus_csv_parser!(sys::System, data::PowerSystemTableData)
    for bus in iterate_rows(data, BUS::InputCategory)
        bus_type = get_enum_value(BusTypes.BusType, bus.bus_type)
        number = bus.bus_id
        voltage_limits = (min = 0.95, max = 1.05)

        area_name = string(bus.area)
        area = get_component(Area, sys, area_name)
        if isnothing(area)
            area = Area(area_name)
            add_component!(sys, area)
        end

        ps_bus = Bus(;
            number = number,
            name = bus.name,
            bustype = bus_type,
            angle = bus.angle,
            voltage = bus.voltage,
            voltagelimits = voltage_limits,
            basevoltage = bus.base_voltage,
            area = area,
            load_zone = get_component(LoadZone, sys, string(bus.zone)),
        )

        add_component!(sys, ps_bus)
    end
end

"""
    branch_csv_parser!(sys::System, data::PowerSystemTableData)

Add branches to the System from the raw data.

"""
function branch_csv_parser!(sys::System, data::PowerSystemTableData)
    available = true

    for branch in iterate_rows(data, BRANCH::InputCategory)
        bus_from = get_bus(sys, branch.connection_points_from)
        bus_to = get_bus(sys, branch.connection_points_to)
        connection_points = Arc(bus_from, bus_to)
        pf = get(branch, :pf, 0.0)
        qf = get(branch, :qf, 0.0)

        #TODO: noop math...Phase-Shifting Transformer angle
        alpha = (branch.primary_shunt / 2) - (branch.primary_shunt / 2)

        branch_type = get_branch_type(branch.tap, alpha)

        if branch_type == Line
            b = branch.primary_shunt / 2
            anglelimits = (min = -π / 2, max = π / 2) #TODO: add field in CSV
            value = Line(
                name = branch.name,
                available = available,
                activepower_flow = pf,
                reactivepower_flow = qf,
                arc = connection_points,
                r = branch.r,
                x = branch.x,
                b = (from = b, to = b),
                rate = branch.rate,
                anglelimits = anglelimits,
            )
        elseif branch_type == Transformer2W
            value = Transformer2W(
                name = branch.name,
                available = available,
                activepower_flow = pf,
                reactivepower_flow = qf,
                arc = connection_points,
                r = branch.r,
                x = branch.x,
                primaryshunt = branch.primary_shunt,
                rate = branch.rate,
            )
        elseif branch_type == TapTransformer
            value = TapTransformer(
                name = branch.name,
                available = available,
                activepower_flow = pf,
                reactivepower_flow = qf,
                arc = connection_points,
                r = branch.r,
                x = branch.x,
                primaryshunt = branch.primary_shunt,
                tap = branch.tap,
                rate = branch.rate,
            )
        elseif branch_type == PhaseShiftingTransformer
            # TODO create PhaseShiftingTransformer
            error("Unsupported branch type $branch_type")
        else
            error("Unsupported branch type $branch_type")
        end

        add_component!(sys, value)
    end
end

"""
    dc_branch_csv_parser!(sys::System, data::PowerSystemTableData)

Add DC branches to the System from raw data.

"""
function dc_branch_csv_parser!(sys::System, data::PowerSystemTableData)
    for dc_branch in iterate_rows(data, DC_BRANCH::InputCategory)
        available = true
        bus_from = get_bus(sys, dc_branch.connection_points_from)
        bus_to = get_bus(sys, dc_branch.connection_points_to)
        connection_points = Arc(bus_from, bus_to)
        pf = get(dc_branch, :pf, 0.0)

        if dc_branch.control_mode == "Power"
            mw_load = dc_branch.mw_load / data.basepower

            #TODO: is there a better way to calculate these?,
            activepowerlimits_from = (min = -1 * mw_load, max = mw_load)
            activepowerlimits_to = (min = -1 * mw_load, max = mw_load)
            reactivepowerlimits_from = (min = 0.0, max = 0.0)
            reactivepowerlimits_to = (min = 0.0, max = 0.0)
            loss = (l0 = 0.0, l1 = dc_branch.loss) #TODO: Can we infer this from the other data?,

            value = HVDCLine(
                name = dc_branch.name,
                available = available,
                activepower_flow = pf,
                arc = connection_points,
                activepowerlimits_from = activepowerlimits_from,
                activepowerlimits_to = activepowerlimits_to,
                reactivepowerlimits_from = reactivepowerlimits_from,
                reactivepowerlimits_to = reactivepowerlimits_to,
                loss = loss,
            )
        else
            rectifier_taplimits = (
                min = dc_branch.rectifier_tap_limits_min,
                max = dc_branch.rectifier_tap_limits_max,
            )
            rectifier_xrc = dc_branch.rectifier_xrc #TODO: What is this?,
            rectifier_firingangle = dc_branch.rectifier_firingangle
            inverter_taplimits = (
                min = dc_branch.inverter_tap_limits_min,
                max = dc_branch.inverter_tap_limits_max,
            )
            inverter_xrc = dc_branch.inverter_xrc #TODO: What is this?
            inverter_firingangle = (
                min = dc_branch.inverter_firing_angle_min,
                max = dc_branch.inverter_firing_angle_max,
            )
            value = VSCDCLine(
                name = dc_branch.name,
                available = true,
                activepower_flow = pf,
                arc = connection_points,
                rectifier_taplimits = rectifier_taplimits,
                rectifier_xrc = rectifier_xrc,
                rectifier_firingangle = rectifier_firingangle,
                inverter_taplimits = inverter_taplimits,
                inverter_xrc = inverter_xrc,
                inverter_firingangle = inverter_firingangle,
            )
        end

        add_component!(sys, value)
    end
end

"""
    gen_csv_parser!(sys::System, data::PowerSystemTableData)

Add generators to the System from the raw data.

"""
function gen_csv_parser!(sys::System, data::PowerSystemTableData)
    output_percent_fields = Vector{Symbol}()
    heat_rate_fields = Vector{Symbol}()
    fields = get_user_fields(data, GENERATOR::InputCategory)
    for field in fields
        if occursin("output_percent", field)
            push!(output_percent_fields, Symbol(field))
        elseif occursin("heat_rate_", field)
            push!(heat_rate_fields, Symbol(field))
        end
    end

    @assert length(output_percent_fields) > 0
    cost_colnames = zip(heat_rate_fields, output_percent_fields)

    for gen in iterate_rows(data, GENERATOR::InputCategory)
        bus = get_bus(sys, gen.bus_id)
        if isnothing(bus)
            throw(DataFormatError("could not find $(gen.bus_id)"))
        end

        generator = make_generator(data, gen, cost_colnames, bus)
        if !isnothing(generator)
            add_component!(sys, generator)
        end
    end
end

"""
    load_csv_parser!(sys::System, data::PowerSystemTableData)

Add loads to the System from the raw data.

"""
function load_csv_parser!(sys::System, data::PowerSystemTableData)
    for ps_bus in get_components(Bus, sys)
        max_active_power = 0.0
        max_reactive_power = 0.0
        active_power = 0.0
        reactive_power = 0.0
        found = false
        for bus in iterate_rows(data, BUS::InputCategory)
            if bus.bus_id == ps_bus.number
                max_active_power = bus.max_active_power
                max_reactive_power = bus.max_reactive_power
                active_power = get(bus, :active_power, max_active_power)
                reactive_power = get(bus, :reactive_power, max_reactive_power)
                found = true
                break
            end
        end

        if !found
            throw(DataFormatError("Did not find bus index in Load data $(ps_bus.name)"))
        end

        if (max_active_power != 0.0) || (max_reactive_power != 0.0)
            load = PowerLoad(
                name = ps_bus.name,
                available = true,
                bus = ps_bus,
                model = LoadModels.ConstantPower,
                activepower = active_power,
                reactivepower = reactive_power,
                maxactivepower = max_active_power,
                maxreactivepower = max_reactive_power,
            )
            add_component!(sys, load)
        end
    end
end

"""
    loadzone_csv_parser!(sys::System, data::PowerSystemTableData)

Add branches to the System from the raw data.

"""
function loadzone_csv_parser!(sys::System, data::PowerSystemTableData)
    buses = get_dataframe(data, BUS::InputCategory)
    zone_column = get_user_field(data, BUS::InputCategory, "zone")
    if !in(zone_column, names(buses))
        @warn "Missing Data : no 'zone' information for buses, cannot create loads based on zones"
        return
    end

    zones = unique(buses[!, zone_column])
    for zone in zones
        bus_numbers = Set{Int}()
        active_powers = Vector{Float64}()
        reactive_powers = Vector{Float64}()
        for bus in iterate_rows(data, BUS::InputCategory)
            if bus.zone == zone
                bus_number = bus.bus_id
                push!(bus_numbers, bus_number)

                active_power = bus.max_active_power
                push!(active_powers, active_power)

                reactive_power = bus.max_reactive_power
                push!(reactive_powers, reactive_power)
            end
        end

        name = string(zone)
        load_zone = LoadZone(name, sum(active_powers), sum(reactive_powers))
        add_component!(sys, load_zone)
    end
end

"""
    services_csv_parser!(sys::System, data::PowerSystemTableData)

Add services to the System from the raw data.

"""
function services_csv_parser!(sys::System, data::PowerSystemTableData)
    bus_id_column = get_user_field(data, BUS::InputCategory, "bus_id")
    bus_area_column = get_user_field(data, BUS::InputCategory, "area")

    # Shortcut for data that looks like "(val1,val2,val3)"
    make_array(x) = isnothing(x) ? x : split(strip(x, ['(', ')']), ",")

    function _add_device!(contributing_devices, device_categories, name)
        component = []
        for dev_category in device_categories
            component_type = _get_component_type_from_category(dev_category)
            components = get_components_by_name(component_type, sys, name)
            if length(components) == 0
                # There multiple categories, so we might not find a match in some.
                continue
            elseif length(components) == 1
                push!(component, components[1])
            else
                msg = "Found duplicate names type=$component_type name=$name"
                throw(DataFormatError(msg))
            end
        end
        if length(component) > 1
            msg = "Found duplicate components with name=$name"
            throw(DataFormatError(msg))
        elseif length(component) == 1
            push!(contributing_devices, component[1])
        end
    end

    for reserve in iterate_rows(data, RESERVE::InputCategory)
        device_categories = make_array(reserve.eligible_device_categories)
        device_subcategories =
            make_array(get(reserve, :eligible_device_subcategories, nothing))
        devices = make_array(get(reserve, :contributing_devices, nothing))
        regions = make_array(reserve.eligible_regions)
        requirement = get(reserve, :requirement, nothing)
        contributing_devices = Vector{Device}()

        if isnothing(device_subcategories)
            @info("Adding contributing components for $(reserve.name) by component name")
            for device in devices
                _add_device!(contributing_devices, device_categories, device)
            end
        else
            @info("Adding contributing generators for $(reserve.name) by category")
            @warn "Adding contributing components by category only supports generators" maxlog =
                1
            for gen in iterate_rows(data, GENERATOR::InputCategory)
                buses = get_dataframe(data, BUS::InputCategory)
                bus_ids = buses[!, bus_id_column]
                sys_gen = get_components_by_name(Generator, sys, gen.name)
                if length(sys_gen) == 1
                    sys_gen = sys_gen[1]
                    area = string(buses[
                        bus_ids .== get_number(get_bus(sys_gen)),
                        bus_area_column,
                    ][1])
                    if gen.category in device_subcategories && area in regions
                        _add_device!(contributing_devices, device_categories, gen.name)
                    end
                else
                    @warn "Found $(length(sys_gen)) Generators with name=$(gen.name)"
                end
            end
        end

        if length(contributing_devices) == 0
            throw(DataFormatError("did not find contributing devices for service $(reserve.name)"))
        end

        direction = get_reserve_direction(reserve.direction)
        if isnothing(requirement)
            service = StaticReserve{direction}(reserve.name, reserve.timeframe, 0.0)
        else
            service =
                VariableReserve{direction}(reserve.name, reserve.timeframe, requirement)
        end

        add_service!(sys, service, contributing_devices)
    end
end

function get_reserve_direction(direction::AbstractString)
    if direction == "Up"
        return ReserveUp
    elseif direction == "Down"
        return ReserveDown
    else
        throw(DataFormatError("invalid reserve direction $direction"))
    end
end

"""Creates a generator of any type."""
function make_generator(data::PowerSystemTableData, gen, cost_colnames, bus)
    generator = nothing
    gen_type =
        get_generator_type(gen.fuel, get(gen, :unit_type, nothing), data.generator_mapping)

    if gen_type == ThermalStandard
        generator = make_thermal_generator(data, gen, cost_colnames, bus)
    elseif gen_type <: HydroGen
        generator = make_hydro_generator(gen_type, data, gen, bus)
    elseif gen_type <: RenewableGen
        generator = make_renewable_generator(gen_type, data, gen, bus)
    elseif gen_type == GenericBattery
        generator = make_storage(data, gen, bus)
    else
        @error "Skipping unsupported generator" gen_type
    end

    return generator
end

function make_thermal_generator(data::PowerSystemTableData, gen, cost_colnames, bus)
    fuel_cost = gen.fuel_price / 1000

    var_cost = [(getfield(gen, hr), getfield(gen, mw)) for (hr, mw) in cost_colnames]
    var_cost = [
        (tryparse(Float64, string(c[1])), tryparse(Float64, string(c[2])))
        for c in var_cost if !in(nothing, c)
    ]
    if length(unique(var_cost)) > 1
        var_cost[2:end] = [
            (
                var_cost[i][1] *
                (var_cost[i][2] - var_cost[i - 1][2]) *
                fuel_cost *
                data.basepower,
                var_cost[i][2],
            ) .* gen.active_power_limits_max
            for i in 2:length(var_cost)
        ]
        var_cost[1] =
            (
                var_cost[1][1] * var_cost[1][2] * fuel_cost * data.basepower,
                var_cost[1][2],
            ) .* gen.active_power_limits_max

        fixed = max(
            0.0,
            var_cost[1][1] -
                (var_cost[2][1] / (var_cost[2][2] - var_cost[1][2]) * var_cost[1][2]),
        )
        var_cost[1] = (var_cost[1][1] - fixed, var_cost[1][2])
        for i in 2:length(var_cost)
            var_cost[i] = (var_cost[i - 1][1] + var_cost[i][1], var_cost[i][2])
        end
    else
        var_cost = [(0.0, var_cost[1][2]), (1.0, var_cost[1][2])]
        fixed = 0.0
    end

    available = true
    rating = sqrt(gen.active_power_limits_max^2 + gen.reactive_power_limits_max^2)
    active_power_limits =
        (min = gen.active_power_limits_min, max = gen.active_power_limits_max)
    reactive_power_limits =
        (min = gen.reactive_power_limits_min, max = gen.reactive_power_limits_max)
    tech = TechThermal(
        rating = rating,
        primemover = convert(PrimeMovers.PrimeMover, gen.unit_type),
        fuel = convert(ThermalFuels.ThermalFuel, gen.fuel),
        activepowerlimits = active_power_limits,
        reactivepowerlimits = reactive_power_limits,
        ramplimits = (up = gen.ramp_limits, down = gen.ramp_limits),
        timelimits = (up = gen.min_up_time, down = gen.min_down_time),
    )

    capacity = gen.active_power_limits_max
    startup_cost = get(gen, :startup_cost, nothing)
    if isnothing(startup_cost) && hasfield(typeof(gen), :startup_heat_cold_cost)
        startup_cost = gen.startup_heat_cold_cost * fuel_cost * 1000
    else
        @warn("No startup_cost defined for $(gen.name), setting to 0.0")
        startup_cost = 0.0
    end

    shutdown_cost = get(gen, :shutdown_cost, nothing)
    if isnothing(shutdown_cost)
        @warn("No shutdown_cost defined for $(gen.name), setting to 0.0")
        shutdown_cost = 0.0
    end
    op_cost = ThreePartCost(var_cost, fixed, startup_cost, shutdown_cost)

    return ThermalStandard(
        gen.name,
        available,
        bus,
        gen.active_power,
        gen.reactive_power,
        tech,
        op_cost,
    )
end

function make_hydro_generator(gen_type, data::PowerSystemTableData, gen, bus)
    available = true

    rating = calculate_rating(gen.active_power_limits_max, gen.reactive_power_limits_max)
    active_power_limits =
        (min = gen.active_power_limits_min, max = gen.active_power_limits_max)
    reactive_power_limits =
        (min = gen.reactive_power_limits_min, max = gen.reactive_power_limits_max)
    tech = TechHydro(
        rating,
        convert(PrimeMovers.PrimeMover, gen.unit_type),
        active_power_limits,
        reactive_power_limits,
        (up = gen.ramp_limits, down = gen.ramp_limits),
        (up = gen.min_down_time, down = gen.min_down_time),
    )

    if gen_type == HydroEnergyReservoir
        if !haskey(data.category_to_df, STORAGE)
            throw(DataFormatError("Storage information must defined in storage.csv"))
        end
        @debug("Creating $(gen.name) as HydroEnergyReservoir")
        storage = get_storage_by_generator(data, gen.name)
        curtailcost = 0.0
        hydro_gen = HydroEnergyReservoir(
            name = gen.name,
            available = available,
            bus = bus,
            activepower = gen.active_power,
            reactivepower = gen.reactive_power,
            tech = tech,
            op_cost = TwoPartCost(curtailcost, 0.0),
            storage_capacity = storage.storage_capacity,
            inflow = storage.inflow_limit,
            initial_storage = storage.initial_storage,
        )
    elseif gen_type == HydroDispatch
        @debug("Creating $(gen.name) as HydroDispatch")
        hydro_gen = HydroDispatch(
            name = gen.name,
            available = available,
            bus = bus,
            activepower = gen.active_power,
            reactivepower = gen.reactive_power,
            tech = tech,
        )
    else
        error("Tabular data parser does not currently support $gen_type creation")
    end
    return hydro_gen
end

function get_storage_by_generator(data::PowerSystemTableData, gen_name::AbstractString)
    for storage in iterate_rows(data, STORAGE::InputCategory)
        if storage.generator_name == gen_name
            return storage
        end
    end

    throw(DataFormatError("no storage exists with generator $gen_name"))
end

function make_renewable_generator(gen_type, data::PowerSystemTableData, gen, bus)
    generator = nothing
    available = true
    rating = gen.active_power_limits_max

    tech = TechRenewable(
        rating,
        convert(PrimeMovers.PrimeMover, gen.unit_type),
        (min = gen.reactive_power_limits_min, max = gen.reactive_power_limits_max),
        1.0,
    )
    if gen_type == RenewableDispatch
        generator = RenewableDispatch(
            gen.name,
            available,
            bus,
            gen.active_power,
            gen.reactive_power,
            tech,
            TwoPartCost(0.0, 0.0),
        )
    elseif gen_type == RenewableFix
        generator = RenewableFix(
            gen.name,
            available,
            bus,
            gen.active_power,
            gen.reactive_power,
            tech,
        )
    else
        error("Unsupported type $gen_type")
    end

    return generator
end

function make_storage(data::PowerSystemTableData, gen, bus)
    available = true
    energy = 0.0
    capacity = (min = gen.active_power_limits_min, max = gen.active_power_limits_max)
    rating = gen.active_power_limits_max
    input_active_power_limits = (min = 0.0, max = gen.active_power_limits_max)
    output_active_power_limits = (min = 0.0, max = gen.active_power_limits_max)
    efficiency = (in = 0.9, out = 0.9)
    reactive_power_limits = (min = 0.0, max = 0.0)

    battery = GenericBattery(
        name = gen.name,
        available = available,
        bus = bus,
        primemover = convert(PrimeMovers.PrimeMover, gen.unit_type),
        energy = energy,
        capacity = capacity,
        rating = rating,
        activepower = gen.active_power,
        inputactivepowerlimits = input_active_power_limits,
        outputactivepowerlimits = output_active_power_limits,
        efficiency = efficiency,
        reactivepower = gen.reactive_power,
        reactivepowerlimits = reactive_power_limits,
    )

    return battery
end

const CATEGORY_STR_TO_COMPONENT = Dict{String, DataType}(
    "Bus" => Bus,
    "Generator" => Generator,
    "Reserve" => Service,
    "LoadZone" => LoadZone,
    "ElectricLoad" => ElectricLoad,
)

function _get_component_type_from_category(category::AbstractString)
    component_type = get(CATEGORY_STR_TO_COMPONENT, category, nothing)
    if isnothing(component_type)
        throw(DataFormatError("unsupported category=$category"))
    end

    return component_type
end

function _read_config_file(file_path::String)
    return open(file_path) do io
        data = YAML.load(io)
        # Replace keys with enums.
        config_data = Dict{InputCategory, Vector}()
        for (key, val) in data
            # TODO: need to change user_descriptors.yaml to use reserve instead.
            if key == "reserves"
                key = "reserve"
            end
            config_data[get_enum_value(InputCategory, key)] = val
        end
        return config_data
    end
end

"""Stores user-customized information for required dataframe columns."""
struct _FieldInfo
    name::String
    custom_name::Symbol
    needs_per_unit_conversion::Bool
    unit_conversion::Union{NamedTuple{(:From, :To), Tuple{String, String}}, Nothing}
    # TODO unit, value ranges and options
end

function _get_field_infos(data::PowerSystemTableData, category::InputCategory, df_names)
    if !haskey(data.user_descriptors, category)
        throw(DataFormatError("Invalid category=$category"))
    end

    if !haskey(data.descriptors, category)
        throw(DataFormatError("Invalid category=$category"))
    end

    # Cache whether PowerSystems uses a column's values as system-per-unit.
    # The user's descriptors indicate that the raw data is already system-per-unit or not.
    per_unit = Dict{String, Bool}()
    unit = Dict{String, Union{String, Nothing}}()
    descriptor_names = Vector{String}()
    for descriptor in data.descriptors[category]
        per_unit[descriptor["name"]] = get(descriptor, "system_per_unit", false)
        unit[descriptor["name"]] = get(descriptor, "unit", nothing)
        push!(descriptor_names, descriptor["name"])
    end

    fields = Vector{_FieldInfo}()
    try
        for item in data.user_descriptors[category]
            custom_name = Symbol(item["custom_name"])
            name = item["name"]
            if custom_name in df_names
                if !(name in descriptor_names)
                    if occursin("heat_rate_", name) || occursin("output_percent_", name)
                        base = name[(findlast("_", name)[end] + 1):end]
                        d_name = descriptor_names[occursin.(base, descriptor_names)][end]
                        per_unit[name] = per_unit[d_name]
                        unit[name] = unit[d_name]
                    else
                        throw(DataFormatError("$name is not defined in $POWER_SYSTEM_DESCRIPTOR_FILE"))
                    end
                end

                if !per_unit[name] && get(item, "system_per_unit", false)
                    throw(DataFormatError("$name cannot be defined as system_per_unit"))
                end

                needs_pu_conversion =
                    per_unit[name] &&
                    haskey(item, "system_per_unit") && !item["system_per_unit"]

                custom_unit = get(item, "unit", nothing)
                if !isnothing(unit[name]) &&
                   !isnothing(custom_unit) && custom_unit != unit[name]
                    unit_conversion = (From = custom_unit, To = unit[name])
                else
                    unit_conversion = nothing
                end

                push!(
                    fields,
                    _FieldInfo(name, custom_name, needs_pu_conversion, unit_conversion),
                )
            else
                # TODO: This should probably be a fatal error. However, the parsing code
                # doesn't use all the descriptor fields, so skip for now.
                @warn "User-defined column name $custom_name is not in dataframe."
            end
        end
        return fields
    catch err
        if err == KeyError
            msg = "Failed to find category=$category field=$field in input descriptors $err"
            throw(DataFormatError(msg))
        else
            throw(err)
        end
    end

    msg = "Failed to find category=$category field=$field in input descriptors"
    throw(DataFormatError(msg))
end

"""Reads values from dataframe row and performs necessary conversions."""
function _read_data_row(data::PowerSystemTableData, row, field_infos; na_to_nothing = true)
    fields = Vector{String}()
    vals = Vector()
    for field_info in field_infos
        value = row[field_info.custom_name]
        if ismissing(value)
            throw(DataFormatError("$(field_info.custom_name) value missing"))
        end
        if na_to_nothing && value == "NA"
            value = nothing
        end

        if field_info.needs_per_unit_conversion
            @debug "convert to system_per_unit" field_info.custom_name
            value /= data.basepower
        end

        # TODO: need special handling for units
        if !isnothing(field_info.unit_conversion)
            @debug "convert units" field_info.custom_name
            value = convert_units!(value, field_info.unit_conversion)
        end
        # TODO: validate ranges and option lists

        push!(fields, field_info.name)
        push!(vals, value)
    end

    return NamedTuple{Tuple(Symbol.(fields))}(vals)
end


text:
# Additions to Base's searchsorted functionality
"""
    searchsortednearest(vec::AbstractVector, x)

Like `searchsortedfirst` or `searchsortedlast`, this returns the the index of
the element in the sorted vector `vec` whose value is closest to `x`, rounding
up. If there are multiple elements that are equally close to `x`, this will
return the first index if `x` is less than or equal to those in the vector or
the last index if `x` is greater.
"""
function searchsortednearest(vec::AbstractVector, x)
    idx = searchsortedfirst(vec, x) # Returns the first idx | vec[idx] >= x
    if idx > 1 && (idx > length(vec) || (vec[idx] - x) > (x - vec[idx-1]))
        idx -= 1 # The previous element is closer
    end
    return idx
end

# We depend upon extrapolative behaviors in searching ranges to shift axes.
# This can be done by stealing Base's implementations and removing the bounds-
# correcting min/max.

# TODO: This could plug into the sorting system better, but it's fine for now
# TODO: This needs to support Dates.
"""
    unsafe_searchsorted(a::Range, I::ClosedInterval)

Return the indices of the range that fall within an interval without checking
bounds, possibly extrapolating outside the range if needed.
"""
function unsafe_searchsorted(a::Range, I::ClosedInterval)
    unsafe_searchsortedfirst(a, I.left):unsafe_searchsortedlast(a, I.right)
end
# Base only specializes searching ranges by Numbers; so optimize for Intervals
function Base.searchsorted(a::Range, I::ClosedInterval)
    searchsortedfirst(a, I.left):searchsortedlast(a, I.right)
end

# When running with "--check-bounds=yes" (like on Travis), the bounds-check isn't elided
@inline function inbounds_getindex{T}(v::Range{T}, i::Integer)
    convert(T, first(v) + (i-1)*step(v))
end
@inline function inbounds_getindex{T<:Integer}(r::Range, s::Range{T})
    f = first(r)
    st = oftype(f, f + (first(s)-1)*step(r))
    range(st, step(r)*step(s), length(s))
end
if VERSION < v"0.6.0-dev.2390"
    include_string("""
    @inline function inbounds_getindex{T}(r::FloatRange{T}, i::Integer)
        convert(T, (r.start + (i-1)*r.step)/r.divisor)
    end
    @inline function inbounds_getindex(r::FloatRange, s::OrdinalRange)
        FloatRange(r.start + (first(s)-1)*r.step, step(s)*r.step, length(s), r.divisor)
    end
    """)
else
    include_string("""
    @inline inbounds_getindex(r::StepRangeLen, i::Integer) = Base.unsafe_getindex(r, i)
    @inline function inbounds_getindex(r::StepRangeLen, s::OrdinalRange)
        vfirst = Base.unsafe_getindex(r, first(s))
        StepRangeLen(vfirst, step(r)*step(s), length(s))
    end
    """)
end

function unsafe_searchsortedlast{T<:Number}(a::Range{T}, x::Number)
    step(a) == 0 && throw(ArgumentError("ranges with a zero step are unsupported"))
    n = round(Integer,(x-first(a))/step(a))+1
    isless(x, inbounds_getindex(a, n)) ? n-1 : n
end
function unsafe_searchsortedfirst{T<:Number}(a::Range{T}, x::Number)
    step(a) == 0 && throw(ArgumentError("ranges with a zero step are unsupported"))
    n = round(Integer,(x-first(a))/step(a))+1
    isless(inbounds_getindex(a, n), x) ? n+1 : n
end
function unsafe_searchsortedlast{T<:Integer}(a::Range{T}, x::Number)
    step(a) == 0 && throw(ArgumentError("ranges with a zero step are unsupported"))
    fld(floor(Integer,x)-first(a),step(a))+1
end
function unsafe_searchsortedfirst{T<:Integer}(a::Range{T}, x::Number)
    step(a) == 0 && throw(ArgumentError("ranges with a zero step are unsupported"))
    -fld(floor(Integer,-x)+first(a),step(a))+1
end
function unsafe_searchsortedfirst{T<:Integer}(a::Range{T}, x::Unsigned)
    step(a) == 0 && throw(ArgumentError("ranges with a zero step are unsupported"))
    -fld(first(a)-signed(x),step(a))+1
end
function unsafe_searchsortedlast{T<:Integer}(a::Range{T}, x::Unsigned)
    step(a) == 0 && throw(ArgumentError("ranges with a zero step are unsupported"))
    fld(signed(x)-first(a),step(a))+1
end


text:
mutable struct SwitchGenerateState{T}
    score::Float64
    noise::Float64
    weight::Float64
    index::Int
    subtrace::Trace
    retval::T
    SwitchGenerateState{T}(score::Float64, noise::Float64, weight::Float64) where T = new{T}(score, noise, weight)
end

function process!(gen_fn::Switch{C, N, K, T},
                  index::Int,
                  args::Tuple,
                  choices::ChoiceMap,
                  state::SwitchGenerateState{T}) where {C, N, K, T}

    (subtrace, weight) = generate(getindex(gen_fn.branches, index), args, choices)
    state.index = index
    state.subtrace = subtrace
    state.weight += weight
    state.retval = get_retval(subtrace)
end

@inline process!(gen_fn::Switch{C, N, K, T}, index::C, args::Tuple, choices::ChoiceMap, state::SwitchGenerateState{T}) where {C, N, K, T} = process!(gen_fn, getindex(gen_fn.cases, index), args, choices, state)

function generate(gen_fn::Switch{C, N, K, T},
                  args::Tuple,
                  choices::ChoiceMap) where {C, N, K, T}

    index = args[1]
    state = SwitchGenerateState{T}(0.0, 0.0, 0.0)
    process!(gen_fn, index, args[2 : end], choices, state)
    return SwitchTrace(gen_fn, state.subtrace,
                       state.retval, args,
                       state.score, state.noise), state.weight
end


text:
using Memoize
using Base.Test

@doc """
    ncases_swb(minblocklen::Integer, rowlen::Integer)

Compute the number of cases starting with a block.

# Arguments
* `minblocklen::Integer`: mininal block length
* `rowlen::Integer`: row length
""" ncases_swb

@memoize function ncases_swb(minblocklen::Integer, rowlen::Integer)
    if rowlen < minblocklen
        return 0
    end

    result = 0
    for blocklen in minblocklen:rowlen
        result += 1
        for nspace = 1:(rowlen-blocklen)
            result += ncases_swb(minblocklen, rowlen-blocklen-nspace)
        end
    end

    result
end

function ncases(minblocklen::Integer, rowlen::Integer)
    result = 1 # no block at all
    for nheadspace = 0:rowlen-minblocklen
        result += ncases_swb(minblocklen, rowlen-nheadspace)
    end

    result
end

@test ncases(3, 7) == 17
@test ncases(3, 50) == 16475640049
@test ncases(3, 29) == 673135
@test ncases(3, 30) == 1089155
@test ncases(10, 56) == 880711
@test ncases(10, 57) == 1148904

rowlen = 51
while ncases(50, rowlen) < 1_000_000
    rowlen += 1
end

@printf("Euler 0115: F(50, %d) = %d\n", rowlen-1, ncases(50, rowlen-1))
@printf("Euler 0115: F(50, %d) = %d\n", rowlen  , ncases(50, rowlen  ))

max_stars_repo_name: dpieroux/euler

max_issues_repo_name: dpieroux/euler

max_forks_repo_name: dpieroux/euler


text:
struct Centerer{T} <: Number
    s::T
end

Centerer(x::Vector{T}; robust::Bool=false) where T = (m = robust ? median(x) : mean(x); Centerer(m))

(s::Centerer)(x::Vector) = x .- s.s

mutable struct StepCenterer{S,F} <: AbstractStep
    selections::S
    params::Union{Nothing, Dict{Symbol,Centerer}}
    robust::Bool
    trained::Bool
    skip::Bool
    prehook::F
end

function transform!(s::StepCenterer, df)
    for (k,f) in s.params
        df[k] = f(df[k])
    end
    df
end
function transform(s::StepCenterer, df)
    transform!(s, copy(df))
end
function fit!(s::StepCenterer, df)
    s.params = Dict(col => Centerer(s.prehook(df[col])) for col in getselectionkeys(df, s.selections))
    s.trained = true
end
function step_center!(r::Recipe, s...; robust::Bool=false, abs::Bool=false, skip::Bool=false, prehook=identity)
    push!(r.steps,
          StepCenterer([s...],
                       nothing,
                       robust,
                       false,
                       skip,
                       prehook))
end


text:
function Poisson_solve( pw_grid::PWGrid, rho )
    #
    Ω  = pw_grid.Ω
    G2 = pw_grid.G2
    Ns = pw_grid.Ns
    Npoints = pw_grid.Npoints
    #
    ctmp = 4.0*pi*R_to_G( Ns, rho )
    for ip = 2:Npoints
        ctmp[ip] = ctmp[ip] / G2[ip]
    end
    ctmp[1] = 0.0
    phi = real( G_to_R( Ns, ctmp ) )
    return phi
end

# The same as `solve_poisson`, but using FFT that calls to C directly
function solve_poisson_new( pw_grid::PWGrid, rho )
    #
    Ω  = pw_grid.Ω
    G2 = pw_grid.G2
    Ns = pw_grid.Ns
    Npoints = pw_grid.Npoints
    #
    ctmp = 4.0*pi*c_R_to_G( Ns, convert(Array{Complex128},rho) )
    for ip = 2:Npoints
        ctmp[ip] = ctmp[ip] / G2[ip]
    end
    ctmp[1] = 0.0
    phi = real( c_G_to_R( Ns, ctmp ) )
    return phi
end


text:
using InteractiveUtils: versioninfo

versioninfo()

@info("Sys.CPU_THREADS is $(Sys.CPU_THREADS)")
@info("Threads.nthreads() is $(Threads.nthreads()) threads")

function is_coverage()
    return !iszero(Base.JLOptions().code_coverage)
end

const coverage = is_coverage()

@info("Code coverage is $(coverage ? "enabled" : "disabled")")


text:
using NLopt
using DataFrames

include(joinpath("..", "src", "problems", "single_objective_base_functions.jl"))

abstract Evaluator

type NonGradientEvaluator <: Evaluator
  func::Function
  fevals::Int
  NonGradientEvaluator(func) = new(func, 0)
end
import Base.reset
reset(e::Evaluator) = (e.fevals = 0)
fevals(e::Evaluator) = e.fevals

function makeobjfunc(e::Evaluator)
  ofunc = (x::Vector, grad::Vector) -> begin
    e.fevals::Int += 1
    e.func(x)
  end
  ofunc
end

LocalNLoptAlgs = [:LN_COBYLA, :LN_BOBYQA, :LN_NEWUOA, :LN_PRAXIS, :LN_NELDERMEAD, :LN_SBPLX]
GlobalNLoptAlgs = [:GN_DIRECT, :GN_DIRECT_L, :GN_ORIG_DIRECT, :GN_ORIG_DIRECT_L,
  :GN_CRS2_LM, :GN_ISRES, :GN_ESCH]

# Randomly select np pairs of algs where the first is global and the second is local.
function random_pairs_of_algs(np)
  algs = Symbol[]
  for i in 1:np
    push!(algs, shuffle(GlobalNLoptAlgs)[1])
    push!(algs, shuffle(LocalNLoptAlgs)[1])
  end
  algs
end

function compare_nlopt_algs(funcToOpt, dim;
  LowBound = -10.0,
  HighBound = 10.0,
  NumReps = 3,
  MaxTimes = [1.0, 3.0, 10.0],
  Algs = [:LN_COBYLA, :LN_BOBYQA, :LN_NEWUOA, :LN_PRAXIS, :LN_NELDERMEAD, :LN_SBPLX,
    :GN_DIRECT, :GN_DIRECT_L, :GN_CRS2_LM, :GN_ISRES, :GN_ESCH]
  )

  function objfunc(x::Vector, grad::Vector)
    # We don't have gradient...

    global numfevals
    numfevals::Int += 1

    funcToOpt(x)
  end

  res = DataFrame()

  for maxt in MaxTimes
    for alg in Algs
      for rep in 1:NumReps

        opt = Opt(alg, dim)
        lower_bounds!(opt, (LowBound  * ones(dim)))
        upper_bounds!(opt, (HighBound * ones(dim)))
        maxtime!(opt, maxt)
        min_objective!(opt, objfunc)

        println("Solving...")

        global numfevals = 0
        tic()
        (fmin, solution, ret) = optimize(opt, randn(dim))
        t = toq()

        println("Dims = ", dim)
        println("Alg = ", alg)
        println("time given = ", maxt, " seconds")
        println("time taken = ", t, " seconds (", round(t/maxt*100.0, 2), "%)")
        println("fmin = ", fmin)
        println("Fevals = ", numfevals)
        println("Return value = ", ret, "\n\n")

        res = rbind(res, DataFrame(Alg = alg, D = dim, Rep = rep, Fevals = numfevals,
          TimeGiven = maxt, TimeTaken = t, TimeUsed = round(100.0*t/maxt, 2), fmin = fmin))
      end
    end
  end
  res
end

# Given an array of NLopt algs and a vector of relative times to apply them
# we run them in sequence on the value of the previous one.
function combine_nlopt_algs(func, dim;
  Algs = nothing, RelTimes = nothing,
  LowBound = -100.0,
  HighBound = 100.0,
  MaxTime = nothing,
  ResetToBestProbability = 0.80,
  FTOL_REL = -Inf,
  FTOL_ABS = -Inf
  )

  if Algs == nothing
    numpairs = iceil(log10(dim))
    Algs = vcat([:GN_DIRECT_L, :LN_PRAXIS], random_pairs_of_algs(numpairs),
      shuffle(LocalNLoptAlgs), [:LN_PRAXIS])
    RelTimes = vcat([1, 3], rand(1:3, 2*numpairs), 0.30 * rand(length(LocalNLoptAlgs)), 2)
  end

  if RelTimes == nothing
    RelTimes = rand(length(Algs))
  end

  if MaxTime == nothing
    MaxTime = 3*sqrt(dim)
  end

  numsteps = length(Algs)
  @assert numsteps == length(RelTimes)

  # Ensure percentages
  percenttimes = RelTimes / sum(RelTimes)

  fmin = ret = 0
  global numfevals = 0

  evaluator = NonGradientEvaluator(func)
  objfunc = makeobjfunc(evaluator)

  xbest = x = randn(dim)
  fbest = objfunc(xbest, zeros(dim))
  numtolreached = fevalsbest = 0

  start_time = time()

  println("\nStarting solve sequence...")
  for i in 1:numsteps
    opt = Opt(Algs[i], dim)
    lower_bounds!(opt, (LowBound  * ones(dim)))
    upper_bounds!(opt, (HighBound * ones(dim)))
    maxt = MaxTime * percenttimes[i]
    maxtime!(opt, maxt)
    ftol_rel!(opt, FTOL_REL)
    ftol_abs!(opt, FTOL_ABS)
    min_objective!(opt, objfunc)
    try
      tic()
      (fmin, x, ret) = optimize(opt, x)
      t = toq()
      println(Algs[i], " for ", round(t, 3), " s (", round(t/maxt*100.0, 2), "%), fmin = ", fmin, ", ret = ", ret)
    catch err
      println(Algs[i], " failed. Skipping...")
    end
    nowtime = time()
    if fmin < fbest
      fbest = fmin
      xbest = x
      fevalsbest = fevals(evaluator)
      timebest = nowtime
    else
      #if rand() < ResetToBestProbability
      #  x = xbest
      #end
    end
    elapsed = nowtime - start_time
    if (nowtime - timebest) > 0.20 * MaxTime
      if in(ret, [:SUCCESS, :FTOL_REACHED, :XTOL_REACHED, :ROUNDOFF_LIMITED])
        if numtolreached < 1
          # If first time to converge we try once more from a new random position
          numtolreached += 1
          x = randn(dim)
          println("Restart since ret = ", ret)
        else
          println("Return since ret = ", ret, " and numtolreached = ", numtolreached)
          return (fbest, xbest, :CONVERGED, fevalsbest, fevalsbest/dim, fevals(evaluator), fevals(evaluator)/dim, time() - start_time)
        end
      end
    end
  end

  return (fbest, xbest, ret, fevalsbest, fevalsbest/dim, fevals(evaluator), fevals(evaluator)/dim, time() - start_time)
end

dfmapper(r) = begin
  DataFrame(nfevalsbest = r[4], nfevals = r[6], nfbpd = r[5], fmin = r[1])
end

SOPS = [sphere, ellipsoid, elliptic, cigar, cigtab, ackley, schwefel1_2,
  schwefel2_22, schwefel2_21, rastrigin, rosenbrock]

SOPSNonzero = [step, griewank, schwefel2_26, shekel10]

function xrotatedandshifted(n, f, shiftAmplitude = 1.0, rotateAmplitude = 1.0)
  shift = shiftAmplitude * randn(n, 1)
  rotmatrix = rotateAmplitude * rand(n, n)
  transformed_f(x) = f(rotmatrix * (x .- shift))
end

#res8 = map((sop) -> combine_nlopt_algs(sop, 8), SOPS)
#resrot8 = map((sop) -> combine_nlopt_algs(xrotatedandshifted(8, sop), 8), SOPS)

#res64 = map((sop) -> combine_nlopt_algs(sop, 64), SOPS)
#resrot64 = map((sop) -> combine_nlopt_algs(xrotatedandshifted(64, sop), 64), SOPS)

#res256 = map((sop) -> combine_nlopt_algs(sop, 256), SOPS)
#resrot256 = map((sop) -> combine_nlopt_algs(xrotatedandshifted(256, sop), 256), SOPS)

function compare_optimizers(optimizers, problems, dim, numreps = 10;
  dataframe_mapper = identity)

  timestamp = strftime("%Y%m%d_%H%M%S", time())
  filename = "results_nlopt_$(timestamp).csv"
  res = DataFrame()
  its = 0
  for rep in 1:numreps
    for (probname, origproblem) in problems
      problem = xrotatedandshifted(dim, origproblem)
      for (optname, opt) in optimizers
        its += 1
        println("Optimization ", its)
        tic()
        r = opt(problem, dim)
        t = toq()
        res = rbind(res, cbind(
          DataFrame(Opt = optname, Prob = probname, Dim = dim, Time = round(t, 2)),
          dataframe_mapper(r)
          ))
        writetable(filename, res)
        println("Wrote results to file: ", filename)
      end
    end
  end
  return res, filename
end

Problems = {
  "Sphere      " => sphere,
  "Rosenbrock  " => rosenbrock,
  "Schwefel2_22" => schwefel2_22,
  "Schwefel2_21" => schwefel2_21,
  "Rastrigin   " => rastrigin,
}

DP = [:GN_DIRECT_L, :LN_PRAXIS]
DPN = [:GN_DIRECT_L, :LN_PRAXIS, :LN_NELDERMEAD]

Optimizers = {
  "DirL1+Prax3       " => (pr, d) -> combine_nlopt_algs(pr, d; Algs = DP, RelTimes = [1,3]),
  "DirL1+Prax2+NelM2 " => (pr, d) -> combine_nlopt_algs(pr, d; Algs = DPN, RelTimes = [1,2,2]),
  "RandNLopt non-grad" => (pr, d) -> combine_nlopt_algs(pr, d),
}

function compare_and_save(opts, probs, dims, reps, mapper)
  for dim in dims
    res, filename = compare_optimizers(opts, probs, dim, reps; dataframe_mapper = mapper)
    println("Wrote results to file: ", filename)
  end
end

compare_and_save(Optimizers, Problems, [256, 1024], 1, dfmapper)


text:
include("timeop.jl")
include("Cooling.jl")

max_stars_repo_name: Roger-luo/QuCmp

max_issues_repo_name: Roger-luo/QuCmp

max_forks_repo_name: Roger-luo/QuCmp


text:
@testset "Req Body" begin
    @testset "JSON Parser" begin
        app = App()
        req_body = Dict("a" => 1, "b" => 3, "c" => -4)

        req = HTTP.Request("POST", "", [], JSON3.write(req_body))
        req_empty = HTTP.Request("POST", "", [], "")

        @test app.reqparser isa JSONParser

        @test WebAPI.parsebody(app.reqparser, req) == JSON3.read(req.body)
        @test WebAPI.parsebody(app.reqparser, req_empty) == JSON3.Object()
    end

    @testset "Text Parser" begin
        app = App(reqparser=TextParser())
        req_body = "Test String"
        req = HTTP.Request("POST", "", [], req_body)

        @test app.reqparser isa TextParser

        @test WebAPI.parsebody(app.reqparser, req) == req_body
    end
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "gperftools"
version = v"2.7.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/gperftools/gperftools/releases/download/gperftools-2.7/gperftools-2.7.tar.gz", "1ee8c8699a0eff6b6a203e59b43330536b22bbcbe6448f54c7091e5efb0763c9")
]

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

# Remove misleading libtool files
rm -f ${prefix}/lib/*.la
rm -f /opt/${target}/${target}/lib*/*.la
rm -f /opt/${MACHTYPE}/${MACHTYPE}/lib*/*.la

if [[ "${target}" == *-linux-* ]]; then
    # https://github.com/gperftools/gperftools/blob/e5f77d6485bd2f6ce43862e3e57118b1bb97d30a/README
    export CXXFLAGS="-fno-builtin-malloc -fno-builtin-calloc -fno-builtin-realloc -fno-builtin-free"
elif [[ "${target}" == *-apple-* ]]; then
    # https://github.com/Homebrew/homebrew-core/blob/5e8ab5f092bd4dfe9658bab6d86150e26a326de3/Formula/gperftools.rb#L28
    export CXXFLAGS=-D_XOPEN_SOURCE
elif [[ "${target}" == *-freebsd* ]]; then
    # Fix the error
    #   undefined reference to `backtrace_symbols'
    export LDFLAGS="-lexecinfo"
fi

./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target}
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:i686, libc=:glibc),
    Linux(:x86_64, libc=:glibc),
    Linux(:aarch64, libc=:glibc),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf),
    Linux(:powerpc64le, libc=:glibc),
    MacOS(:x86_64),
    FreeBSD(:x86_64)
]
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("libtcmalloc_debug", :libtcmalloc_debug),
    LibraryProduct("libtcmalloc", :libtcmalloc),
    LibraryProduct("libtcmalloc_and_profiler", :libtcmalloc_and_profiler),
    LibraryProduct("libtcmalloc_minimal", :libtcmalloc_minimal),
    LibraryProduct("libtcmalloc_minimal_debug", :libtcmalloc_minimal_debug),
    LibraryProduct("libprofiler", :libprofiler)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("LibUnwind_jll"),
]

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


text:
module TryTests

include("test_core.jl")
include("test_base.jl")
include("test_errortrace.jl")
include("test_tools.jl")
include("test_inferrability.jl")
include("test_show.jl")
include("test_aqua.jl")
include("test_doctest.jl")

end  # module TryTests


text:
using DSGE, JLD, DifferentialEquations, StateSpaceRoutines, Statistics
include("../../../../src/estimate/ct_filters/ct_kalman_filter.jl")

# This file will test how well this Kalman Filter finds the
# true outcome when we generate new data many different times

# Solve model
m = KrusellSmithCT()
global T, R, C, inverse_basis, basis = DSGE.solve(m)
# Get quarterly data using the truth to change from reduced to full basis
trials = load("many_fine_data_red_basis.jld", "trials")
ntrials = length(keys(trials))
global states = trials[1]
global data = zeros(convert(Dims, (1, (length(states) - 1) / 90)))

# created 50 years worth of daily data -> transform into quarterly data for discrete observations
for i = 1:length(data)
    full_state = inverse_basis * states[1+ (i-1)*90]
    data[i] = full_state[m.endogenous_states[:output][1]]
end

# Compute logliks
σ_vec = collect(range(0.005, stop=.009, length=20))
σ_vec = sort([σ_vec; .007])
loglik_mat = zeros(length(σ_vec), ntrials)

# Run once to compile
# Define measurement equation
global Δ_t = 1. # length of time traversed b/n observations
global dt = 1/.90 # discretization fineness in simulated data
global D = [0.]
global Z = zeros(1, size(inverse_basis, 1))
global Z[m.endogenous_states[:output][1]] = 1
global Z = Z * inverse_basis # change to reduced basis
global E = zeros(1, 1) * Δ_t # Shock is Browian motion -> variance is Δ_t b/c this is quarterly observation data
global Q = Matrix{Float64}(I, 1, 1) * dt    # Shock is Brownian motion -> variance is dt b/c this is state data

global out = ct_kalman_filter(data, T, R, C, Q, Z, D, E, 1.)
true_lik = sum(out[1])

# Start timing
max_σ = zeros(ntrials)
max_lik = zeros(ntrials)
for n = 1:ntrials
    global states = trials[n] # get trial data
    global data = zeros(convert(Dims, (1, (length(states) - 1) / 90)))
    # created 2 years worth of daily data -> transform into quarterly data for discrete observations
    for i = 1:length(data)
        full_state = inverse_basis * states[1+ (i-1)*90]
        data[i] = full_state[m.endogenous_states[:output][1]]
    end

    for i = 1:length(σ_vec)
        # Recompute steady state solution
        m[:σ_tfp] = σ_vec[i]
        steadystate!(m)
        global T, R, C, inverse_basis, basis = DSGE.solve(m)

        # Define measurement equation
        global Δ_t = 1. # length of time traversed b/n observations
        global dt = 1/.90 # discretization fineness in simulated data
        global  D = [0.]
        global Z = zeros(1, size(inverse_basis, 1))
        global Z[m.endogenous_states[:output][1]] = 1
        global Z = Z * inverse_basis # change to reduced basis
        global E = zeros(1, 1) * Δ_t # Shock is Browian motion -> variance is Δ_t b/c this is quarterly observation data
        global Q = Matrix{Float64}(I, 1, 1) * dt    # Shock is Brownian motion -> variance is dt b/c this is state data
        global out = ct_kalman_filter(data, T, R, C, Q, Z, D, E, 1.)
        loglik_mat[i, n] = sum(out[1])
    end
    max_ind = findall(maximum(loglik_mat[:, n]) .== loglik_mat[:, n])[1]
    max_lik[n] = maximum(loglik_mat[:, n])
    max_σ[n] = σ_vec[max_ind]
end

# Report mean abs distance (true lik sigma- estimated max lik sigma), st dev, and max
abs_dist = max_σ .- .007
st_dev = std(max_σ .- .007)
max_abs_dist = maximum(abs.(max_σ .- .007))

println("Mean absolute distance")
println(mean(abs_dist))

println("\n")

println("Standard deviation of error")
println(st_dev)
println("\n")

println("Max absolute error")
println(max_abs_dist)

save("ct_results.jld", "abs_dist", abs_dist, "st_dev", st_dev, "max_abs_dist", max_abs_dist, "max_σ", max_σ, "true_σ", .007, "σ_vec", σ_vec, "max_lik", max_lik, "true_lik", true_lik)