text:
#
# methods.jl -
#
# Methods for OI-FITS package.
#
#------------------------------------------------------------------------------

"""
    OIFITS.is_same(a, b)

yields whether `a` and `b` are sufficiently identical for OI-FITS data.  If `a`
and `b` are strings, trailing spaces and case of letters are irrelevant.  If
`a` and `b` are arrays, they are compared element-wise.  This operator is used
in OI-FITS to merge contents.  For an `OITargetEntry`, the target identifier is
ignored.

"""
is_same(a, b) = (a === a)
is_same(a::Integer, b::Integer) = (a == b)
is_same(a::AbstractFloat, b::AbstractFloat) = (a == b)
is_same(a::Complex{<:AbstractFloat}, b::Complex{<:AbstractFloat}) = (a == b)

function is_same(A::AbstractArray{<:Any,N},
                 B::AbstractArray{<:Any,N}) where {N}
    axes(A) == axes(B) || return false
    @inbounds for i in eachindex(A, B)
        is_same(A[i], B[i]) || return false
    end
    return true
end

function is_same(A::AbstractString, B::AbstractString)
    i_last = lastindex(A)
    j_last = lastindex(B)
    i = firstindex(A)
    j = firstindex(B)
    @inbounds begin
        while ((i ≤ i_last)&(j ≤ j_last)) && (uppercase(A[i]) == uppercase(B[j]))
            i = nextind(A, i)
            j = nextind(B, j)
        end
        while i ≤ i_last
            if !isspace(A[i])
                return false
            end
            i = nextind(A, i)
        end
        while j ≤ j_last
            if !isspace(B[j])
                return false
            end
            j = nextind(B, j)
        end
    end
    return true
end

function is_same(A::OITargetEntry, B::OITargetEntry)
    # Private helper.
    @inline _is_same(A::OITargetEntry, B::OITargetEntry, sym::Symbol) =
        is_same(getfield(A, sym), getfield(B, sym))

    # Purposely ignore target_id.
    return (_is_same(A, B, :target   ) &&
            _is_same(A, B, :raep0    ) &&
            _is_same(A, B, :decep0   ) &&
            _is_same(A, B, :equinox  ) &&
            _is_same(A, B, :ra_err   ) &&
            _is_same(A, B, :dec_err  ) &&
            _is_same(A, B, :sysvel   ) &&
            _is_same(A, B, :veltyp   ) &&
            _is_same(A, B, :veldef   ) &&
            _is_same(A, B, :pmra     ) &&
            _is_same(A, B, :pmdec    ) &&
            _is_same(A, B, :pmra_err ) &&
            _is_same(A, B, :pmdec_err) &&
            _is_same(A, B, :parallax ) &&
            _is_same(A, B, :para_err ) &&
            _is_same(A, B, :spectyp  ) &&
            _is_same(A, B, :category ))
end

"""
    OIFITS.is_same(A, B, sym)

yields whether field `sym` of data-blocks `A` and `B` are sufficiently
identical for OI-FITS data.

"""
@inline is_same(A::T, B::T, sym::Symbol) where {T<:OIDataBlock} =
    is_same(A, B, Val(sym))

is_same(A::T, B::T, ::Val{S}) where {T<:OIDataBlock,S} = begin
    if (defined = isdefined(A, S)) != isdefined(B, S)
        return false
    elseif defined
        return is_same(getfield(A, S), getfield(B, S))
    else
        return true
    end
end

"""
    OIFITS.fix_name(str)

yields string `str` converted to uppercase letters and with trailing spaces
removed.  This is to follow FITS conventions that names comparisons should be
done ignoring the case of letters and trailing spaces.

"""
fix_name(str::AbstractString) = uppercase(rstrip(str))

# Extend A[name] syntax to get an OI_ARRAY, OI_WAVELENGTH, or OI_CORR by its
# name (according to FITS conventions).  The syntax keys(String,A) can be used
# to get an iterator over the string keys of A.
for (type, name) in ((:OI_ARRAY,      :arrname),
                     (:OI_WAVELENGTH, :insname),
                     (:OI_CORR,       :corrname))
    @eval begin
        getindex(A::AbstractArray{<:$type}, key::AbstractString) = begin
            @inbounds for i in eachindex(A)
                if is_same(A[i].$name, key)
                    return A[i]
                end
            end
            throw(KeyError(fix_name(key)))
        end
        keys(::Type{String}, A::AbstractArray{<:$type}) =
            Iterators.map(x -> fix_name(x.$name), A)
        haskey(A::AbstractArray{<:$type}, key::AbstractString) = begin
            @inbounds for i in eachindex(A)
                if is_same(A[i].$name, key)
                    return true
                end
            end
            return false
        end
        get(A::AbstractArray{<:$type}, key::AbstractString, def) = begin
            @inbounds for i in eachindex(A)
                if is_same(A[i].$name, key)
                    return A[i]
                end
            end
            return def
        end
    end
end

"""
    OIFITS.extname(arg)

yields the name of the FITS extension for argument `arg` which can be an
OI-FITS datablock instance or type or a FITS table HDU.  In this latter case,
an empty string is returned if the `EXTNAME` keyword is not found and the
extension name is converted to upper case letters and trailing spaces discarded
otherwise.  This method is not exported, the same result for a data-block
instance `db` is obtained by `db.extname`.

For a data-block instance or type `db`, an optional first argument, `Symbol` or
`String`, may be used to specify the type of the result:

    OIFITS.extname(Symbol, db) -> ext::Symbol
    OIFITS.extname(String, db) -> ext::String

"""
extname(hdu::TableHDU) = fix_name(read_keyword(String, hdu, "EXTNAME", ""))
extname(db::OIDataBlock) = extname(typeof(db))
extname(T::Type{<:OIDataBlock}) = extname(String, T)
extname(S::Type{<:Union{String,Symbol}}, db::OIDataBlock) =
    extname(S, typeof(db))

for sym in (:OI_TARGET,
            :OI_ARRAY,
            :OI_WAVELENGTH,
            :OI_CORR,
            :OI_VIS,
            :OI_VIS2,
            :OI_T3,
            :OI_FLUX,
            :OI_INSPOL)
    @eval begin
        extname(::Type{String}, ::Type{<:$sym}) = $(String(sym))
        extname(::Type{Symbol}, ::Type{<:$sym}) = $(QuoteNode(sym))
    end
end

"""
    OIFITS.get_format(db,  revn=db.revn; throw_errors=false)
    OIFITS.get_format(T,   revn;         throw_errors=false)
    OIFITS.get_format(ext, revn;         throw_errors=false)

all yield OI-FITS definitions for data-block `db`, of data-block type `T`, or
of OI-FITS extension `ext` (a string or a symbol).

"""
get_format(::Type{T}, revn::Integer; kwds...) where {T<:OIDataBlock} =
    get_format(extname(Symbol, T), revn; kwds...)
get_format(db::OIDataBlock, revn::Integer=db.revn; kwds...) =
    get_format(extname(Symbol, db), revn; kwds...)
get_format(extname::Union{Symbol,AbstractString}, revn::Integer; kwds...) =
    get_format(Symbol(extname), revn; kwds...)

#------------------------------------------------------------------------------
# PROPERTIES OF DATA-BLOCKS
#
# We use the `db.key` syntax to give access to other properties (e.g.,
# `extname`, `name`, `eff_wave`, or `eff_band`) than the fields of the
# data-block `db`.
#
# To allow for efficient handling of the `db.key` syntax, we wrap the `key`
# symbol as a `Val(key)` so that each `getproperty` or `setproperty!` method
# can be specific to the value of `key`.

# _properties(T) yields public properties for type `T`.
for T in (OI_TARGET, OI_ARRAY, OI_WAVELENGTH, OI_CORR,
          OI_VIS, OI_VIS2, OI_T3, OI_FLUX, OI_INSPOL)
    isconcretetype(T) || continue
    extra_fields = if T <: Union{OI_ARRAY,OI_WAVELENGTH,OI_CORR}
        (:extname, :name,)
    elseif T <: Union{OI_VIS,OI_VIS2,OI_T3,OI_FLUX}
        (:extname, :eff_wave, :eff_band,)
    else
        (:extname,)
    end
    @eval _properties(::Type{$T}) = $((fieldnames(T)..., extra_fields...))
end

propertynames(db::OIDataBlock) = _properties(typeof(db))

getproperty(db::OIDataBlock, sym::Symbol) = getproperty(db, Val(sym))
getproperty(db::OIDataBlock, ::Val{S}) where {S} = getfield(db, S)
getproperty(db::OIDataBlock, ::Val{:extname}) = extname(typeof(db))

setproperty!(db::OIDataBlock, sym::Symbol, val) =
    setproperty!(db, Val(sym), val)

setproperty!(db::T, ::Val{S}, val) where {S,T<:OIDataBlock} = begin
    try
        setfield!(db, S, convert(fieldtype(T, S), val))
    catch ex
        rethrow_convert_field_error(ex, T, sym, val)
    end
end

for (type, name) in ((:OI_ARRAY,      :arrname),
                     (:OI_WAVELENGTH, :insname),
                     (:OI_CORR,       :corrname))
    @eval begin
        getproperty(db::$type, ::Val{:name}) = db.$name
        setproperty!(db::$type, ::Val{:name}, val) = (db.$name = val)
    end
end

# NOTE: Directly calling `getfield` is necessary to optimize the indirection
#       (speedup: 1.4ns instead of 370ns).
_instr(db::Union{OI_VIS,OI_VIS2,OI_T3,OI_FLUX}) = getfield(db, :instr)
_instr(db::OI_WAVELENGTH) = db
for type in (:OI_VIS, :OI_VIS2, :OI_T3, :OI_FLUX)
    @eval begin
        getproperty(db::$type, ::Val{:eff_wave}) =
            _instr(db).eff_wave
        getproperty(db::$type, ::Val{:eff_band}) =
            _instr(db).eff_band
        setproperty!(db::$type, ::Val{:eff_wave}, val) =
            _instr(db).eff_wave = val
        setproperty!(db::$type, ::Val{:eff_band}, val) =
            _instr(db).eff_band = val
    end
end

@noinline function rethrow_convert_field_error(ex::Exception, T::Type,
                                               sym::Symbol, val)
    if isa(ex, MethodError) && ex.f === convert
        rethrow(ArgumentError(string(
            "Cannot `convert` an object of type `", typeof(val),
            "` to an object of type `", fieldtype(T, sym), "` for field `",
            sym, "` of structure of type `", T, "`")))
    else
        rethrow(ex)
    end
end

#------------------------------------------------------------------------------
# PROPERTIES OF DATA-SETS
#
# Here, we want to make the "private" fields of data-sets (the dictionaries)
# "hidden" to the user.  They can only be retrieved by getfield(ds,sym).

# Public properties for data-sets.
_properties(::Type{OIDataSet}) = (:target, :array, :instr, :correl,
                                  :vis, :vis2, :t3, :flux, :inspol)

# Only the "public" fields of data-sets are part of the properties.
propertynames(ds::OIDataSet) = _properties(OIDataSet)

# An OIDataSet is immutable so explicitely forbid writing a field is not
# necessary but the following is to provide a more explicit message.
@noinline setproperty!(db::OIDataSet, sym::Symbol, val) =
    error("attempt to ", (sym ∈ propertynames(db) ? "set read-only" :
                          "access non-existing or private"),
          " field `", sym, "` in ", nameof(typeof(db)), " instance")

getproperty(db::OIDataSet, sym::Symbol) = getproperty(db, Val(sym))

# By default, accessing any field by the dot notation is forbidden.
getproperty(db::OIDataSet, ::Val{S}) where {S} =
    error("attempt to access non-existing or private field `", S,
          "` in ", nameof(typeof(db)), " instance")

# Each public field is explicitely allowed.
for sym in _properties(OIDataSet)
    quoted_sym = QuoteNode(sym)
    @eval begin
        getproperty(db::OIDataSet, ::Val{$quoted_sym}) =
            getfield(db, $quoted_sym)
    end
end

#------------------------------------------------------------------------------
# COPY DATA-BLOCKS AND DATA-SETS

# Copy the outer structure of an OI-FITS data-block.
function copy(A::T) where {T<:OIDataBlock}
    B = T(undef)
    for sym in fieldnames(T)
        if isdefined(A, sym)
            setfield!(B, sym, getfield(A, sym))
        end
    end
    return B
end

function copy(A::OIDataSet)
    # Create new instance and copy contents.  Including the dictionary for
    # re-writing target identifiers.  The copy is thus in the same "state" as
    # the original.  Note that using `map` is ~ 3.6 times slower here.
    B = OIDataSet()
    _copy!(B, A, :target)
    _copy!(B, A, :target_dict)
    _copy!(B, A, :target_id_map)
    _copy!(B, A, :array)
    _copy!(B, A, :array_dict)
    _copy!(B, A, :instr)
    _copy!(B, A, :instr_dict)
    _copy!(B, A, :correl)
    _copy!(B, A, :correl_dict)
    _copy!(B, A, :vis)
    _copy!(B, A, :vis2)
    _copy!(B, A, :t3)
    _copy!(B, A, :flux)
    _copy!(B, A, :inspol)
    return B
end

function _copy!(dest::Vector{T}, src::Vector{T}) where {T}
    if dest !== src && length(src) > 0
        copyto!(resize!(dest, length(src)), src)
    end
    return dest
end

function _copy!(dest::Dict{K,V}, src::Dict{K,V}) where {K,V}
    empty!(dest)
    for (key, val) in src
        dest[key] = val
    end
    return dest
end

function _copy!(dest::OI_TARGET, src::OI_TARGET)
    dest.revn = src.revn
    _copy!(dest.list, src.list)
    return dest
end

# @inline is needed here to make this as fast as possible.
# An alternative is to use Val(sym) (without @inline).
@inline _copy!(dest::T, src::T, sym::Symbol) where {T} =
    _copy!(getfield(dest, sym), getfield(src, sym))
#_copy!(dest::T, src::T, ::Val{S}) where {T,S} =
#    _copy!(getfield(dest, S), getfield(src, S))

#------------------------------------------------------------------------------
# PUSH DATA-BLOCKS IN DATA-SETS

# Create empty instance, then push all data-blocks which must be ordered.
OIDataSet(args::OIDataBlock...) = push!(OIDataSet(), args...)

function empty!(ds::OIDataSet)
    @inline _empty!(A, sym::Symbol) = empty!(getfield(A, sym))
    empty!(ds.target.list)
    ds.target.revn = 0
    _empty!(ds, :target_dict)
    _empty!(ds, :target_id_map)
    _empty!(ds, :array)
    _empty!(ds, :array_dict)
    _empty!(ds, :instr)
    _empty!(ds, :instr_dict)
    _empty!(ds, :correl)
    _empty!(ds, :correl_dict)
    _empty!(ds, :vis)
    _empty!(ds, :vis2)
    _empty!(ds, :t3)
    _empty!(ds, :flux)
    _empty!(ds, :inspol)
    return ds
end

for (type, pass) in ((:OI_TARGET,     1),
                     (:OI_ARRAY,      1),
                     (:OI_WAVELENGTH, 1),
                     (:OI_CORR,       1),
                     (:OI_VIS,        2),
                     (:OI_VIS2,       2),
                     (:OI_T3,         2),
                     (:OI_FLUX,       2),
                     (:OI_INSPOL,     2))
    @eval function _push!(ds::OIDataSet, db::$type, pass::Integer)
        if pass == $pass
            push!(ds, db)
        end
        nothing
    end
end

"""
    push!(ds::OIDataSet, db) -> ds

adds an OI-FITS data-block `db` in the data-set `ds` and returns `ds`.  After
the operation, the contents of `db` is left unchanged but may be partially
shared by `ds`.  Depending on type of `db`, different things can happen:

- If `db` is an `OI_TARGET` instance, the targets from `db` not already in the
  data-set `ds` are added to the list of targets in `ds` and the internal
  dictionary in `ds` mapping target identifiers in `db` to those in `ds` is
  reinitialized for subsequent data-blocks.

- If `db` is an `OI_ARRAY`, `OI_WAVELENGTH`, or `OI_CORR` instance, it is added
  to the corresponding list of data-blocks in the data-set `ds` unless an entry
  with a name matching that of `db` already exists in `ds`.  In this latter
  case, nothing is done except checking that the two data-blocks with matching
  names have the same contents.  This is need to ensure the consistency of the
  data-set.

- If `db` is an `OI_VIS`, `OI_VIS2`, `OI_T3`, `OI_FLUX`, or `OI_INSPOL`
  instance, it is added to the corresponding list of data-blocks in the
  data-set `ds` after having rewritten its target identifiers according to the
  mapping set by the last `OI_TARGET` pushed into the data-set `ds`.  Add
  keyword `rewrite_target_id=false` to avoid rewritting target identifiers.

"""
function push!(ds::OIDataSet, db::OI_TARGET)
    # Check data-block.
    check(db)

    # Get current list of targets in data-set and corresponding mapping of
    # target names to identifiers.
    target_list = ds.target.list
    target_dict = getfield(ds, :target_dict)
    length(target_list) == length(target_dict) || error(
        "inconsistent list and dictionary of targets")

    # Get and reset dictionary to reindex target identifiers in this and
    # subsequent data-blocks.
    target_id_map = getfield(ds, :target_id_map)
    empty!(target_id_map)

    # Insert all targets from the data-block to the data-set.
    for tgt in db.list
        name = fix_name(tgt.target)
        id = Int(tgt.target_id)
        if haskey(target_id_map, id)
            is_same(target_list[target_id_map[id]], tgt) || error(
                "duplicate target identifier ", id, " in ", db.extname,
                " but for different targets")
        end
        if haskey(target_dict, name)
            # The same target already existed in destination.  Check that they
            # have the same parameters.
            index = target_dict[name]
            is_same(target_list[index], tgt) || error(
                "target \"", tgt.target, "\" in ", db.extname, " has ",
                "different parameters than a previous one matching this name")
        else
            # This target is new.
            index = length(target_list) + 1
            push!(target_list, OITargetEntry(tgt; target_id = index))
            target_dict[name] = index
        end
        target_id_map[id] = index
    end

    # Finally fix revision number.
    ds.target.revn = max(ds.target.revn, db.revn)
    return ds
end

function push!(ds::OIDataSet,
               db::T) where {T<:Union{OI_ARRAY,OI_WAVELENGTH,OI_CORR}}
    check(db)
    dict = _dict(T, ds)
    list = _list(T, ds)
    name = fix_name(db.name)
    if haskey(dict, name)
        # Make sure the two data-blocks with the same name are
        # identical.
        assert_identical(db, list[dict[name]])
    else
        # Push new data-block in destination.
        dict[name] = length(push!(list, db))
    end
    return ds
end

function push!(ds::OIDataSet, db::T;
               rewrite_target_id::Bool = true) where {
                   T<:Union{OI_VIS,OI_VIS2,OI_T3,OI_FLUX,OI_INSPOL}}
    # Define private helper.
    function find_depencency(dest_list::Vector{T},
                             dest_dict::Dict{String,Int},
                             dep::AbstractString,
                             optional::Bool = false) where {S,T<:OIDataBlock}
        name = fix_name(dep)
        if haskey(dest_dict, name)
            return dest_list[dest_dict[name]]
        else
            optional || error(
                "no extensions ", extname(T), " match name \"", name, "\"")
            return nothing
        end
    end

    check(db)
    new_db = copy(db) # make a copy to avoid side effects
    if rewrite_target_id
        new_db.target_id = rewrite_indices(db.target_id,
                                           getfield(ds, :target_id_map))
    else
        id_min, id_max = extrema(db.target_id)
        if id_min < 1 || id_max > length(db.target)
            error("out of range target identifier(s) (", id_min, ":", id_max,
                  " ⊈ 1:", length(db.target), ")")
        end
    end
    if isdefined(new_db, :arrname)
        let dep = find_depencency(ds.array, getfield(ds, :array_dict),
                                  new_db.arrname, true)
            if dep !== nothing
                new_db.array = dep
            end
        end
    end
    if isdefined(new_db, :insname)
        new_db.instr = find_depencency(
            ds.instr, getfield(ds, :instr_dict), new_db.insname)
    end
    if isdefined(new_db, :corrname)
        new_db.correl = find_depencency(
            ds.correl, getfield(ds, :correl_dict), new_db.corrname)
    end
    push!(_list(T, ds), new_db)
    return ds
end

# `_list(T,ds)` yields the list of entries of type `T` stored in `ds`.
_list(::Type{OI_ARRAY},      ds::OIDataSet) = ds.array
_list(::Type{OI_WAVELENGTH}, ds::OIDataSet) = ds.instr
_list(::Type{OI_CORR},       ds::OIDataSet) = ds.correl
_list(::Type{OI_VIS},        ds::OIDataSet) = ds.vis
_list(::Type{OI_VIS2},       ds::OIDataSet) = ds.vis2
_list(::Type{OI_T3},         ds::OIDataSet) = ds.t3
_list(::Type{OI_FLUX},       ds::OIDataSet) = ds.flux
_list(::Type{OI_INSPOL},     ds::OIDataSet) = ds.inspol
#_list(::Type{OITargetEntry},ds::OIDataSet) = ds.target.list

# `_dict(T,ds)` yields the dictionary associated to entries of type `T` stored
# by the OI-FITS data-set `ds`.
_dict(::Type{OI_ARRAY},      ds::OIDataSet) = getfield(ds, :array_dict)
_dict(::Type{OI_WAVELENGTH}, ds::OIDataSet) = getfield(ds, :instr_dict)
_dict(::Type{OI_CORR},       ds::OIDataSet) = getfield(ds, :correl_dict)
#_dict(::Type{OITargetEntry},ds::OIDataSet) = getfield(ds, :target_dict)

"""
    OIFITS.rewrite_indices(inds, dict) -> inds′

yields array of indices similar to `inds` and rewritten according to the
dictionary `dict` used as a mapping from integers to integers.  If the indices
are left unchanged, the input array is returned.

"""
function rewrite_indices(inds::AbstractArray{<:Integer,N},
                         dict::Dict{<:Integer,<:Integer}) where {N}
    new_inds = similar(inds)
    unchanged = true
    for i in eachindex(new_inds, inds)
        idx = inds[i]
        new_idx = dict[idx]
        new_inds[i] = new_idx
        unchanged &= (new_idx == idx)
    end
    return (unchanged ? inds : new_inds)
end

"""
    OIFITS.assert_identical(A, B)

throws an exception if *named* data-blocks `A` and `B` are not identical.  A
*named* data-block is an `OI_ARRAY`, `OI_WAVELENGTH`, or `OI_CORR` extension
that can be uniquely identifed by its name.

The variant

    OIFITS.assert_identical(A, B, sym)

is to assert that fields `A.sym` and `B.sym` are identical, it shall only be
called if `A` and `B` have the same type and the same name.

""" assert_identical

function assert_identical(A::T, B::T) where {T<:Union{OI_ARRAY,OI_WAVELENGTH,OI_CORR}}
    if !(A === B)
        is_same(A.name, B.name) || throw_assertion_error(
            "two ", extname(T), " data-blocks named \"", A.name,
            "\" and \"", B.name, "\" cannot be identical")
        assert_identical(A, B, *)
    end
end

function assert_identical(A::OI_ARRAY, B::OI_ARRAY, ::typeof(*))
    assert_identical(A, B, :revn)
    assert_identical(A, B, :frame)
    assert_identical(A, B, :arrayx)
    assert_identical(A, B, :arrayy)
    assert_identical(A, B, :arrayz)
    assert_identical(A, B, :tel_name)
    assert_identical(A, B, :sta_name)
    assert_identical(A, B, :sta_index)
    assert_identical(A, B, :diameter)
    assert_identical(A, B, :staxyz)
    assert_identical(A, B, :fov)
    assert_identical(A, B, :fovtype)
end

function assert_identical(A::OI_WAVELENGTH, B::OI_WAVELENGTH, ::typeof(*))
    assert_identical(A, B, :revn)
    assert_identical(A, B, :eff_wave)
    assert_identical(A, B, :eff_band)
end

function assert_identical(A::OI_CORR, B::OI_CORR, ::typeof(*))
    assert_identical(A, B, :revn)
    assert_identical(A, B, :ndata)
    assert_identical(A, B, :iindx)
    assert_identical(A, B, :jindx)
    assert_identical(A, B, :corr)
end

@inline assert_identical(A::T, B::T, sym::Symbol) where {T<:NamedDataBlock} =
    assert_identical(A, B, Val(sym))

assert_identical(A::T, B::T, val::Val{sym}) where {T<:NamedDataBlock,sym} =
    is_same(A, B, val) || throw_assertion_error(
        "two ", extname(T), " extensions both named \"", fix_name(A.name),
        " have different field `", sym, "`")

throw_assertion_error(msg::AbstractString) = throw(AssertionError(msg))
@noinline throw_assertion_error(args...) =
    throw_assertion_error(string(args...))

"""
    OIFITS.check(db) -> (nchns, nfrms)

checks that mandatory fields of data-block `db` are defined and that defined
fields have correct sizes.  The number of spectral channels `nchns` and of
temporal frames `nfrms` are returned (with a value of `-1` if not relevant for
the type of `db`).

    OIFITS.check(ds)

checks all data-blocks of data-set `ds`.

""" check

function check(ds::OIDataSet)
    check(ds.target)
    for db in ds.array
        check(db)
    end
    for db in ds.instr
        check(db)
    end
    for db in ds.correl
        check(db)
    end
    for db in ds.vis
        check(db)
    end
    for db in ds.vis2
        check(db)
    end
    for db in ds.t3
        check(db)
    end
    for db in ds.flux
        check(db)
    end
    for db in ds.inspol
        check(db)
    end
    nothing
end

function check(db::OIDataBlock)
    nchns = if isdefined(db, :instr)
        length(getfield(db, :instr).eff_wave)
    else
        -1 # number of spectral channels not yet known
    end
    nfrms = -1 # number of temporal frames not yet known
    for spec in get_format(db; throw_errors=true)
        nchns, nfrms = _check(db, Val(spec.symb), spec, nchns, nfrms)
    end
    return nchns, nfrms
end

function check(db::OI_TARGET)
    isdefined(db, :list) || throw_undefined_field(db, :list)
    return -1, -1
end

function _check(db::OIDataBlock, ::Val{S}, spec::FieldDefinition,
                nchns::Int, nfrms::Int) where {S}
    if !isdefined(db, S)
        spec.optional || throw_undefined_field(db, S)
    else
        val = getfield(db, S)
        dims = (isa(val, AbstractString) ? () : size(val))
        ndims_spec = ndims(spec)
        length(dims) == ndims_spec || error(
            "field `", S, "` of ", db.extname, " must be a ",
            (ndims_spec == 0 ? "scalar" : "$ndims_spec-D array"))
        rank = spec.rank
        if ndims_spec ≥ 1
            # Last dimension is number of temporal frames.
            val_nfrms = dims[end]
            if nfrms == -1
                nfrms = val_nfrms
            else
                nfrms == val_nfrms || error(
                    "incompatible number of temporal frames in field `", S,
                    "` of ", db.extname, " (got ", val_nfrms, ", should be ",
                    nfrms, ")")
            end
        end
        if rank > 1 &&  ndims_spec == 2 && spec.type !== :A
            dims[1] == rank || error(
                "invalid 1st dimension in field `", S,
                "` of ", db.extname, " (got ", dims[1], ", should be ",
                rank, ")")
        elseif rank < 0
            # All leading dimensions are equal to the number of spectral
            # channels.
            val_nchns = dims[1]
            for i in 2:length(dims)-1
                dims[i] == val_nchns || error(
                    "leading dimensions in field `", S, "` of ", db.extname,
                    " must be all equal to the numeber of spectral channels,",
                    " size is ", dims)
            end
            if nchns == -1
                nchns = val_nchns
            else
                nchns == val_nchns || error(
                    "incompatible number of spectral channels in field `", S,
                    "` of ", db.extname, " (got ", val_nchns, ", should be ",
                    nchns, ")")
            end
        end
    end
    return nchns, nfrms
end

const FieldName = Union{AbstractString,Symbol}

@noinline throw_undefined_field(name::FieldName) =
    throw(ErrorException("undefined field `$name`"))

@noinline throw_undefined_field(obj::T, name::FieldName) where {T} =
    throw(ErrorException(string(
        "undefined field `", name, "` in object of type `", nameof(T), "`")))

@noinline throw_undefined_field(obj::OIDataBlock, name::FieldName) =
    throw(ErrorException(string(
        "mandatory field `", name, "` undefined in `", obj.extname,
        "` data-block")))

#------------------------------------------------------------------------------
# MERGE DATA-SETS

merge(A::OIDataSet) = copy(A)
merge(A::OIDataSet, args::OIDataSet...) = merge!(OIDataSet(), A, args...)

merge!(A::OIDataSet) = A
merge!(A::OIDataSet, others::OIDataSet...) = begin
    for B in others
        merge!(A, B)
    end
    return A
end

function merge!(A::OIDataSet, B::OIDataSet)
    # First push targets and dependencies.
    push!(A, B.target)
    for db in B.array
        push!(A, db)
    end
    for db in B.instr
        push!(A, db)
    end
    for db in B.correl
        push!(A, db)
    end

    # Then push all other data-blocks.
    for db in B.vis
        push!(A, db)
    end
    for db in B.vis2
        push!(A, db)
    end
    for db in B.t3
        push!(A, db)
    end
    for db in B.flux
        push!(A, db)
    end
    for db in B.inspol
        push!(A, db)
    end

    # Return updated data-set.
    return A
end

#------------------------------------------------------------------------------
# OI_TARGET DATA-BLOCKS

"""
    OITargetEntry([def::OITargetEntry]; kwds...)

yields an entry of the `OI_TARGET` table in an OI-FITS file.  All fields are
specified by keywords.  If template entry `def` is specified, keywords have
default values taken from `def`; otherwise all keywords are mandatory but
`category` which is assumed to be `""` if unspecified.

"""
function OITargetEntry(;
                       target_id ::Integer,
                       target    ::AbstractString,
                       raep0     ::AbstractFloat,
                       decep0    ::AbstractFloat,
                       equinox   ::AbstractFloat,
                       ra_err    ::AbstractFloat,
                       dec_err   ::AbstractFloat,
                       sysvel    ::AbstractFloat,
                       veltyp    ::AbstractString,
                       veldef    ::AbstractString,
                       pmra      ::AbstractFloat,
                       pmdec     ::AbstractFloat,
                       pmra_err  ::AbstractFloat,
                       pmdec_err ::AbstractFloat,
                       parallax  ::AbstractFloat,
                       para_err  ::AbstractFloat,
                       spectyp   ::AbstractString,
                       category  ::AbstractString = empty_string)
    return OITargetEntry(target_id,
                         target,
                         raep0,
                         decep0,
                         equinox,
                         ra_err,
                         dec_err,
                         sysvel,
                         veltyp,
                         veldef,
                         pmra,
                         pmdec,
                         pmra_err,
                         pmdec_err,
                         parallax,
                         para_err,
                         spectyp,
                         category)
end

function OITargetEntry(def::OITargetEntry;
                       target_id ::Integer        = def.target_id,
                       target    ::AbstractString = def.target,
                       raep0     ::AbstractFloat  = def.raep0,
                       decep0    ::AbstractFloat  = def.decep0,
                       equinox   ::AbstractFloat  = def.equinox,
                       ra_err    ::AbstractFloat  = def.ra_err,
                       dec_err   ::AbstractFloat  = def.dec_err,
                       sysvel    ::AbstractFloat  = def.sysvel,
                       veltyp    ::AbstractString = def.veltyp,
                       veldef    ::AbstractString = def.veldef,
                       pmra      ::AbstractFloat  = def.pmra,
                       pmdec     ::AbstractFloat  = def.pmdec,
                       pmra_err  ::AbstractFloat  = def.pmra_err,
                       pmdec_err ::AbstractFloat  = def.pmdec_err,
                       parallax  ::AbstractFloat  = def.parallax,
                       para_err  ::AbstractFloat  = def.para_err,
                       spectyp   ::AbstractString = def.spectyp,
                       category  ::AbstractString = def.category)
    return OITargetEntry(target_id,
                         target,
                         raep0,
                         decep0,
                         equinox,
                         ra_err,
                         dec_err,
                         sysvel,
                         veltyp,
                         veldef,
                         pmra,
                         pmdec,
                         pmra_err,
                         pmdec_err,
                         parallax,
                         para_err,
                         spectyp,
                         category)
end

"""
    OI_TARGET(lst=OITargetEntry[]; revn=0)

yields an `OI_TARGET` data-block.  Optional argument `lst` is a vector of
`OITargetEntry` specifying the targets (none by default).  Keyword `revn`
specifies the revision number.

An instance, say `db`, of `OI_TARGET` has the following properties:

    db.extname # OI-FITS extension name
    db.list    # list of targets
    db.revn    # revision number

and can be used as an iterable or as an array to access the targets
individually by their index or by their name:

    length(db) # the number of targets
    db[i]      # the i-th target
    db[key]    # the target whose name matches string `key`

    for tgt in db; ...; end # loop over all targets

Note that, when an `OI_TARGET` instance is pushed in a data-set, target
identifiers (field `target_id`) are automatically rewritten to be identical to
the index in the list of targets of the data-set.

Standard methods `get` and `haskey` work as expected and according to the type
(integer or string) of the key.  For the `keys` method, the default is to
return an iterator over the target names, but the type of the expected keys
can be specified:

    keys(db)          # iterator over target names
    keys(String, db)  # idem
    keys(Integer, db) # iterator over target indices
    keys(Int, db)     # idem

Call [`OIFITS.get_column`](@ref) to retrieve a given target field for all
targets of `OI_TARGET` data-block in the form of a vector.

""" OI_TARGET

# Make OI_TARGET instances iterable and behave more or less like vectors
# (with properties).
ndims(db::OI_TARGET)   =  ndims(db.list)
size(db::OI_TARGET)    =   size(db.list)
size(db::OI_TARGET, i) =   size(db.list, i)
axes(db::OI_TARGET)    =   axes(db.list)
axes(db::OI_TARGET, i) =   axes(db.list, i)
length(db::OI_TARGET)  = length(db.list)

IndexStyle(db::OI_TARGET) = IndexStyle(OI_TARGET)
IndexStyle(::Type{OI_TARGET}) = IndexStyle(fieldtype(OI_TARGET, :list))

eltype(db::OI_TARGET) = eltype(OI_TARGET)
eltype(::Type{OI_TARGET}) = OITargetEntry

@inline @propagate_inbounds getindex(db::OI_TARGET, i::Integer) = db.list[i]

getindex(A::AbstractVector{OITargetEntry}, name::AbstractString) = begin
    x = get(A, name, undef)
    x === undef && error("no targets named \"", fix_name(name), "\" found")
    return x
end

getindex(db::OI_TARGET, name::AbstractString) = begin
    x = get(db.list, name, undef)
    x === undef && error( "no targets named \"", fix_name(name),
                          "\" found in ", db.extname, "data-block")
    return x
end

firstindex(db::OI_TARGET) = firstindex(db.list)
lastindex(db::OI_TARGET) = lastindex(db.list)
eachindex(db::OI_TARGET) = eachindex(db.list)

function iterate(db::OI_TARGET,
                 state::Tuple{Int,Int} = (firstindex(db),
                                          lastindex(db)))
    i, n = state
    if i > n
        return nothing
    else
        return db.list[i], (i + 1, n)
    end
end

haskey(db::OI_TARGET, i::Integer) = (1 ≤ i ≤ length(db))

haskey(db::OI_TARGET, key::AbstractString) = haskey(db.list, key)

haskey(A::AbstractVector{OITargetEntry}, key::AbstractString) = begin
    for tgt in A
        if is_same(tgt.target, key)
            return true
        end
    end
    return false
end

get(db::OI_TARGET, i::Integer, def) =
    if 1 ≤ i ≤ length(db)
        return @inbounds db[i]
    else
        return def
    end

get(db::OI_TARGET, key::AbstractString, def) = get(db.list, key, def)

get(A::AbstractVector{OITargetEntry}, key::AbstractString, def) = begin
    for tgt in A
        if is_same(tgt.target, key)
            return tgt
        end
    end
    return def
end

# Note: this is the default behavior.
values(db::OI_TARGET) = db

keys(db::OI_TARGET) = Iterators.map(x -> fix_name(x.target), db)
keys(::Type{String}, db::OI_TARGET) = keys(db)
keys(::Type{Int}, db::OI_TARGET) = eachindex(db)
keys(::Type{Integer}, db::OI_TARGET) = eachindex(db)

"""
    OIFITS.get_column([T,] db::OI_TARGET, col)

yields the column `col` of an OI-FITS data-block `db`.  Column is identified by
`col` which is either `sym` or `Val(sym)` where `sym` is the symbolic name of
the corresponding field in `OITargetEntry`.

Optional argument `T` is to specify the element type of the returned array.

""" get_column

@inline get_column(db::OI_TARGET, sym::Symbol) =
    get_column(db, Val(sym))

@inline get_column(T::Type, db::OI_TARGET, sym::Symbol) =
    get_column(T, db, Val(sym))

get_column(db::OI_TARGET, ::Val{sym}) where {sym} =
    get_column(fieldtype(OITargetEntry, sym), db, sym)

function get_column(::Type{T}, db::OI_TARGET, ::Val{sym}) where {T,sym}
    len = length(db)
    col = Vector{T}(undef, len)
    for i in 1:len
        col[i] = getfield(db[i], sym)
    end
    return col
end


text:
#### UpdraftVars

export Cloud, SurfaceBC, UpdraftVar, compute_cloud_base_top_cover!

mutable struct Cloud{FT}
    base::FT
    top::FT
    cover::FT
end

mutable struct SurfaceBC{FT}
    a::FT
    w::FT
    θ_liq::FT
    q_tot::FT
end

struct UpdraftVar{FT}
    cloud::Cloud{FT}
    surface_bc::SurfaceBC{FT}
    surface_scalar_coeff::FT
    function UpdraftVar(i::I, surface_area::FT, n_updrafts::I) where {FT, I}
        a_s = surface_area / FT(n_updrafts)
        c = Cloud{FT}(0, 0, 0)
        s_bc = SurfaceBC{FT}(0, 0, 0, 0)
        surface_scalar_coeff = 1.7072226094205676
        # TODO: FIXME:
        # surface_scalar_coeff = percentile_bounds_mean_norm(1-surface_area + i    *a_s,
        #                                                    1-surface_area + (i+1)*a_s , 1000)
        return new{FT}(c, s_bc, surface_scalar_coeff)
    end
end

function compute_cloud_base_top_cover!(
    UpdVar,
    grid::Grid{FT},
    q::StateVec,
    aux::StateVec,
) where {FT}
    gm, en, ud, sd, al = allcombinations(q)
    k_2 = first_interior(grid, Zmax())
    @inbounds for i in ud
        UpdVar[i].cloud.base = grid.zc[k_2]
        UpdVar[i].cloud.top = FT(0)
        UpdVar[i].cloud.cover = FT(0)
        @inbounds for k in over_elems_real(grid)
            a_ik = q[:a, k, i]
            z_k = grid.zc[k]
            if aux[:q_liq, k, i] > FT(1e-8) && a_ik > FT(1e-3)
                UpdVar[i].cloud.base = min(UpdVar[i].cloud.base, z_k)
                UpdVar[i].cloud.top = max(UpdVar[i].cloud.top, z_k)
                UpdVar[i].cloud.cover = max(UpdVar[i].cloud.cover, a_ik)
            end
        end
    end
end

ext: jl


text:
using Distributions
using Test, LinearAlgebra

v = 7.0
S = Matrix(1.0I, 2, 2)
S[1, 2] = S[2, 1] = 0.5

W = Wishart(v,S)
IW = InverseWishart(v,S)

for d in [W,IW]
    local d
    @test size(d) == size(rand(d))
    @test length(d) == length(rand(d))
    @test typeof(d)(params(d)...) == d
    @test partype(d) == Float64
end

@test partype(Wishart(7, Matrix{Float32}(I, 2, 2))) == Float32
@test partype(InverseWishart(7, Matrix{Float32}(I, 2, 2))) == Float32

@test isapprox(mean(rand(W,100000)) , mean(W) , atol=0.1)
@test isapprox(mean(rand(IW,100000)), mean(IW), atol=0.1)

v = 3.0

@test isapprox(pdf(Wishart(v,S), S)         , 0.04507168, atol=1e-8)
@test isapprox(pdf(Wishart(v,S), inv(S))    , 0.01327698, atol=1e-8)
@test isapprox(pdf(Wishart(v,inv(S)),S)     , 0.0148086 , atol=1e-8)
@test isapprox(pdf(Wishart(v,inv(S)),inv(S)), 0.01901462, atol=1e-8)

@test isapprox(pdf(Wishart(v,S), [S, S])               , [0.04507168, 0.04507168], atol=1e-6)
@test isapprox(pdf(Wishart(v,S), [inv(S), inv(S)])     , [0.01327698, 0.01327698], atol=1e-6)
@test isapprox(pdf(Wishart(v,inv(S)), [S, S])          , [0.0148086, 0.0148086]  , atol=1e-6)
@test isapprox(pdf(Wishart(v,inv(S)), [inv(S), inv(S)]), [0.01901462, 0.01901462], atol=1e-6)

@test logpdf(Wishart(v,S), S)           ≈ log.(pdf(Wishart(v,S), S))
@test logpdf(Wishart(v,S), inv(S))      ≈ log.(pdf(Wishart(v,S), inv(S)))
@test logpdf(Wishart(v,inv(S)), S)      ≈ log.(pdf(Wishart(v,inv(S)), S))
@test logpdf(Wishart(v,inv(S)), inv(S)) ≈ log.(pdf(Wishart(v,inv(S)), inv(S)))

@test logpdf(Wishart(v,S), [S, S])                ≈ log.(pdf(Wishart(v,S), [S, S]))
@test logpdf(Wishart(v,S), [inv(S), inv(S)])      ≈ log.(pdf(Wishart(v,S), [inv(S), inv(S)]))
@test logpdf(Wishart(v,inv(S)), [S, S])           ≈ log.(pdf(Wishart(v,inv(S)), [S, S]))
@test logpdf(Wishart(v,inv(S)), [inv(S), inv(S)]) ≈ log.(pdf(Wishart(v,inv(S)), [inv(S), inv(S)]))

@test isapprox(pdf(InverseWishart(v,S), S)         , 0.04507168 , atol=1e-8)
@test isapprox(pdf(InverseWishart(v,S), inv(S))    , 0.006247377, atol=1e-8)
@test isapprox(pdf(InverseWishart(v,inv(S)),S)     , 0.03147137 , atol=1e-8)
@test isapprox(pdf(InverseWishart(v,inv(S)),inv(S)), 0.01901462 , atol=1e-8)

@test isapprox(pdf(InverseWishart(v,S), [S, S])               , [0.04507168,  0.04507168] , atol=1e-6)
@test isapprox(pdf(InverseWishart(v,S), [inv(S), inv(S)])     , [0.006247377, 0.006247377], atol=1e-6)
@test isapprox(pdf(InverseWishart(v,inv(S)), [S, S])          , [0.03147137,  0.03147137] , atol=1e-6)
@test isapprox(pdf(InverseWishart(v,inv(S)), [inv(S), inv(S)]), [0.01901462,  0.01901462] , atol=1e-6)

@test logpdf(InverseWishart(v,S), S)           ≈ log(pdf(InverseWishart(v,S), S))
@test logpdf(InverseWishart(v,S), inv(S))      ≈ log(pdf(InverseWishart(v,S), inv(S)))
@test logpdf(InverseWishart(v,inv(S)), S)      ≈ log(pdf(InverseWishart(v,inv(S)), S))
@test logpdf(InverseWishart(v,inv(S)), inv(S)) ≈ log(pdf(InverseWishart(v,inv(S)), inv(S)))

@test logpdf(InverseWishart(v,S), [S, S])                ≈ log.(pdf(InverseWishart(v,S), [S, S]))
@test logpdf(InverseWishart(v,S), [inv(S), inv(S)])      ≈ log.(pdf(InverseWishart(v,S), [inv(S), inv(S)]))
@test logpdf(InverseWishart(v,inv(S)), [S, S])           ≈ log.(pdf(InverseWishart(v,inv(S)), [S, S]))
@test logpdf(InverseWishart(v,inv(S)), [inv(S), inv(S)]) ≈ log.(pdf(InverseWishart(v,inv(S)), [inv(S), inv(S)]))


text:
using FMIndexes
using FactCheck

srand(12345)

# simple DNA sequence type
@enum Nuc A C G T
type DNASeq
    data::Vector{Nuc}
end
Base.getindex(seq::DNASeq, i::Integer) = seq.data[i]
Base.length(seq::DNASeq) = length(seq.data)
Base.endof(seq::DNASeq)  = length(seq.data)

facts("construct") do
    context("one") do
        σ = 2
        seq = [0x00]
        index = FMIndex(seq, σ)
        @fact typeof(index) --> FMIndex{1,UInt8}
    end

    context("short") do
        σ = 4
        seq = [0x00, 0x01, 0x02, 0x03]
        index = FMIndex(seq, σ)
        @fact typeof(index) --> FMIndex{2,UInt8}
    end

    context("long") do
        σ = 2
        seq = rand(0x00:0x01, 2^24)
        index = FMIndex(seq, σ)
        @fact typeof(index) --> FMIndex{1,UInt32}
    end

    context("dna") do
        σ = 4
        seq = DNASeq([A, C, G, T, A, T])
        index = FMIndex(seq, σ)
        @fact typeof(index) --> FMIndex{2,UInt8}
    end

    context("boundaries") do
        σ = 4
        # 8 bits
        n = 2^8
        seq = rand(0x00:0x03, n)
        index = FMIndex(seq, σ)
        @fact typeof(index) --> FMIndex{2,UInt8}
        # 16 bits
        n = 2^16
        seq = rand(0x00:0x03, n)
        index = FMIndex(seq, σ)
        @fact typeof(index) --> FMIndex{2,UInt16}
        # 32 bits
        # too long to test!
    end

    context("mmap") do
        σ = 4
        seq = rand(0x00:0x03, 2^10)
        index = FMIndex(seq, σ, mmap=true)
        @fact typeof(index) --> FMIndex{2,UInt16}
    end

    context("use pSAscan") do
        σ = 4
        seq = rand(0x00:0x03, 2^24)
        if haskey(ENV, "PSASCAN")
            index = FMIndex(seq, σ, program=:psascan, psascan=ENV["PSASCAN"], mmap=true)
            @fact typeof(index) --> FMIndex{2,UInt32}
        else
            @pending typeof(index) --> FMIndex{2,UInt32}
        end
    end
end

facts("restore") do
    context("examples") do
        σ = 2
        seq = [0x01]
        index = FMIndex(seq, σ)
        @fact restore(index) --> seq

        σ = 2
        seq = [0x00, 0x00, 0x01, 0x01]
        index = FMIndex(seq, σ)
        @fact restore(index) --> seq

        seq = [0x01, 0x01, 0x00, 0x00]
        index = FMIndex(seq, σ)
        @fact restore(index) --> seq

        σ = 128
        seq = "abracadabra".data
        index = FMIndex(seq, σ)
        @fact restore(index) --> seq
    end

    context("dna") do
        σ = 4
        seq = DNASeq([A, C, G, T, A, T])
        index = FMIndex(seq, σ)
        @fact restore(index) --> [0x00, 0x01, 0x02, 0x03, 0x0, 0x03]
    end

    context("random") do
        σ = 128
        for n in [2, 15, 100], _ in 1:100
            seq = randstring(n)
            index = FMIndex(seq.data, σ)
            @fact restore(index) --> seq.data
        end
    end
end

facts("count") do
    context("[0x00]") do
        σ = 2
        index = FMIndex([0x00], σ)
        @fact count([0x00], index) --> 1
        @fact count([0x01], index) --> 0
    end

    context("\"abracadabra\"") do
        σ = 128
        seq = "abracadabra".data
        index = FMIndex(seq, σ)

        @fact count("a", index) --> 5
        @fact count("ab", index) --> 2
        @fact count("abr", index) --> 2
        @fact count("abra", index) --> 2
        @fact count("abrac", index) --> 1

        @fact count("d", index) --> 1
        @fact count("da", index) --> 1
        @fact count("dab", index) --> 1

        @fact count("r", index) --> 2
        @fact count("ra", index) --> 2
        @fact count("rac", index) --> 1

        # not existing
        @fact count("x", index) --> 0
        @fact count("ax", index) --> 0
        @fact count("aa", index) --> 0
        @fact count("aaa", index) --> 0
        @fact count("braa", index) --> 0
    end

    context("[0x00, 0x00, 0x01, 0x01]") do
        σ = 2
        seq = [0x00, 0x00, 0x01, 0x01]
        index = FMIndex(seq, σ)

        @fact count([0x00], index) --> 2
        @fact count([0x00, 0x00], index) --> 1
        @fact count([0x00, 0x00, 0x01], index) --> 1
        @fact count([0x00, 0x01, 0x01], index) --> 1
        @fact count([0x00, 0x00, 0x01, 0x01], index) --> 1

        @fact count([0x01], index) --> 2
        @fact count([0x01, 0x01], index) --> 1

        @fact count([0x01, 0x00], index) --> 0
        @fact count([0x01, 0x01, 0x00], index) --> 0
        @fact count([0x00, 0x01, 0x01, 0x00], index) --> 0
        @fact count([0x01, 0x01, 0x00, 0x00], index) --> 0
    end
end

function linear_search(query, seq)
    locs = Int[]
    for i in 1:endof(seq)-length(query)+1
        j = 1
        while j ≤ length(query) && query[j] == seq[i+j-1]
            j += 1
        end
        if j > length(query)
            push!(locs, i)
        end
    end
    return locs
end

facts("locate/locateall") do
    context("[0x00]") do
        σ = 2
        index = FMIndex([0x00], σ)
        @fact locateall([0x00], index) --> [1]
        @fact locateall([0x01], index) --> isempty
    end

    context("\"abracadabra\"") do
        σ = 128
        seq = "abracadabra".data
        index = FMIndex(seq, σ)

        @fact locate("a", index) |> collect |> sort --> [1, 4, 6, 8, 11]

        @fact locateall("a", index) |> sort --> [1, 4, 6, 8, 11]
        @fact locateall("ab", index) |> sort --> [1, 8]
        @fact locateall("abr", index) |> sort --> [1, 8]
        @fact locateall("abra", index) |> sort --> [1, 8]
        @fact locateall("abrac", index) |> sort --> [1]

        @fact locateall("d", index) |> sort --> [7]
        @fact locateall("da", index) |> sort --> [7]
        @fact locateall("dab", index) |> sort --> [7]

        @fact locateall("r", index) |> sort --> [3, 10]
        @fact locateall("ra", index) |> sort --> [3, 10]
        @fact locateall("rac", index) |> sort --> [3]

        @fact locate("x", index) --> isempty

        @fact locateall("x", index) --> isempty
        @fact locateall("ax", index) --> isempty
        @fact locateall("aa", index) --> isempty
        @fact locateall("aaa", index) --> isempty
        @fact locateall("braa", index) --> isempty
    end

    context("random") do
        σ = 4
        seq = rand(0x00:0x03, 10_000)
        index = FMIndex(seq, σ)
        for m in [1, 2, 3, 5, 10]
            query = rand(0x00:0x03, m)
            @fact locateall(query, index) |> sort --> linear_search(query, seq)
        end
    end

    context("boundaries") do
        σ = 4
        for n in [2^8, 2^16]
            seq = rand(0x00:0x03, n)
            index = FMIndex(seq, σ)
            for set in ([0x00],
                        [0x03],
                        [0x00, 0x00],
                        [0x00, 0x01],
                        [0x02, 0x03],
                        [0x00, 0x01, 0x02],
                        [0x01, 0x02, 0x03],
                        [0x00, 0x01, 0x02, 0x03])
                for query in permutations(set)
                    @fact locateall(query, index) |> sort --> linear_search(query, seq)
                end
            end
        end
    end
end

facts("full-text search") do
    function linear_search(query)
        locs = Int[]
        loc = 0
        while (loc = searchindex(text, query, loc + 1)) > 0
            push!(locs, loc)
        end
        return locs
    end

    text = open(readall, Pkg.dir("FMIndexes", "test", "lorem_ipsum.txt"))
    index = FMIndex(text, r=2)

    @fact count("Lorem", index) --> 1
    @fact locateall("Lorem", index) --> [1]
    @fact count("hoge", index) --> 0
    @fact locateall("hoge", index) --> isempty
    for query in ["a", ".",
                  "ex", "In",
                  "non", "Sed",
                  "odio", "Cras",
                  "sollicitudin"]
        locs = linear_search(query)
        @fact count(query, index) --> length(locs)
        @fact locateall(query, index) |> sort --> locs
    end

    @fact bytestring(restore(index)) --> text
end


text:
export substellarlatitude, hourangle
export diurnalfluxfactor, diurnalfluxfactors
export annualfluxfactor, annualfluxfactors

"""
    substellarlatitude(f, γ)

Compute the latitude of the substellar point for a given solar longitude `f` (true anomaly) and obliquity `γ`
"""
substellarlatitude(f, γ) = asin(cos(f)*sin(γ))

"""
    hourangle(θ, θₛ)

Compute the [hour angle](https://en.wikipedia.org/wiki/Hour_angle)
"""
function hourangle(θ, θₛ)
    x = -sin(θ)*sin(θₛ)/(cos(θ)*cos(θₛ))
    if x <= -1
        return π
    elseif x >= 1
        return 0.0
    end
    return acos(x)
end

#θ - latitude
#θₛ - substellar latitude
"""
    diurnalfluxfactor(θ, θₛ)

Compute the diurnally averaged fraction of incoming stellar flux received by a point at latitude `θ` when the substellar latitude is `θₛ`
"""
function diurnalfluxfactor(θ, θₛ)
    h = hourangle(θ, θₛ)
    return (sin(h)*cos(θ)*cos(θₛ) + h*sin(θ)*sin(θₛ))/π
end

#θ - latitude
#f - solar longitude
#γ - obliquity
"""
    diurnalfluxfactor(θ, f, γ)

Compute the diurnally averaged fraction of incoming stellar flux received by a point at latitude `θ` when the planet is at solar longitude (true anomaly) `f`, with obliquity `γ`
"""
diurnalfluxfactor(θ, f, γ) = diurnalfluxfactor(θ, substellarlatitude(f, γ))

"""
    diurnalfluxfactor(t, a, m, e, θ, γ, p)

Compute the diurnally averaged fraction of incoming stellar flux received by a point at latitude `θ` for a general elliptical orbit
"""
function diurnalfluxfactor(t, a, m, e, θ, γ, p)
    f = trueanomaly(t, a, m, e)
    r = orbitaldistance(a, f, e)
    return diurnalfluxfactor(θ, f - p, γ)*(a/r)^2
end

"""
    diurnalfluxfactors(γ; nf=251, nθ=181)

Compute a grid of diurnally averaged fractions of incoming stellar flux received by a point at latitude `θ` for a planet with obliquity `γ` in a circular orbit. Returns  a solar longitude vector (column values), latitude vector (row values), and the grid of flux factors. `nf` indicates the number of points around the orbit and `nθ` indicates the number of latitudes.
"""
function diurnalfluxfactors(γ; nf::Int=251, nθ::Int=181)
    θ = LinRange(-π/2, π/2, nθ)
    f = LinRange(0, 2π, nf)
    F, Θ = meshgrid(f, θ)
    return (f, θ, diurnalfluxfactor.(Θ, F, γ))
end

"""
    diurnalfluxfactors(a, m, e, γ, p; nt=251, nθ=181)

Compute a grid of diurnally averaged fractions of incoming stellar flux for a planet in a general elliptical orbit. Returns a time vector (column values) over one orbital period, latitude vector (row values), and the grid of flux factors. `nt` indicates the number of time samples around the orbit and `nθ` indicates the number of latitudes.
"""
function diurnalfluxfactors(a, m, e, γ, p; nt::Int=251, nθ::Int=181)
    t = LinRange(0, orbitalperiod(a, m), nt)
    θ = LinRange(-π/2, π/2, nθ)
    T, Θ = meshgrid(t, θ)
    return (t, θ, diurnalfluxfactor.(T, a, m, e, Θ, γ, p))
end

"""
    annualfluxfactor(e, θ, γ, p)

Compute the annually averaged flux factor for a latitude `θ` on a planet in a general elliptical orbit.
"""
function annualfluxfactor(e, θ, γ, p; tol::Float64=1e-4)
    T = orbitalperiod(1.0, 1.0)
    f(t) = diurnalfluxfactor(t, 1.0, 1.0, e, θ, γ, p)
    F, _ = hquadrature(f, 0, T, reltol=tol, abstol=tol)
    return F/T
end

"""
    annualfluxfactors(e, γ, p; nθ=181)

Compute a range of annually averaged flux factors for a planet in a general elliptical orbit. Returns a latitude vector (row values) and a vector of flux factors. `nθ` indicates the number of latitude samples.
"""
function annualfluxfactors(e, γ, p; nθ::Int=181)
    θ = LinRange(-π/2, π/2, nθ)
    F = annualfluxfactor.(e, θ, γ, p)
    return θ, F
end


text:
module LDPCStorage

include("utils.jl")

include("alist.jl")
export load_alist, save_to_alist

include("cscmat.jl")
export save_to_cscmat, load_cscmat, load_matrix_from_qc_cscmat_file, CSCMAT_FORMAT_VERSION

end # module


text:
module FIFOBuffers

import ..bytesavailable
import ..bytes
using ..IOExtras
import Base.==

export FIFOBuffer

"""
    FIFOBuffer([max::Integer])
    FIFOBuffer(string_or_bytes_vector)
    FIFOBuffer(io::IO)

A `FIFOBuffer` is a first-in, first-out, in-memory, async-friendly IO buffer type.

`FIFOBuffer([max])`: creates a "open" `FIFOBuffer` with a maximum size of `max`; this means that bytes can be written
up until `max` number of bytes have been written (with none being read). At this point, the `FIFOBuffer` is full
and will return 0 for all subsequent writes. If no `max` (`FIFOBuffer()`) argument is given, then a default size of `typemax(Int32)^2` is used;
this essentially allows all writes every time. Note that providing a string or byte vector argument mirrors the behavior of `Base.IOBuffer`
in that the `max` size of the `FIFOBuffer` is the length of the string/byte vector; it is also not writeable.

Reading is supported via `readavailable(f)` and `read(f, nb)`, which returns all or `nb` bytes, respectively, starting at the earliest bytes written.
All read functions will return an empty byte vector, even if the buffer has been closed. Checking `eof` will correctly reflect when the buffer has
been closed and no more bytes will be available for reading.

You may call `String(f::FIFOBuffer)` to view the current contents in the buffer without consuming them.

A `FIFOBuffer` is built to be used asynchronously to allow buffered reading and writing. In particular, a `FIFOBuffer`
detects if it is being read from/written to the main task, or asynchronously, and will behave slightly differently depending on which.

Specifically, when reading from a `FIFOBuffer`, if accessed from the main task, it will not block if there are no bytes available to read, instead returning an empty `UInt8[]`.
If being read from asynchronously, however, reading will block until additional bytes have been written. An example of this in action is:

```julia
f = HTTP.FIFOBuffer(5) # create a FIFOBuffer that will hold at most 5 bytes, currently empty
f2 = HTTP.FIFOBuffer(5) # a 2nd buffer that we'll write to asynchronously

# start an asynchronous writing task with the 2nd buffer
tsk = @async begin
    while !eof(f)
        write(f2, readavailable(f))
    end
end

# now write some bytes to the first buffer
# writing triggers our async task to wake up and read the bytes we just wrote
# leaving the first buffer empty again and blocking again until more bytes have been written
write(f, [0x01, 0x02, 0x03, 0x04, 0x05])

# we can see that `f2` now holds the bytes we wrote to `f`
String(readavailable(f2))

# our async task will continue until `f` is closed
close(f)

istaskdone(tsk) # true
```
"""
mutable struct FIFOBuffer <: IO
    len::Int64 # length of buffer in bytes
    max::Int64 # the max size buffer is allowed to grow to
    nb::Int64  # number of bytes available to read in buffer
    f::Int64   # buffer index that should be read next, unless nb == 0, then buffer is empty
    l::Int64   # buffer index that should be written to next, unless nb == len, then buffer is full
    buffer::Vector{UInt8}
    cond::Condition
    task::Task
    eof::Bool
end

const DEFAULT_MAX = Int64(typemax(Int32))^Int64(2)

FIFOBuffer(f::FIFOBuffer) = f
FIFOBuffer(max) = FIFOBuffer(0, max, 0, 1, 1, UInt8[], Condition(), current_task(), false)
FIFOBuffer() = FIFOBuffer(DEFAULT_MAX)

const EMPTYBODY = FIFOBuffer()

FIFOBuffer(str::String) = FIFOBuffer(Vector{UInt8}(str))
function FIFOBuffer(bytes::Vector{UInt8})
    len = length(bytes)
    return FIFOBuffer(len, len, len, 1, 1, bytes, Condition(), current_task(), true)
end
FIFOBuffer(io::IOStream) = FIFOBuffer(read(io))
FIFOBuffer(io::IO) = FIFOBuffer(readavailable(io))

==(a::FIFOBuffer, b::FIFOBuffer) = String(a) == String(b)
Base.length(f::FIFOBuffer) = f.nb
bytesavailable(f::FIFOBuffer) = f.nb
bytes(f::FIFOBuffer) = readavailable(f)
Base.wait(f::FIFOBuffer) = wait(f.cond)
Base.read(f::FIFOBuffer) = readavailable(f)
Base.flush(f::FIFOBuffer) = nothing
Base.position(f::FIFOBuffer) = f.f, f.l, f.nb
function Base.seek(f::FIFOBuffer, pos::Tuple{Int64, Int64, Int64})
    f.f = pos[1]
    f.l = pos[2]
    f.nb = pos[3]
    return
end

Base.eof(f::FIFOBuffer) = f.eof && f.nb == 0
Base.isopen(f::FIFOBuffer) = !f.eof
function Base.close(f::FIFOBuffer)
    f.eof = true
    notify(f.cond)
    return
end

# 0 | 1 | 2 | 3 | 4 | 5 |
#---|---|---|---|---|---|
#   |f/l| _ | _ | _ | _ | empty, f == l, nb = 0, can't read, can write from l to l-1, don't need to change f, l = l, nb = len
#   | _ | _ |f/l| _ | _ | empty, f == l, nb = 0, can't read, can write from l:end, 1:l-1, don't need to change f, l = l, nb = len
#   | _ | f | x | l | _ | where f < l, can read f:l-1, then set f = l, can write l:end, 1:f-1, then set l = f, nb = len
#   | l | _ | _ | f | x | where l < f, can read f:end, 1:l-1, can write l:f-1, then set l = f
#   |f/l| x | x | x | x | full l == f, nb = len, can read f:l-1, can't write
#   | x | x |f/l| x | x | full l == f, nb = len, can read f:end, 1:l-1, can't write
function Base.readavailable(f::FIFOBuffer)
    # no data to read
    if f.nb == 0
        if current_task() == f.task || f.eof
            return UInt8[]
        else # async + still open: block till there's data to read
            wait(f.cond)
            f.nb == 0 && return UInt8[]
        end
    end
    if f.f < f.l
        @inbounds bytes = f.buffer[f.f:f.l-1]
    else
        # we've wrapped around
        @inbounds bytes = f.buffer[f.f:end]
        @inbounds append!(bytes, view(f.buffer, 1:f.l-1))
    end
    f.f = f.l
    f.nb = 0
    notify(f.cond)
    return bytes
end

# read at most `nb` bytes
function Base.read(f::FIFOBuffer, nb::Int)
    # no data to read
    if f.nb == 0
        if current_task() == f.task || f.eof
            return UInt8[]
        else # async: block till there's data to read
            wait(f.cond)
            f.nb == 0 && return UInt8[]
        end
    end
    if f.f < f.l
        l = (f.l - f.f) <= nb ? (f.l - 1) : (f.f + nb - 1)
        @inbounds bytes = f.buffer[f.f:l]
        f.f = mod1(l + 1, f.max)
    else
        # we've wrapped around
        if nb <= (f.len - f.f + 1)
            # we can read all we need between f.f and f.len
            @inbounds bytes = f.buffer[f.f:(f.f + nb - 1)]
            f.f = mod1(f.f + nb, f.max)
        else
            @inbounds bytes = f.buffer[f.f:f.len]
            l = min(f.l - 1, nb - length(bytes))
            @inbounds append!(bytes, view(f.buffer, 1:l))
            f.f = mod1(l + 1, f.max)
        end
    end
    f.nb -= length(bytes)
    notify(f.cond)
    return bytes
end

function Base.read(f::FIFOBuffer, ::Type{Tuple{UInt8,Bool}})
    # no data to read
    if f.nb == 0
        if current_task() == f.task || f.eof
            return 0x00, false
        else # async: block till there's data to read
            f.eof && return 0x00, false
            wait(f.cond)
            f.nb == 0 && return 0x00, false
        end
    end
    # data to read
    @inbounds b = f.buffer[f.f]
    f.f = mod1(f.f + 1, f.max)
    f.nb -= 1
    notify(f.cond)
    return b, true
end

function Base.read(f::FIFOBuffer, ::Type{UInt8})
    byte, valid = read(f, Tuple{UInt8,Bool})
    valid || throw(EOFError())
    return byte
end

function Base.String(f::FIFOBuffer)
    f.nb == 0 && return ""
    if f.f < f.l
        return String(f.buffer[f.f:f.l-1])
    else
        bytes = f.buffer[f.f:end]
        append!(bytes, view(f.buffer, 1:f.l-1))
        return String(bytes)
    end
end

function Base.write(f::FIFOBuffer, b::UInt8)
    # buffer full, check if we can grow it
    if f.nb == f.len || f.len < f.l
        if f.len < f.max
            push!(f.buffer, 0x00)
            f.len += 1
        else
            if current_task() == f.task || f.eof
                return 0
            else # async: block until there's room to write
                wait(f.cond)
                f.nb == f.len && return 0
            end
        end
    end
    # write our byte
    @inbounds f.buffer[f.l] = b
    f.l = mod1(f.l + 1, f.max)
    f.nb += 1
    notify(f.cond)
    return 1
end

function Base.write(f::FIFOBuffer, bytes::Vector{UInt8}, i, j)
    len = j - i + 1
    if f.nb == f.len || f.len < f.l
        # buffer full, check if we can grow it
        if f.len < f.max
            append!(f.buffer, zeros(UInt8, min(len, f.max - f.len)))
            f.len = length(f.buffer)
        else
            if current_task() == f.task || f.eof
                return 0
            else # async: block until there's room to write
                wait(f.cond)
                f.nb == f.len && return 0
            end
        end
    end
    if f.f <= f.l
        # non-wraparound
        avail = f.len - f.l + 1
        if len > avail
            # need to wrap around, and check if there's enough room to write full bytes
            # write `avail` # of bytes to end of buffer
            unsafe_copy!(f.buffer, f.l, bytes, i, avail)
            if len - avail < f.f
                # there's enough room to write the rest of bytes
                unsafe_copy!(f.buffer, 1, bytes, avail + 1, len - avail)
                f.l = len - avail + 1
            else
                # not able to write all of bytes
                unsafe_copy!(f.buffer, 1, bytes, avail + 1, f.f - 1)
                f.l = f.f
                f.nb += avail + f.f - 1
                notify(f.cond)
                return avail + f.f - 1
            end
        else
            # there's enough room to write bytes through the end of the buffer
            unsafe_copy!(f.buffer, f.l, bytes, i, len)
            f.l = mod1(f.l + len, f.max)
        end
    else
        # already in wrap-around state
        if len > mod1(f.f - f.l, f.max)
            # not able to write all of bytes
            nb = f.f - f.l
            unsafe_copy!(f.buffer, f.l, bytes, i, nb)
            f.l = f.f
            f.nb += nb
            notify(f.cond)
            return nb
        else
            # there's enough room to write bytes
            unsafe_copy!(f.buffer, f.l, bytes, i, len)
            f.l  = mod1(f.l + len, f.max)
        end
    end
    f.nb += len
    notify(f.cond)
    return len
end

if isdefined(Base, :CodeUnits)
Base.write(f::FIFOBuffer, bytes::Base.CodeUnits) = write(f, Vector{UInt8}(String(bytes)))
end
Base.write(f::FIFOBuffer, bytes::Vector{UInt8}) = write(f, bytes, 1, length(bytes))
Base.write(f::FIFOBuffer, str::String) = write(f, IOExtras.bytes(str))

end # module


text:
# Example which demonstrates adding an edge to a Graph

using LightGraphs
using NetworkViz
using ThreeJS
using Escher

g = Graph(10)

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

  inp = Signal(Dict{Any, Any}(:node1=>"3",:node2=>"3"))

  s = Escher.sampler()
  sampler = Escher.sampler()
  form = vbox(
      h1("Add New Edge"),
      watch!(s, :node1, textinput("", label="Source Vertex")),
      watch!(s, :node2, textinput("", label="Destination Vertex")),
      trigger!(s, :submit, button("Add Edge"))
  ) |> maxwidth(400px)

  map(inp) do dict
      node1 = parse(Int,dict[:node1])
      node2 = parse(Int,dict[:node2])
      vbox(
          intent(s, form) >>> inp,
          vskip(2em),
          addEdge(g, node1, node2, 0)
      ) |> Escher.pad(2em)
  end
end


text:
# This function calculates the largest eigenvalue
# (absolute value wise) by power iteration.

function maxeig!(integrator, cache::OrdinaryDiffEqConstantCache)
  isfirst = integrator.iter == 1 || integrator.u_modified
  @unpack t, dt, uprev, u, f, p, fsalfirst = integrator
  maxiter = 50
  safe = 1.2
  # Initial guess for eigenvector `z`
  if isfirst
    fz = fsalfirst
    z = f(fz, p, t)
  else
    z = cache.zprev
  end
  # Perturbation
  u_norm = integrator.opts.internalnorm(uprev)
  z_norm = integrator.opts.internalnorm(z)
  pert   = eps(u_norm)
  sqrt_pert = sqrt(pert)
  is_u_zero = u_norm == zero(u_norm)
  is_z_zero = z_norm == zero(z_norm)
  # Normalize `z` such that z-u lie in a circle
  if ( !is_u_zero && !is_z_zero )
    dz_u = u_norm * sqrt_pert
    quot = dz_u/z_norm
    z = uprev + quot*z
  elseif !is_u_zero
    dz_u = u_norm * sqrt_pert
    z = uprev + uprev*dz_u
  elseif !is_z_zero
    dz_u = pert
    quot = dz_u/z_norm
    z *= quot
  else
    dz_u = pert
    z = dz_u
  end # endif
  # Start power iteration
  integrator.eigen_est = 0
  for iter in 1:maxiter
    fz = f(z, p, t)
    tmp = fz - fsalfirst
    Δ  = integrator.opts.internalnorm(tmp)
    eig_prev = integrator.eigen_est
    integrator.eigen_est = Δ/dz_u * safe
    # Convergence
    if iter >= 2 && abs(eig_prev - integrator.eigen_est) < integrator.eigen_est*0.05
      # Store the eigenvector
      cache.zprev = z
      return true
    end
    # Next `z`
    if Δ != zero(Δ)
      quot = dz_u/Δ
      z = uprev + quot*tmp
    else
      # An arbitrary change on `z`
      cache.z = shuffle(cache.z)
    end
  end
  return false
end

function maxeig!(integrator, cache::OrdinaryDiffEqMutableCache)
  isfirst = integrator.iter == 1 || integrator.u_modified
  @unpack t, dt, uprev, u, f, p, fsalfirst = integrator
  fz, z, atmp = cache.k, cache.tmp, cache.atmp
  ccache = cache.constantcache
  maxiter = 50
  safe = 1.2
  # Initial guess for eigenvector `z`
  if isfirst
    @. fz = u
    f(z, fz, p, t)
  else
    @. z = ccache.zprev
  end
  # Perturbation
  u_norm = integrator.opts.internalnorm(uprev)
  z_norm = integrator.opts.internalnorm(z)
  pert   = eps(u_norm)
  sqrt_pert = sqrt(pert)
  is_u_zero = u_norm == zero(u_norm)
  is_z_zero = z_norm == zero(z_norm)
  # Normalize `z` such that z-u lie in a circle
  if ( !is_u_zero && !is_z_zero )
    dz_u = u_norm * sqrt_pert
    quot = dz_u/z_norm
    @. z = uprev + quot*z
  elseif !is_u_zero
    dz_u = u_norm * sqrt_pert
    @. z = uprev + uprev*dz_u
  elseif !is_z_zero
    dz_u = pert
    quot = dz_u/z_norm
    @. z *= quot
  else
    dz_u = pert
    @. z = dz_u
  end # endif
  # Start power iteration
  integrator.eigen_est = 0
  for iter in 1:maxiter
    f(fz, z, p, t)
    @. atmp = fz - fsalfirst
    Δ  = integrator.opts.internalnorm(atmp)
    eig_prev = integrator.eigen_est
    integrator.eigen_est = Δ/dz_u * safe
    # Convergence
    if iter >= 2 && abs(eig_prev - integrator.eigen_est) < integrator.eigen_est*0.05
      # Store the eigenvector
      @. ccache.zprev = z
      return true
    end
    # Next `z`
    if Δ != zero(Δ)
      quot = dz_u/Δ
      @. z = uprev + quot*atmp
    else
      # An arbitrary change on `z`
      nind = length(uprev)
      ind = 1 + iter % nind
      z[ind] = uprev[ind] - (z[ind] - uprev[ind])
    end
  end
  return false
end

"""
    choosedeg!(cache) -> nothing

Calculate `ms[mdeg]` (the degree of the Chebyshev polynomial)
and `cache.recind` (the index of recurrence parameters for that
degree), where `recf[recind:(recind+ms[mdeg]-2)]` are the `μ,κ` pairs
for the `mdeg` degree method.
"""
function choosedeg!(cache::T) where T
  isconst = T <: OrdinaryDiffEqConstantCache
  isconst || ( cache = cache.constantcache )
  @unpack ms, fp1, fp2, recf, zprev = cache
  recind = 0
  @inbounds for i in 1:46
    recind += ms[i]
    if ms[i] > cache.mdeg
      cache.mdeg = i
      cache.recind = recind
      break
    end # end if
  end # end for
  return nothing
end


text:
import BasisFunctions: DualType

"Trait type to select a dual dictionary suitable for the AZ algorithm."
struct AZDual <: DualType end

dual(::AZDual, dict::Dictionary, measure; options...) =
    azdual(dict, measure; options...)

"""
The dual that is used to create an AZ `Z` matrix.
"""
azdual(ss::DiscreteStyle, ap::ApproximationProblem; options...) =
    dualdictionary(platform(ap), platformparameter(ap), discretemeasure(ap); samplingstyle=ss, options...)
azdual(ss::SamplingStyle, ap::ApproximationProblem; options...) =
    dualdictionary(platform(ap), platformparameter(ap), measure(ap); samplingstyle=ss, options...)

export azdual
@ap_interface azdual
@ap_sampling_property azdual

## The AZ algorithm
export AZ_A
@ap_interface AZ_A
AZ_A(ap::ApproximationProblem; options...) = discretization(ap; options...)

export AZ_Z
@ap_interface AZ_Z
AZ_Z(ap::ApproximationProblem; options...) = dualdiscretization(ap; options...)

export AZ_Zt
@ap_interface AZ_Zt
AZ_Zt(ap::ApproximationProblem; options...) = AZ_Z(ap; options...)'

export plungeoperator
@ap_interface plungeoperator
function plungeoperator(ap::ApproximationProblem; options...)
    A = AZ_A(ap; options...)
    Zt = AZ_Zt(ap; options...)
    I - A*Zt
end

export plungematrix, firstAZstepoperator
@ap_interface plungematrix
function plungematrix(ap::ApproximationProblem; options...)
    A = AZ_A(ap; options...)
    P = plungeoperator(ap; options...)
    P * A
end
const firstAZstepoperator = plungematrix

export plungerank
@ap_interface plungerank
function plungerank(ap::ApproximationProblem;
            rankestimate = 40,
            options...)
    C = plungematrix(ap; options...)
    Q = rSVD_solver(C; rankestimate = rankestimate, options...)
    length(Q.Sinv)
end


text:
using ProfileSVG

"""
    ProfileMetrics

Tracks compute profile traces.
"""
struct ProfileMetrics end
function (::ProfileMetrics)(ev)
    pr = ev.profiler_samples
    samples, lidict = pr.samples, pr.lineinfo
    function demote_frames(samples)
        new_samples = UInt[]
        startidx = 1
        for i in 1:length(samples)
            if samples[i] == 0
                append!(new_samples, samples[startidx:(i-1)]) # FIXME: i-5
                push!(new_samples, 0)
                startidx = i+1
            end
        end
        new_samples
    end
    @static if VERSION >= v"1.8-"
        samples = demote_frames(samples)
    end
    if ev.category != :compute || length(samples) == 0
        return ""
    end
    path, io = mktemp(; cleanup=false)
    ProfileSVG.save(io, samples; lidict)
    close(io)
    path = replace(path, tempdir()*"/" => "")
    ExpiringFile(path)
end

"An automatically-deleting file."
mutable struct ExpiringFile
    path::String
    function ExpiringFile(path)
        ef = new(path)
        finalizer(ef) do ef
            try
                rm(tempdir()*"/"*path)
            catch
                @warn "Failed to remove /tmp/$path"
            end
        end
        ef
    end
end
StructTypes.StructType(::Type{ExpiringFile}) = StructTypes.CustomStruct()
StructTypes.lower(ef::ExpiringFile) = ef.path


text:
# operators
export joPAdistribute, joPAdistributedLinOp, joPAgather
global joPAdistribute       = joSAdistribute
global joPAdistributedLinOp = joSAdistributedLinOp
global joPAgather           = joSAgather
# utilities
export palloc, pzeros, pones, prand, prandn
global palloc   = salloc
global pzeros   = szeros
global pones    = sones
global prand    = srand
global prandn   = srandn
global PAmode = :SA

export jo_PAmode
"""
jo_PAmode - set and query parallel mode

jo_PAmode(verbose=false) - query parallel mode

jo_PAmode([ :SA | :DA ]; verbose=false) - set parallel mode to respectively shared or distributed arrays

"""
function jo_PAmode(;verbose::Bool=false)
    if PAmode==:SA
        verbose && println("Parallel mode is set to :$PAmode (shared arrays)")
    elseif PAmode==:DA
        verbose && println("Parallel mode is set to :$PAmode (distributed arrays)")
    else
        error("Unknown parallel mode $m")
    end
    return PAmode
end
function jo_PAmode(m::Symbol;verbose::Bool=false)
    if m==:SA
        global joPAdistribute       = joSAdistribute
        global joPAdistributedLinOp = joSAdistributedLinOp
        global joPAgather           = joSAgather
        global palloc   = salloc
        global pzeros   = szeros
        global pones    = sones
        global prand    = srand
        global prandn   = srandn
        global PAmode = :SA
        verbose && println("Parallel mode is set to :$PAmode (shared arrays)")
    elseif m==:DA
        global joPAdistribute       = joDAdistribute
        global joPAdistributedLinOp = joDAdistributedLinOp
        global joPAgather           = joDAgather
        global palloc   = dalloc
        global pzeros   = dzeros
        global pones    = dones
        global prand    = drand
        global prandn   = drandn
        global PAmode = :DA
        verbose && println("Parallel mode is set to :$PAmode (distributed arrays)")
    else
        error("Unknown parallel mode $m")
    end
    return nothing
end


text:
# Ref: http://blog.otoro.net/2016/03/25/generating-abstract-patterns-with-tensorflow/
using Images
using Flux
using Parameters: @with_kw

# set parameters
@with_kw mutable struct Args
    z_dim::Int = 2		# Dim of Latent Vector
    x_dim::Int = 512	# X-Dimension of Image
    y_dim::Int = 512	# Y-Dimension of Image
    N::Int = 14			#
    hidden::Int = 9		# Number of hidden layers in the image
    batch_size::Int = 1024	# Batch Size for prediction
end

# cast 0:x-1 to -0.5:0.5
cast(x) = [range(-0.5, stop=0.5, step=1/(x - 1))...]

function getdata(args)
    xs, ys = cast(args.x_dim), cast(args.y_dim)
    xs = repeat(xs, inner=(args.y_dim))
    ys = repeat(ys, outer=(args.x_dim))
	# Radius term for each point of input
    rs = sqrt.(xs.^2 + ys.^2)
    return xs,ys,rs
end

#Definition for each individual layer
# sample weigths from a gaussian distribution
unit(args, in=args.N, out=args.N, f=tanh) = Dense(in, out, f, initW=randn)

function build_model(args)
    # input -> [x, y, r, z...]
    layers = Any[unit(args, 3 + args.z_dim)]
    for i=1:args.hidden
        push!(layers, unit(args))
    end
    push!(layers, unit(args,args.N, 1, σ))
    model = Chain(layers...)
    return model
end

function batch(arr, s)
    batches = []
    l = size(arr, 2)
    for i=1:s:l
        push!(batches, arr[:, i:min(i+s-1, l)])
    end
    batches
end

function get_image(z, model, args)
    n = args.x_dim * args.y_dim
    z = repeat(reshape(z, 1, args.z_dim), outer=(n, 1))
    xs, ys, rs = getdata(args)
    coords = hcat(xs, ys, rs, z)'

    coords = batch(coords, args.batch_size)

    # Pixel value at a position x is defined as output of model at that point
    getColorAt(x) = model(x)

    pixels = [Gray.(hcat(getColorAt.(coords)...))...]
    reshape(pixels, args.y_dim, args.x_dim)
end

function save_img(z, model, args, image_path=joinpath(dirname(@__FILE__),"sample.png"))
    imgg = get_image(z, model, args)
    save(image_path, imgg)
    imgg
end

function generate_img(; kws...)
    args = Args(; kws...)

    model = build_model(args)

    #Saving image as "sample.png"
    save_img(rand(args.z_dim), model, args)
end

cd(@__DIR__)
generate_img()

max_stars_repo_name: thomascherickal/model-zoo

max_issues_repo_name: thomascherickal/model-zoo

max_forks_repo_name: thomascherickal/model-zoo


text:
# Tests for inters.jl

test_pdb_filepath = testfile("4AKE_H.pdb")
test_dssp_filepath = testfile("4AKE.dssp")

@testset "Interactions" begin
    atom_one = Atom("CA", "ALA", 'A', 20, [1.0, 0.0, 0.0], "C")
    atom_two = Atom("CA", "ALA", 'A', 21, [0.0, 1.0, 0.0], "C")
    atom_thr = Atom("CA", "ALA", 'A', 22, [0.0, 0.0, 1.0], "C")
    atoms = Atom[atom_one, atom_two, atom_thr]
    dists_test = distances(atoms)
    dists_real = [
        0.0   sqrt(2)   sqrt(2);
        sqrt(2)   0.0   sqrt(2);
        sqrt(2)   sqrt(2)   0.0;
    ]
    @test isapprox(dists_test, dists_real)

    bonded_pairs = BondedPair[
        BondedPair("ALA", "N", "CA"),
        BondedPair("ALA", "C", "O"),
        BondedPair("ALA", "C", "N")
    ]
    atom_one = Atom("N", "ALA", 'A', 20, [0.0, 0.0, 0.0], "N")
    atom_two = Atom("CA", "ALA", 'A', 20, [0.0, 0.0, 0.0], "C")
    atom_thr = Atom("O", "ALA", 'A', 20, [0.0, 0.0, 0.0], "O")
    atoms = Atom[atom_one, atom_two, atom_thr]
    bonds_test = bonds(atoms, bonded_pairs)
    bonds_real = [
        false true false;
        true false false;
        false false false;
    ]
    @test bonds_test == bonds_real
    # Test if the amide bond is found correctly
    atom_one = Atom("C", "ALA", 'A', 20, [0.0, 0.0, 0.0], "C")
    atom_two = Atom("N", "ALA", 'A', 21, [0.0, 0.0, 0.0], "N")
    atom_thr = Atom("N", "ALA", 'A', 22, [0.0, 0.0, 0.0], "N")
    atoms = Atom[atom_one, atom_two, atom_thr]
    bonds_test = bonds(atoms, bonded_pairs)
    bonds_real = [
        false true false;
        true false false;
        false false false;
    ]
    @test bonds_test == bonds_real

    bonded_pairs = BondedPair[
        BondedPair("ALA", "N", "CA"),
        BondedPair("ALA", "CA", "C"),
    ]
    atom_one = Atom("N", "ALA", 'A', 20, [0.0, 0.0, 0.0], "N")
    atom_two = Atom("CA", "ALA", 'A', 20, [0.0, 0.0, 0.0], "C")
    atom_thr = Atom("C", "ALA", 'A', 20, [0.0, 0.0, 0.0], "C")
    atoms = Atom[atom_one, atom_two, atom_thr]
    angles_test = angles(atoms, bonds(atoms, bonded_pairs))
    angles_real = [
        false false true;
        false false false;
        true false false;
    ]
    @test angles_test == angles_real
    # Test angles are detected across residues
    atom_one = Atom("CA", "ALA", 'A', 20, [0.0, 0.0, 0.0], "C")
    atom_two = Atom("C", "ALA", 'A', 20, [0.0, 0.0, 0.0], "C")
    atom_thr = Atom("N", "ALA", 'A', 21, [0.0, 0.0, 0.0], "N")
    atoms = Atom[atom_one, atom_two, atom_thr]
    angles_test = angles(atoms, bonds(atoms, bonded_pairs))
    angles_real = [
        false false true;
        false false false;
        true false false;
    ]
    @test angles_test == angles_real

    bonded_pairs = BondedPair[
        BondedPair("ALA", "N", "CA"),
        BondedPair("ALA", "CA", "C"),
    ]
    atom_one = Atom("N", "ALA", 'A', 20, [0.0, 0.0, 0.0], "N")
    atom_two = Atom("CA", "ALA", 'A', 20, [0.0, 0.0, 0.0], "C")
    atom_three = Atom("C", "ALA", 'A', 20, [0.0, 0.0, 0.0], "C")
    atom_four = Atom("N", "ALA", 'A', 21, [0.0, 0.0, 0.0], "N")
    atoms = [atom_one, atom_two, atom_three, atom_four]
    bs = bonds(atoms, bonded_pairs)
    divalents_test = divalents(atoms, bs, angles(atoms, bs))
    divalents_real = [
        false false false true;
        false false false false;
        false false false false;
        true false false false;
    ]
    @test divalents_test == divalents_real

    # It is assumed that atoms are different
    # Test if atoms in different residues return false
    atom_one = Atom("CD2", "HIS", 'A', 20, [0.0, 0.0, 0.0], "C")
    atom_two = Atom("CG", "HIS", 'A', 30, [0.0, 0.0, 0.0], "C")
    @test isring(atom_one, atom_two) == false
    # Test if atoms on same residue but not both in ring return false
    atom_one = Atom("CD2", "HIS", 'A', 20, [0.0, 0.0, 0.0], "C")
    atom_two = Atom("CA", "HIS", 'A', 20, [0.0, 0.0, 0.0], "C")
    @test isring(atom_one, atom_two) == false
    # Test if atoms on same residue and both in ring return true
    atom_one = Atom("CD2", "HIS", 'A', 20, [0.0, 0.0, 0.0], "C")
    atom_two = Atom("CG", "HIS", 'A', 20, [0.0, 0.0, 0.0], "C")
    @test isring(atom_one, atom_two) == true

    # It is assumed that the atoms are a 1-4 pair
    # Test if atoms in a double bond as required return true
    atom_one = Atom("CB", "ASN", 'A', 20, [0.0, 0.0, 0.0], "C")
    atom_two = Atom("1HD1", "ASN", 'A', 20, [0.0, 0.0, 0.0], "H")
    @test isdbonefour(atom_one, atom_two) == true
    # Test if atoms not in a double bond as required return false
    atom_one = Atom("CG", "ASN", 'A', 20, [0.0, 0.0, 0.0], "C")
    atom_two = Atom("C", "ASN", 'A', 20, [0.0, 0.0, 0.0], "C")
    @test isdbonefour(atom_one, atom_two) == false

    # Assumes atoms are a 1-4 pair
    # Test that omega 1-4 bonds return true
    atom_one = Atom("O", "ALA", 'A', 20, [0.0, 0.0, 0.0], "O")
    atom_two = Atom("CA", "ALA", 'A', 21, [0.0, 0.0, 0.0], "CA")
    @test isomonefour(atom_one, atom_two) == true
    # Test that non omega 1-4 bonds return false
    atom_one = Atom("N", "ALA", 'A', 20, [0.0, 0.0, 0.0], "N")
    atom_two = Atom("O", "ALA", 'A', 21, [0.0, 0.0, 0.0], "O")
    @test isomonefour(atom_one, atom_two) == false

    # Assumes atoms are a 1-4 pair
    # Test that phi/psi bonds return true
    atom_one = Atom("N", "ALA", 'A', 20, [0.0, 0.0, 0.0], "N")
    atom_two = Atom("O", "ALA", 'A', 21, [0.0, 0.0, 0.0], "O")
    @test isphipsi(atom_one, atom_two) == true
    # Test that non phi/psi bonds return false
    atom_one = Atom("O", "ALA", 'A', 20, [0.0, 0.0, 0.0], "O")
    atom_two = Atom("CA", "ALA", 'A', 21, [0.0, 0.0, 0.0], "CA")
    @test isphipsi(atom_one, atom_two) == false

    # Assumes atoms are a 1-4 pair around the phi/psi angle
    # Test that atoms in the same helix return true
    atom_one = Atom("N", "ALA", 'A', 20, [0.0, 0.0, 0.0], "N")
    atom_two = Atom("O", "ALA", 'A', 21, [0.0, 0.0, 0.0], "O")
    dssps = Dict{String, Char}("20A"=> 'H', "21A"=> 'H')
    @test istightphipsi(atom_one, atom_two, dssps) == true
    # Test that atoms where one residue is proline return true
    atom_one = Atom("N", "PRO", 'A', 20, [0.0, 0.0, 0.0], "N")
    atom_two = Atom("O", "ALA", 'A', 21, [0.0, 0.0, 0.0], "O")
    dssps = Dict{String, Char}("20A"=> ' ', "21A"=> ' ')
    @test istightphipsi(atom_one, atom_two, dssps) == true
    # Test that atoms not in a tight phi/psi bond return false
    atom_one = Atom("N", "ALA", 'A', 20, [0.0, 0.0, 0.0], "N")
    atom_two = Atom("O", "ALA", 'A', 21, [0.0, 0.0, 0.0], "O")
    dssps = Dict{String, Char}("20A"=> ' ', "21A"=> ' ')
    @test istightphipsi(atom_one, atom_two, dssps) == false

    # Assumes atoms are a 1-4 pair around the phi/psi angle
    # Test that atoms in a loop region return true
    atom_one = Atom("N", "ALA", 'A', 20, [0.0, 0.0, 0.0], "N")
    atom_two = Atom("O", "ALA", 'A', 21, [0.0, 0.0, 0.0], "O")
    dssps = Dict{String, Char}("20A"=> ' ', "21A"=> ' ')
    @test isloosephipsi(atom_one, atom_two, dssps) == true
    # Test that atoms where one residue is glycine return true
    atom_one = Atom("N", "GLY", 'A', 20, [0.0, 0.0, 0.0], "N")
    atom_two = Atom("O", "ALA", 'A', 21, [0.0, 0.0, 0.0], "O")
    dssps = Dict{String, Char}("20A"=> 'H', "21A"=> 'H')
    @test isloosephipsi(atom_one, atom_two, dssps) == true
    # Test that atoms not in a loose phi/psi bond return false
    atom_one = Atom("N", "ALA", 'A', 20, [0.0, 0.0, 0.0], "N")
    atom_two = Atom("O", "ALA", 'A', 21, [0.0, 0.0, 0.0], "O")
    dssps = Dict{String, Char}("20A"=> 'H', "21A"=> 'H')
    @test isloosephipsi(atom_one, atom_two, dssps) == false

    # Test two backbone atoms in the same helix return true
    atoms = readpdb(test_pdb_filepath)
    dssps = readdssp(test_dssp_filepath, atoms)
    @test issecstr(atoms[145], atoms[178], atoms, dssps) == true
    # Test two backbone atoms in the same helix but more than 4 residues apart return false
    @test issecstr(atoms[145], atoms[210], atoms, dssps) == false
    # Test non-backbone atoms return false
    @test issecstr(atoms[145], atoms[181], atoms, dssps) == false
    # Test two backbone atoms not in the same secondary structural element return false
    @test issecstr(atoms[145], atoms[127], atoms, dssps) == false

    # Test if close groups of the same charge return false
    atom_one = Atom("CZ", "ARG", 'A', 20, [0.0, 0.0, 0.0], "C")
    atom_two = Atom("NZ", "LYS", 'A', 30, [3.0, 0.0, 0.0], "N")
    atoms = [atom_one, atom_two]
    @test issaltbridge(atom_one, atom_two, atoms, distances(atoms)) == false
    # Test if close groups of opposite charge return true
    atom_one = Atom("CZ", "ARG", 'A', 20, [0.0, 0.0, 0.0], "C")
    atom_two = Atom("CG", "ASP", 'A', 30, [3.0, 0.0, 0.0], "C")
    atoms = [atom_one, atom_two]
    @test issaltbridge(atom_one, atom_two, atoms, distances(atoms)) == true
    # Test if groups of opposite charge with other close atoms return true
    atom_one = Atom("CZ", "ARG", 'A', 20, [0.0, 0.0, 0.0], "C")
    atom_two = Atom("CG", "ASP", 'A', 30, [5.0, 0.0, 0.0], "C")
    atom_three = Atom("NH1", "ARG", 'A', 20, [1.0, 0.0, 0.0], "N")
    atom_four = Atom("OD1", "ASP", 'A', 30, [4.0, 0.0, 0.0], "O")
    atoms = [atom_one, atom_two, atom_three, atom_four]
    @test issaltbridge(atom_one, atom_two, atoms, distances(atoms)) == true

    # Test donor and acceptor atoms far apart return false
    bonded_pairs = BondedPair[BondedPair("ALA", "N", "H")]
    atom_one = Atom("N", "ALA", 'A', 20, [0.0, 0.0, 0.0], "N")
    atom_two = Atom("H", "ALA", 'A', 20, [1.0, 0.0, 0.0], "H")
    atom_three = Atom("O", "ALA", 'A', 30, [4.0, 0.0, 0.0], "O")
    atoms = [atom_one, atom_two, atom_three]
    @test ishbond(1, 3, atoms, distances(atoms), bonds(atoms, bonded_pairs)) == false
    # Test donor and acceptor atoms return true if H is at correct angle
    atom_one = Atom("N", "ALA", 'A', 20, [0.0, 0.0, 0.0], "N")
    atom_two = Atom("H", "ALA", 'A', 20, [1.0, 0.0, 0.0], "H")
    atom_three = Atom("O", "ALA", 'A', 30, [3.0, 0.0, 0.0], "O")
    atoms = [atom_one, atom_two, atom_three]
    @test ishbond(1, 3, atoms, distances(atoms), bonds(atoms, bonded_pairs)) == true
    # Test donor and acceptor atoms return false if H is at wrong angle
    atom_one = Atom("N", "ALA", 'A', 20, [0.0, 0.0, 0.0], "N")
    atom_two = Atom("H", "ALA", 'A', 20, [0.0, 1.0, 0.0], "H")
    atom_three = Atom("O", "ALA", 'A', 30, [3.0, 0.0, 0.0], "O")
    atoms = [atom_one, atom_two, atom_three]
    @test ishbond(1, 3, atoms, distances(atoms), bonds(atoms, bonded_pairs)) == false
    # Test acceptor and H atoms return true if they are close enough
    atom_one = Atom("N", "ALA", 'A', 20, [0.0, 0.0, 0.0], "N")
    atom_two = Atom("H", "ALA", 'A', 20, [1.0, 0.0, 0.0], "H")
    atom_three = Atom("O", "ALA", 'A', 30, [3.0, 0.0, 0.0], "O")
    atoms = [atom_one, atom_two, atom_three]
    @test ishbond(2, 3, atoms, distances(atoms), bonds(atoms, bonded_pairs)) == true

    # Test if hydrophobic atoms within distance return true
    atom_one = Atom("CA", "ALA", 'A', 20, [0.0, 0.0, 0.0], "C")
    atom_two = Atom("H", "ALA", 'A', 30, [3.5, 0.0, 0.0], "H")
    @test ishydrophobic(atom_one, atom_two, 3.5, 1.0) == true
    # Test if hydrophobic atoms out of distance return false
    atom_one = Atom("CA", "ALA", 'A', 20, [0.0, 0.0, 0.0], "C")
    atom_two = Atom("H", "ALA", 'A', 30, [4.5, 0.0, 0.0], "H")
    @test ishydrophobic(atom_one, atom_two, 4.5, 1.0) == false
    # Test if non-hydrophobic atoms in distance return false
    atom_one = Atom("CA", "ALA", 'A', 20, [0.0, 0.0, 0.0], "C")
    atom_two = Atom("O", "ALA", 'A', 30, [3.5, 0.0, 0.0], "O")
    @test ishydrophobic(atom_one, atom_two, 3.5, 1.0) == false

    atom_one = Atom("O", "ALA", 'A', 20, [0.0, 0.0, 0.0], "O")
    atom_two = Atom("C", "ALA", 'A', 20, [1.2, 0.0, 0.0], "C")
    atom_thr = Atom("H", "ALA", 'A', 30, [-1.0, -1.0, -1.0], "H")
    atoms = [atom_one, atom_two, atom_thr]
    dists = distances(atoms)
    inters = [
        0 1 12;
        1 0 15;
        12 15 0;
    ]
    constraints = Constraints(atoms, dists, inters, other_ratio=0.0, tolerance_weight=1.0)
    @test constraints.pres_inds == [2 1; 3 1]
    inter_types, tolerances = interactioninfo(1.0)
    lower_real = [
        dists[2, 1] - tolerances[1],
        dists[3, 1] - tolerances[12],
    ]
    upper_real = [
        dists[2, 1] + tolerances[1],
        dists[3, 1] + tolerances[12],
    ]
    @test isapprox(constraints.lower, lower_real)
    @test isapprox(constraints.upper, upper_real)

    atom_a_one = Atom("O", "ALA", 'A', 20, [0.0, 0.0, 0.0], "O")
    atom_a_two = Atom("C", "ALA", 'A', 20, [1.2, 0.0, 0.0], "C")
    atom_a_thr = Atom("H", "ALA", 'A', 30, [-1.0, -1.0, -1.0], "H")
    atoms_a = [atom_a_one, atom_a_two, atom_a_thr]
    atom_b_one = Atom("O", "ALA", 'A', 20, [0.0, 0.0, 0.0], "O")
    atom_b_two = Atom("C", "ALA", 'A', 20, [1.25, 0.0, 0.0], "C")
    atom_b_thr = Atom("H", "ALA", 'A', 30, [-1.1, -1.1, -1.1], "H")
    atoms_b = [atom_b_one, atom_b_two, atom_b_thr]
    dists_a = distances(atoms_a)
    dists_b = distances(atoms_b)
    inters_a = [
        0 1 12;
        1 0 15;
        12 15 0;
    ]
    inters_b = [
        0 1 10;
        1 0 15;
        10 15 0;
    ]
    constraints = Constraints(atoms_a, dists_a, dists_b, inters_a, inters_b, other_ratio=0.0, tolerance_weight=1.0)
    @test constraints.pres_inds == [2 1; 3 1]
    inter_types, tolerances = interactioninfo(1.0)
    lower_real = [
        min(dists_a[2, 1] - tolerances[1], dists_b[2, 1] - tolerances[1]),
        min(dists_a[3, 1] - tolerances[12], dists_b[3, 1] - tolerances[10]),
    ]
    upper_real = [
        max(dists_a[2, 1] + tolerances[1], dists_b[2, 1] + tolerances[1]),
        max(dists_a[3, 1] + tolerances[12], dists_b[3, 1] + tolerances[10]),
    ]
    @test isapprox(constraints.lower, lower_real)
    @test isapprox(constraints.upper, upper_real)
end


text:
#=
test_control:
- Julia version: 1.5.0
- Author: shisa
- Date: 2020-08-08
=#

module TestControl
    # packages
    using Test

    # external modules
    include("if.jl")
    include("for_while.jl")
    include("try_catch_finally.jl")

    # methods
    function test()
        @testset "Control" begin
            @testset "If" begin
                @test_nowarn If.main()
                @test_throws ErrorException If.check_negative_value(40)
                @test_throws ErrorException If.check_non_negative_value(-40)
            end
            @testset "ForWhile" begin
                @test_nowarn ForWhile.main()
            end
            @testset "TryCatchFinally" begin
                @test_nowarn TryCatchFinally.main()
            end
        end
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    using.TestControl
    TestControl.test()
end

max_stars_repo_name: ShisatoYano/JuliaPractice

max_issues_repo_name: ShisatoYano/JuliaPractice

max_forks_repo_name: ShisatoYano/JuliaPractice


text:
using Dates
using DelimitedFiles
using Statistics
using TimerOutputs

const to = TimerOutput();

# Logging file
io_log = open(
    string(
        ".\\log\\perform_results_",
        Dates.format(now(), "yyyy-mm-dd_HH-MM-SS"),
        ".txt",
    ),
    "w",
)

function randmsq_timed(rows::Int64, cols::Int64)
 @timeit to "randmsq" begin
    x = @timeit to "rand" rand(rows, cols)
    y = @timeit to "mean" mean(x.^2, dims=1)
    return y
    end
end

for i in 0:9
    n_rows = 5_000 + 5_000*i
    n_cols = 5_000 + 5_000*i
    size = n_rows*n_cols

    write(io_log, "Size:\t$size, Rows:\t$n_rows, Columns:\t$n_cols\n")
    randmsq_timed(n_rows, n_cols);
    print_timer(to)
    write(io_log, "$to\n\n\n")

end

close(io_log);

max_stars_repo_name: power-systems-op/grace-arpae

max_issues_repo_name: power-systems-op/grace-arpae

max_forks_repo_name: power-systems-op/grace-arpae


text:
using FromDigits
using Test

# we use `===` instead of `==` to check also the type

@testset "FromDigits.jl" begin
    @testset "empty" begin
        @test fromdigits(Int[]) === 0
        @test fromdigits(Int32[]) === Int32(0)
        @test fromdigits(Int64, Int32[]) === Int64(0)
    end
    @testset "valid" begin
        @test fromdigits([0]) === 0
        @test fromdigits([0,0]) === 0
        @test fromdigits([1,2,3]) === 321
        @test fromdigits(0:9) === 9876543210
        @test fromdigits(UInt16, [1,2,3]) === 0x0141
        @test fromdigits(fill(1,19)) === 1111111111111111111
        @test fromdigits([false, true, false, true, false, true]; base=2) === 42
    end
    @testset "invalid" begin
        @testset "invalid digit" begin
            @test_throws ArgumentError fromdigits(1:19)
        end
        @testset "overflow" begin
            @test_throws OverflowError fromdigits(fill(1,20))
        end
    end
end