text:
"""
`module  JAC.BiOrthogonal`
    ... a submodel of JAC that contains all methods for a bi-orthogonal transformation of two non-orthogonal orbitals set.
"""
module BiOrthogonal

    using Printf, JAC, ..Basics,  ..Defaults, ..ManyElectron, ..Radial

    """
    `BiOrthogonal.computeTransformation(leftMp::Multiplet, rightMp::Multiplet)`
        ... computes the bi-orthogonal transformation of the two given multiplets by rotating the radial functions
            by counter-rotating the corresponding CI coefficients. A tuple of two multiplets
            tpl(newLeftMp::::Multiplet, newRightMp::Multiplet) is returned.
    """
    function  computeTransformation(leftMp::Multiplet, rightMp::Multiplet)
        return( 1.0 )
    end

    """
    `BiOrthogonal.computeTransformationMatrices()`
        ... computes the transformation matrices for rotating the radial functions and for counter-rotating CI coefficients.
    """
    function  computeTransformationMatrices()
        return( 1.0 )
    end

    """
    `BiOrthogonal.generateBiorthogonalShellMatrices()`
        ... generates the biorthogonal forms of initial and final shells.
    """
    function  generateBiorthogonalShellMatrices()
        return( 1.0 )
    end

    """
    `BiOrthogonal.generateCounterRotatingCiMatrices()`
        ... generates the counter-rotating matrices for the CI coefficients of the initial and final states.
    """
    function  generateCounterRotatingCiMatrices()
        return( 1.0 )
    end

end # module

#========================================================
Comments & thoughts on the implementation of this module
--------------------------------------------------------

++  I looked briefly through your rbiotransform implementation in GRASP2018 to get a very first impression which basic ingredients
    are needed for such an implementation. A large fraction of code refers to 'internal data handling of quantum numbers and
    functions' which might become obsolete or can be coded much more compact in Julia.
++  Moreover, for the manipulation of matrices (multiplication, inversion, ...), etc., we shall use the internal Julia features
    as much as possible. This will likely reduced the size of the code considerably since some of the procedures in rbiotransform
    can be written as a single line.

++  Overall, therefore, I belief that the code could be made more compact by a factor 3-7, and also more easy to check.

++  A first goal should be to write-down the basic formulas for and in the notation of the UserGuide of JAC.
    This will enable us to bring the notation and internal use much closer to each other; I noticed that large parts in the
    rbiotransform implementation in GRASP2018 are highly technical with very little use of the underlying physics language
    which is required to formulate the transformation itself.

++  I will be happy (and grateful) if we can agree about proper names of functions and if we make use of the
    standard Julia features.

++  Please, use ? AngularMomentum.Wigner_6j(a, b, c, d, e, f) and similar for Wigner_3j and _9j symbols to get the Wigner symbols
    if needed.

++  Let's use CamelCase notation an start all variables with a lowercase letter; arrays are usually given the same name
    with an additional s:  coeff::Coefficient1p  vs.  coeffs::Array{Coefficient1p,1}
++  All angular momenta should be of type j::AngularJ64 (for j-values) and of type m::AngularM64 (for m-values,
    and if they appear at all in your implementation), please.
++  For a given subhsell sh::Subshell, get n and kappa simply by: sh.n, sh.kappa
++  For a given subhsell sh::Subshell, get l and j simply by:     Basics.subshell_l(sh), Basics.subshell_j(sh)

++  Look for ?CsfR ... to get the underlying definition of a CSF in JAC; it allows a direct access to J, parity,
    the occupation, seniority and all angular momenta from the coupling.

++  if csf.useStandardSubshells  subshells = Defaults.GBL_STANDARD_SUBSHELL_LIST
    else                         subshells = csf.subshells
    end
    In the latter case, these subshells could be different for the leftCsf and rightCsf.

========================================================#


text:
@testset "TestData" begin

    force = get(ENV, "PS_FORCE_DOWNLOAD", "false")
    if force == "true"
        force = true
    elseif force == "false"
        force = false
    else
        error("PS_FORCE_DOWNLOAD must be 'true' or 'false'")
    end

    base = abspath(joinpath(@__DIR__, ".."))
    directory = download(PowerSystems.TestData; folder=base, force=force) |> abspath
    @test directory == joinpath(base, "data")

end # testset

ext: jl


text:
using Plots
using Statistics
using ProgressMeter
using FileIO
using JLD2

# These functions are for grid searches.

"""
    sensitivity

plots a sensitivity curve over sweep arg with all settings producted according to product_args
"""
function sensitivity(exp_loc,
                     sweep_arg::String,
                     product_args::Vector{String};
                     results_file="results.jld2",
                     clean_func=identity,
                     ci_const = 1.96,
                     sweep_args_clean=identity,
                     save_dir="sensitivity",
                     ylim=nothing)

    gr()

    if exp_loc[end] == '/'
        exp_loc = exp_loc[1:end-1]
    end
    head_dir = dirname(exp_loc)

    ic = ItemCollection(exp_loc)
    diff_dict = diff(ic.items)
    args = Iterators.product([diff_dict[arg] for arg in product_args]...)

    p1 = ProgressMeter.Progress(length(args), 0.1, "Args: ", offset=0)

    for arg in args

        plt=nothing
        μ = zeros(length(diff_dict[sweep_arg]))
        σ = zeros(length(diff_dict[sweep_arg]))

        p2 = ProgressMeter.Progress(length(diff_dict[sweep_arg]), 0.1, "$(sweep_arg): ", offset=1)
        for (idx, s_a) in enumerate(diff_dict[sweep_arg])
            search_dict = Dict(sweep_arg=>s_a, [product_args[idx]=>key for (idx, key) in enumerate(arg)]...)
            _, hashes, _ = search(ic, search_dict)
            # println(search_dict)
            # println(length(hashes))
            μ_runs = zeros(length(hashes))
            for (idx_d, d) in enumerate(hashes)

                if isfile(joinpath(head_dir, d, results_file))
                    results = load(joinpath(head_dir, d, results_file))
                    μ_runs[idx_d] = clean_func(results)
                # catch e
                else
                    # println(joinpath(head_dir, d, results_file))
                    μ_runs[idx_d] = Inf
                end

            end
            μ[idx] = mean(μ_runs)
            # println(μ)
            σ[idx] = ci_const * std(μ_runs)/sqrt(length(μ_runs))
            next!(p2)
        end

        if plt == nothing
            plt = plot(sweep_args_clean(diff_dict[sweep_arg]), μ, yerror=σ, ylim=ylim)
        else
            plot!(plt, sweep_args_clean(diff_dict[sweep_arg]), μ, yerror=σ)
        end

        if !isdir(joinpath(exp_loc, save_dir))
            mkdir(joinpath(exp_loc, save_dir))
        end

        save_file_name = join(["$(key)_$(arg[idx])" for (idx, key) in enumerate(product_args)], "_")

        savefig(plt, joinpath(exp_loc, save_dir, "$(save_file_name).pdf"))
        next!(p1)
    end

end

"""
    sensitivity_multiline

plots a sensitivity curve over sweep arg with all settings producted according to product_args with lines with args according to line_arg
"""
function sensitivity_multiline(exp_loc, sweep_arg::String, line_arg::String, product_args::Vector{String};
                               results_file="results.jld2", clean_func=identity,
                               sweep_args_clean=identity, save_dir="sensitivity_line",
                               ylim=nothing, ci_const = 1.96, kwargs...)

    gr()

    if exp_loc[end] == '/'
        exp_loc = exp_loc[1:end-1]
    end
    head_dir = dirname(exp_loc)

    ic = ItemCollection(exp_loc)
    diff_dict = diff(ic.items)
    args = Iterators.product([diff_dict[arg] for arg in product_args]...)

    p1 = ProgressMeter.Progress(length(args), 0.1, "Args: ", offset=0)

    for arg in args

        plt=nothing

        p2 = ProgressMeter.Progress(length(diff_dict[line_arg]), 0.1, "$(line_arg): ", offset=1)

        for (idx_line, l_a) in enumerate(diff_dict[line_arg])

            μ = zeros(length(diff_dict[sweep_arg]))
            σ = zeros(length(diff_dict[sweep_arg]))

            p3 = ProgressMeter.Progress(length(diff_dict[sweep_arg]), 0.1, "$(sweep_arg): ", offset=2)
            for (idx, s_a) in enumerate(diff_dict[sweep_arg])
                search_dict = Dict(sweep_arg=>s_a, line_arg=>l_a, [product_args[idx]=>key for (idx, key) in enumerate(arg)]...)
                _, hashes, _ = search(ic, search_dict)
                μ_runs = zeros(length(hashes))
                for (idx_d, d) in enumerate(hashes)
                    if isfile(joinpath(head_dir, d, results_file))
                        results = load(joinpath(head_dir, d, results_file))
                        μ_runs[idx_d] = clean_func(results)
                        # catch e
                    else
                        # println(joinpath(head_dir, d, results_file))
                        μ_runs[idx_d] = Inf
                    end
                end
                μ[idx] = mean(μ_runs)
                σ[idx] = ci_const * std(μ_runs)/sqrt(length(μ_runs))
                next!(p3)
            end

            if plt == nothing
                plt = plot(sweep_args_clean(diff_dict[sweep_arg]), μ, yerror=σ, ylim=ylim, label="$(line_arg)=$(l_a)"; kwargs...)
            else
                plot!(plt, sweep_args_clean(diff_dict[sweep_arg]), μ, yerror=σ, label="$(line_arg)=$(l_a)"; kwargs...)
            end
            next!(p2)
        end

        if !isdir(joinpath(exp_loc, save_dir))
            mkdir(joinpath(exp_loc, save_dir))
        end

        save_file_name = join(["$(key)_$(arg[idx])" for (idx, key) in enumerate(product_args)], "_")

        savefig(plt, joinpath(exp_loc, save_dir, "$(save_file_name).pdf"))
        next!(p1)
    end

end

"""
    sensitivity_best_arg

plots a sensitivity curve over sweep arg with all settings producted according to product_args selecting the best over best_arg
"""
function sensitivity_best_arg(exp_loc,
                              sweep_arg::String,
                              best_arg::String,
                              product_args::Vector{String};
                              results_file="results.jld2",
                              clean_func=identity,
                              sweep_args_clean=identity,
                              compare=(new, old)->new<old,
                              save_dir="sensitivity_best",
                              ylim=nothing, ci_const = 1.96, kwargs...)

    gr()

    if exp_loc[end] == '/'
        exp_loc = exp_loc[1:end-1]
    end
    head_dir = dirname(exp_loc)

    ic = ItemCollection(exp_loc)
    diff_dict = diff(ic.items)
    args = Iterators.product([diff_dict[arg] for arg in product_args]...)

    p1 = ProgressMeter.Progress(length(args), 0.1, "Args: ", offset=0)

    for arg in args

        plt=nothing

        p2 = ProgressMeter.Progress(length(diff_dict[best_arg]), 0.1, "$(best_arg): ", offset=1)

        μ = zeros(length(diff_dict[sweep_arg]))
        fill!(μ, Inf)
        σ = zeros(length(diff_dict[sweep_arg]))

        for (idx_line, b_a) in enumerate(diff_dict[best_arg])

            p3 = ProgressMeter.Progress(length(diff_dict[sweep_arg]), 0.1, "$(sweep_arg): ", offset=2)
            for (idx, s_a) in enumerate(diff_dict[sweep_arg])
                search_dict = Dict(sweep_arg=>s_a, best_arg=>b_a, [product_args[idx]=>key for (idx, key) in enumerate(arg)]...)
                _, hashes, _ = search(ic, search_dict)
                μ_runs = zeros(length(hashes))
                for (idx_d, d) in enumerate(hashes)
                    if isfile(joinpath(head_dir, d, results_file))
                        results = load(joinpath(head_dir, d, results_file))
                        μ_runs[idx_d] = clean_func(results)
                        # catch e
                    else
                        # println(joinpath(head_dir, d, results_file))
                        μ_runs[idx_d] = Inf
                    end
                end
                if compare(mean(μ_runs), μ[idx])
                    μ[idx] = mean(μ_runs)
                    σ[idx] = ci_const * std(μ_runs)/sqrt(length(μ_runs))
                end
                next!(p3)
            end

            next!(p2)
        end

        if plt == nothing
            plt = plot(sweep_args_clean(diff_dict[sweep_arg]), μ, yerror=σ, ylim=ylim; kwargs...)
        else
            plot!(plt, sweep_args_clean(diff_dict[sweep_arg]), μ, yerror=σ; kwargs...)
        end

        if !isdir(joinpath(exp_loc, save_dir))
            mkdir(joinpath(exp_loc, save_dir))
        end

        save_file_name = join(["$(key)_$(arg[idx])" for (idx, key) in enumerate(product_args)], "_")

        savefig(plt, joinpath(exp_loc, save_dir, "$(save_file_name).pdf"))
        next!(p1)
    end

end

function plot_sens_files(file_list, line_settings_list, save_file="tmp.pdf", ci = 1.97; plot_back=gr, kwargs...)

    plot_back()

    plt = nothing

    for (idx, f) in enumerate(file_list)

        ret = load(f)
        println(ret)

        if plt == nothing
            plt = plot(ret["sens"], ret["avg"], ribbon=ci.*ret["std_err"]; line_settings_list[idx]..., kwargs...)
        else
            plot!(plt, ret["sens"], ret["avg"], ribbon=ci.*ret["std_err"]; line_settings_list[idx]..., kwargs...)
        end
    end

    savefig(plt, save_file)

end

function plot_lc_files(file_list, line_settings_list; save_file="tmp.pdf", ci=1.97, n=1, clean_func=identity, plot_back=gr, ignore_nans=false, kwargs...)

    plot_back()

    plt = nothing

    for (idx, f) in enumerate(file_list)

        ret = load(f)
        l = length(clean_func(ret["results"][1]))

        filtered = ret["results"]
        if ignore_nans
            filtered = filter(x->mean(x)!=NaN, ret["results"])
        end
        avg = mean([mean(reshape(clean_func(v), n, Int64(l/n)); dims=1) for v in filtered])'
        std_err = (std([mean(reshape(clean_func(v), n, Int64(l/n)); dims=1) for v in filtered])./sqrt(length(filtered)))'

        x = 0:n:l

        if plt == nothing
            plt = plot(avg, ribbon=ci.*std_err; line_settings_list[idx]..., kwargs...)
        else
            plot!(plt, avg, ribbon=ci.*std_err; line_settings_list[idx]..., kwargs...)
        end
    end

    savefig(plt, save_file)

end


text:
using Chemostat_Folsom2014
using Test

@testset "Chemostat_Folsom2014.jl" begin
    # Write your tests here.
end


text:
export  AbstractMeasure,
        AbstractCanonicalMeasure,
        AbstractQuad,
        AbstractOrthoPoly,
        AbstractCanonicalOrthoPoly,
        AbstractTensor

abstract type AbstractQuad end

abstract type AbstractMeasure end
abstract type AbstractCanonicalMeasure <: AbstractMeasure end

abstract type AbstractOrthoPoly end
abstract type AbstractCanonicalOrthoPoly <: AbstractOrthoPoly end

abstract type AbstractTensor end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Algorithm test for SGP4. All tests are based on [1].
#
#   Notice that only the tests related to the SGP4 in [1] are considered,
#   because the SDP4 (orbit propagator for deep space) are not available yet in
#   SatelliteToolbox.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
#
#   [1] Vallado, D. A., Crawford, P., Hujsak, R., Kelso, T. S (2006). Revisiting
#       Spacetrack Report #3: Rev1. AIAA.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

################################################################################
#                                   Test 01
################################################################################

@testset "Tests from the paper AIAA 2006-6753" begin
    # Read all TLEs that will be used to test.
    tles = read_tle("./sgp4_tests/sgp4_tests.tle")

    st_sgp4_result = []
    SGP4_results = []

    for tle in tles
        filename = @sprintf("./sgp4_tests/aiaa-2006-6753/sgp4_tle_%d_result.txt",
                            tle.sat_num)
        SGP4_results = readdlm(filename; comments=true)

        # Initialize the orbit propagator.
        sgp4d = sgp4_init(tle, sgp4_gc_wgs72)

        t = SGP4_results[:,1]
        @inbounds for k = 1:length(t)

            # Propagate the orbit.
            r_TEME, v_TEME = sgp4!(sgp4d, t[k])

            # Assemble the result vector.
            st_sgp4_result = vcat(t[k], r_TEME, v_TEME)

            # Compare the values.
            @test st_sgp4_result[1] == SGP4_results[k,1]
            @test st_sgp4_result[2] ≈  SGP4_results[k,2] atol=1e-8
            @test st_sgp4_result[3] ≈  SGP4_results[k,3] atol=1e-8
            @test st_sgp4_result[4] ≈  SGP4_results[k,4] atol=1e-8
            @test st_sgp4_result[5] ≈  SGP4_results[k,5] atol=1e-9
            @test st_sgp4_result[6] ≈  SGP4_results[k,6] atol=1e-9
            @test st_sgp4_result[7] ≈  SGP4_results[k,7] atol=1e-9
        end
    end
end


text:
using BigCleaner
using Documenter

DocMeta.setdocmeta!(BigCleaner, :DocTestSetup, :(using BigCleaner); recursive=true)

makedocs(;
    modules=[BigCleaner],
    authors="Ronald Gamez",
    repo="https://github.com/TheRoniOne/BigCleaner.jl/blob/{commit}{path}#{line}",
    sitename="BigCleaner.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://TheRoniOne.github.io/BigCleaner.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/TheRoniOne/BigCleaner.jl",
    devbranch="main",
)


text:
# FITSHeader methods

# -----------------------------------------------------------------------------
# Helper functions
#
# Used here and in other files. Functions that operate on FITSFile
# start with `fits_`.

function try_parse_hdrval(::Type{Bool}, s::Compat.ASCIIString)
    if length(s) == 1
        if s[1] == 'T'
            return Nullable(true)
        elseif s[1] == 'F'
            return Nullable(false)
        end
    end
    return Nullable{Bool}()
end

# Note that trailing whitespace is not significant in FITS header
# keywords, but *leading* whitespace is, so "'    '" parses as " " (a
# single space).  See CFITSIO manual section 4.5 for details.
#
# TODO: parse '' within the string as a single '.
function try_parse_hdrval(::Type{Compat.ASCIIString}, s::Compat.ASCIIString)
    if length(s) < 2 || s[1] != '\'' || s[end] != '\''
        return Nullable{Compat.ASCIIString}()
    end

    i = endof(s) - 1
    while i > 2
        if s[i] != ' '
            return Nullable(s[2:i])
        end
        i -= 1
    end
    return Nullable(s[2:i])
end

try_parse_hdrval(::Type{Float64}, s::Compat.ASCIIString) = tryparse(Float64, s)

# hack for integers in Julia v0.3: tryparse(Int, s) not available in Compat.
if VERSION > v"0.4.0-dev+3864"
    try_parse_hdrval(::Type{Int}, s::Compat.ASCIIString) = tryparse(Int, s)
else
    try_parse_hdrval(::Type{Int}, s::Compat.ASCIIString) = try
        Nullable(parseint(s))
    catch e
        Nullable{Int}()
    end
end

# Try to parse the header value as any type
function try_parse_hdrval(s::Compat.ASCIIString)
    length(s) == 0 && return Nullable(nothing)

    nb = try_parse_hdrval(Bool, s)
    isnull(nb) || return nb

    ns = try_parse_hdrval(Compat.ASCIIString, s)
    isnull(ns) || return ns

    ni = try_parse_hdrval(Int, s)
    isnull(ni) || return ni

    nf = try_parse_hdrval(Float64, s)
    isnull(nf) || return nf

    return Nullable{Any}()
end

# functions for displaying header values in show(io, header)
hdrval_repr(v::Bool) = v ? "T" : "F"
hdrval_repr(v::Compat.ASCIIString) = @sprintf "'%s'" v
hdrval_repr(v::@compat(Union{AbstractFloat, Integer})) = string(v)

# returns one of: Compat.ASCIIString, Bool, Int, Float64, nothing
# (never error)
function parse_header_val(s::Compat.ASCIIString)
    nval = try_parse_hdrval(s)
    return isnull(nval) ? s : get(nval)
end

# Try to read the raw keys in order given; returns Nullable.
# (null if no key exists or if parsing an existing key is unsuccessful.)
function fits_try_read_keys{T}(f::FITSFile, ::Type{T}, keys)
    status = Cint[0]
    value = Vector{UInt8}(71)
    for key in keys
        ccall((:ffgkey, libcfitsio), Cint,
              (Ptr{Void},Ptr{UInt8},Ptr{UInt8},Ptr{UInt8},Ptr{Cint}),
              f.ptr, key, value, C_NULL, status)

        # If the key is found, return it. If there was some other error
        # besides key not found, throw an error.
        if status[1] == 0
            return try_parse_hdrval(T, unsafe_string(pointer(value)))
        elseif status[1] != 202
            error(fits_get_errstatus(status[1]))
        end
    end
    return Nullable{T}()
end

# Build a string with extension keywords, if present.
# This is a helper function for show(::HDU).
const EXTNAME_KEYS = ["EXTNAME", "HDUNAME"]
const EXTVER_KEYS = ["EXTVER", "HDUVER"]
fits_try_read_extname(f::FITSFile) =
    fits_try_read_keys(f, Compat.ASCIIString, EXTNAME_KEYS)
fits_try_read_extver(f::FITSFile) = fits_try_read_keys(f, Int, EXTVER_KEYS)

function fits_get_ext_info_string(f::FITSFile)
    extname = fits_try_read_extname(f)
    extver = fits_try_read_extver(f)
    if !isnull(extname) && !isnull(extver)
        return " (name=$(repr(get(extname))), ver=$(get(extver)))"
    elseif !isnull(extname)
        return " (name=$(repr(get(extname))))"
    end
    return ""
end

# Return indicies of reserved keys in a header.
# This is more complex than you would think because some reserved keys
# are only reserved when other keys are present. Also, in general a key
# may appear more than once in a header.
const RESERVED_KEYS = ["SIMPLE","EXTEND","XTENSION","BITPIX","PCOUNT","GCOUNT",
                       "THEAP","EXTNAME","BSCALE","BZERO","BLANK",
                       "ZQUANTIZ","ZDITHER0","ZIMAGE","ZCMPTYPE","ZSIMPLE",
                       "ZTENSION","ZPCOUNT","ZGCOUNT","ZBITPIX","ZEXTEND",
                       "CHECKSUM","DATASUM"]
function reserved_key_indicies(hdr::FITSHeader)
    nhdr = length(hdr)
    indicies = Int[]
    for i=1:nhdr
        if in(hdr.keys[i], RESERVED_KEYS)
            push!(indicies, i)
        end
    end

    # Note that this removes anything matching NAXIS\d regardless of # of axes.
    if in("NAXIS", hdr.keys)
        for i=1:nhdr
            if ismatch(r"^NAXIS\d*$", hdr.keys[i])
                push!(indicies, i)
            end
        end
    end

    if in("ZNAXIS", hdr.keys)
        for i=1:nhdr
            if (ismatch(r"^ZNAXIS\d*$", hdr.keys[i]) ||
                ismatch(r"^ZTILE\d*$", hdr.keys[i]) ||
                ismatch(r"^ZNAME\d*$", hdr.keys[i]) ||
                ismatch(r"^ZVAL\d*$", hdr.keys[i]))
                push!(indicies, i)
            end
        end
    end

    if in("TFIELDS", hdr.keys)
        for i=1:nhdr
            for re in [r"^TFORM\d*$", r"^TTYPE\d*$", r"^TDIM\d*$",
                       r"^TUNIT\d*$", r"^TSCAL\d*$", r"^TZERO\d*$",
                       r"^TNULL\d*$", r"^TDISP\d*$", r"^TDMIN\d*$",
                       r"^TDMAX\d*$", r"^TDESC\d*$", r"^TROTA\d*$",
                       r"^TRPIX\d*$", r"^TRVAL\d*$", r"^TDELT\d*$",
                       r"^TCUNI\d*$", r"^TFIELDS$"]
                if ismatch(re, hdr.keys[i])
                    push!(indicies, i)
                end
            end
        end
    end

    return indicies
end

# Write header to CHDU.
# If `clean` is true, skip writing reserved header keywords.
function fits_write_header(f::FITSFile, hdr::FITSHeader, clean::Bool=true)
    indicies = clean? reserved_key_indicies(hdr): Int[]
    for i=1:length(hdr)
        if clean && in(i, indicies)
            continue
        end
        if hdr.keys[i] == "COMMENT"
            fits_write_comment(f, hdr.comments[i])
        elseif hdr.keys[i] == "HISTORY"
            fits_write_history(f, hdr.comments[i])
        elseif hdr.comments[i] == ""
            fits_update_key(f, hdr.keys[i], hdr.values[i])
        else
            fits_update_key(f, hdr.keys[i], hdr.values[i], hdr.comments[i])
        end
    end
end

# -----------------------------------------------------------------------------
# Public API

function read_key(hdu::HDU, key::Integer)
    fits_assert_open(hdu.fitsfile)
    fits_movabs_hdu(hdu.fitsfile, hdu.ext)
    keyout, value, comment = fits_read_keyn(hdu.fitsfile, key)
    keyout, parse_header_val(value), comment
end

function read_key(hdu::HDU, key::Compat.ASCIIString)
    fits_assert_open(hdu.fitsfile)
    fits_movabs_hdu(hdu.fitsfile, hdu.ext)
    value, comment = fits_read_keyword(hdu.fitsfile, key)
    parse_header_val(value), comment
end

function read_header(hdu::HDU)
    fits_assert_open(hdu.fitsfile)
    fits_movabs_hdu(hdu.fitsfile, hdu.ext)

    # Below, we use a direct call to ffgkyn so that we can keep reusing the
    # same buffers.
    key = Vector{UInt8}(81)
    value = Vector{UInt8}(81)
    comment = Vector{UInt8}(81)
    status = Cint[0]

    nkeys, morekeys = fits_get_hdrspace(hdu.fitsfile)

    # Initialize output arrays
    keys = Vector{Compat.ASCIIString}(nkeys)
    values = Vector{Any}(nkeys)
    comments = Vector{Compat.ASCIIString}(nkeys)
    for i=1:nkeys
        ccall((:ffgkyn,libcfitsio), Cint,
              (Ptr{Void},Cint,Ptr{UInt8},Ptr{UInt8},Ptr{UInt8},Ptr{Cint}),
              hdu.fitsfile.ptr, i, key, value, comment, status)
        keys[i] = unsafe_string(pointer(key))
        values[i] = parse_header_val(unsafe_string(pointer(value)))
        comments[i] = unsafe_string(pointer(comment))
    end
    fits_assert_ok(status[1])
    FITSHeader(keys, values, comments)
end

function read_header(hdu::HDU, ::Type{Compat.ASCIIString})
    # Return the header as a raw string.

    fits_assert_open(hdu.fitsfile)
    fits_movabs_hdu(hdu.fitsfile, hdu.ext)

    fits_hdr2str(hdu.fitsfile)
end

length(hdr::FITSHeader) = length(hdr.keys)
haskey(hdr::FITSHeader, key::Compat.ASCIIString) = in(key, hdr.keys)
keys(hdr::FITSHeader) = hdr.keys
values(hdr::FITSHeader) = hdr.values
getindex(hdr::FITSHeader, key::Compat.ASCIIString) = hdr.values[hdr.map[key]]
getindex(hdr::FITSHeader, i::Integer) = hdr.values[i]

function setindex!(hdr::FITSHeader, value::Any, key::Compat.ASCIIString)
    if in(key, hdr.keys)
        hdr.values[hdr.map[key]] = value
    else
        push!(hdr.keys, key)
        push!(hdr.values, value)
        push!(hdr.comments, "")
        hdr.map[key] = length(hdr.keys)
    end
end

function setindex!(hdr::FITSHeader, value::Any, i::Integer)
    hdr.values[i] = value
end

# Comments
get_comment(hdr::FITSHeader, key::Compat.ASCIIString) = hdr.comments[hdr.map[key]]
get_comment(hdr::FITSHeader, i::Integer) = hdr.comments[i]
function set_comment!(hdr::FITSHeader, key::Compat.ASCIIString, comment::Compat.ASCIIString)
    hdr.comments[hdr.map[key]] = comment
end
function set_comment!(hdr::FITSHeader, i::Integer, comment::Compat.ASCIIString)
    hdr.comments[i] = comment
end

# Display the header
function show(io::IO, hdr::FITSHeader)
    for i=1:length(hdr)
        cl = length(hdr.comments[i])
        if hdr.keys[i] == "COMMENT" || hdr.keys[i] == "HISTORY"
            if cl > 71
                @printf io "%s %s\n" hdr.keys[i] hdr.comments[i][1:71]
            else
                @printf io "%s %s\n" hdr.keys[i] hdr.comments[i]
            end
        else
            @printf io "%-8s" hdr.keys[i]
            if hdr.values[i] === nothing
                print(io, "                      ")
                rc = 50  # remaining characters on line
            else
                val = hdrval_repr(hdr.values[i])
                @printf io "= %20s" val
                rc = length(val) <= 20 ? 50 : 70 - length(val)
            end

            if cl > 0
                if cl > rc - 3
                    @printf io " / %s" hdr.comments[i][1:rc-3]
                else
                    @printf io " / %s" hdr.comments[i]
                end
            end
            print(io, "\n")
        end
    end
end


text:
using Gadfly
using Distributions

plot_beta(α, β) =
    plot(x -> pdf(Beta(α, β), x), 0, 1)

main(window) = begin
    push!(window.assets, "widgets")

    αᵗ = Signal(1.0)
    βᵗ = Signal(1.0)

    vbox(md"## Static Plot",
        drawing(4Gadfly.inch, 2Gadfly.inch, plot(sin, 0, 25)),
        md"## Dynamic plot",
        hbox("Alpha: " |>
            width(4em), slider(1:100) >>> αᵗ) |>
            packacross(center),
        hbox("Beta: "  |>
            width(4em), slider(1:100) >>> βᵗ) |>
            packacross(center),
        map(αᵗ, βᵗ) do α, β
            plot_beta(α,β) |> drawing(4Gadfly.inch, 3Gadfly.inch)
        end
    ) |> Escher.pad(2em)
end


text:
using GeoData, Test
using GeoData: cleankeys, shiftindexloci

@testset "cleankeys" begin
    @test cleankeys(["a", "b", "c"]) == (:a, :b, :c)
    @test cleankeys(("a", "b", "c")) == (:a, :b, :c)
end

@testset "shiftindexloci" begin
    dim = X(1.0:3.0; mode=Sampled(Ordered(), Regular(1.0), Intervals(Center())))
    @test val(shiftindexloci(Start(), dim)) == 0.5:1.0:2.5
    @test val(shiftindexloci(End(), dim)) == 1.5:1.0:3.5
    @test val(shiftindexloci(Center(), dim)) == 1.0:1.0:3.0
    dim = X([3, 4, 5]; mode=Sampled(Ordered(), Regular(1), Intervals(Start())))
    @test val(shiftindexloci(End(), dim)) == [4, 5, 6]
    @test val(shiftindexloci(Center(), dim)) == [3.5, 4.5, 5.5]
    @test val(shiftindexloci(Start(), dim)) == [3, 4, 5]
    dim = X([3, 4, 5]; mode=Sampled(Ordered(), Regular(1), Intervals(End())))
    @test val(shiftindexloci(End(), dim)) == [3, 4, 5]
    @test val(shiftindexloci(Center(), dim)) == [2.5, 3.5, 4.5]
    @test val(shiftindexloci(Start(), dim)) == [2, 3, 4]
end


text:
const STORE_CONTAINER_DUALS = :duals
const STORE_CONTAINER_PARAMETERS = :parameters
const STORE_CONTAINER_VARIABLES = :variables
const STORE_CONTAINERS =
    Set((STORE_CONTAINER_DUALS, STORE_CONTAINER_PARAMETERS, STORE_CONTAINER_VARIABLES))

"""
Provides storage of simulation data
"""
abstract type SimulationStore end

# Required methods:
# - open_store
# - Base.isopen(store::SimulationStore)
# - Base.close(store::SimulationStore)
# - Base.flush(store::SimulationStore)
# - get_params(store::SimulationStore)
# - initialize_problem_storage!
# - list_fields(store::SimulationStore, problem::Symbol, container_type::Symbol)
# - list_problems(store::SimulationStore)
# - log_cache_hit_percentages(store::SimulationStore)
# - write_result!
# - read_result!
# - write_optimizer_stats!
# - read_problem_optimizer_stats

struct SimulationStoreProblemRequirements
    duals::Dict{Symbol, Dict{String, Any}}
    parameters::Dict{Symbol, Dict{String, Any}}
    variables::Dict{Symbol, Dict{String, Any}}
end

function SimulationStoreProblemRequirements()
    return SimulationStoreProblemRequirements(
        Dict{Symbol, Dict{String, Any}}(),
        Dict{Symbol, Dict{String, Any}}(),
        Dict{Symbol, Dict{String, Any}}(),
    )
end

struct SimulationStoreProblemParams
    num_executions::Int
    horizon::Int
    interval::Dates.Period
    resolution::Dates.Period
    end_of_interval_step::Int
    base_power::Float64
    system_uuid::Base.UUID

    function SimulationStoreProblemParams(
        num_executions,
        horizon,
        interval,
        resolution,
        end_of_interval_step,
        base_power,
        system_uuid,
    )
        new(
            num_executions,
            horizon,
            Dates.Millisecond(interval),
            Dates.Millisecond(resolution),
            end_of_interval_step,
            base_power,
            system_uuid,
        )
    end
end

get_num_executions(params::SimulationStoreProblemParams) = params.num_executions
get_horizon(params::SimulationStoreProblemParams) = params.horizon
get_interval(params::SimulationStoreProblemParams) = params.interval
get_resolution(params::SimulationStoreProblemParams) = params.resolution
get_end_of_interval_step(params::SimulationStoreProblemParams) = params.end_of_interval_step
get_base_power(params::SimulationStoreProblemParams) = params.base_power
get_system_uuid(params::SimulationStoreProblemParams) = params.system_uuid

struct SimulationStoreParams
    initial_time::Dates.DateTime
    step_resolution::Dates.Period
    num_steps::Int
    # The key order is the problem execution order.
    problems::OrderedDict{Symbol, SimulationStoreProblemParams}

    function SimulationStoreParams(initial_time, step_resolution, num_steps, problems)
        new(initial_time, Dates.Millisecond(step_resolution), num_steps, problems)
    end
end

function SimulationStoreParams(initial_time, step_resolution, num_steps)
    return SimulationStoreParams(
        initial_time,
        step_resolution,
        num_steps,
        OrderedDict{Symbol, SimulationStoreProblemParams}(),
    )
end

function SimulationStoreParams()
    return SimulationStoreParams(
        Dates.DateTime("1970-01-01T00:00:00"),
        Dates.Millisecond(0),
        0,
        OrderedDict{Symbol, SimulationStoreProblemParams}(),
    )
end

get_initial_time(store_params::SimulationStoreParams) = store_params.initial_time
get_problems(store_params::SimulationStoreParams) = store_params.problems