text:
# Autogenerated wrapper script for LLVMLibUnwind_jll for x86_64-apple-darwin
export libunwind

JLLWrappers.@generate_wrapper_header("LLVMLibUnwind")
JLLWrappers.@declare_library_product(libunwind, "@rpath/libunwind.1.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libunwind,
        "lib/libunwind.1.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using PyCall
using Distributions

@pyimport matplotlib.pyplot as plt
@pyimport matplotlib.figure as pyfigure
@pyimport matplotlib.pylab as pylab
@pyimport matplotlib.colors as colors

function main(keepinput=true)

  #
  # 1 means white, 0 means black
  #

  # number of lines in filter
  n = 500
  # number of columns in filter
  m = 500
  # radius of the filter circle
  r = 60
  # filter off if nothing, keep input
  x = keepinput ? nothing : 0
  # filter on
  y = 1

  filterfunc = f_cria_func_transf(n,m,r,x,y)

  # 1 is white, 0 is black
  inputmatrix = rand(n,m)

  filterfunc(inputmatrix)

end

function f_cria_func_transf(n,m,r,x,y)

  function applyfilter(inputmatrix)

    # input matrix must be at least as large as the requested filter
    @assert size(inputmatrix,1) >= n
    @assert size(inputmatrix,2) >= m
    # only square matrices
    @assert m == n

    totalsize = n*m

    outputmatrix = reshape(deepcopy(inputmatrix),1,totalsize)

    # add indices to elements so that we know where we are
    outputmatrix = collect(enumerate(outputmatrix))

    # apply the filter per se
    # each (index,element) tuple becomes a single filtered element

    filtered = map( tuple -> mapper(tuple[1],tuple[2],x,y,r,m,n),outputmatrix)

    # turn it back into a matrix
    filtered = reshape(filtered,n,m)

    # display
    plt.matshow(filtered,cmap = plt.get_cmap("hot") )
    plt.show()

  end
end

# returns a possibly filtered element
function mapper(index,elem,off,on,radius,numrows,numcols)

  if withinradius(numrows,numcols,radius,index)
    on
  else
    if off == nothing
      elem
    else
      off
    end
  end

end

function withinradius(numrows,numcols,radius,elemindex)

  # otherwise we can't find the center of the filter circle
  @assert numrows == numcols

  center = if numrows % 2 == 0
    (numrows/2,numcols/2)
  else
    ((numrows+1)/2,(numcols+1)/2)
  end

  # this is the position of the element we want to evaluate
  element_i = div(elemindex,numcols) +1
  element_j = elemindex % numcols +1

  # do you remember pythagoras' theorem?
  dist1 = abs(center[1]-element_i)
  dist2 = abs(center[2]-element_j)
  hyp = sqrt( (dist1^2) + (dist2^2) )

  if hyp <= radius
    true
  else
    false
  end

end

# main(true) sets the off state of the filter to be the original input
# main(false) sets the off state of the filter be 0 (i.e. black)
main(true)

ext: jl


text:
#
# utils.jl --
#
# Implement utility methods for dealing with FITS files and for converting
# values.
#
#------------------------------------------------------------------------------

"""
   OIFITS.get_hdu_type(arg)

yields the FITS Header Data Unit (HDU) type of `arg` as one of the following
symbols: `:binary_table` for a binary table, `:image_hdu` for an image,
`:ascii_table` for an ASCII table, or `:unknown` otherwise.  Argument `arg` can
be a FITS header, a FITS HDU or the value of the XTENSION keyword.

The returned symbol should match the result of the low level method
`FITSIO.Libcfitsio.fits_get_hdu_type`.

"""
get_hdu_type(xtension::AbstractString) =
    (xtension == "BINTABLE"   ? :binary_table :
     xtension == "IMAGE"      ? :image_hdu    :
     xtension == "TABLE"      ? :ascii_table  : :unknown)
get_hdu_type(::Union{T,Type{T}}) where {T<:TableHDU} = :binary_table
get_hdu_type(::Union{T,Type{T}}) where {T<:ImageHDU} = :image_hdu
get_hdu_type(::Union{T,Type{T}}) where {T<:ASCIITableHDU} = :ascii_table
get_hdu_type(::Any) = :unknown
function get_hdu_type(hdr::FITSHeader)
    i = get_key_index(hdr, "XTENSION", 0)
    if i > 0
        let val = hdr[i]
            if isa(val, AbstractString)
                return get_hdu_type(fix_name(val))
            end
        end
    end
    i = get_key_index(hdr, "SIMPLE", 0)
    if i > 0 && hdr[i] == true
        return :image_hdu
    else
        return :unknown
    end
end

"""
   OIFITS.get_key_index(hdr, key[, def])

yields the index of the keyword `key` in FITS header `hdr`.  If `key` is not
found in `hdr`, the default value `def` is returned if it is specified,
otherwise a `KeyError` exception is thrown.

"""
get_key_index(hdr::FITSHeader, key::String) = getindex(hdr.map, key)
get_key_index(hdr::FITSHeader, key::String, def) = get(hdr.map, key, def)
# FIXME: This method is a hack, it should be part of FITSIO.

"""
    OIFITS.get_value(hdr, key[, def])

yields the value of keyword `key` in FITS header `hdr`.  If `key` is not found
in `hdr`, the default value `def` is returned if it is specified, otherwise an
error is thrown.

"""
get_value(hdr::FITSHeader, key::String) =
    ((i = get_key_index(hdr, key, 0)) > 0 ? hdr[i] : missing_fits_keyword(key))
get_value(hdr::FITSHeader, key::String, def) =
    ((i = get_key_index(hdr, key, 0)) > 0 ? hdr[i] : def)

@noinline missing_fits_keyword(key) =
    error("missing FITS keyword \"", key, "\"")

"""
    OIFITS.get_comment(hdr, key[, def])

yields the comment of keyword `key` in FITS header `hdr`.  If `key` is not
found in `hdr`, the default value `def` is returned if it is specified,
otherwise an error is thrown.

"""
get_comment(hdr::FITSHeader, key::String) =
    ((i = get_key_index(hdr, key, 0)) > 0 ? FITSIO.get_comment(hdr, i) :
     missing_fits_keyword(key))
get_comment(hdr::FITSHeader, key::String, def) =
    ((i = get_key_index(hdr, key, 0)) > 0 ? FITSIO.get_comment(hdr, i) : def)

"""
    OIFITS.get_integer(hdr, key[, def])

yields the value of keyword `key` in FITS header `hdr` as an integer.  If `key`
is not found in `hdr`, the default value `def` is returned if it is specified,
otherwise an error is thrown.  If `key` is found, it is checked that it can be
converted into an `Int` which is returned.  The possible types of the result
are thus `Int` (if `key` is found) or `typeof(def)` (if `key` not found and
`def` specified).

"""
function get_integer(hdr::FITSHeader, key::String, def)
    i = get_key_index(hdr, key, 0)
    i == 0 && return def
    val = hdr[i]
    isa(val, Integer) || noninteger_fits_keyword(key)
    to_integer(val)
end

function get_integer(hdr::FITSHeader, key::String)
    i = get_key_index(hdr, key, 0)
    i == 0 && missing_fits_keyword(key)
    val = hdr[i]
    isa(val, Integer) || noninteger_fits_keyword(key)
    to_integer(val)
end

@noinline noninteger_fits_keyword(key) =
    error("non-integer value for FITS keyword \"", key, "\"")

"""
    OIFITS.get_float(hdr, key[, def])

yields the value of keyword `key` in FITS header `hdr` as a floating-point.  If
`key` is not found in `hdr`, the default value `def` is returned if it is
specified, otherwise an error is thrown.  If `key` is found, it is checked that
it can be converted into a `Float64` which is returned.  The possible types of
the result are thus `Float64` (if `key` is found) or `typeof(def)` (if `key`
not found and `def` specified).

"""
function get_float(hdr::FITSHeader, key::String, def)
    i = get_key_index(hdr, key, 0)
    i == 0 && return def
    val = hdr[i]
    isa(val, Real) || nonfloat_fits_keyword(key)
    to_float(val)
end

function get_float(hdr::FITSHeader, key::String)
    i = get_key_index(hdr, key, 0)
    i == 0 && missing_fits_keyword(key)
    val = hdr[i]
    isa(val, Real) || nonfloat_fits_keyword(key)
    to_float(val)
end

@noinline nonfloat_fits_keyword(key) =
    error("non-floating point value for FITS keyword \"", key, "\"")

"""
    OIFITS.get_logical(hdr, key[, def])

yields the value of keyword `key` in FITS header `hdr` as a boolean.  If `key`
is not found in `hdr`, the default value `def` is returned if it is specified,
otherwise an error is thrown.  If `key` is found, it is checked that it can be
converted into a `Bool` which is returned.  The possible types of the result
are thus `Bool` (if `key` is found) or `typeof(def)` (if `key` not found and
`def` specified).

""" get_logical

function get_logical(hdr::FITSHeader, key::String, def)
    i = get_key_index(hdr, key, 0)
    i == 0 && return def
    val = hdr[i]
    isa(val, Bool) || nonlogical_fits_keyword(key)
    to_logical(val)
end

function get_logical(hdr::FITSHeader, key::String)
    i = get_key_index(hdr, key, 0)
    i == 0 && missing_fits_keyword(key)
    val = hdr[i]
    isa(val, Bool) || nonlogical_fits_keyword(key)
    to_logical(val)
end

@noinline nonlogical_fits_keyword(key) =
    error("non-logical value for FITS keyword \"", key, "\"")

"""
    OIFITS.get_string(hdr, key[, def]; fix=false)

yields the value of keyword `key` in FITS header `hdr` as a string.  If `key`
is not found in `hdr`, the default value `def` is returned if it is specified,
otherwise an error is thrown.  If `key` is found, it is checked that it can be
converted into a `String` which is returned.  The possible types of the result
are thus `String` (if `key` is found) or `typeof(def)` (if `key` not found and
`def` specified).

If keyword `fix` is true, trailing spaces are removed and characters converted
to upper case letters in the returned value.  This is to follow FITS
conventions that letter case and trailing spaces are insignificant when
comparing names.

"""
function get_string(hdr::FITSHeader, key::String, def; fix::Bool = false)
    i = get_key_index(hdr, key, 0)
    i == 0 && return def
    val = hdr[i]
    isa(val, AbstractString) || nonstring_fits_keyword(key)
    fix ? fix_name(val) : to_string(val)
end

function get_string(hdr::FITSHeader, key::String; fix::Bool = false)
    i = get_key_index(hdr, key, 0)
    i == 0 && missing_fits_keyword(key)
    val = hdr[i]
    isa(val, AbstractString) || nonstring_fits_keyword(key)
    fix ? fix_name(val) : to_string(val)
end

@noinline nonstring_fits_keyword(key) =
    error("non-string value for FITS keyword \"", key, "\"")

"""
   OIFITS.fix_name(str)

converts string `str` to uppercase letters and removes trailing spaces.  This
is useful to compare names according to FITS conventions that letter case and
trailing spaces are insignificant.  Only ASCII characters are converted (see
[`OIFITS.to_upper`](@ref)).

"""
fix_name(name::AbstractString) = map(to_upper, rstrip(c -> c == ' ', name))

"""
    OIFITS.to_upper(c)

yields character `c` converted to uppercase.  Conversion is only performed for
ASCII characters, that is for `c ∈ 'a':'z'`, other characters are returned
unchanged.

"""
to_upper(c::Char) = ((c < 'a')|(c > 'z')) ? c : Char(c - 32)

"""
    OIFITS.same_name(A, B)

yields whether strings `A` and `B` are the same name according to FITS
conventions, that is ignoring trailing spaces and the case of ASCII letters.

"""
function same_name(A::AbstractString, B::AbstractString)
    na = length(A)
    nb = length(B)
    for i in 1:min(na, nb)
        ca = A[i]
        cb = B[i]
        if ca != cb
            to_upper(ca) == to_upper(cb) || return false
        end
    end
    if na > nb
        for i in nb+1:na
            A[i] == ' ' || return false
        end
    elseif na < nb
        for i in na+1:nb
            B[i] == ' ' || return false
        end
    end
    return true
end

"""
    OIFITS.to_logical(arg)

converts argument `arg` to boolean(s) of type `Bool`.

"""
to_logical(x::Bool) = x
to_logical(x::Integer) = Bool(x)
to_logical(x::Array{Bool}) = x
to_logical(x::AbstractArray{<:Integer,N}) where {N} = convert(Array{Bool,N}, x)
to_logical(x::Tuple{Vararg{Bool}}) = x
to_logical(x::Tuple{Vararg{Integer}}) = map(to_logical, x)

"""
    OIFITS.to_integer(arg)

converts argument `arg` to integer(s) of type `Int`.

"""
to_integer(x::Int) = x
to_integer(x::Integer) = Int(x)
to_integer(x::Array{Int}) = x
to_integer(x::AbstractArray{<:Integer,N}) where {N} = convert(Array{Int,N}, x)
to_integer(x::Tuple{Vararg{Int}}) = x
to_integer(x::Tuple{Vararg{Integer}}) = map(to_integer, x)

"""
    OIFITS.to_float(arg)

converts argument `arg` to float(s) of type `Float64`.

"""
to_float(x::Float64) = x
to_float(x::Real) = Float64(x)
to_float(x::Array{Float64}) = x
to_float(x::AbstractArray{<:Real,N}) where {N} = convert(Array{Float64,N}, x)
to_float(x::Tuple{Vararg{Float64}}) = x
to_float(x::Tuple{Vararg{Real}}) = map(to_float, x)

"""
    OIFITS.to_complex(arg)

converts argument `arg` to complex(es) of type `Complex{Float64}`.

"""
to_complex(x::Complex{Float64}) = x
to_complex(x::Union{Real,Complex{<:Real}}) = convert(Complex{Float64}, x)
to_complex(x::Array{Complex{Float64}}) = x
to_complex(x::AbstractArray{<:Union{Real,Complex{<:Real}},N}) where {N} =
    convert(Array{Complex{Float64},N}, x)
to_complex(x::Tuple{Vararg{Complex{Float64}}}) = x
to_complex(x::Tuple{Vararg{Union{Real,Complex{<:Real}}}}) = map(to_complex, x)

"""
    OIFITS.to_string(arg)

converts argument `arg` to string(s) of type `String`.  Argument can be a
scalar, an array or a tuple.

""" to_string
# `string(x::SubString)` calls `String(x)` so the two are as
# fast. `string(x::Symbol)` is 2-3 tiles slower than `String(x::Symbol)`.
# Hence we call `String`, not `string` in all other cases.
to_string(x::String) = x
to_string(x::AbstractString) = String(x)
to_string(x::Symbol) = String(x)
to_string(x::Array{String}) = x
to_string(x::AbstractArray{T}) where {T<:Union{AbstractString,Symbol}} =
    map(to_string, x)
to_string(x::Tuple{Vararg{String}}) = x
to_string(x::Tuple{Vararg{Union{AbstractString,Symbol}}}) = map(to_string, x)

"""
    OIFITS.is_logical(arg)

yields whether argument `arg` is a logical value, or an array or a tuple of
logical values.

"""
is_logical(::Any) = false
is_logical(::Bool) = true
is_logical(::Array{<:Bool}) = true
is_logical(::Tuple{Vararg{Bool}}) = true

"""
    OIFITS.is_string(arg)

yields whether argument `arg` is a string value, or an array or a tuple of
string values.

"""
is_string(::Any) = false
is_string(::AbstractString) = true
is_string(::Array{<:AbstractString}) = true
is_string(::Tuple{Vararg{AbstractString}}) = true

"""
    OIFITS.is_integer(arg)

yields whether argument `arg` is an integer value, or an array or a tuple of
integers.

"""
is_integer(::Any) = false
is_integer(::Integer) = true
is_integer(::Array{<:Integer}) = true
is_integer(::Tuple{Vararg{Integer}}) = true

"""
    OIFITS.is_float(arg)

yields whether argument `arg` is a real value, or an array or a tuple of real
values.  The name means that `to_float(arg)` can be applied to `arg`.

"""
is_float(::Any) = false
is_float(::Real) = true
is_float(::Array{<:Real}) = true
is_float(::Tuple{Vararg{Real}}) = true

"""
    OIFITS.is_complex(arg)

yields whether argument `arg` is a complex value, or an array or a tuple of
complex values.  The name means that `to_complex(arg)` can be applied to `arg`.

"""
is_complex(::Any) = false
is_complex(::Complex) = true
is_complex(::Array{<:Complex{<:Real}}) = true
is_complex(::Tuple{Vararg{Complex{<:Real}}}) = true

"""
    OIFITS.to_fieldname(s)

converts OI-FITS column name `s` into a `Symbol` that can be used as a field
name in OI-FITS structures.  Letters are converted to lowercase and
non-alphanumeric characters replaced by underscores.

"""
to_fieldname(sym::Symbol) = sym

function to_fieldname(name::AbstractString)
    key = lowercase(name)
    key == "oi_revn" ? :revn : Symbol(replace(key, r"[^a-z0-9_]" => '_'))
end

"""
    OIFITS.convert_key_type(K, dict)

yields a dictionary with the same values as `dict` but whose keys have type
`K`.  If `keytype(dict) == K`, `dict` itself is returned.

"""
convert_key_type(::Type{K}, dict::AbstractDict{K}) where {K} = dict

function convert_key_type(::Type{K}, dict::AbstractDict{<:Any,V}) where {K,V}
    newdict = Dict{K,V}()
    for (key, val) in dict
        newdict[K(key)] = val
    end
    return newdict
end

"""
    OIFITS.warn([io=stderr,] args...)

prints a warning message made of `args...` to the stream `io`, the standard
error output by default.

"""
warn(args...) = warn(stderr, args...)
warn(io::IO, args...) = message(io, args...; head="WARNING", color=:yellow)

"""
    OIFITS.inform([io=stdout,] args...)

prints an informational message made of `args...` to the stream `io`, the
standard output by default.

"""
inform(args...) = inform(stdout, args...)
inform(io::IO, args...) = message(io, args...; head="INFO", color=:blue)

"""
    OIFITS.message([io=stdout,] args...; head=nothing, color=:normal)

prints a message made of `args...` to the stream `io` and using specified
`color`.  If `head` is not `nothing`, it is printed first in bold and followed
by a space.

"""
message(args...; kwds...) = message(stdout, args...; kwds...)
@noinline function message(io::IO, args...; head = nothing, color::Symbol = :normal)
    head === nothing || printstyled(io, head, " "; bold=true, color=color)
    printstyled(io, args..., "\n"; bold=false, color=color)
    nothing
end

#
# The functions `type_to_eqcoltype`, `eqcoltype_to_type`, `bitpix_to_type` and
# `type_to_bitpix` deal with conversion between CFITSIO type code or BITPIX
# value and actual Julia data types.  They can be used as follows (assuming `T`
# is a Julia data type, while `eqcoltype` and `bitpix` are integers):
#
#     type_to_eqcoltype(T) -----------> eqcoltype (e.g., TBYTE, TFLOAT, etc.)
#     eqcoltype_to_type(eqcoltype) ---> T
#
#     type_to_bitpix(T) --------------> bitpix (e.g., BYTE_IMG, FLOAT_IMG, etc.)
#     bitpix_to_type(bitpix) ---------> T
#
# The following table gives the correspondances between CFITSIO "types",
# the BITPIX keyword and Julia types.
#
#     -------------------------------------------------
#     BITPIX   CFISTIO         Julia     Comments
#     -------------------------------------------------
#              int             Cint
#              long            Clong
#              LONGLONG        Int64     64-bit integer
#     -------------------------------------------------
#        8     BYTE_IMG        UInt8
#       16     SHORT_IMG       Int16
#       32     LONG_IMG        Int32
#       64     LONGLONG_IMG    Int64
#      -32     FLOAT_IMG       Float32
#      -64     DOUBLE_IMG      Float64
#     -------------------------------------------------
#              TBIT
#              TBYTE           Cuchar = UInt8
#              TSBYTE          Cchar = Int8
#              TLOGICAL        Bool
#              TSHORT          Cshort
#              TUSHORT         Cushort
#              TINT            Cint
#              TUINT           Cuint
#              TLONG           Clong
#              TLONGLONG       Int64
#              TULONG          Culong
#              TFLOAT          Cfloat
#              TDOUBLE         Cdouble
#              TCOMPLEX        Complex{Cfloat}
#              TDBLCOMPLEX     Complex{Cdouble}
#     -------------------------------------------------

"""
    OIFITS.type_to_bitpix(T) -> val

yields the FITS BITPIX value `val` which corresponds to Julia type `T`.

"""
type_to_bitpix(::Type{UInt8})   =   8
type_to_bitpix(::Type{Int16})   =  16
type_to_bitpix(::Type{Int32})   =  32
type_to_bitpix(::Type{Int64})   =  64
type_to_bitpix(::Type{Float32}) = -32
type_to_bitpix(::Type{Float64}) = -64

"""
    OIFITS.bitpix_to_type(val) -> T

yields the Julia type `T` which corresponds to the BITPIX value `val`.

"""
bitpix_to_type(val::Integer) =
    (val ==   8 ? UInt8   :
     val ==  16 ? Int16   :
     val ==  32 ? Int32   :
     val ==  64 ? Int64   :
     val == -32 ? Float32 :
     val == -64 ? Float64 : bad_bitpix(val))

@noinline bad_bitpix(val) =
    error(string("bad BITPIX value ", val))

"""
    OIFITS.eqcoltype_to_type(val) -> T

yields the Julia type `T` which best corresponds to the FITSIO equivalent
column type `val` or `Nothing` if there is no correspondence.

"""
eqcoltype_to_type(val::Integer) =
    (#   ==   1 ? ...              : # TBIT
     val ==  11 ? UInt8            : # TBYTE
     val ==  12 ? Int8             : # TSBYTE
     val ==  14 ? Bool             : # TLOGICAL
     val ==  16 ? String           : # TSTRING
     val ==  20 ? Cushort          : # TUSHORT      (Uint16?)
     val ==  21 ? Cshort           : # TSHORT       (Int16?)
     val ==  30 ? Cuint            : # TUINT        (Uint32?)
     val ==  31 ? Cint             : # TINT         (Int32?)
     val ==  40 ? Culong           : # TULONG
     val ==  41 ? Clong            : # TLONG
     #   ==  42 ? ...              : # TINT32BIT
     val ==  42 ? Cfloat           : # TFLOAT      (Float32?)
     val ==  81 ? Int64            : # TLONGLONG
     val ==  82 ? Cdouble          : # TDOUBLE     (Float64?)
     val ==  83 ? Complex{Cfloat}  : # TCOMPLEX    (ComplexF32?)
     val == 163 ? Complex{Cdouble} : # TDBLCOMPLEX (ComplexF64?)
     Nothing)

"""
    OIFITS.type_to_eqcoltype(T) -> val

yields the FITSIO equivalent column type `val` which best corresponds to the
Julia type `T`.

"""
type_to_eqcoltype(::Type{UInt8})            =  11
type_to_eqcoltype(::Type{Int8})             =  12
type_to_eqcoltype(::Type{Bool})             =  14
type_to_eqcoltype(::Type{String})           =  16
type_to_eqcoltype(::Type{Cushort})          =  20
type_to_eqcoltype(::Type{Cshort})           =  21
type_to_eqcoltype(::Type{Cuint})            =  30
type_to_eqcoltype(::Type{Cint})             =  31
@static if Culong !== Cuint
    type_to_eqcoltype(::Type{Culong})       =  40
end
@static if Clong !== Cint && Clong !== Int64
    type_to_eqcoltype(::Type{Clong})        =  41
end
type_to_eqcoltype(::Type{Cfloat})           =  42
type_to_eqcoltype(::Type{Int64})            =  81
type_to_eqcoltype(::Type{Cdouble})          =  82
type_to_eqcoltype(::Type{Complex{Cfloat} }) =  83
type_to_eqcoltype(::Type{Complex{Cdouble}}) = 163


text:
""" Computing with permutations: the computer algebra of the symmetric group.
"""
module Permutations
export cycles, permutation_to_expr, adjacent_transpositions_by_bubble_sort!,
  adjacent_transpositions_by_insertion_sort!

using ...Syntax
using ...Theories: dom, codom, compose, id, otimes, munit, braid

# Decomposition
###############

""" Decompose a permutation into its cycles.

Returns a vector of vectors, the cycles of the permutation.
"""
cycles(σ::AbstractVector{Int}) = cycles(i -> σ[i], length(σ))

function cycles(σ, n::Int)
  cycles = Vector{Int}[]
  used = falses(n)
  for i in 1:n
    if used[i]; continue end
    cycle, j = Int[], i
    while true
      push!(cycle, j)
      used[j] = true
      j = σ(j)
      @assert 1 <= j <= n
      if j == i; break end
    end
    push!(cycles, cycle)
  end
  cycles
end

""" Decompose permutation into adjacent transpositions using bubble sort.

An *adjacent transposition*, also known as a *simple transposition*, is a
transposition of form (i i+1), represented here as simply the number i.

This algorithm appears as Algorithm 2.7 in the PhD thesis of Jonathan Huang,
"Probabilistic reasonsing and learning on permutations: Exploiting structural
decompositions of the symmetric group". As Huang notes, the algorithm is
very similar to the well-known bubble sort. It has quadratic complexity.

See also: [`adjacent_transpositions_by_insertion_sort!`](@ref).
"""
function adjacent_transpositions_by_bubble_sort!(σ::AbstractVector{Int})
  n = length(σ)
  transpositions = Int[]
  for i in 1:n-1
    for j = n-1:-1:i
      if σ[j+1] < σ[j]
        σ[j], σ[j+1] = σ[j+1], σ[j]
        push!(transpositions, j)
      end
    end
  end
  transpositions
end

""" Decompose permutation into adjacent transpositions using insertion sort.

An *adjacent transposition*, also known as a *simple transposition*, is a
transposition of form (i i+1), represented here as simply the number i.

Bubble sort and insertion sort are, in a sense, dual algorithms (Knuth, TAOCP,
Vol 3: Searching and Sort, Sec 5.3.4: Networks for sorting, Figures 45 & 46). A
minimal example on which they give different decompositions is the permutation:

  [1,2,3] ↦ [3,2,1]

See also: [`adjacent_transpositions_by_bubble_sort!`](@ref).
"""
function adjacent_transpositions_by_insertion_sort!(σ::AbstractVector{Int})
  n = length(σ)
  transpositions = Int[]
  for i in 2:n
    for j in i:-1:2
      if σ[j-1] > σ[j]
        σ[j-1], σ[j] = σ[j], σ[j-1]
        push!(transpositions, j-1)
      end
    end
  end
  transpositions
end

# Conversion to expression
##########################

""" Convert a typed permutation into a morphism expression.

Warning: The morphism expression is not simplified.
"""
function permutation_to_expr(σ::AbstractVector{Int}, xs::AbstractVector;
                             sort::Symbol=:insertion)
  permutation_to_expr!(copy(σ), copy(xs); sort=sort)
end
function permutation_to_expr!(σ::AbstractVector{Int}, xs::AbstractVector;
                              sort::Symbol=:insertion)
  n = length(σ)
  @assert length(xs) == n

  transpositions = if sort == :bubble
    adjacent_transpositions_by_bubble_sort!(σ)
  elseif sort == :insertion
    adjacent_transpositions_by_insertion_sort!(σ)
  else
    error("Sorting algorithm not supported: $sort")
  end
  if isempty(transpositions)
    return id(otimes(xs))
  end

  layers = map(transpositions) do τ
    layer = [
      τ > 1 ? id(otimes(xs[1:τ-1])) : nothing,
      braid(xs[τ], xs[τ+1]),
      τ+1 < n ? id(otimes(xs[τ+2:n])) : nothing,
    ]
    xs[τ], xs[τ+1] = xs[τ+1], xs[τ]
    foldl(otimes, filter(!isnothing, layer))
  end
  foldl(compose, layers)
end

end


text:
using CSV, GLM, StatsBase, Random, Distributions, Plots, Measures; pyplot()
Random.seed!(0)

n, N = 10^3, 10^6

function waldWolfowitz(data)
    n = length(data)
    sgns = data .>= mean(data);
    nPlus, nMinus = sum(sgns), n - sum(sgns)
    wwMu = 2*nPlus*nMinus/n + 1
    wwVar = (wwMu-1)*(wwMu-2)/(n-1)

    R = 1
    for i in 1:n-1
        R += sgns[i] != sgns[i+1]
    end

    zStat = abs((R-wwMu)/sqrt(wwVar))
    2*ccdf(Normal(),zStat)
end

experimentPvals = [waldWolfowitz(rand(Normal(),n)) for _ in 1:N]
for alpha in [ 0.001, 0.005, 0.01, 0.05, 0.1]
    pva = sum(experimentPvals .< alpha)/N
    println("For alpha = $(alpha), p-value area = $(pva)")
end

p1 = histogram(experimentPvals,bins = 5n, legend = false,
    xlabel = "p-value", ylabel = "Frequency")

Fhat = ecdf(experimentPvals)

pGrid = 0:0.001:1
p2 = plot(pGrid,Fhat.(pGrid), legend =false, xlabel = "p-value", ylabel = "ECDF")
plot(p1, p2, size = (1000,400), margin = 5mm)

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
module BondiToy

export Param, NestedIBVP, CoupledIBVP, run_toy

using Parameters
using DifferentialEquations
using Interpolations
using HDF5
using Printf
using RecursiveArrayTools

@with_kw struct Param
    # discretization parameters
    NX                :: Int
    Nz                :: Int
    cX                :: Float64 = 10.0
    rmin              :: Float64 = 2.0
    umax              :: Float64
    # directory to save data
    out_dir           :: String

    # how often to save data
    out_every               :: Int
    compute_L2_norm_every   :: Int = 1 # set to negative values to turn off
    compute_Lop_norm_every  :: Int = 1 # set to negative values to turn off
end

"""
Extend this type for the initial and boundary condition parameters. The IBVP can
be nested (each intrinsic equation can be solved one after the other) or coupled
(the intrinsic equations are coupled).
"""
abstract type IBVP end

"""
The NestedIBVP structure encodes all the models with nested intrinsic
equations, both SH and WH, and for sources B0, B1 and B2.  It needs
the paremeters az_21, b11, b13, b21, b22, b33, which are the elements
of the angular principal matrix Az and the lower order matrix B,
respectively.
"""
abstract type NestedIBVP  <: IBVP end
abstract type CoupledIBVP <: IBVP end

"""
Function that defines the BCs of ϕ at r=rmin

Needs signature ϕ0_of_uz(u::T, z::T, ibvp::IBVP) where {T<:Real}
"""
function ϕ0_of_uz end

"""
Function that defines the BCs of ψv at r=rmin

Needs signature ψv0_of_uz(u::T, z::T, ibvp::IBVP) where {T<:Real}
"""
function ψv0_of_uz end

"""
Function that defines the ingoing mode of ψ at u=0, i.e. ψ(u=0,X,z)

Needs signature ψ0_of_Xz(X::T, z::T, ibvp::IBVP) where {T<:Real}
"""
function ψ0_of_Xz end

rtoX(r, cX, rmin) = (r .- rmin) ./ sqrt.(cX*cX .+ (r .- rmin).^2)
Xtor(X, cX, rmin) = rmin .+ cX * X ./ sqrt.(1.0 .- X.^2)
dX_dr(X, cX) = (1.0 .- X.^2).^1.5 ./ cX

struct System{T}
    X    :: Vector{T}
    r    :: Vector{T}
    dXdr :: Vector{T}
    hX   :: T
    z    :: Vector{T}
    hz   :: T
end
function System(p::Param)
    hX = 1.0 / (p.NX-1)
    X  = range(0, 1, length=p.NX)

    hz = 2.0*π / (p.Nz) # remove last point; 0=2π (periodic)
    z  = range(0, 2.0*π-hz, length=p.Nz)

    rr    = Xtor(X, p.cX, p.rmin)
    dXdr  = dX_dr(X, p.cX)

    System{typeof(hX)}(X, rr, dXdr, hX, z, hz)
end

# 2nd order accurate finite differences for 1st order derivatives along a
# periodic direction
function Dz!(f_z::Matrix, f::Matrix, sys::System)
    NX, Nz = size(f)
    odz2 = 0.5 / sys.hz

    @inbounds for j in 2:Nz-1
        @inbounds for i in 1:NX
            f_z[i,j] = (f[i,j+1] - f[i,j-1]) * odz2
        end
    end

    @inbounds for i in 1:NX
        f_z[i,1]   = (f[i,2] - f[i,end]) * odz2
        f_z[i,end] = (f[i,1] - f[i,end-1]) * odz2
    end

    f_z
end

# now acting on vectors
function Dz!(f_z::Vector, f::Vector, sys::System)
    Nz = length(f)
    odz2 = 0.5 / sys.hz

    @inbounds for j in 2:Nz-1
        f_z[j] = (f[j+1] - f[j-1]) * odz2
    end

    f_z[1]   = (f[2] - f[end]) * odz2
    f_z[end] = (f[1] - f[end-1]) * odz2

    f_z
end

function Dz(f, sys::System)
    f_z = similar(f)
    Dz!(f_z, f, sys)
end

# 2nd order accurate finite difference operator for 1st order derivatives along
# a non-periodic direction. forward and backward 2nd order accurate finite
# differences are used on the first and last points of the grid respectively
function DX!(f_X::Matrix, f::Matrix, sys::System)
    NX, Nz = size(f)
    odX2 = 0.5 / sys.hX

    @inbounds for j in 1:Nz
        @inbounds for i in 2:NX-1
            f_X[i,j] = (f[i+1,j] - f[i-1,j]) * odX2
        end
    end

    @inbounds for j in 1:Nz
        f_X[1,j]   = (-3.0* f[1,j]   + 4.0*f[2,j]     - f[3,j]) * odX2
        f_X[end,j] = ( 3.0* f[end,j] - 4.0*f[end-1,j] + f[end-2,j]) * odX2
    end

    f_X
end
function DX(f, sys::System)
    f_X = similar(f)
    DX!(f_X, f, sys)
end

function init_ψ(sys::System{T}, ibvp::IBVP) where {T}
    NX = length(sys.X)
    Nz = length(sys.z)

    ψ0 = zeros(T, NX, Nz)

    @inbounds for j in 1:Nz
        @inbounds for i in 1:NX
            ψ0[i,j] = ψ0_of_Xz(sys.X[i], sys.z[j], ibvp)
        end
    end
    ψ0
end

# integrate in the null hypersurface for nested intrinsic eqs
function get_ϕψv!(ϕ, ψv, ψ, u, sys::System, ibvp::NestedIBVP)
    NX, Nz = size(ψ)

    S_ϕ  = ibvp.b13 * ψ

    Threads.@threads for j in 1:Nz
        itp  = interpolate( S_ϕ[:,j], BSpline(Cubic(Flat(OnGrid()))) )
        sitp = scale(itp, sys.X[1]:sys.hX:sys.X[end])
        rhs_ϕ(f, p, x) = sitp(x) + ibvp.b11 * f

        Xspan = (sys.X[1], sys.X[end])

        # this defines the outgoing mode i.e. boundary condition at each
        # timestep
        ϕ0 = ϕ0_of_uz(u, sys.z[j], ibvp)
        # define the PDE problem
        prob_ϕ = ODEProblem(rhs_ϕ, ϕ0, Xspan)
        # solve the PDE problem
        sol_ϕ  = solve(prob_ϕ, SSPRK22(), dt=sys.hX, adaptive=false)
        # pass the above solution
        ϕ[:, j] = sol_ϕ.(sys.X)
    end

    ϕ_z = Dz(ϕ, sys)
    S_ψv  = ibvp.az_21 * ϕ_z .+ ibvp.b21 * ϕ .+ ibvp.b23 * ψ

    Threads.@threads for j in 1:Nz
        itp  = interpolate( S_ψv[:,j], BSpline(Cubic(Flat(OnGrid()))) )
        sitp = scale(itp, sys.X[1]:sys.hX:sys.X[end])
        rhs_ψv(f, p, x) = sitp(x) + ibvp.b22 * f

        Xspan = (sys.X[1], sys.X[end])

        # this defines the outgoing mode i.e. boundary condition at each
        # timestep
        ψv0 = ψv0_of_uz(u, sys.z[j], ibvp)
        # define the PDE problem
        prob_ψv = ODEProblem(rhs_ψv, ψv0, Xspan)
        # solve the PDE problem
        sol_ψv  = solve(prob_ψv, SSPRK22(), dt=sys.hX, adaptive=false)
        # pass the above solution
        ψv[:, j] = sol_ψv.(sys.X)
    end

    ϕ, ψv
end

# integrate in the null hypersurface for non-nested intrinsic eqs
function intrinsic_rhs!(dv, v, sys, x)
    dϕ  = dv[1]
    dψv = dv[2]
    ϕ   = v[1]
    ψv  = v[2]

    # the rhs for the intrinsic coupled PDE of ϕ and ψv
    dϕ .= ψv
    Dz!(dψv, ϕ, sys)  # dψv = ∂_z ϕ
    dv
end
function get_ϕψv!(ϕ, ψv, ψ, u, sys::System{T}, ibvp::CoupledIBVP) where {T}
    NX, Nz = size(ψ)
    ϕ0  = zeros(T, Nz)
    ψv0 = zeros(T, Nz)

    # this defines the outgoing mode i.e. boundary condition at each timestep
    @inbounds for j in 1:Nz
        ϕ0[j]  = ϕ0_of_uz(u, sys.z[j], ibvp)
        ψv0[j] = ψv0_of_uz(u, sys.z[j], ibvp)
    end

    # save the outgoing mode at X=X0 to ϕ and ψv of (x,z)
    ϕ[1,:]  = ϕ0
    ψv[1,:] = ψv0

    # write the ID as a vector for the coupled PDE system
    v0 = VectorOfArray([ϕ0, ψv0])

    # define the x-span to solve the PDE for
    Xspan = (sys.X[1], sys.X[end])

    # write the PDE as an ODE prob
    intrinsic_prob = ODEProblem(intrinsic_rhs!, v0, Xspan, sys)

    # def the x-integrator
    x_integrator = init(intrinsic_prob, SSPRK22(),
                        save_everytimestep=false, dt=sys.hX, adaptive=false)

    # the following assumes that the integrator is stepping *exactly* at the X
    # grid points. this should be true with the above options of dt=sys.hX,
    # adaptive=false, but we must remember to change things if using an adaptive
    # integrator or something similar.
    iter = 1
    for (f,t) in tuples(x_integrator)
        iter += 1

        ϕ[iter,:]  = f[1]
        ψv[iter,:] = f[2]
    end

    ϕ, ψv
end
function get_ϕψv(ψ, u, sys, ibvp)
    ϕ    = similar(ψ)
    ψv   = similar(ψ)
    get_ϕψv!(ϕ, ψv, ψ, u, sys, ibvp)
end

# the rhs for the time evolution equation
function get_ψ_u!(ψ_u, ψ, ϕ, ψv, sys, ibvp)
    ψ_X = DX(ψ, sys)
    ψ_z = Dz(ψ, sys)
    ψ_r = sys.dXdr .* ψ_X

    ψ_u .= 0.5 * ψ_r .+ ψ_z .+ ibvp.b31 * ϕ .+ ibvp.b32 * ψv .+ ibvp.b33 * ψ
end
function rhs_ψ!(dψ, ψ, (sys, ibvp), u)
    ϕ, ψv = get_ϕψv(ψ, u, sys, ibvp)
    get_ψ_u!(dψ, ψ, ϕ, ψv, sys, ibvp)
    nothing
end

function write_2D(it::Int, t, data_dir::String, ψ::Array, ψv::Array, ϕ::Array)
    it_str  = lpad(it, 4, "0")
    outfile = joinpath(data_dir, "data_$(it_str).h5")
    h5open(outfile, "w") do file
        write(file, "ψ",  ψ)
        write(file, "ψv", ψv)
        write(file, "ϕ",  ϕ)
        attrs(file)["time"] = t
    end
    nothing
end

# function that performs the time evolution
function run_toy(p::Param, ibvp::IBVP)

    # pass the parameters of the system
    sys = System(p)

    # create the folders where data are saved
    data_dir = mkpath(p.out_dir)

    # initialize ψ
    ψ = init_ψ(sys, ibvp)

    # time span of the simulation
    tspan = (0.0, p.umax)

    # timestep. careful with CFL condition
    dt0   = 0.25 * minimum([sys.hX, sys.hz])

    # define the PDE problem for time integration
    prob  = ODEProblem(rhs_ψ!, ψ, tspan, (sys, ibvp))
    # http://docs.juliadiffeq.org/latest/basics/integrator.html
    integrator = init(prob, RK4(), save_everystep=false, dt=dt0, adaptive=false)

    # write the coordinates. if there's a problem with this for windows users,
    # it can be done differently (with joinpath, or similar)
    h5write(data_dir*"/x.h5", "x", sys.X)
    h5write(data_dir*"/z.h5", "z", sys.z)
    h5write(data_dir*"/r.h5", "r", sys.r)

    it = 0
    t  = 0.0
    ϕ, ψv = get_ϕψv(ψ, t, sys, ibvp)

    # save initial data
    write_2D(it, t, data_dir, ψ, ψv, ϕ)

    if p.compute_L2_norm_every > 0
        # initiate the grid function to be used for the sum in the time domain
        # here we sum over the z domain for the initial time
        l2_out_x  = dt0*sys.hz*sum(ψv.*ψv + ϕ.*ϕ, dims=2)
        L2_norm = sqrt(sys.hz*sys.hX*sum(ψ.*ψ)) + maximum(sqrt.(l2_out_x))
        outfile = joinpath(data_dir, "L2_norm.dat")
        open(outfile, "w") do io
            println(io, "# t        |      norm_L2")
            println(io, "$t  \t  $L2_norm")
        end
    end

    if p.compute_Lop_norm_every > 0
        ϕz = Dz(ϕ, sys)
        # initiate the grid function to be used for the sum in the time domain
        # here we sum over the z domain for the initial time
        lop_out_x = dt0*sys.hz*sum(ψv.*ψv + ϕ.*ϕ + ϕz.*ϕz, dims=2)
        Lopside_norm = sqrt(sys.hz*sys.hX*sum(ψ.*ψ)) + maximum(sqrt.(lop_out_x))
        outfile = joinpath(data_dir, "Lopside_norm.dat")
        open(outfile, "w") do io
            println(io, "# t        |      norm_Lop")
            println(io, "$t  \t  $Lopside_norm")
        end
    end

    println("-------------------------------------------------------------")
    println("Iteration      Time |            ψ ")
    println("                    |    minimum      maximum")
    println("-------------------------------------------------------------")
    @printf "%9d %9.3f |  %9.4g    %9.4g\n" it t minimum(ψ) maximum(ψ)

    # start time evolution
    for (f,t) in tuples(integrator)
        it += 1

        ψ = f
        get_ϕψv!(ϕ, ψv, ψ, t, sys, ibvp)

        @printf "%9d %9.3f |  %9.4g    %9.4g\n" it t minimum(ψ) maximum(ψ)

        if p.compute_L2_norm_every > 0 && it % p.compute_L2_norm_every == 0
            # add on the l2_out_x grid function the z sums for this time instant
            l2_out_x + = dt0*sys.hz*sum(ψv.*ψv + ϕ.*ϕ, dims=2)
            L2_norm = sqrt(sys.hz*sys.hX*sum(ψ.*ψ)) + maximum(sqrt.(l2_out_x))
            outfile = joinpath(data_dir, "L2_norm.dat")
            open(outfile, "a") do io
                println(io, "$t  \t  $L2_norm")
            end
        end

        if p.compute_Lop_norm_every > 0 && it % p.compute_Lop_norm_every == 0
            Dz!(ϕz, ϕ, sys)
            # add on the lop_out_x grid function the z sums for this time instant
            lop_out_x += dt0*sys.hz*sum(ψv.*ψv + ϕ.*ϕ + ϕz.*ϕz, dims=2)
            Lopside_norm = sqrt(sys.hz*sys.hX*sum(ψ.*ψ)) + maximum(sqrt.(lop_out_x))
            outfile = joinpath(data_dir, "Lopside_norm.dat")
            open(outfile, "a") do io
                println(io, "$t  \t  $Lopside_norm")
            end
        end

        if it % p.out_every == 0
            write_2D(it, t, data_dir, ψ, ψv, ϕ)
        end
    end

    println("-------------------------------------------------------------")
    println("Done.")
end

end # module

max_stars_repo_name: ThanasisGiannakopoulos/BondiToy

max_issues_repo_name: ThanasisGiannakopoulos/BondiToy

max_forks_repo_name: ThanasisGiannakopoulos/BondiToy


text:
samples = hcat([rand(predictor_joint) for sample = 1:1_000]...)
modes = [mode(samples[i, :]) for i = 1:size(samples, 1)]
overall_accuracy = Statistics.mean(modes .== iris[!, label_column])
@test overall_accuracy > 0.9 # make sure the overall accuracy is greater than 90%
classes = sort(unique(iris[!, label_column]))
@test length(classes) == 3
@test classes == ["setosa", "versicolor", "virginica"]
setosa_rows = findall(iris[!, label_column] .== "setosa")
versicolor_rows = findall(iris[!, label_column] .== "versicolor")
virginica_rows = findall(iris[!, label_column] .== "virginica")
@test length(setosa_rows) == 50
@test length(versicolor_rows) == 50
@test length(virginica_rows) == 50
@test length(iris[setosa_rows, label_column]) == 50
@test length(iris[versicolor_rows, label_column]) == 50
@test length(iris[virginica_rows, label_column]) == 50
setosa_accuracy  = Statistics.mean(modes[setosa_rows] .== iris[setosa_rows, label_column])
versicolor_accuracy  = Statistics.mean(modes[versicolor_rows] .== iris[versicolor_rows, label_column])
virginica_accuracy  = Statistics.mean(modes[virginica_rows] .== iris[virginica_rows, label_column])
@test setosa_accuracy > 0.9 # make sure the accuracy on the "setosa" class is greater than 90%
@test versicolor_accuracy > 0.9 # make sure the accuracy on the "setosa" class is greater than 90%
@test virginica_accuracy > 0.9 # make sure the accuracy on the "setosa" class is greater than 90%

predictor_marginal = MLJBase.predict(mach, iris[1:size(iris, 1), feature_columns])
@test predictor_marginal isa MLJBase.UnivariateFiniteVector
@test size(predictor_marginal) == (size(iris, 1),)

# Now, we make sure that the cross-validated accuracy is greater than 90%
nfolds = 4
evaluation_results = evaluate!(mach, resampling=CV(; nfolds = nfolds, shuffle = true), measure=MLJBase.accuracy, operation=MLJBase.predict_mode)
@test evaluation_results.measurement isa Vector
@test length(evaluation_results.measurement) == 1
@test evaluation_results.measurement[1] > 0.9 # make sure that the overall cross-validated accuracy is greater than 90%
@test evaluation_results.per_fold isa Vector
@test length(evaluation_results.per_fold) == 1
@test evaluation_results.per_fold[1] isa Vector
@test length(evaluation_results.per_fold[1]) == nfolds # make sure that the accuracy is greater than 80% in each fold
@test all(evaluation_results.per_fold[1] .> 0.8)


text:
const lookup = [i / 255 for i in 0:510]
@inline function unsafe_load_UInt64(v::Vector{UInt8}, i::Integer)
    p = convert(Ptr{UInt64}, pointer(v, i))
    unsafe_load(p)
end

"""
    Genotypes{T}(p::Preamble, d::Vector{UInt8}) where T <: AbstractFloat
Create `Genotypes` struct from the preamble and decompressed data string.
"""
function Genotypes{T}(p::Preamble, d::Vector{UInt8}) where T <: AbstractFloat
    Genotypes{T}(p, d, nothing, 0, nothing, false, false)
end

const zlib = ZlibDecompressor()

@inline function zstd_uncompress!(input::Vector{UInt8}, output::Vector{UInt8})
    r = ccall((:ZSTD_decompress, CodecZstd.libzstd),
        Csize_t, (Ptr{Cchar}, Cint, Ptr{Cchar}, Cint),
        pointer(output), length(output), pointer(input), length(input))
    @assert r == length(output) "zstd decompression returned data of wrong length"
end

@inline function check_decompressed_length(io, v, h)
    seek(io, v.geno_offset)
    decompressed_field = 0
    if h.compression != 0
        if h.layout == 1
            decompressed_length = 6 * h.n_samples
        elseif h.layout == 2
            decompressed_field = 4
            decompressed_length = read(io, UInt32)
        end
    end
    return decompressed_length, decompressed_field
end

"""
    decompress(io, v, h; decompressed=nothing)
Decompress the compressed byte string for genotypes.
"""
function decompress(io::IOStream, v::Variant, h::Header;
    decompressed::Union{Nothing, AbstractVector{UInt8}}=nothing
    )
    compression = h.compression
    seek(io, v.geno_offset)
    decompressed_length, decompressed_field = check_decompressed_length(io, v, h)
    if decompressed !== nothing
        @assert length(decompressed) ==
            decompressed_length "decompressed length mismatch"
    else
        decompressed = Vector{UInt8}(undef, decompressed_length)
    end
    compressed_length = v.next_var_offset - v.geno_offset - decompressed_field
    buffer_compressed = read(io, compressed_length)
    if compression == 0
        decompressed .= buffer_compressed
    else
        if compression == 1
            codec = zlib
        elseif compression == 2
            codec = nothing
        else
            @error "invalid compression"
        end

        if compression == 1
            input = Buffer(buffer_compressed)
            output = Buffer(decompressed)
            error = Error()
            initialize(codec)
            _, _, e = process(codec, buffermem(input), buffermem(output), error)
            if e === :error
                throw(error[])
            end
            finalize(codec)
        elseif compression == 2
            zstd_uncompress!(buffer_compressed, decompressed)
        end
    end
    return decompressed
end

"""
    parse_ploidy(ploidy, d, idx, n_samples)
Parse ploidy part of the preamble.
"""
function parse_ploidy(ploidy::Union{UInt8,AbstractVector{UInt8}}, d::AbstractVector{UInt8},
    n_samples::Integer)
    missings = Int[]
    mask = 0x3f # 63 in UInt8
    mask_8 = 0x8080808080808080 # UInt64, mask for missingness
    idx1 = 9
    if typeof(ploidy) == UInt8 # if constant ploidy, just scan for missingness
        # check eight samples at a time

        if n_samples >= 8
            @inbounds for i in 0:8:(n_samples - (n_samples % 8) - 1)
                if mask_8 & unsafe_load_UInt64(d, idx1 + i) != 0
                    for j in (i+1):(i+8)
                        if d[idx1 + j - 1] & 0x80 != 0
                            push!(missings, j)
                        end
                    end
                end

            end
        end
        # remainder not in multiple of 8
        @inbounds for j in (n_samples - (n_samples % 8) + 1):n_samples
            if d[idx1 + j - 1] & 0x80 != 0
                push!(missings, j)
            end
        end
    else
        @inbounds for j in 1:n_samples
            ploidy[j] = mask & d[idx1 + j - 1]
            if d[idx1 + j - 1] & 0x80 != 0
                push!(missings, j)
            end
        end
    end
    return missings
end

@inline function get_max_probs(max_ploidy, n_alleles, phased)
    phased == 1 ?
        n_alleles : binomial(max_ploidy + n_alleles - 1, n_alleles - 1)
end

"""
    parse_preamble(d, idx, h, v)
Parse preamble of genotypes.
"""
function parse_preamble(d::AbstractVector{UInt8}, h::Header, v::Variant)
    startidx = 1
    if h.layout == 1
        n_samples = h.n_samples
        n_alleles = 2
        phased = false
        min_ploidy = 0x02
        max_ploidy = 0x02
        bit_depth = 16
    elseif h.layout == 2
        n_samples = reinterpret(UInt32, @view(d[startidx:startidx+3]))[1]
        startidx += 4
        @assert n_samples == h.n_samples "invalid number of samples"
        n_alleles = reinterpret(UInt16, @view(d[startidx:startidx+1]))[1]
        startidx += 2
        @assert n_alleles == v.n_alleles "invalid number of alleles"
        min_ploidy = d[startidx]
        startidx += 1
        max_ploidy = d[startidx]
        startidx += 1
    else
        @error "invalid layout"
    end
    constant_ploidy = (min_ploidy == max_ploidy)

    if constant_ploidy
        ploidy = max_ploidy
    else
        ploidy = Vector{UInt8}(undef, n_samples)
        fill!(ploidy, 0)
    end
    missings = []
    if h.layout == 2
        # this function also parses missingness.
        missings = parse_ploidy(ploidy, d, n_samples)
        startidx += n_samples
        phased = d[startidx]
        startidx += 1
        bit_depth = d[startidx]
        startidx += 1
    end
    max_probs = get_max_probs(max_ploidy, n_alleles, phased)
    Preamble(n_samples, n_alleles, phased, min_ploidy, max_ploidy, ploidy,
        bit_depth, max_probs, missings)
end

"""
    parse_layout1!(data, p, d, startidx)
Parse probabilities from layout 1.
"""
function parse_layout1!(data::AbstractArray{<:AbstractFloat},
    p::Preamble, d::AbstractArray{UInt8}, startidx::Integer
    )
    @assert length(data) == p.n_samples * p.max_probs
    factor = 1.0 / 32768
    idx = startidx
    @fastmath @inbounds for i in 1:p.max_probs:(p.n_samples * p.max_probs)
        j = idx
        data[i]     = reinterpret(UInt16, @view(d[j     : j + 1]))[1] * factor
        data[i + 1] = reinterpret(UInt16, @view(d[j + 2 : j + 3]))[1] * factor
        data[i + 2] = reinterpret(UInt16, @view(d[j + 4 : j + 5]))[1] * factor
        idx += 6

        # triple zero denotes missing for layout1
        if data[i] == 0.0 && data[i+1] == 0.0 && data[i+2] == 0.0
            data[i:i+2] .= NaN
            push!(p.missings, (i-1) ÷ p.max_probs + 1)
        end
    end
    return data
end

"""
   parse_layout2!(data, p, d, startidx)
Parse probabilities from layout 2.
"""
function parse_layout2!(data::AbstractArray{<:AbstractFloat},
    p::Preamble, d::AbstractArray{UInt8}, startidx::Integer
    )
    constant_ploidy = p.max_ploidy == p.min_ploidy
    if p.phased == 0
        nrows = p.n_samples
    else
        if constant_ploidy
            nrows = p.n_samples * p.max_ploidy
        else
            nrows = sum(p.ploidy)
        end
    end
    @assert length(data) == p.max_probs * nrows

    max_less_1 = p.max_probs - 1
    prob = 0.0
    factor = 1.0 / (2 ^ p.bit_depth - 1)
    # mask for depth not multiple of 8
    probs_mask = 0xFFFFFFFFFFFFFFFF >> (64 - p.bit_depth)
    bit_idx = 0

    if constant_ploidy && p.max_probs == 3 && p.bit_depth == 8
        # fast path for unphased, ploidy==2, 8 bits per prob.
        idx1 = startidx
        @inbounds for offset in 1:3:(3 * nrows)
            idx2 = 2 * ((offset-1) ÷ 3)
            first = d[idx1 + idx2]
            second = d[idx1 + idx2 + 1]
            data[offset] = lookup[first + 1]
            data[offset + 1] = lookup[second + 1]
            data[offset + 2] = lookup[256 - first - second]
        end
    else
        idx1 = startidx
        @inbounds for offset in 1:p.max_probs:(nrows * p.max_probs)
            # number of probabilities to be read per row
            if constant_ploidy
                n_probs = max_less_1
            elseif p.phased == 1
                n_probs = p.n_alleles - 1
            elseif p.ploidy[offset ÷ p.max_probs + 1] == 2 && p.n_alleles == 2
                n_probs = 2
            else
                n_probs = binomial(p.ploidy[offset ÷ p.max_probs + 1] +
                    p.n_alleles - 1, p.n_alleles - 1) - 1
            end

            remainder = 1.0
            @inbounds for i in 1:n_probs
                j = idx1 + bit_idx ÷ 8
                @inbounds prob = (unsafe_load_UInt64(d, j) >> (bit_idx % 8) &
                    probs_mask) * factor
                bit_idx += p.bit_depth
                remainder -= prob
                data[offset + i - 1] = prob
            end
            data[offset + n_probs] = remainder
            if n_probs + 1 < p.max_probs
                data[(offset + n_probs + 1):(offset + p.max_probs - 1)] .= NaN
            end
        end
    end
    for m in p.missings
        offset = p.max_probs * (m - 1) + 1
        data[offset:(offset + p.max_probs - 1)] .= NaN
    end
    return data
end

"""
    ref_dosage_fast!(data, p, d, idx, layout)
Dosage retrieval for 8-bit biallele case, no floating-point operations!
"""
const one_255th = 1.0f0 / 255.0f0
const mask_odd = reinterpret(Vec{16, UInt16}, Vec{32, UInt8}(
    tuple(repeat([0xff, 0x00], 16)...)))
const mask_even = reinterpret(Vec{16, UInt16}, Vec{32, UInt8}(
    tuple(repeat([0x00, 0xff], 16)...)))
function ref_dosage_fast!(data::Vector{T}, p::Preamble,
    d::Vector{UInt8}, startidx::Integer, layout::UInt8
    ) where {T <:AbstractFloat}
    @assert length(data) == p.n_samples
    @assert layout == 2
    @assert p.bit_depth == 8 && p.max_probs == 3 && p.max_ploidy == p.min_ploidy
    idx1 = startidx

    if p.n_samples >= 16
        @inbounds for n in 1:16:(p.n_samples - p.n_samples % 16)
            idx_base = idx1 + ((n-1) >> 1) << 2
            r = reinterpret(Vec{16, UInt16}, vload(Vec{32, UInt8}, d, idx_base))
            second = (r & mask_even) >> 8
            first  = (r & mask_odd) << 1
            dosage_level = first + second
            dosage_level_float = one_255th * convert(
                Vec{16, T}, dosage_level)
            vstore(dosage_level_float, data, n)
        end
    end
    rem = p.n_samples % 16
    if rem != 0
        @inbounds for n in ((p.n_samples - rem) + 1) : p.n_samples
            idx_base = idx1 + ((n - 1) << 1)
            data[n] = lookup[d[idx_base] * 2 +
                                d[idx_base + 1] + 1]
        end
    end
    return data
end

"""
    ref_dosage_slow!(data, p, d, idx, layout)
Dosage computation for general case.
"""
function ref_dosage_slow!(data::Vector{<:AbstractFloat}, p::Preamble,
    d::Vector{UInt8}, startidx::Integer, layout::UInt8
    )
    @assert length(data) == p.n_samples
    ploidy = p.max_ploidy
    half_ploidy = ploidy / 2

    maxval = 2 ^ p.bit_depth - 1
    factor = layout == 2 ? 1.0 / maxval : 1.0 / 32768
    probs_mask = 0xFFFFFFFFFFFFFFFF >> (64 - p.bit_depth)
    bit_idx = 0
    for n = 1:p.n_samples
        if p.max_ploidy != p.min_ploidy
            ploidy = ploidy[n]
            half_ploidy = ploidy ÷ 2
        end
        j = startidx + bit_idx ÷ 8
        hom = (unsafe_load_UInt64(d, j) >> (bit_idx % 8)) & probs_mask
        bit_idx += p.bit_depth
        j = startidx + bit_idx ÷ 8
        het = (unsafe_load_UInt64(d, j) >> (bit_idx % 8)) & probs_mask
        bit_idx += p.bit_depth
        data[n] = ((hom * ploidy) + (het * half_ploidy)) * factor
        if layout == 1
            # layout 1 also stores hom_alt probability, and it indicates missing by
            # triple zero.
            j = startidx + bit_idx ÷ 8
            hom_alt = (unsafe_load_UInt64(d, j) >> (bit_idx % 8)) & probs_mask
            bit_idx += p.bit_depth
            if hom == 0 && het == 0 && hom_alt == 0
                push!(p.missings, n)
            end
        end
    end
    return data
end

"""
    ref_dosage_phased!(data, p, d, idx, layout)
Dosage computation for phased genotypes.
"""
function ref_dosage_phased!(data::Vector{<:AbstractFloat}, p::Preamble,
    d::Vector{UInt8}, startidx::Integer, layout::UInt8
    )
    @assert length(data) == p.n_samples
    @assert layout == 2 "Phased genotypes not supported for Layout 1"
    ploidy = p.max_ploidy
    half_ploidy = ploidy / 2

    maxval = 2 ^ p.bit_depth - 1
    factor = 1.0 / maxval
    probs_mask = 0xFFFFFFFFFFFFFFFF >> (64 - p.bit_depth)
    bit_idx = 0
    for n = 1:p.n_samples
        if p.max_ploidy != p.min_ploidy
            ploidy = ploidy[n]
            half_ploidy = ploidy ÷ 2
        end
        ref_level = 0
        for _ in 1:ploidy
            j = startidx + bit_idx ÷ 8
            ref_level += (unsafe_load_UInt64(d, j) >> (bit_idx % 8)) & probs_mask
            bit_idx += p.bit_depth
        end
        data[n] = ref_level * factor
    end
    return data
end

"""
    alt_dosage(data, p)
Switch ref allele dosage `data` to alt allele dosage.
"""
function alt_dosage!(data::Vector{<:AbstractFloat}, p::Preamble)
    if p.n_samples >= 8
        @inbounds for n in 1:8:(p.n_samples - p.n_samples % 8)
            data[n]     = 2.0 - data[n]
            data[n + 1] = 2.0 - data[n + 1]
            data[n + 2] = 2.0 - data[n + 2]
            data[n + 3] = 2.0 - data[n + 3]
            data[n + 4] = 2.0 - data[n + 4]
            data[n + 5] = 2.0 - data[n + 5]
            data[n + 6] = 2.0 - data[n + 6]
            data[n + 7] = 2.0 - data[n + 7]
        end
    end
    @inbounds for n in (p.n_samples - p.n_samples % 8 + 1):p.n_samples
        data[n] = 2.0 - data[n]
    end
end

"""
    find_minor_allele(data, p)
Find minor allele index, returns 1 (ref) or 2 (alt)
"""
function find_minor_allele(data::Vector{<:AbstractFloat}, p::Preamble)
    batchsize = 100
    increment = max(p.n_samples ÷ batchsize, 1)
    total = 0.0
    freq = 0.0
    cnt = 0
    for idx2 in 1:increment
        for n in idx2:increment:p.n_samples
            cnt += 1
            total += data[n]
        end
        freq = total / (cnt * 2)
        @assert 0 <= freq <= 1
        if minor_certain(freq, batchsize * idx2, 5.0)
            break
        end
    end
    if freq <= 0.5
        return 1
    else
        return 2
    end
end

@inline function get_data_size(p::Preamble, layout::Integer)
    if layout == 1
        return p.n_samples * p.max_probs
    else
        constant_ploidy = p.max_ploidy == p.min_ploidy
        if p.phased == 0
            nrows = p.n_samples
        else
            if constant_ploidy
                nrows = p.n_samples * p.max_ploidy
            else
                nrows = sum(p.ploidy)
            end
        end
        return p.max_probs * nrows
    end
end

function _get_prob_matrix(d::Vector{T}, p::Preamble) where {T <: AbstractFloat}
    reshaped = reshape(d, p.max_probs, :)
    if p.phased == 1
        current = 1
        ragged = Matrix{T}(undef, p.max_ploidy * p.max_probs, p.n_samples)
        fill!(ragged, NaN)
        if p.max_ploidy == p.min_ploidy
            for i in 1:p.n_samples
                for j in 1:p.max_ploidy
                    first = (j-1) * p.max_probs + 1
                    last = j * p.max_probs
                    ragged[first:last, i] = reshaped[:, current]
                    current += 1
                end
            end
        else
            for (i, v) in enumerate(p.ploidy)
                for j in 1:v
                    first = (j-1) * p.max_probs + 1
                    last = j * p.max_probs
                    ragged[first:last, i] = reshaped[:, current]
                    current += 1
                end
            end
        end
        return ragged
    else
        return reshaped
    end
end

"""
    probabilities!(b::Bgen, v::Variant; T=Float32, clear_decompressed=false)
Given a `Bgen` struct and a `Variant`, compute probabilities.
The result is stored inside `v.genotypes.probs`, which can be cleared using
`clear!(v)`.

- T: type for the resutls
- `clear_decompressed`: clears decompressed byte string after execution if set `true`
"""
function probabilities!(b::Bgen, v::Variant;
        T=Float32, clear_decompressed=false, data=nothing, decompressed=nothing, is_decompressed=false)
    io, h = b.io, b.header
    if (decompressed !== nothing && !is_decompressed) ||
        (decompressed === nothing && (v.genotypes === nothing ||
        v.genotypes.decompressed === nothing))
        decompressed = decompress(io, v, h; decompressed=decompressed)
    else
        decompressed = v.genotypes.decompressed
    end
    if v.genotypes === nothing
        p = parse_preamble(decompressed, h, v)
        v.genotypes = Genotypes{T}(p, decompressed)
    else
        p = v.genotypes.preamble
    end
    startidx = 1
    if h.layout == 2
        startidx += 10 + h.n_samples
    end

    genotypes = v.genotypes
    data_size = get_data_size(p, h.layout)

    # skip parsing if already parsed
    if genotypes.probs !== nothing && length(genotypes.probs) >= data_size
        _get_prob_matrix(genotypes.probs, p)
    end

    if data !== nothing
        @assert length(data) == data_size
        genotypes.probs = data
    else
        genotypes.probs = Vector{T}(undef, data_size)
    end
    if h.layout == 1
        parse_layout1!(genotypes.probs, p, decompressed, startidx)
    elseif h.layout == 2
        parse_layout2!(genotypes.probs, p, decompressed, startidx)
    end
    if clear_decompressed
        clear_decompressed!(genotypes)
    end
    return _get_prob_matrix(genotypes.probs, p)
end

function ref_allele_dosage!(b::Bgen, v::Variant;
        T=Float32, mean_impute=false, clear_decompressed=false,
        data=nothing, decompressed=nothing, is_decompressed=false)
    io, h = b.io, b.header
    # just return it if already computed
    if v.genotypes !== nothing && v.genotypes.dose !== nothing
        genotypes = v.genotypes
        p = genotypes.preamble
        if genotypes.dose_mean_imputed && !mean_impute
            genotypes.dose[p.missings] .= NaN
        elseif !genotypes.dose_mean_imputed && mean_impute
            genotypes.dose[p.missings] .= mean(filter(!isnan, genotypes.dose))
        end
        if genotypes.minor_allele_dosage && genotypes.minor_idx != 1
            alt_dosage!(genotypes.dose, p)
        end
        return v.genotypes.dose
    end
    if (decompressed !== nothing && !is_decompressed) ||
        (decompressed === nothing && (v.genotypes === nothing ||
        v.genotypes.decompressed === nothing))
        decompressed = decompress(io, v, h; decompressed=decompressed)
    else
        decompressed = v.genotypes.decompressed
    end
    startidx = 1
    if v.genotypes === nothing
        p = parse_preamble(decompressed, h, v)
        v.genotypes = Genotypes{T}(p, decompressed)
    else
        p = v.genotypes.preamble
    end
    if h.layout == 2
        startidx += 10 + h.n_samples
    end

    @assert p.n_alleles == 2 "allele dosages are available for biallelic variants"
    #@assert p.phased == 0

    genotypes = v.genotypes
    if data !== nothing
        @assert length(data) == h.n_samples
    else
        genotypes.dose = Vector{T}(undef, h.n_samples)
        data = genotypes.dose
    end

    if p.phased == 0
        if p.max_ploidy == p.min_ploidy && p.max_probs == 3 && p.bit_depth == 8 &&
                b.header.layout == 2
            ref_dosage_fast!(data, p, decompressed, startidx, h.layout)
        else
            ref_dosage_slow!(data, p, decompressed, startidx, h.layout)
        end
    else # phased
        ref_dosage_phased!(data, p, decompressed, startidx, h.layout)
    end
    genotypes.minor_idx = find_minor_allele(data, p)
    data[p.missings] .= NaN
    if mean_impute
        data[p.missings] .= mean(filter(!isnan, data))
    end
    if clear_decompressed
        clear_decompressed!(genotypes)
    end
    return data
end

"""
    minor_allele_dosage!(b::Bgen, v::Variant; T=Float32,
    mean_impute=false, clear_decompressed=false)
Given a `Bgen` struct and a `Variant`, compute minor allele dosage.
The result is stored inside `v.genotypes.dose`, which can be cleared using
`clear!(v)`.

- `T`: type for the results
- `mean_impute`: impute missing values with the mean of nonmissing values
- `clear_decompressed`: clears decompressed byte string after execution if set `true`
"""
function minor_allele_dosage!(b::Bgen, v::Variant;
        T=Float32, mean_impute=false, clear_decompressed=false,
        data=nothing, decompressed=nothing, is_decompressed=false)
    # just return it if already computed
    io, h = b.io, b.header
    if v.genotypes !== nothing && v.genotypes.dose !== nothing
        genotypes = v.genotypes
        p = genotypes.preamble
        if genotypes.dose_mean_imputed && !mean_impute
            genotypes.dose[p.missings] .= NaN
        elseif !genotypes.dose_mean_imputed && mean_impute
            genotypes.dose[p.missings] .= mean(filter(!isnan, genotypes.dose))
        end
        if !genotypes.minor_allele_dosage && genotypes.minor_idx != 1
            alt_dosage!(genotypes.dose, p)
        end
        return v.genotypes.dose
    end
    ref_allele_dosage!(b, v; T=T, mean_impute=mean_impute,
        clear_decompressed=clear_decompressed, data=data, decompressed=decompressed,
        is_decompressed=is_decompressed)
    genotypes = v.genotypes
    if data === nothing
        data = genotypes.dose
    end
    if genotypes.minor_idx != 1
        alt_dosage!(data, genotypes.preamble)
    end
    return data
end

"""
    clear!(g::Genotypes)
    clear!(v::Variant)
Clears cached decompressed byte representation, probabilities, and dose.
If `Variant` is given, it removes the corresponding `.genotypes` altogether.
"""
function clear!(g::Genotypes)
    g.decompressed = nothing
    g.probs = nothing
    g.dose = nothing
    return
end

"""
    clear_decompressed!(g::Genotypes)
Clears cached decompressed byte representation.
"""
function clear_decompressed!(g::Genotypes)
    g.decompressed = nothing
    return
end


text:
function menu_export_dialog(data)
    fname = Gtk.save_dialog("Export as...",
                            Gtk.Null(),
                            (Gtk.GtkFileFilter("*.png, *.jpg",
                                               name="All supported formats"),
                             "*.png",
                             "*.jpg"))
    if(length(fname) < 4)
        return
    end

    println(fname)

    if(fname[end-2:end] == "png")
        println(fname)
        FileIO.save(fname, data)
    end
end

function menu_signals!(signals)
    menu_open_impl = map(signals[:file_open], init=false) do sig
        if(sig)
            Gtk.save_dialog(
                "Load an image or video",
                Gtk.GtkNullContainer(),
                ("*.png",
                 "*.jpg",
                 Gtk.GtkFileFilter("*.png, *.jpg",
                                   name="All supported formats")))
        end
        ""
    end
    signals[:file_open] = filter(fname-> fname != "", menu_open_impl)
end

function add_item!(menu,
                   signals::Dict,
                   menu_name::String,
                   signal_name::Symbol)
    signal = GtkReactive.Signal(false)
    item   = Gtk.MenuItem(menu_name)
    push!(menu, item)
    Gtk.signal_connect(item, :activate) do widget
        push!(signal, true)
    end
    signals[signal_name] = signal
end

function create_menubar(window)
    signals = Dict{Symbol, GtkReactive.Signal}()

    file     = Gtk.MenuItem("_File")
    filemenu = Gtk.Menu(file)
    add_item!(filemenu, signals, "Open",                    :file_open)
    add_item!(filemenu, signals, "Export current image",    :file_export_image)
    add_item!(filemenu, signals, "Export best image",       :file_export_best)
    add_item!(filemenu, signals, "Export best parameter",   :file_export_params)
    add_item!(filemenu, signals, "Export Gaussian process", :file_export_gp)

    action     = Gtk.MenuItem("_Action")
    actionmenu = Gtk.Menu(action)
    add_item!(actionmenu, signals, "Start",                :action_start)
    #add_item!(actionmenu, signals, "Show Best",            :action_best)
    add_item!(actionmenu, signals, "Restart", :action_restart)

    menubar = Gtk.MenuBar()
    push!(menubar, file)
    push!(menubar, action)

    menu_signals!(signals)

    menubar, signals
end


text:
"""
Table Schema field
https://github.com/frictionlessdata/TableSchema.jl#field
"""
mutable struct Field
    descriptor::Dict
    name::String
    typed::String
    format::String
    constraints::Constraints
    # required::Bool

    function Field(d::Dict)
        name = haskey(d, "name") ? d["name"] : ""
        typed = haskey(d, "type") ? d["type"] : DEFAULT_TYPE
        format = haskey(d, "format") ? d["format"] : DEFAULT_FORMAT
        constraints = haskey(d, "constraints") ?
            Constraints(d["constraints"]) : Constraints()
        # required = cons.required
        new(d, name, typed, format, constraints)
    end

    Field(name::String) = Field(Dict( "name" => name ))
end

OPTION_KEYS = ["decimalChar", "groupChar", "bareNumber", "trueValues", "falseValues"]

_TRUE_VALUES = ["true", "True", "TRUE", "1"]
_FALSE_VALUES = ["false", "False", "FALSE", "0"]
_DEFAULT_BARE_NUMBER = true
_DEFAULT_GROUP_CHAR = ""
_DEFAULT_DECIMAL_CHAR = "."

function cast_by_type(value, typed::String, format::String, options::Dict)
    if typed == "any"
        return value

    elseif typed == "boolean"
        isa(value, Bool) && return value
        if !isa(value, String); return CastError(); end
        value = strip(value)
        value in _TRUE_VALUES && return true
        value in _FALSE_VALUES && return false

    elseif typed == "integer"
        isa(value, Integer) && return value
        if !isa(value, String); return CastError(); end
        if !get(options, "bareNumber", _DEFAULT_BARE_NUMBER)
            value = replace(value, r"((^\D*)|(\D*$))" => "")
        end
        try; return parse(Int64, value); catch; return CastError(); end

    elseif typed == "number"
        isa(value, AbstractFloat) && return value
        isa(value, Integer) && return Float64(value)
        if !isa(value, String); return CastError(); end
        group_char = get(options, "groupChar", _DEFAULT_GROUP_CHAR)
        decimal_char = get(options, "decimalChar", _DEFAULT_DECIMAL_CHAR)
        value = replace(value, r"\s" => "")
        value = replace(value, decimal_char => ".")
        value = replace(value, group_char => "")
        if !get(options, "bareNumber", _DEFAULT_BARE_NUMBER)
            value = replace(value, r"((^\D*)|(\D*$))" => "")
        end
        try; return parse(Float64, value); catch; return CastError(); end

    elseif typed == "object"
        typeof(value) == Dict && return value
        if !isa(value, String); return CastError(); end
        try; value = JSON.parse(value); catch; return CastError(); end
        isa(value, Dict) && return value

    elseif typed == "string"
        isa(value, AbstractString) && return value
        isa(value, Number) && return repr(value)
        try; return parse(String, value); catch; return CastError(); end

    else
        throw(ErrorException("Cast by type $typed not implemented"))

    end
    CastError()
end

function build(f::Field)
    d = Dict()
    d["name"] = f.name
    d["type"] = f.typed
    d["format"] = f.format
    d["constraints"] = build(f.constraints)
    f.descriptor = d
    d
end

function cast_value(f::Field, value, constraints=true)

    # TODO: Ignore Missing values
    cast_value = value
    if !isempty(value)

        # Collect options
        d = f.descriptor
        options = Dict(k => d[k] for k in OPTION_KEYS if haskey(d, k) && !isempty(d[k]))
        cast_value = cast_by_type(value, f.typed, f.format, options)
        if isa(cast_value, CastError)
            throw(CastError(
                "Field $(f.name) cannot cast value \"$(value)\" for type $(f.typed) with format $(f.format)"
            ))
        end

    end

    if constraints
        throw(ErrorException("Not implemented"))
    end

    cast_value
end

# test_value = throw(ErrorException("Not implemented"))

function validate(f::Field)
    isempty(f.descriptor) &&
        throw(FieldError("Missing Descriptor"))
    f.name == "" &&
        throw(FieldError("Name is empty"))
    f.typed == "" &&
        throw(FieldError("Type is empty"))
    return true
end

max_stars_repo_name: loleg/tableschema-jl

max_issues_repo_name: loleg/tableschema-jl

max_forks_repo_name: loleg/tableschema-jl


text:
##=================####==============##
## TEST SIMULATION ABM NETWORK-EPIDEMICS
##=================####==============##

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

using ArgParse      # For CLI arguments
using StatsBase     # For sample function
using SparseArrays  # For Adyacency Matrix
using Distributions # For probability distributions
using DelimitedFiles # For data output
using LightGraphs   # For Netwrork stuff

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

include("NetABM.jl") # Include file with functions and definitions

### ============== ### ============== ### ============== ###
### SYSTEM'S PARAMETERS FROM COMMAND LINE
### ============== ### ============== ### ============== ###

args_settings = ArgParseSettings()

@add_arg_table! args_settings begin
    "--numAgents", "-n"
        help = "Number of Agents"
        arg_type = Int
        default = 65536
    "--initialI", "-i"
        help = "Initial fraction of Infected Agents"
        arg_type = Float64
        default = 0.01
    "--attackR", "-a"
        help = "Attack Rate"
        arg_type = Float64
        default = 0.25
    "--meanRecT", "-t"
        help = "Mean recovery time"
        arg_type = Int64
        default = 10
    "--meanCop", "-c"
        help = "Mean cooperation probability"
        arg_type = Float64
        default = 0.5
    "--stdCop", "-s"
        help = "Std Deviation cooperation probability"
        arg_type = Float64
        default = 0.1
    "--rep", "-r"
        help = "Repetition ID (for ensemble statistics)"
        arg_type = Int
        default = 1
end

args = parse_args(args_settings)

##=================####==============##
## ASIGN SYSTEM'S PARAMETERS FROM CLI

params = Params()

params.num_agents      = args["numAgents"] # Number of agents (32768)
params.p_infected_t0   = args["initialI"] # Initial fraction of infected agents
params.μ_recovT_agents = args["meanRecT"] # Std deviation of cooperation probability
params.μ_cop_agents    = args["meanCop"] # Mean of cooperation probability
params.σ_cop_agents    = args["stdCop"] # Std deviation of cooperation probability
params.attack_rate     = args["attackR"] # Virus attack rate
params.repetition      = args["rep"] # Repetition ID for ensemble statistics

##=================####==============##
# TEST PARAMETERS

# params.num_agents    = 2^16 # Number of agents (32768)
# params.p_infected_t0 = 0.01 # Initial fraction of infected agents
# params.attack_rate   = 0.25 # Virus attack rate
# params.repetition    = 1 # Repetition ID for ensemble statistics

# params.μ_cop_agents  = 0.15 # Mean of cooperation probability
# params.σ_cop_agents  = 0.1 # Std deviation of cooperation probability
# params.μ_recovT_agents = 7 # Std deviation of cooperation probability

##=================####==============##
# PATH TO SAVE DATA (MODIFY WITH YOUR OWN!)

repo_path = joinpath(homedir(), "Code", "Collaborations", "networkEpidemiology_MX")

data_out_path = joinpath(homedir(),
    "Code", "Collaborations", "networkEpidemiology_MX", "julia_ABM", "data", "LFR")

##=================####==============##
# INITIALIZE ADJACENCY MATRIX (LFR Model)

# the_root = joinpath(repo_path, "src", "binary_networks")
#
# println("BUILDING NETWORK...")
#
# the_net, true_com = get_network(
#     the_root,
#     "julia_ABM/sims_scripts",
#     N    = params.num_agents,
#     t2   = 0.25,
#     minc = 100,
#     maxc = 1000);

the_net = readdlm(joinpath(repo_path, "julia_ABM", "sims_scripts", "network.dat"))

g, Nodes = lectura_uw(the_net)

adj_mat = adjacency_matrix(g)
rows = rowvals(adj_mat)

println("NETWORK DONE")
##=================####==============##
# INITIALIZE MEETS, COOPERATION AND RECOVERY TIMES DISTRIBUTION

rt_μ = params.μ_recovT_agents # MEAN RECOVERY TIME
pc_μ = params.μ_cop_agents # MEAN COOPERATION PROBABILITY
pc_σ = params.σ_cop_agents # STD_DEV COOPERATION PROBABILITY

# meets_dist = Exponential()
meets_dist  = Geometric()
coop_dist   = truncated(Normal(pc_μ, pc_σ), 0.0, 1.0)
recovt_dist = truncated(Poisson(rt_μ), 7, 35)

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

# DEFINE AGENTS POPULATION (DEFAULT VALUES)
agents = [Agent(id) for id in 1:params.num_agents]

# MODIFY AGENT'S CHARACTERISTICS ACCORDING TO ITS DEMOGRAPHICS
initialize_demographics!(agents, params, coop_dist, meets_dist, recovt_dist)

# ASSGINS CONTACTS TO agent FROM adj_mat
for ag in agents
    assign_contacts!(ag, agents, adj_mat, rows)
end

##=================####==============##
# SYSTEM EVOLUTION

# TIME AT WHICH LOCKDOWN IS ANNOUNCED
time_lockdown = 15

tau = []
num_S = []
num_I = []
num_R = []

# INITIALIZE TIME AND POPULATIONS
params.now_t = 0
populations = get_populations(agents, params)

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

println("SIM STARTS")

while populations[2] > 0.0 # WHILE THERE ARE INFECTED AGENTS

    # RECORD DATA
    push!(tau, params.now_t)
    push!(num_S, populations[1])
    push!(num_I, populations[2])
    push!(num_R, populations[3])

    # EACH AGENTS DECIDES TO STAY AT HOME/ NOT TO STAY AT HOME
    if params.now_t >= time_lockdown
        set_coop_agents!(agents, params)
    end

    # FINDS NEXT STATE AND UPDATES ALL AGENTS
    update_all_agents!(agents, params)

    # GET CURRENT POPULATIONS AND UPDATE TIME
    global populations = get_populations(agents, params);
    params.now_t += 1
    println(params.now_t, populations)
end

println("NO MORE INFECTED!")

##=================####==============##
# EXPORT PARAMETERS
export_parameters(params, data_out_path)

##=================####==============##
# EXPORT POPULATIONS DATA

params_filename =
    "_I0_$(params.p_infected_t0)_AR_$(params.attack_rate)_CA_$(params.μ_cop_agents)_$(params.repetition).csv"

file_name = "populations_N_$(params.num_agents)"*params_filename

open(joinpath(data_out_path, file_name), "w") do io
    println(io, "t,S,I,R")
    writedlm(io, hcat(tau, num_S, num_I, num_R), ',')
end

##=================####==============##
# EXPORT INFECTION TIME OF AGENTS

file_name = "infTime_N_$(params.num_agents)"*params_filename

open(joinpath(data_out_path, file_name), "w") do io
    println(io, "agent_ID,time_of_infection")
    for ag in agents
        println(io, "$(ag.id),$(ag.infection_t)")
    end
end

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

println("TAN TAN...")

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


text:
#  Due to @andreasnoack in PR #7 (modified)
macro atleastversion(x)
    if eval(x) >= v"0.4-"
        return :(qrfact(Q*X, Val{true}))
    else
        return :(qrfact(Q*X, pivot = true))
    end
end

# A Frisch-Newton algorithm as described in http://projecteuclid.org/euclid.ss/1030037960
# qreg_ip_coef is a modified lp_fnm of Daniel Morillo & Roger Koenker
# Translated from Ox to Matlab by Paul Eilers 1999
# Modified by Roger Koenker 2000
# More changes by Paul Eilers 2004
# QuantileRegression.jl developers 2015
#
# The functions in this file is released under the BSD-3 license, although the original
# versions were not. Permissing was given to Patrick Kofod Mogensen by Roger Koenker and
# Paul Eilers on April 24 2015 via E-mail. Original mail correspondences can be provided
# upon request.

function bound(x, dx)
# Fill vector with allowed step lengths
# Replace with -x/dx for negative dx
    b = 1e20 + 0.0 * x
    for i = 1:length(dx)
        if dx[i] < 0.0
            @inbounds b[i] = -x[i] / dx[i]
        end
    end
    return b
end

function qreg_coef(Y::Vector, X::Matrix, p, s::IP)
# input: X is an n x k matrix of exogenous regressors,
#        Y is an n x 1 vector of outcome variables
#        p \in (0,1) is the quantile of interest
# Output: p^{th} regression quantiles.
# Construct the dual problem of quantile regression
c = -Y

m = size(X, 1)
u = ones(m)
x = (1 - p) .* u
b = X'x

# Solve a linear program by the interior point method:
# min(c * u), s.t. A * x = b and 0 < x < u
# An initial feasible solution has to be provided as x

# Set some constants
beta = 0.99995
small = 1e-6
max_it = 50
n, m = size(X)

# Generate inital feasible point
s = u - x
y = -X\Y
r = c - X*y
BLAS.axpy!(0.001, (r .== 0.0).*1.0, r)
z = r .* (r .> 0.0)
w = z - r
gap = Base.LinAlg.BLAS.dot(c, x) - Base.LinAlg.BLAS.dot(y, b) + Base.LinAlg.BLAS.dot(w, u)

# Start iterations
it = 0
for it = 1:max_it
    #   Compute affine step
    q = 1 ./ (z ./ x + w ./ s)
    r = z - w
    Q = Diagonal(sqrt.(q)) # Very efficient to do since Q diagonal
    # AQtF = @atleastversion(VERSION)
    AQtF = qrfact(Q*X, Val{true}) # PE 2004
    rhs = Q*r        # "
    dy = AQtF\rhs   # "
    dx = q.*(X*dy - r)
    ds = -dx
    dz = -z .* (1 + dx ./ x)
    dw = -w .* (1 + ds ./ s)

    # Compute maximum allowable step lengths
    fx = bound(x, dx)
    fs = bound(s, ds)
    fw = bound(w, dw)
    fz = bound(z, dz)
    fpv = min.(fx, fs)
    fdv = min.(fw, fz)
    fp = min.(minimum(beta * fpv), 1)
    fd = min.(minimum(beta * fdv), 1)

    # If full step is feasible, take it. Otherwise modify it
    if min.(fp, fd) < 1.0

        # Update mu
        mu = Base.LinAlg.BLAS.dot(z, x) + Base.LinAlg.BLAS.dot(w, s)
        g = Base.LinAlg.BLAS.dot(z + fd*dz, x + fp*dx) + Base.LinAlg.BLAS.dot(w + fd*dw, s + fp*ds)
        mu = mu * (g / mu)^3 / (2 * n)

        # Compute modified step
        dxdz = dx .* dz
        dsdw = ds .* dw
        xinv = 1 ./ x
        sinv = 1 ./ s
        xi = mu .* (xinv - sinv)
        #rhs = rhs + Q * (dxdz - dsdw - xi)
        BLAS.axpy!(1.0, Q * (dxdz - dsdw - xi), rhs) # no gemv-wrapper gemv(Q, (dxdz - dsdw - xi), rhs,1,1,n)?

        dy = AQtF\rhs
        dx = q .* (X*dy + xi - r - dxdz + dsdw)
        ds = -dx
        for i = 1:length(dz)
            dz[i] = mu * xinv[i] - z[i] - xinv[i] * z[i] * dx[i] - dxdz[i]
            dw[i] = mu * sinv[i] - w[i] - sinv[i] * w[i] * ds[i] - dsdw[i]
        end

        # Compute maximum allowable step lengths
        fx = bound(x, dx)
        fs = bound(s, ds)
        fw = bound(w, dw)
        fz = bound(z, dz)
        fp = min.(fx, fs)
        fd = min.(fw, fz)
        fp = min.(minimum(beta .* fp), 1)
        fd = min.(minimum(beta .* fd), 1)

    end

    # Take the steps
    BLAS.axpy!(fp, dx, x)
    BLAS.axpy!(fp, ds, s)
    BLAS.axpy!(fd, dy, y)
    BLAS.axpy!(fd, dw, w)
    BLAS.axpy!(fd, dz, z)

    gap = Base.LinAlg.BLAS.dot(c, x) - Base.LinAlg.BLAS.dot(y, b) + Base.LinAlg.BLAS.dot(w, u)

    if gap < small
        break
    end
end

return -y
end


text:
# Computes the M_{p,q} coefficients defined in (D37) for
# the r0 >> 1 limit in which k2 << 1 and r0-1 < b < r0+1
# (this is a huge occultor).

using PyPlot

function eps3_of(k2::T,Kofk::T,Eofk::T) where {T <: Real}
eps3 = Kofk-Eofk-0.5*k2*Kofk
return eps3
end

function eps3_series(k2::Real)
# Series evaluation of eps3 for k2 < 1:
@assert(k2 < 1)
tol = eps(k2)
eps3 = zero(k2)
# Compute n=1 term
term = pi/32
eps3 += term
n=2
while abs(term) > tol*abs(eps3)
  term *= (2*n-1)^2*k2/4/(n^2-1)
  eps3 += term
  n +=1
end
return eps3*k2*k2
end

function mpq_init(k2::Real)
# Initialize the M_{p,q} computation:
k4 = k2*k2
m001 = (8.-16.*k2)/3.  # eps3 term
m002 = 4k4/3  # K(m) term
m021 = (8.-28k2-12k4)/15.
m022 = (22.-6k2)*k4/15
m201 = (32.-52k2+12k4)/15.
m202 = (-2.+6k2)*k4/15
m221 = (32.-76k2+36k4-24k4*k2)/105.
m222 = (-2.+30k2-12k4)*k4/105.
return m001,m002,m021,m022,m201,m202,m221,m222
end

function mpq_of(k2::T,pmax::Int64,qmax::Int64,Kofk::T,Eofk::T,alt::Bool) where {T <: Real}
mpq = zeros(typeof(k2),pmax+1,qmax+1)
if k2 < 1 && alt
  mpq1 = zeros(typeof(k2),pmax+1,qmax+1)
  mpq2 = zeros(typeof(k2),pmax+1,qmax+1)
  m001 = 0.; m021=0.; m201=0.; m221=0.
  m002 = 0.; m022=0.; m202=0.; m222=0.
  eps3 = zero(typeof(k2))
# Alternate form:
  if k2 < 0.1
    eps3 = eps3_series(k2)
  else
    eps3=eps3_of(k2,Kofk,Eofk)
  end
  m001,m002,m021,m022,m201,m202,m221,m222 = mpq_init(k2)
  mpq[1,1]=m001*eps3+m002*Kofk; mpq[1,3]=m021*eps3+m022*Kofk
  mpq[3,1]=m201*eps3+m202*Kofk; mpq[3,3]=m221*eps3+m222*Kofk
#  println("m221: ",m221," m222: ",m222," eps3: ",eps3," K(m): ",Kofk)
else
  k4  = k2^2
  if k2 < 1
    eps1 = (1-k2)*Kofk; eps2 = Eofk
  elseif k2 > 1
  # In this case, no special handling is required as there
  # is not cancellation as there is in k2 < 1 case.
  # However, for k2 >> 1, are the recursion relations stable?
    k2inv = 1.0/k2; k = sqrt(k2)
    eps1 = (1.0-k2)/k*Kofk; eps2 = k*Eofk+eps1
  else
    eps1 = zero(k2); eps2= one(k2)
  end
  mpq[1,1] = (8-12k2)*eps1/3+(-8+16*k2)*eps2/3
  mpq[1,3] = (8-24k2)*eps1/15+(-8+28*k2+12*k4)*eps2/15
  mpq[3,1] = (32-36k2)*eps1/15+(-32+52*k2-12k4)*eps2/15
  mpq[3,3] = (32-60k2+12k4)*eps1/105+(-32+76k2-36k4+24k4*k2)*eps2/105
end
#println("p,q: ",0," ",0," M_{pq}: ",mpq[1,1])
#println("p,q: ",0," ",2," M_{pq}: ",mpq[1,3])
#println("p,q: ",2," ",0," M_{pq}: ",mpq[3,1])
#println("p,q: ",2," ",2," M_{pq}: ",mpq[3,3])

# Now, use recursion relation to compute M_{p,q}:
for p=0:2:pmax
  if p > 2
    for q = 0:2:2
      d3 = 2p+q-(p+q-2)*(1.-k2)
      d4 = (3-p)*k2
      mpq[p+1,q+1] = (d3*mpq[p-1,q+1]+d4*mpq[p-3,q+1])/(p+q+3)
#      println("p,q: ",p," ",q," M_{pq}: ",mpq[p+1,q+1])
    end
  end
  for q=4:2:qmax
    d1 = q+2+(p+q-2)*(1.-k2)
    d2 = (3-q)*(1.-k2)
    mpq[p+1,q+1]=(d1*mpq[p+1,q-1]+d2*mpq[p+1,q-3])/(p+q+3)
#    println("p,q: ",p," ",q," M_{pq}: ",mpq[p+1,q+1])
  end
end
return mpq
end

max_stars_repo_name: tagordon/limbdark

max_issues_repo_name: tagordon/limbdark

max_forks_repo_name: tagordon/limbdark


text:
using Test, Statistics, Distributions, Distances, DensityRatioEstimation
using DensityRatioEstimation: pairwise_sqd, gaussian_gram, gaussian_gram_by_pairwise_sqd
using JuMP, Ipopt

@testset "Correctness of `pairwise_sqd` and `gaussian_gram_by_pairwise_sqd`" begin
    pairwise_sqd_golden(x) = pairwise(SqEuclidean(), x; dims=2)
    pairwise_sqd_golden(x, y) = pairwise(SqEuclidean(), x, y; dims=2)

    xtest = [
        1.0 2.0 4.0;
        1.0 2.0 4.0
    ]

    ytest = [
        0.0 2.0 18.0;
        2.0 0.0  8.0;
        18.0 8.0  0.0
    ]

    @testset "`pairwise_sqd(x)`" begin
        n_randtests = 10
        @test pairwise_sqd(xtest) == ytest
        for _ = 1:n_randtests
            xrand = randn(784, 100)
            @test pairwise_sqd(xrand) ≈ pairwise_sqd_golden(xrand)
        end
    end

    @testset "`pairwise_sqd(x, y)`" begin
        n_randtests = 10
        @test pairwise_sqd(xtest, xtest) == ytest
        for _ = 1:n_randtests
            xrand = randn(784, 100)
            yrand = randn(784, 200)
            @test pairwise_sqd(xrand, yrand) ≈ pairwise_sqd_golden(xrand, yrand)
        end
    end

    @testset "`gaussian_gram`" begin
        for σ in sqrt.([1, 2, 4, 8, 16])
            K1 = gaussian_gram(xtest, σ)
            K2 = gaussian_gram(xtest, xtest, σ)
            K3 = gaussian_gram_by_pairwise_sqd(ytest, σ)
            @test K1 == K2 == K3 == exp.(-ytest / (2 * σ ^ 2))
        end
    end
end

@testset "MMDNumerical`" begin
    @testset "Positivity and normalisation constraints" begin
        dist_nu = Normal(1, 1)
        dist_de = Normal(0, 2)

        n_nu = 100
        x_nu = rand(dist_nu, 1, n_nu)
        n_de = 200
        x_de = rand(dist_de, 1, n_de)

        r_numerical = estimate_ratio(MMDNumerical(), x_de, x_nu)

        @test mean(r_numerical) ≈ 1
        @test all(r_numerical .> 0)
    end
end

@testset "MMDAnalytical`" begin
    @testset "Consistency between solve and inv methods" begin
        dist_nu = Normal(1, 1)
        dist_de = Normal(0, 2)

        n_nu = 100
        x_nu = rand(dist_nu, 1, n_nu)
        n_de = 200
        x_de = rand(dist_de, 1, n_de)

        r_solve = estimate_ratio(MMDAnalytical(method=:solve), x_de, x_nu)
        r_inv = estimate_ratio(MMDAnalytical(method=:inv), x_de, x_nu)

        @test r_solve ≈ r_inv

        r_solve = estimate_ratio(MMDAnalytical(0), x_de, x_nu)
        r_inv = estimate_ratio(MMDAnalytical(0), x_de, x_nu)

        @test r_solve ≈ r_inv
    end
end


text:
using HTTP

URL = "https://adventofcode.com/"

function get_aoc_input(day::Int, year::Int=2018)
	url = URL * string(year) * "/day/" * string(day) * "/input"
	session = ENV["AOC_SESSION"]
	header = ("Cookie", "session=" * session)
	response = HTTP.get(url, headers=[header])
	return String(response.body)
end

export get_aoc_input


text:
using Printf, Serialization, BSON
function default_cb(step, load_time, calculation_time, training_loss)
	println(@sprintf("%d: training loss = %.4g load time= %.0fs (%.2fs per step)  compute time %.0fs (%.2fs per step) ",step , training_loss, load_time, load_time / step, calculation_time, calculation_time / step))
end

function wronggrad(model)
	ps = params(model)
	a = any(map(p -> any(isnan.(Flux.Tracker.grad(p))), ps))
	a && @warn("nan in a gradient")
	b = any(map(p -> any(isinf.(Flux.Tracker.grad(p))), ps))
	b && @warn("inf in a gradient")
	a || b
end
correctgrad(model) = !wronggrad(model)

"""
		serializestate(filename, model, loss, data)

		serialize the model and data to filename. Used mainly for debugging nans

"""
function serializestate(filename, model, loss, data)
	open(filename, "w") do fid
		serialize(fid, (cpu(model), loss, cpu(data)))
	end
end

"""
		deserializestate(filename)

		serialize the model and data to filename. Used mainly for debugging nans

"""
function deserializestate(filename)
	open(filename, "r") do fid
		deserialize(fid)
	end
end

"""
	meminfo()

	runs the garbage collector and output content of `/proc/meminfo` to track the memory usage. Use only on linux

"""
function meminfo()
	gc();gc();gc();gc();gc();gc();
	run(pipeline(`/bin/cat /proc/meminfo`, `/bin/grep MemFree`))
end

"""
		learn(model, loss, opt, data_provider, max_steps;cb = default_cb, breaks = 100, state_filename = nothing)

		learns `model` with `loss` function using `opt` optimalization algorithm

		`model` is a model accepting data provided by data provider, which can be as simple as

"""
function learn(model, loss, opt, data_provider, max_steps; cb = default_cb, breaks = 100, state_filename = nothing)
	step, load_time, calculation_time, training_loss =0, 0.0, 0.0, 0.0
	state = nothing
	while step < max_steps
	  # load_time += @elapsed  data = state == nothing ?  Base.iterate(data_provider) : Base.iterate(data_provider, state)
	  load_time += @elapsed  data = data_provider()
	  data == nothing && break
	  # data, state = data
	  calculation_time += @elapsed training_loss +=  ∇loss!(loss, model, data, state_filename)
	  calculation_time += @elapsed opt()
	  step+=1
	  if mod(step,breaks) == 0
	  	cb(step, load_time, calculation_time, training_loss / breaks)
	  	serializestate("state.ser", model, loss, data)
	  	GC.gc();GC.gc();GC.gc();GC.gc();GC.gc()
	    training_loss = 0.0
	  end
	end
end

function learn(model, loss, opt::Tuple, data_provider, max_steps; cb = default_cb, breaks = 100, state_filename = nothing)
	step, load_time, calculation_time, training_loss =0, 0.0, 0.0, 0.0
	state = nothing
	while step < max_steps
	  # load_time += @elapsed  data = state == nothing ?  Base.iterate(data_provider) : Base.iterate(data_provider, state)
	  load_time += @elapsed  data = data_provider()
	  data == nothing && break
	  # data, state = data
	  calculation_time += @elapsed training_loss +=  ∇loss!(loss, model, data, state_filename)
	  calculation_time += @elapsed Flux.Optimise._update_params!(opt...)
	  step+=1
	  if mod(step,breaks) == 0
	  	cb(step, load_time, calculation_time, training_loss / breaks)
	  	BSON.@save state_filename model loss data
	  	GC.gc();GC.gc();GC.gc();GC.gc();GC.gc()
	    training_loss = 0.0
	  end
	end
end

function ∇loss!(loss, model, data, filename = nothing)
	fVal = loss(model, data)
	if filename != nothing && (isnan(Flux.data(fVal)) || isinf(Flux.data(fVal)))
		serializestate(filename, model, loss, data)
		error("nan or inf in function evaluation")
  end
	Flux.back!(fVal)
  if filename != nothing && wronggrad(model)
  	serializestate(filename, model, loss, data)
  	error("nan or inf in the functionloc(f::Function, types)ction gradient")
  end
  Flux.data(fVal)
end

∇loss!(loss, models::NTuple{N, A}, datas::NTuple{N, B}, filename)  where {N, A, B} = ∇loss!(loss, models, datas,  map(params, models),  filename)
function ∇loss!(loss, models::NTuple{N, A}, datas::NTuple{N, B}, mparams, filename) where {N, A, B}
	@assert length(models) == length(datas)
	n = length(models)
	fVals = fill(0.0, n)
	Threads.@threads for i in 1:n
		i != 1 && copyvalues!(mparams[i], mparams[1])
		zerograds!(mparams[i])
		fVals[i] = ∇loss!(loss, models[i], datas[i], filename)
	end
	for i in 2:n
		addgrads!(mparams[1], mparams[i])
	end
	scalegrads!(mparams[1], 1/length(models))
	mean(fVals)
end

max_stars_repo_name: pevnak/FluxExtensions

max_issues_repo_name: pevnak/FluxExtensions

max_forks_repo_name: pevnak/FluxExtensions


text:
prob, opts = Problems.DoubleIntegrator()
solver = Altro.iLQRSolver(prob, verbose=true, cost_tolerance=1e-1)
@test solver.opts.verbose == true
@test solver.opts.cost_tolerance == 1e-1

solver = Altro.AugmentedLagrangianSolver(prob, verbose=true, cost_tolerance=1e-1, square_root=true)
@test solver.opts.verbose == true
@test solver.solver_uncon.opts.verbose == true
@test solver.opts.cost_tolerance == 1e-1
@test solver.solver_uncon.opts.cost_tolerance == 1e-1
@test solver.solver_uncon.opts.square_root == true
solver = Altro.AugmentedLagrangianSolver(prob, verbose=1, cost_tolerance=1e-1, square_root=true)
@test solver.opts.verbose == true
@test solver.solver_uncon.opts.verbose == 1
Altro.set_options!(solver, verbose=0)
@test solver.opts.verbose == false
@test solver.solver_uncon.opts.verbose == false

# Test ALTRO solver options
solver = ALTROSolver(prob, cost_tolerance=1e-1, square_root=true,
    iterations=210, iterations_outer=60, ρ_chol=1e-1, constraint_tolerance=1e-2)
@test solver.opts.cost_tolerance == 1e-1
@test solver.solver_al.opts.cost_tolerance == 1e-1
@test solver.solver_al.solver_uncon.opts.cost_tolerance == 1e-1
@test solver.solver_al.solver_uncon.opts.iterations == 210
@test solver.solver_al.solver_uncon.opts.square_root == true
@test solver.solver_al.opts.iterations_outer == 60
@test solver.solver_pn.opts.ρ_chol == 1e-1
@test solver.solver_pn.opts.constraint_tolerance == 1e-2
@test solver.solver_al.opts.constraint_tolerance == 1e-2
@test solver.opts.constraint_tolerance == 1e-2


text:
module Metrics

export xmetrics

include("xmetrics.jl")

end


text:
# LDU factorization for unsymmetric systems

function ldu_factorization_acyclic!(diagonal_v, offdiagonal_l, diagonal_c, offdiagonal_u, diagonal_inverse_c)
    if diagonal_inverse_c.isinverted
        invdiagonal_c = diagonal_inverse_c.value
    else
        invdiagonal_c = inv(diagonal_c.value)
        diagonal_inverse_c.value = invdiagonal_c
        diagonal_inverse_c.isinverted = true
    end
    offdiagonal_l.value = offdiagonal_l.value * invdiagonal_c
    offdiagonal_u.value = invdiagonal_c * offdiagonal_u.value
    diagonal_v.value -= offdiagonal_l.value*diagonal_c.value*offdiagonal_u.value
    return
end
function ldu_factorization_cyclic!(entry_lu, offdiagonal_lu, diagonal_c, offdiagonal_ul)
    entry_lu.value -= offdiagonal_lu.value*diagonal_c.value*offdiagonal_ul.value
    return
end

function ldu_factorization!(system)
    matrix_entries = system.matrix_entries
    diagonal_inverses = system.diagonal_inverses
    acyclic_children = system.acyclic_children
    cyclic_children = system.cyclic_children

    for v in system.dfs_list
        for c in cyclic_children[v]
            for cc in cyclic_children[v]
                cc == c && break
                (cc ∉ children(system,c) && cc ∉ cyclic_children[c]) && continue
                ldu_factorization_cyclic!(matrix_entries[v,c], matrix_entries[v,cc], matrix_entries[cc,cc], matrix_entries[cc,c])
                ldu_factorization_cyclic!(matrix_entries[c,v], matrix_entries[c,cc], matrix_entries[cc,cc], matrix_entries[cc,v])
            end
            ldu_factorization_acyclic!(matrix_entries[v,v], matrix_entries[v,c], matrix_entries[c,c], matrix_entries[c,v], diagonal_inverses[c])
        end
        for c in acyclic_children[v]
            ldu_factorization_acyclic!(matrix_entries[v,v], matrix_entries[v,c], matrix_entries[c,c], matrix_entries[c,v], diagonal_inverses[c])
        end
    end
    return
end

function ldu_backsubstitution_l!(vector_v, offdiagonal, vector_c)
    vector_v.value -= offdiagonal.value*vector_c.value
    return
end
function ldu_backsubstitution_u!(vector_v, offdiagonal, vector_p)
    vector_v.value -= offdiagonal.value*vector_p.value
    return
end
function ldu_backsubstitution_d!(vector, diagonal, diagonal_inverse)
    if diagonal_inverse.isinverted
        vector.value = diagonal_inverse.value*vector.value
    else
        vector.value = diagonal.value\vector.value
    end
    diagonal_inverse.isinverted = false
    return
end

function ldu_backsubstitution!(system)
    matrix_entries = system.matrix_entries
    diagonal_inverses = system.diagonal_inverses
    vector_entries = system.vector_entries
    acyclic_children = system.acyclic_children
    cyclic_children = system.cyclic_children
    parents = system.parents
    dfs_list = system.dfs_list

    for v in dfs_list
        for c in cyclic_children[v]
            ldu_backsubstitution_l!(vector_entries[v], matrix_entries[v,c], vector_entries[c])
        end
        for c in acyclic_children[v]
            ldu_backsubstitution_l!(vector_entries[v], matrix_entries[v,c], vector_entries[c])
        end
    end
    for v in reverse(dfs_list)
        ldu_backsubstitution_d!(vector_entries[v], matrix_entries[v,v], diagonal_inverses[v])
        for p in parents[v]
            ldu_backsubstitution_u!(vector_entries[v], matrix_entries[v,p], vector_entries[p])
        end
    end
end

function ldu_solve!(system)
    ldu_factorization!(system)
    ldu_backsubstitution!(system)
    return
end


text:
module CamiXon

using Dates
using IntervalSets
#using LaTeXStrings
using LinearAlgebra

export sup
export sub
export frac

export mendeleev
export Atom
export createAtom
export Term
export createTerm
export Grid
export gridfunction
export createGrid
export matG

export FITS_HDU
export FITS_header
export FITS_data
export FITS_table
export parse_FITS_TABLE
export FITS_name
export cast_FITS_name

export f_diff_weight
export f_diff_weights
export f_diff_weights_array
export f_diff_expansion_coeffs_lagrange
export f_diff_expansion_weights
export summation_range
export f_diff_function_sequences
export lagrange_interpolation
export lagrange_extrapolation
export f_diff_expansion_coeffs_differentiation
export create_lagrange_differentiation_weights
export create_lagrange_differentiation_matrix
export lagrange_differentiation
export f_diff_expansion_coeffs_adams_moulton
export create_adams_moulton_weights
export f_diff_expansion_coeffs_adams_bashford
export trapezoidal_weights
export trapezoidal_integration

export fits_create
export fits_read
export fits_extend
export fits_info
export fits_copy
export fits_combine
export fits_add_key
export fits_edit_key
export fits_delete_key
export fits_rename_key

#export plot_matrices
#export plot!
export step125
export select125
export edges
export steps
export stepcenters
export stepedges
#export centerticks
#export edgeticks
#export centers
#export edges

export FORTRAN_format
export cast_FORTRAN_format
export cast_FORTRAN_datatype

export find_all
export find_first
export find_last

export bernoulli_numbers
export VectorRational
export normalize_VectorRational
export canonical_partitions
export faulhaber_polynom
export faulhaber_summation
export harmonic_number
export integer_partitions
export log10_characteristic_power
export log10_mantissa
export pascal_triangle
export pascal_next
export polynom
export polynom_derivative
export polynom_derivatives
export polynom_primitive
export polynom_product
export polynom_product_expansion
export polynom_power
export polynom_powers
export permutations_unique_count

include("strings.jl")
include("adams.jl")
include("atom.jl")
include("fits_pointers.jl")
include("read_io.jl")
include("write_io.jl")
include("finite_differences.jl")
include("fits_objects.jl")
include("fits_private_sector.jl")
include("fits_public_sector.jl")
include("plot_private_sector.jl")
include("plot_public_sector.jl")
include("Header-and-Data-Input.jl")
include("FORTRAN.jl")
include("search_algorithms.jl")
include("mathematics.jl")

end


text:
"""There is the possibility of the created vector existing outside of the bounds specified at initialization.
 Because of this, we created a separate function i.e BoundEnsure that checks and corrects for this. In the event
 that one of these rogue points are found, we’ll simply move it to the nearest boundary """

function bound_ensure(m::SVector{N, T}, u::SVector{N, T}, X::IntervalBox{N, T}) where {N, T}

    for i in 1:length(X)
      if m[i] < X[i].lo
        m = setindex(m, u[i]+rand()*(X[i].lo-u[i]), i)
      elseif X[i].hi < m[i]
        m = setindex(m, u[i]+rand()*(X[i].hi-u[i]), i)
      end
    end
    return m
end