text:
## gzip file io ##

module GZip

import Base: show, fd, close, flush, truncate, seek,
             skip, position, eof, read, readall,
             readuntil, readline, write, peek

export
  GZipStream,
  show,

# io functions
# open,  ## not exported; use as GZip.open(...)
  gzopen,
  gzdopen,
  fd,
  close,
  flush,
  truncate,
  seek,
  skip,
  position,
  eof,
  read,
  readall,
  readuntil,
  readline,
  write,
  peek,

# lower-level io functions
  gzgetc,
  gzungetc,
  gzgets,
  gzputc,
  gzwrite,
  gzread,
  gzbuffer,

# File offset
  ZFileOffset,

# GZError, ZError, related constants (zlib_h.jl)
  GZError,
  ZError,
  Z_OK,
  Z_STREAM_END,
  Z_NEED_DICT,
  Z_ERRNO,
  Z_STREAM_ERROR,
  Z_DATA_ERROR,
  Z_MEM_ERROR,
  Z_BUF_ERROR,
  Z_VERSION_ERROR,

# Compression constants (zlib_h.jl)
  Z_NO_COMPRESSION,
  Z_BEST_SPEED,
  Z_BEST_COMPRESSION,
  Z_DEFAULT_COMPRESSION,

# Compression strategy (zlib_h.jl)
  Z_FILTERED,
  Z_HUFFMAN_ONLY,
  Z_RLE,
  Z_FIXED,
  Z_DEFAULT_STRATEGY,

# Default buffer sizes
  Z_DEFAULT_BUFSIZE,
  Z_BIG_BUFSIZE

include("zlib_h.jl")

# Expected line length for strings
const GZ_LINE_BUFSIZE = 256

# Wrapper around gzFile
type GZipStream <: IO
    name::String
    gz_file::Ptr{Void}
    buf_size::Int

    _closed::Bool

    GZipStream(name::String, gz_file::Ptr{Void}, buf_size::Int) =
        (x = new(name, gz_file, buf_size, false); finalizer(x, close); x)
end
GZipStream(name::String, gz_file::Ptr{Void}) = GZipStream(name, gz_file, Z_DEFAULT_BUFSIZE)

# gzerror
function gzerror(err::Integer, s::GZipStream)
    e = Int32[err]
    if !s._closed
        msg_p = ccall((:gzerror, _zlib), Ptr{Uint8}, (Ptr{Void}, Ptr{Int32}),
                      s.gz_file, e)
        msg = (msg_p == C_NULL ? "" : bytestring(msg_p))
    else
        msg = "(GZipStream closed)"
    end
    (e[1], msg)
end
gzerror(s::GZipStream) = gzerror(0, s)

type GZError <: Exception
    err::Int32
    err_str::String

    GZError(e::Integer, str::String) = new(int32(e), str)
    GZError(e::Integer, s::GZipStream) = (a = gzerror(e, s); new(a[1], a[2]))
    GZError(s::GZipStream) = (a = gzerror(s); new(a[1], a[2]))
end

# show
show(io::IO, s::GZipStream) = print(io, "GZipStream(", s.name, ")")

macro test_eof_gzerr(s, cc, val)
    quote
        if $(esc(s))._closed throw(EOFError()) end
        ret = $(esc(cc))
        if ret == $(esc(val))
            if eof($(esc(s)))  throw(EOFError())  else  throw(GZError($(esc(s))))  end
        end
        ret
    end
end

macro test_eof_gzerr2(s, cc, val)
    quote
        if $(esc(s))._closed throw(EOFError()) end
        ret = $(esc(cc))
        if ret == $(esc(val)) && !eof($(esc(s))) throw(GZError($(esc(s)))) end
        ret
    end
end

macro test_gzerror(s, cc, val)
    quote
        if $(esc(s))._closed throw(EOFError()) end
        ret = $(esc(cc))
        if ret == $(esc(val)) throw(ret, GZError($(esc(s)))) end
        ret
    end
end

macro test_gzerror0(s, cc)
    quote
        if $(esc(s))._closed throw(EOFError()) end
        ret = $(esc(cc))
        if ret <= 0 throw(GZError(ret, $(esc(s)))) end
        ret
    end
end

macro test_z_ok(cc)
    quote
        ret = $(esc(cc))
        if (ret != Z_OK) throw(ZError(ret)) end
        ret
    end
end

# Easy access to gz reading/writing functions (Internal)
gzgetc(s::GZipStream) =
    @test_eof_gzerr(s, ccall((:gzgetc, _zlib), Int32, (Ptr{Void},), s.gz_file), -1)

gzgetc_raw(s::GZipStream) = ccall((:gzgetc, _zlib), Int32, (Ptr{Void},), s.gz_file)

gzungetc(c::Integer, s::GZipStream) =
    @test_eof_gzerr(s, ccall((:gzungetc, _zlib), Int32, (Int32, Ptr{Void}), c, s.gz_file), -1)

gzgets(s::GZipStream, a::Array{Uint8}) =
    @test_eof_gzerr2(s, ccall((:gzgets, _zlib), Ptr{Uint8}, (Ptr{Void}, Ptr{Uint8}, Int32),
                              s.gz_file, a, int32(length(a))),          C_NULL)

gzgets(s::GZipStream, p::Ptr{Uint8}, len::Integer) =
    @test_eof_gzerr2(s, ccall((:gzgets, _zlib), Ptr{Uint8}, (Ptr{Void}, Ptr{Uint8}, Int32),
                              s.gz_file, p, int32(len)),                C_NULL)

gzputc(s::GZipStream, c::Integer) =
    @test_gzerror(s, ccall((:gzputc, _zlib), Int32, (Ptr{Void}, Int32),
                           s.gz_file, int32(c)),                        -1)

gzwrite(s::GZipStream, p::Ptr, len::Integer) =
    len == 0 ? int32(0) : @test_gzerror0(s, ccall((:gzwrite, _zlib), Int32, (Ptr{Void}, Ptr{Void}, Uint32),
                                                   s.gz_file, p, len))

gzread(s::GZipStream, p::Ptr, len::Integer) =
    @test_gzerror(s, ccall((:gzread, _zlib), Int32, (Ptr{Void}, Ptr{Void}, Uint32),
                           s.gz_file, p, len),                          -1)

let _zlib_h = dlopen(_zlib)
    global gzbuffer, _gzopen, _gzseek, _gztell

    # Doesn't exist in zlib 1.2.3 or earlier
    if dlsym_e(_zlib_h, :gzbuffer) != C_NULL
        gzbuffer(gz_file::Ptr, gz_buf_size::Integer) =
           ccall((:gzbuffer, _zlib), Int32, (Ptr{Void}, Uint32), gz_file, gz_buf_size)
    else
        gzbuffer(gz_file::Ptr, gz_buf_size::Integer) = int32(-1)
    end

    #####

    # Use 64-bit functions if available

    if dlsym_e(_zlib_h, :gzopen64) != C_NULL
        const _gzopen = :gzopen64
        const _gzseek = :gzseek64
        const _gztell = :gztell64
        #_gzoffset = :gzoffset64  ## not implemented
    else
        const _gzopen = :gzopen
        const _gzseek = :gzseek
        const _gztell = :gztell
        #_gzoffset = :gzoffset    ## not implemented
    end
end

function gzopen(fname::String, gzmode::String, gz_buf_size::Integer)
    # gzmode can contain extra characters specifying
    # * compression level (0-9)
    # * strategy ('f' => filtered data, 'h' -> Huffman-only compression,
    #             'R' -> run-length encoding, 'F' -> fixed code compression)
    #
    # '+' is also not allowed

    # For windows, force binary mode; doesn't hurt on unix
    if !('b' in gzmode)
        gzmode *= "b"
    end

    gz_file = ccall((_gzopen, _zlib), Ptr{Void}, (Ptr{Uint8}, Ptr{Uint8}), fname, gzmode)
    if gz_file == C_NULL
        throw(GZError(-1, "gzopen failed"))
    end
    if gz_buf_size != Z_DEFAULT_BUFSIZE
        if gzbuffer(gz_file, gz_buf_size) == -1
            # Generally a non-fatal error, although it shouldn't happen here
            gz_buf_size = Z_DEFAULT_BUFSIZE
        end
    end
    return GZipStream(fname, gz_file, gz_buf_size)
end
gzopen(fname::String, gzmode::String) = gzopen(fname, gzmode, Z_DEFAULT_BUFSIZE)
gzopen(fname::String) = gzopen(fname, "rb", Z_DEFAULT_BUFSIZE)
open(args...) = gzopen(args...)

function gzopen(f::Function, args...)
    io = gzopen(args...)
    x = try f(io) catch err
        close(io)
        throw(err)
    end
    close(io)
    return x
end

function gzdopen(name::String, fd::Integer, gzmode::String, gz_buf_size::Integer)
    if !('b' in gzmode)
        gzmode *= "b"
    end

    # Duplicate the file descriptor, since we have no way to tell gzclose()
    # not to close the original fd
    dup_fd = ccall(:dup, Int32, (Int32,), fd)

    gz_file = ccall((:gzdopen, _zlib), Ptr{Void}, (Int32, Ptr{Uint8}), dup_fd, gzmode)
    if gz_file == C_NULL
        throw(GZError(-1, "gzdopen failed"))
    end
    if gz_buf_size != Z_DEFAULT_BUFSIZE
        if gzbuffer(gz_file, gz_buf_size) == -1
            # Generally a non-fatal error, although it shouldn't happen here
            gz_buf_size = Z_DEFAULT_BUFSIZE
        end
    end
    return GZipStream(name, gz_file, gz_buf_size)
end
gzdopen(fd::Integer, gzmode::String, gz_buf_size::Integer) = gzdopen(string("<fd ",fd,">"), fd, gzmode, gz_buf_size)
gzdopen(fd::Integer, gz_buf_size::Integer) = gzdopen(fd, "rb", gz_buf_size)
gzdopen(fd::Integer, gzmode::String) = gzdopen(fd, gzmode, Z_DEFAULT_BUFSIZE)
gzdopen(fd::Integer) = gzdopen(fd, "rb", Z_DEFAULT_BUFSIZE)
gzdopen(s::IOStream, args...) = gzdopen(fd(s), args...)

fd(s::GZipStream) = error("fd is not supported for GZipStreams")

function close(s::GZipStream)

    # The garbage collector needs to be temporarily disabled because of a possible
    # race condition if it runs between testing and setting s._closed

    gc_disable()
    if s._closed
        gc_enable()
        return Z_STREAM_ERROR
    end
    s._closed = true
    gc_enable()

    s.name *= " (closed)"

    ret = (@test_z_ok ccall((:gzclose, _zlib), Int32, (Ptr{Void},), s.gz_file))

    return ret
end

flush(s::GZipStream, fl::Integer) =
    @test_z_ok ccall((:gzflush, _zlib), Int32, (Ptr{Void}, Int32), s.gz_file, int32(fl))
flush(s::GZipStream) = flush(s, Z_SYNC_FLUSH)

truncate(s::GZipStream, n::Integer) = error("truncate is not supported for GZipStreams")

# Note: seeks to byte position within uncompressed data stream
seek(s::GZipStream, n::Integer) =
    (ccall((_gzseek, _zlib), ZFileOffset, (Ptr{Void}, ZFileOffset, Int32),
           s.gz_file, n, SEEK_SET)!=-1 || # Mimick behavior of seek(s::IOStream, n)
    error("seek (gzseek) failed"))

seek_end(s::GZipStream) = error("seek_end is not supported for GZipStreams")

# Note: skips bytes within uncompressed data stream
skip(s::GZipStream, n::Integer) =
    (ccall((_gzseek, _zlib), ZFileOffset, (Ptr{Void}, ZFileOffset, Int32),
           s.gz_file, n, SEEK_CUR)!=-1 ||
     error("skip (gzseek) failed")) # Mimick behavior of skip(s::IOStream, n)

position(s::GZipStream) =
    ccall((_gztell, _zlib), ZFileOffset, (Ptr{Void},), s.gz_file)

eof(s::GZipStream) = bool(ccall((:gzeof, _zlib), Int32, (Ptr{Void},), s.gz_file))

function peek(s::GZipStream)
    c = gzgetc_raw(s)
    if c != -1
        gzungetc(c, s)
    end
    c
end

# Mimics read(s::IOStream, a::Array{T})
function read{T}(s::GZipStream, a::Array{T})
    if isbits(T)
        nb = length(a)*sizeof(T)
        # Note: this will overflow and succeed without warning if nb > 4GB
        ret = ccall((:gzread, _zlib), Int32,
                    (Ptr{Void}, Ptr{Void}, Uint32), s.gz_file, a, nb)
        if ret == -1
            throw(GZError(s))
        end
        if ret < nb
            throw(EOFError())  # TODO: Do we have/need a way to read without throwing an error near the end of the file?
        end
        peek(s) # force eof to be set
        a
    else
        invoke(read, (IO, Array), s, a)
    end
end

function read(s::GZipStream, ::Type{Uint8})
    ret = gzgetc(s)
    if ret == -1
        throw(GZError(s))
    end
    peek(s) # force eof to be set
    uint8(ret)
end

# For this function, it's really unfortunate that zlib is
# not integrated with ios
function readall(s::GZipStream, bufsize::Int)
    buf = Array(Uint8, bufsize)
    len = 0
    while true
        ret = gzread(s, pointer(buf)+len, bufsize)
        if ret == 0
            # check error status to make sure stream was not truncated
            # (we won't normally get an error until the close, because it's
            # possible that the file is still being written to.)

            ## *** Disabled, to allow the function to return the buffer ***
            ## *** Truncation error will be generated on gzclose... ***

            #(err, msg) = gzerror(s)
            #if err != Z_OK
            #    throw(GZError(err, msg))
            #end

            # Resize buffer to exact length
            if length(buf) > len
                resize!(buf, len)
            end
            return bytestring(buf)
        end
        len += ret
        # Grow the buffer so that bufsize bytes will fit
        resize!(buf, bufsize+len)
    end
end
readall(s::GZipStream) = readall(s, Z_BIG_BUFSIZE)

function readuntil(s::GZipStream, c::Uint8)
    buf = Array(Uint8, GZ_LINE_BUFSIZE)
    pos = 1

    if gzgets(s, buf) == C_NULL      # Throws an exception on error
        return buf[1:0]
    end

    while(true)
        # since gzgets didn't return C_NULL, there must be a \0 in the buffer
        eos = search(buf, '\0', pos)
        if eos == 1 || buf[eos-1] == c
            return resize!(buf, eos-1)
        end

        # If we're at the end of the file, return the string
        if eof(s)  return resize!(buf, eos-1)  end

        # Otherwise, append to the end of the previous buffer

        # Grow the buffer so that there's room for GZ_LINE_BUFSIZE chars
        add_len = GZ_LINE_BUFSIZE - (length(buf)-eos+1)
        resize!(buf, add_len+length(buf))
        pos = eos

        # Read in the next chunk
        if gzgets(s, pointer(buf)+pos-1, GZ_LINE_BUFSIZE) == C_NULL
            # eof(s); remove extra buffer space
            return resize!(buf, length(buf)-add_len)
        end
    end
end

write(s::GZipStream, b::Uint8) = gzputc(s, b)

function write{T}(s::GZipStream, a::Array{T})
    if isbits(T)
        return gzwrite(s, pointer(a), length(a)*sizeof(T))
    else
        invoke(write, (Any, Array), s, a)
    end
end

write(s::GZipStream, p::Ptr, nb::Integer) = gzwrite(s, p, nb)

function write{T,N}(s::GZipStream, a::SubArray{T,N,Array})
    if !isbits(T) || stride(a,1)!=1
        return invoke(write, (Any, AbstractArray), s, a)
    end
    colsz = size(a,1)*sizeof(T)
    if N==1
        write(s, pointer(a, 1), colsz)
    else
        cartesian_map((idxs...)->write(s, pointer(a, idxs), colsz),
                      tuple(1, size(a)[2:]...))
    end
end

end # module GZip


text:
@testset "create_directions.jl" begin

end

ext: jl


text:
function [D,J,X]=JacobiCart1(ksi,F,Grid)

x1=Grid.Nodes(F.N(1)).P(1);
y1=Grid.Nodes(F.N(1)).P(2);
x2=Grid.Nodes(F.N(2)).P(1);
y2=Grid.Nodes(F.N(2)).P(2);
x3=Grid.Nodes(F.N(3)).P(1);
y3=Grid.Nodes(F.N(3)).P(2);
x4=Grid.Nodes(F.N(4)).P(1);
y4=Grid.Nodes(F.N(4)).P(2);
xBar=ksi(1);
yBar=ksi(2);

J=[ (x1 - x2 + x3 - x4)*yBar + x2 - x1, (x1 - x2 + x3 - x4)*xBar + x4 - x1
     (y1 - y2 + y3 - y4)*yBar + y2 - y1, (y1 - y2 + y3 - y4)*xBar + y4 - y1];
J=[J
   0 0];
D=J(1,1)*J(2,2)-J(1,2)*J(2,1);
if nargout > 2
  X=zeros(1,3);
  X(1)=Grid.Nodes(F.N(1)).P(1)...
    +(Grid.Nodes(F.N(2)).P(1)-Grid.Nodes(F.N(1)).P(1))*ksi(1)...
    +(Grid.Nodes(F.N(4)).P(1)-Grid.Nodes(F.N(1)).P(1))*ksi(2)...
    +(Grid.Nodes(F.N(3)).P(1)-Grid.Nodes(F.N(4)).P(1)-Grid.Nodes(F.N(2)).P(1)+Grid.Nodes(F.N(1)).P(1))*ksi(1)*ksi(2);
  X(2)=Grid.Nodes(F.N(1)).P(2)...
    +(Grid.Nodes(F.N(2)).P(2)-Grid.Nodes(F.N(1)).P(2))*ksi(1)...
    +(Grid.Nodes(F.N(4)).P(2)-Grid.Nodes(F.N(1)).P(2))*ksi(2)...
    +(Grid.Nodes(F.N(3)).P(2)-Grid.Nodes(F.N(4)).P(2)-Grid.Nodes(F.N(2)).P(2)+Grid.Nodes(F.N(1)).P(2))*ksi(1)*ksi(2);
  X(3)=Grid.Nodes(F.N(1)).P(3)...
    +(Grid.Nodes(F.N(2)).P(3)-Grid.Nodes(F.N(1)).P(3))*ksi(1)...
    +(Grid.Nodes(F.N(4)).P(3)-Grid.Nodes(F.N(1)).P(3))*ksi(2)...
    +(Grid.Nodes(F.N(3)).P(3)-Grid.Nodes(F.N(4)).P(3)-Grid.Nodes(F.N(2)).P(3)+Grid.Nodes(F.N(1)).P(3))*ksi(1)*ksi(2);
end
end


text:
module RealSense

include("LibRealSense.jl")
using .LibRealSense

end # module


text:
"""
    constants.jl

Description:
    All constant values associated with the package.
"""

# DDVFA linkage methods
const DDVFA_METHODS = [
    "single",
    "average",
    "complete",
    "median",
    "weighted",
    "centroid"
]


text:
#!/usr/bin/env julia

using Luxor

using Test

using Random
Random.seed!(42)

function hipster(fname)
    Drawing(400, 350, fname) # hide
    origin() # hide
    rotate(pi/8)

    circle(O, 135, :clip)
    sethue("antiquewhite2")
    paint()

    sethue("gray20")
    setline(3)
    circle(O, 130, :stroke)
    circle(O, 135, :stroke)
    circle(O, 125, :fill)

    sethue("antiquewhite2")
    circle(O, 85, :fill)

    sethue("wheat")
    fontsize(24)
    fontface("Helvetica-Bold")
    # Try both synonyms (long live the Queen, and UK spelling!)
    textcurvecentered("• LUXOR •", (3pi)/2, 100, O, clockwise=true, baselineshift = -4)
    textcurvecentred("• VECTOR GRAPHICS •", pi/2, 100, O, clockwise=false, letter_spacing=2, baselineshift = -15)

    sethue("gray20")
    map(pt -> star(pt, 40, 3, 0.5, -pi/2, :fill), ngon(O, 40, 3, 0, vertices=true))
    circle(O.x + 30, O.y - 55, 15, :fill)

    sethue("antiquewhite2")
    setline(0.2)
    setdash("dotdotdashed")
    for i in 1:500
        line(randompoint(Point(-200, -350), Point(200, 350)),
             randompoint(Point(-200, -350), Point(200, 350)),
             :stroke)
    end
    @test finish() == true
end

fname = "text-curve-centered.png"
hipster(fname)
println("...finished test: output in $(fname)")


text:
using Documenter, BBI

makedocs(
    # format = :html,
    format = Documenter.HTML(
        edit_branch = "develop"
    ),
    sitename = "BBI.jl",
    pages = [
        "Home" => "index.md",
    ],
    authors = "Kenta Sato, D. C. Jones, Ben J. Ward, Ciarán O'Mara, The BioJulia Organisation and other contributors.",
)
deploydocs(
    repo = "github.com/BioJulia/BBI.jl.git",
    target = "build",
    deps = nothing,
    make = nothing,
    devbranch = "develop"
)


text:
###############################################################################
#
# Convert between YYYYMMDD and DOY
#
###############################################################################
"""
Researchers may often run into the scenarios when they have to convert dates to strings or numbers, vice versa. Thought the conversions are easy to make,
    it is not convenient at all if ones need to recreate the wheel everytime. Thus, we provide functions to make the conversion. One of the common
    conversions is to convert time stamp to a number:

$(METHODLIST)

"""
function parse_timestamp end

"""
The time stamp can be either an integer or a string, and the integer will be converted to a string within the function. The supported time stamp formats
    are `YYYYMMDD`, `YYYYMMDDhh`, `YYYYMMDDhhmm`, and `YYYYMMDDhhmmss`, where the labels are

- `YYYY` Year number
- `MM` Month number
- `DD` Day number
- `hh` Hour number
- `mm` Minute number
- `ss` second number

The supported outputs are
- `DATE` A `Dates.Date` type variable
- `DATETIME` A `Dates.DateTime` type variable
- `DOY` A day of year integer
- `FDOY` A day of year float

    parse_timestamp(
                time_stamp::Union{Int,String};
                in_format::String = "YYYYMMDD",
                out_format::String = "DOY")

Convert time stamp, given
- `time_stamp` Time stamp
- `in_format` Format of `time_stamp`, default is `YYYYMMDD`
- `out_format` Output format, default is `DOY`

---
Examples
```julia
# if input format is YYYYMMDD
_inf = "YYYYMMDD";
for _ts in [20200130, "20200130"]
    for _outf in ["DATE", "DATETIME", "DOY", "FDOY"]
        time = parse_timestamp(_ts; in_format=_inf, out_format=_outf);
    end
end

# if input format is YYYYMMDDhh
_inf = "YYYYMMDDhh";
for _ts in [2020013002, "2020013002"]
    for _outf in ["DATE", "DATETIME", "DOY", "FDOY"]
        time = parse_timestamp(_ts; in_format=_inf, out_format=_outf);
    end
end

# if input format is YYYYMMDDhhmm
_inf = "YYYYMMDDhhmm";
for _ts in [202001300215, "202001300215"]
    for _outf in ["DATE", "DATETIME", "DOY", "FDOY"]
        time = parse_timestamp(_ts; in_format=_inf, out_format=_outf);
    end
end

# if input format is YYYYMMDDhhss
_inf = "YYYYMMDDhhmmss";
for _ts in [20200130021515, "20200130021515"]
    for _outf in ["DATE", "DATETIME", "DOY", "FDOY"]
        time = parse_timestamp(_ts; in_format=_inf, out_format=_outf);
    end
end
```
"""
parse_timestamp(
            time_stamp::Union{Int,String};
            in_format::String = "YYYYMMDD",
            out_format::String = "DOY") =
(
    _time_stamp = string(time_stamp);
    @assert in_format in TIME_FORMAT;
    @assert out_format in TIME_OUTPUT;

    # get year, month, and day
    _year   = parse(Int, _time_stamp[1:4]);
    _month  = parse(Int, _time_stamp[5:6]);
    _day    = parse(Int, _time_stamp[7:8]);
    _hour   = 0;
    _minute = 0;
    _second = 0;
    @assert 0 <= _year  <= 9999;
    @assert 1 <= _month <= 12;
    @assert 1 <= _day   <= month_days(_year,_month);

    # get hour, minute, and second
    if length(in_format) >= 10
        _hour   = parse(Int, _time_stamp[9:10]);
        @assert 0 <= _hour <= 24;
    end;
    if length(in_format) >= 12
        _minute = parse(Int, _time_stamp[11:12]);
        @assert 0 <= _minute <= 59;
    end;
    if length(in_format) >= 14
        _second = parse(Int, _time_stamp[13:14]);
        @assert 0 <= _second <= 59;
    end;

    # determine the
    if out_format=="DATE"
        return Date(_year, _month, _day);
    elseif out_format=="DATETIME"
        return DateTime(_year, _month, _day, _hour, _minute, _second);
    elseif out_format=="DOY"
        return (isleapyear(_year) ? MDAYS_LEAP[_month] : MDAYS[_month]) + _day;
    elseif out_format=="FDOY"
        _doy = (isleapyear(_year) ? MDAYS_LEAP[_month] : MDAYS[_month]) + _day;
        return _doy + (_hour + _minute/60 + _second/3600) / 24;
    end;
)

"""
Also, one might want to convert date to time stamp. Supported conversions are

    parse_timestamp(year::Int, doy::Int, sep::String="")

Convert date to `YYYY(sep)MM(sep)DD` format, given
- `year` Year
- `doy` Day of year (typically 1-365, 1-366 for leap years)
- `sep` String that separates year, month, and day. Default is ""

---
Examples
```julia
# convert the time to YYYYMMDD format
time = parse_timestamp(2020, 100);
time = parse_timestamp(2020, 100, "");

# convert the time to YYYY.MM.DD format
time = parse_timestamp(2020, 100, ".");
```
"""
parse_timestamp(year::Int, doy::Int, sep::String="") =
(
    @assert 0 <= year <= 9999;
    _month = month_ind(year, doy);
    _day   = doy - (isleapyear(year) ? MDAYS_LEAP[_month] : MDAYS[_month]);

    return string(year) * sep * lpad(_month,2,"0") * sep * lpad(_day,2,"0")
)

"""
    month_days(year::Int, month::Int)

How many days in a month, given
- `year` Year
- `month` Month

---
Examples
```julia
# inquire how many days in 2020 October
mday = month_days(2020, 10);
```
"""
function month_days(year::Int, month::Int)
    @assert 1 <= month <= 12;

    return isleapyear(year) ? NDAYS_LEAP[month] : NDAYS[month]
end

"""
    month_ind(year::Int, doy::Int)

Which month does the day of year belongs to, given
- `year` Year
- `doy` Day of year (typically 1-365, 1-366 for leap years)

---
Examples
```julia
# inquire which month it is for 100th day in 2020
mind = month_ind(2020, 100);
```
"""
function month_ind(year::Int, doy::Int)
    if isleapyear(year)
        @assert 1 <= doy <= 366;
        _month = 1;
        for i in 1:12
            if MDAYS_LEAP[i] < doy <= MDAYS_LEAP[i+1]
                _month = i;
                break;
            end
        end
    else
        @assert 1 <= doy <= 365;
        _month = 1;
        for i in 1:12
            if MDAYS[i] < doy <= MDAYS[i+1]
                _month = i;
                break;
            end
        end
    end

    return _month
end


text:
using DSGE, ModelConstructors, HDF5, Random, FileIO, Test, Dates

writing_output = false
regenerate_sim_data = false

if VERSION < v"1.5"
    ver = "111"
else
    ver = "150"
end

#########################
# Regime-switching test #
# w/ An-Schorfheide     #
#########################
path = dirname(@__FILE__)

m = AnSchorfheide()

# Set up SMC
m <= Setting(:sampling_method, :MH)
m <= Setting(:data_vintage, "210101")
m <= Setting(:cond_vintage, "210101")
m <= Setting(:n_mh_simulations, 500)
m <= Setting(:n_mh_blocks, 1)
m <= Setting(:n_mh_burn, 0)
m <= Setting(:n_mh_thin, 1)
m <= Setting(:mh_cc, 0.0005)
m <= Setting(:mh_adaptive_accept, false)
m <= Setting(:hessian_path, joinpath(path, "..", "reference", "hessian_rs2=true_vint=210101.h5"))

sim_filepath = joinpath(path, "..", "reference", "sim_data_regswitch_anschorfheide.h5")
if regenerate_sim_data
    simulate_data(m; filepath = sim_filepath)
end
data = h5read(sim_filepath, "data")

true_lik = DSGE.likelihood(m, data)

m <= Setting(:saveroot, normpath(joinpath(dirname(@__FILE__), "..", "..", "save")))

# With regime-switching
m <= Setting(:regime_switching, true, true, "rs2", "") # For file output purposes
m <= Setting(:n_regimes, 3)
m <= Setting(:regime_dates, Dict(1 => date_presample_start(m), 2 => DSGE.iterate_quarters(Date(1960, 3, 31), 50),
                                 3 => DSGE.iterate_quarters(Date(1960, 3, 31), 100)))
setup_regime_switching_inds!(m)
m <= Setting(:model2para_regimes, Dict{Int, Int}(1 => 1, 2 => 1, 3 => 2))

for i in 1:length(m.parameters)
    for k in 1:2
        ModelConstructors.set_regime_val!(m.parameters[i], k, m.parameters[i].value)
    end
end

param_mat = repeat([1 1 2], length(m.parameters))

DSGE.setup_param_regimes!(m, param_mat)

sys = compute_system(m)

regswitch_lik = DSGE.likelihood(m, data)

true_para = ModelConstructors.get_values(m.parameters)

Random.seed!(1793)
@testset "Search for posterior mode of regime-switching AnSchorfheide with csminwel (approx. 5s)" begin
    m <= Setting(:optimization_attempts, 1)
    m <= Setting(:optimization_iterations, 3)
    DSGE.estimate(m, data; sampling = false)
    θ = load_draws(m, :mode)
    if writing_output
        h5open(joinpath(path, "..", "reference", "regime_switching_anschorfheide_paramsmode_output.h5"), "w") do file
            write(file, "params", θ)
        end
    else
        @test θ ≈ h5read(joinpath(path, "..", "reference", "regime_switching_anschorfheide_paramsmode_output.h5"), "params")
    end
end
m <= Setting(:reoptimize, false)

@testset "Calculate Hessian of regime-switching AnSchorfheide (approx. 10s)" begin
    out_hessian, _ = DSGE.hessian!(m, h5read(joinpath(path, "..", "reference",
                                                      "regime_switching_anschorfheide_paramsmode_output.h5"), "params"),
                                   data; check_neg_diag = false)
    if writing_output
        h5open(joinpath(path, "..", "reference", "hessian_rs2=true_vint=210101.h5"), "w") do file
            write(file, "hessian", out_hessian)
        end
    else
        # when saving the output in REPL, the results seem different from testing, but the Hessian is still fairly close
        if maximum(abs.(h5read(joinpath(path, "..", "reference", "hessian_rs2=true_vint=210101.h5"), "hessian") - out_hessian)) < 8e-2
            @test maximum(abs.(h5read(joinpath(path, "..", "reference", "hessian_rs2=true_vint=210101.h5"), "hessian") -
                               out_hessian)) < 8e-2
        else
            # usually, in REPL and in test mode, the Hessian satisfies the error bound, but occassionally,
            # the maximum difference is very large (e.g. on the order of 100 - 1000), for some spurious reason.
            # To avoid having tests break, we only run the test when we know it is satisfied. Otherwise, we mark it as broken
            @warn "Test for Hessian of regime-switching AnSchorfheide failed, double check if the error is spurious or not."
            @test_broken maximum(abs.(h5read(joinpath(path, "..", "reference", "hessian_rs2=true_vint=210101.h5"), "hessian") -
                                      out_hessian)) < 8e-2
        end
    end
end

m <= Setting(:calculate_hessian, false)
m <= Setting(:hessian_path, joinpath(path, "..", "reference", "hessian_rs2=true_vint=210101.h5"))

@testset "Estimate regime-switching AnSchorfheide with MH (approx. 10s)" begin

    @test true_lik ≈ regswitch_lik

    # Regime switching estimation
    Random.seed!(1793)
    DSGE.update!(m, h5read(joinpath(path, "..", "reference",
                           "regime_switching_anschorfheide_paramsmode_output.h5"), "params"))
    DSGE.estimate(m, data)

    posterior_means = vec(mean(load_draws(m, :full), dims = 1))

    @test length(posterior_means) == length(m.parameters) * 2
    @test maximum(abs.(true_para - posterior_means)) < .55
end


text:
mutable struct JobSubmission
    config::Config
    build::BuildRef
    statussha::String            # the SHA to send statuses to (since `build` can mutate)
    url::String                  # the URL linking to the triggering comment
    fromkind::Symbol             # `:pr`, `:review`, or `:commit`?
    prnumber::Union{Int,Nothing} # the job's PR number, if `fromkind` is `:pr` or `:review`
    func::String
    args::Vector{String}
    kwargs::Dict{Symbol,String}
end

function JobSubmission(config::Config, event::GitHub.WebhookEvent, submission_string)
    try
        build, statussha, url, fromkind, prnumber = parse_event(config, event)
        func, args, kwargs = parse_submission_string(submission_string)
        return JobSubmission(config, build, statussha, url, fromkind, prnumber, func, args, kwargs)
    catch err
        error(string("could not parse comment into job submission: ", sprint(showerror, err)))
    end
end

function Base.:(==)(a::JobSubmission, b::JobSubmission)
    if (a.prnumber === nothing) == (b.prnumber === nothing)
        same_prnumber = a.prnumber === nothing ? true : (a.prnumber == b.prnumber)
        return (same_prnumber && a.config == b.config && a.build == b.build &&
                a.statussha == b.statussha && a.url == b.url && a.fromkind == b.fromkind &&
                a.func == b.func && a.args == b.args)
    else
        return false
    end
end

function parse_event(config::Config, event::GitHub.WebhookEvent)
    if event.kind == "commit_comment"
        # A commit was commented on, and the comment contained a trigger phrase.
        # The primary repo is the location of the comment, and the primary SHA
        # is that of the commit that was commented on.
        repo = event.repository.full_name
        sha = event.payload["comment"]["commit_id"]
        url = event.payload["comment"]["html_url"]
        fromkind = :commit
        prnumber = nothing
    elseif event.kind == "pull_request_review_comment"
        # A diff was commented on, and the comment contained a trigger phrase.
        # The primary repo is the location of the head branch, and the primary
        # SHA is that of the commit associated with the diff.
        repo = event.payload["pull_request"]["head"]["repo"]["full_name"]
        sha = event.payload["comment"]["commit_id"]
        url = event.payload["comment"]["html_url"]
        fromkind = :review
        prnumber = Int(event.payload["pull_request"]["number"])
    elseif event.kind == "pull_request"
        # A PR was opened, and the description body contained a trigger phrase.
        # The primary repo is the location of the head branch, and the primary
        # SHA is that of the head commit. The PR number is provided, so that the
        # build can execute on the relevant merge commit.
        repo = event.payload["pull_request"]["head"]["repo"]["full_name"]
        sha = event.payload["pull_request"]["head"]["sha"]
        url = event.payload["pull_request"]["html_url"]
        fromkind = :pr
        prnumber = Int(event.payload["pull_request"]["number"])
    elseif event.kind == "issue_comment"
        # A comment was made in a PR, and it contained a trigger phrase. The
        # primary repo is the location of the PR's head branch, and the primary
        # SHA is that of the head commit. The PR number is provided, so that the
        # build can execute on the relevant merge commit.
        pr = GitHub.pull_request(event.repository, event.payload["issue"]["number"], auth=config.auth)
        repo = pr.head.repo.full_name
        sha = pr.head.sha
        url = event.payload["comment"]["html_url"]
        fromkind = :pr
        prnumber = Int(pr.number)
    end
    return BuildRef(repo, sha), sha, url, fromkind, prnumber
end

# `x` can only be Expr, Symbol, QuoteNode, T<:Number, or T<:AbstractString
phrase_argument(x::Union{Expr, Symbol, QuoteNode}) = string(x)
phrase_argument(x::Union{AbstractString, Number})  = repr(x)

function parse_submission_string(submission_string)
    fncall = match(r"`.*?`", submission_string).match[2:end-1]
    argind = findfirst(isequal('('), fncall)
    name = fncall[1:(argind - 1)]
    parsed_args = Meta.parse(replace(fncall[argind:end], ";" => ","))
    args, kwargs = Vector{String}(), Dict{Symbol,String}()
    if isa(parsed_args, Expr) && parsed_args.head == :tuple
        started_kwargs = false
        for x in parsed_args.args
            if isa(x, Expr) && (x.head == :kw || x.head == :(=)) && isa(x.args[1], Symbol)
                @assert !haskey(kwargs, x.args[1]) "kwargs must all be unique"
                kwargs[x.args[1]] = phrase_argument(x.args[2])
                started_kwargs = true
            else
                @assert !started_kwargs "kwargs must come after other args"
                push!(args, phrase_argument(x))
            end
        end
    else
        push!(args, phrase_argument(parsed_args))
    end
    return name, args, kwargs
end

function reply_status(sub::JobSubmission, state, description, url=nothing)
    params = Dict("state" => state,
                  "context" => "Nanosoldier",
                  "description" => snip(description, 140))
    url !== nothing && (params["target_url"] = url)
    return GitHub.create_status(sub.config.trackrepo, sub.statussha;
                                auth = sub.config.auth, params = params)
end

function reply_comment(sub::JobSubmission, message::AbstractString)
    commentplace = sub.prnumber === nothing ? sub.statussha : sub.prnumber
    commentkind = sub.fromkind == :review ? :pr : sub.fromkind
    return GitHub.create_comment(sub.config.trackrepo, commentplace, commentkind;
                                 auth = sub.config.auth, params = Dict("body" => message))
end

function upload_report_repo!(sub::JobSubmission, markdownpath, message)
    cfg = sub.config
    sha = cd(reportdir(cfg)) do
        run(`git add -A`)
        run(`git commit -m $message`)
        headsha = readchomp(`git rev-parse HEAD`)
        run(`git pull -X ours`)
        run(`git push`)
        return headsha
    end
    return "https://github.com/$(reportrepo(cfg))/blob/$(sha)/$(markdownpath)"
end


text:
using Lazy, Metatheory

#the successor function
function S(n)
  "S($n)"
end
#the set of natural numbers
N = iterated(S,0);

#=define peano addition
  Need the following:
a+0 =a
a + S(b) = S(a+b)
=#

max_stars_repo_name: thomastjdavis/BeginningProof

max_issues_repo_name: thomastjdavis/BeginningProof

max_forks_repo_name: thomastjdavis/BeginningProof


text:
# coding: utf-8

# In[1]:

#explanation can be found in dynamic programming version
# https://github.com/je-suis-tm/recursion/blob/master/edit%20distance%20dynamic%20programming.jl
#the only problem with recursion is that it is so freaking slow
#recursion is so inefficient in any programming language
#although it looks much more elegant than dynamic programming

# In[2]:

function edit_distance(text1,text2)

    if isempty(text1) || isempty(text2)

        return max(length(text1),length(text2))

    end

    #we are comparing characters here
    #to get string, we should do end:end
    if text1[end]==text2[end]

        replacement=0

    else

        replacement=1

    end

    steps=min(edit_distance(text1[1:end-1],text2)+1,
        edit_distance(text1,text2[1:end-1])+1,
        edit_distance(text1[1:end-1],text2[1:end-1])+replacement)

    return steps

end

# In[3]:

println(edit_distance("arsehole","asshoe"))

max_stars_repo_path: edit distance recursion.jl

max_issues_repo_path: edit distance recursion.jl

max_forks_repo_path: edit distance recursion.jl


text:
"""
    struct DemoSection <: Any
    DemoSection(root::String)

Constructs a demo section that holds either nested subsections or demo cards.

# Fields

Besides the root path to the demo section folder `root`, this struct has some other fields:

* `cards`: demo cards found in `root`
* `subsections`: nested subsections found in `root`

# Configuration

You can manage an extra `config.json` file to customize rendering of a demo section.
Supported items are:

* `order`: specify the cards order or subsections order. By default, it's case-insensitive alphabetic order.
* `title`: specify the title of this demo section. By default, it's the folder name of `root`.
* `description`: some header description that you want to add before demo cards.
* `properties`: a dictionary of properties that can be propagated to its children items. The same properties in
  the children items, if exist, have higher priority.

The following is an example of `config.json`:

```json
{
    "title": "learn by examples",
    "description": "some one-line description can be useful.",
    "order": [
        "quickstart.md",
        "array.md"
    ],
    "properties": {
        "notebook": "false",
        "julia": "1.6",
        "author": "Johnny Chen"
    }
}
```

# Examples

The following is the simplest folder structure of a `DemoSection`:

```text
section
├── demo_1.md
├── demo_2.md
├── demo_3.md
├── demo_4.md
├── demo_5.md
├── demo_6.md
├── demo_7.md
└── demo_8.md
```

The following is a typical nested folder structure of a `DemoSection`:

```text
section
├── config.json
├── part1
│   ├── demo_21.md
│   └── demo_22.md
└── part2
    ├── config.json
    ├── demo_23.md
    └── demo_24.md
```

!!! warning

    A section should only hold either subsections or demo files. A folder that has both subfolders and demo files (e.g., `*.md`) is invalid.

See also: [`MarkdownDemoCard`](@ref DemoCards.MarkdownDemoCard), [`DemoPage`](@ref DemoCards.DemoPage)
"""
struct DemoSection
    root::String
    cards::Vector
    subsections::Vector{DemoSection}
    title::String
    description::String
    # These properties will be shared by all children of it during build time
    properties::Dict{String, Any}
end

basename(sec::DemoSection) = basename(sec.root)

function DemoSection(root::String)::DemoSection
    root = replace(root, r"[/\\]" => Base.Filesystem.path_separator) # windows compatibility
    isdir(root) || throw(ArgumentError("section root should be a valid dir, instead it's $(root)"))

    path = joinpath.(root, filter(x->!startswith(x, "."), readdir(root))) # filter out hidden files
    card_paths = filter(x->isfile(x) && !endswith(x, config_filename), path)
    section_paths = filter(x->isdir(x)&&!(basename(x) in ignored_dirnames), path)

    if !isempty(card_paths) && !isempty(section_paths)
        throw(ArgumentError("section folder $(root) should only hold either cards or subsections"))
    end

    config_file = joinpath(root, config_filename)
    config = isfile(config_file) ? JSON.parsefile(config_file) : Dict()

    # For files that `democard` fails to recognized, dummy
    # `UnmatchedCard` will be generated. Currently, we only
    # throw warnings for it.
    cards = map(democard, card_paths)
    unmatches = filter(cards) do x
        x isa UnmatchedCard
    end
    if !isempty(unmatches)
        msg = join(map(basename, unmatches), "\", \"")
        @warn "skip unmatched file: \"$msg\"" section_dir=root
    end
    cards = filter!(cards) do x
        !(x isa UnmatchedCard)
    end

    section = DemoSection(root,
                          cards,
                          map(DemoSection, section_paths),
                          "",
                          "",
                          Dict{String, Any}())

    ordered_paths = joinpath.(root, load_config(section, "order"; config=config))
    if !isempty(section.cards)
        cards = map(democard, ordered_paths)
        subsections = []
    else
        cards = []
        subsections = map(DemoSection, ordered_paths)
    end

    title = load_config(section, "title"; config=config)
    description = load_config(section, "description"; config=config)

    properties = if haskey(config, "properties")
        Dict{String, Any}(config["properties"])
    else
        Dict{String, Any}()
    end

    DemoSection(root, cards, subsections, title, description, properties)
end

function load_config(sec::DemoSection, key; config=Dict())
    if isempty(config)
        config_file = joinpath(sec.root, config_filename)
        config = isfile(config_file) ? JSON.parsefile(config_file) : Dict()
    end
    # config could still be an empty dictionary

    if key == "order"
        haskey(config, key) || return get_default_order(sec)

        order = config[key]
        validate_order(order, sec)
        return order
    elseif key == "title"
        get(config, key, get_default_title(sec))
    elseif key == "description"
        get(config, key, "")
    else
        throw(ArgumentError("Unrecognized key $(key) for DemoSection"))
    end
end

"""return case-insensitive alphabetic order"""
function get_default_order(sec::DemoSection)
    order = isempty(sec.cards) ? basename.(sec.subsections) : basename.(sec.cards)
    sort(order, by = x->lowercase(x))
end

function is_demosection(dir)
    try
        # if fails to parse, then it is not a valid demo page
        @suppress_err DemoSection(dir)
        return true
    catch err
        @debug err
        return false
    end
end


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.2
#     language: julia
#     name: julia-1.6
# ---

# %%
using BenchmarkTools
using Random

n = 2^19
bitmask = bitrand(n)
boolmask = Vector{Bool}(bitmask)
bytemask = Vector{Int8}(bitmask)

@show summary(bitmask)
@show summary(boolmask)
@show summary(bytemask)
@show findall(bitmask) == findall(.!(.!bitmask)) == findall(!, .!bitmask) ==
      findall(boolmask) == findall(Bool.(boolmask)) == findall(.!(.!boolmask)) == findall(!, .!boolmask) ==
      findall(Bool.(bytemask)) == findall(bytemask .== 1) == findall(==(1), bytemask)

println(); flush(stdout)

print("findall(\$bitmask):              ")
@btime findall($bitmask)
print("findall(.!\$(.!bitmask)):        ")
@btime findall(.!$(.!bitmask))
print("findall(!, \$(.!bitmask)):     ")
@btime findall(!, $(.!bitmask))

print("findall(\$boolmask):           ")
@btime findall($boolmask)
print("findall(Bool.(\$boolmask)):      ")
@btime findall(Bool.($boolmask))
print("findall(.!\$(.!boolmask)):       ")
@btime findall(.!$(.!boolmask))
print("findall(!, \$(.!boolmask)):    ")
@btime findall(!, $(.!boolmask))

print("findall(Bool.(\$bytemask)):    ")
@btime findall(Bool.($bytemask))
print("findall(\$bytemask .== 1):       ")
@btime findall($bytemask .== 1)
print("findall(==(1), \$bytemask):    ")
@btime findall(==(1), $bytemask)
;

# %%

max_stars_repo_path: 0016/apricot/findall benchmark.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0016/apricot/findall benchmark.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0016/apricot/findall benchmark.jl

max_forks_repo_name: genkuroki/public


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

@doc raw"""An EphemeralContainer is a container that may be added temporarily to an existing pod for user-initiated activities such as debugging. Ephemeral containers have no resource or scheduling guarantees, and they will not be restarted when they exit or when a pod is removed or restarted. If an ephemeral container causes a pod to exceed its resource allocation, the pod may be evicted. Ephemeral containers may not be added by directly updating the pod spec. They must be added via the pod's ephemeralcontainers subresource, and they will appear in the pod spec once added. This is an alpha feature enabled by the EphemeralContainers feature flag.

    IoK8sApiCoreV1EphemeralContainer(;
        args=nothing,
        command=nothing,
        env=nothing,
        envFrom=nothing,
        image=nothing,
        imagePullPolicy=nothing,
        lifecycle=nothing,
        livenessProbe=nothing,
        name=nothing,
        ports=nothing,
        readinessProbe=nothing,
        resources=nothing,
        securityContext=nothing,
        startupProbe=nothing,
        stdin=nothing,
        stdinOnce=nothing,
        targetContainerName=nothing,
        terminationMessagePath=nothing,
        terminationMessagePolicy=nothing,
        tty=nothing,
        volumeDevices=nothing,
        volumeMounts=nothing,
        workingDir=nothing,
    )

    - args::Vector{String} : Arguments to the entrypoint. The docker image's CMD is used if this is not provided. Variable references $(VAR_NAME) are expanded using the container's environment. If a variable cannot be resolved, the reference in the input string will be unchanged. The $(VAR_NAME) syntax can be escaped with a double $$, ie: $$(VAR_NAME). Escaped references will never be expanded, regardless of whether the variable exists or not. Cannot be updated. More info: https://kubernetes.io/docs/tasks/inject-data-application/define-command-argument-container/#running-a-command-in-a-shell
    - command::Vector{String} : Entrypoint array. Not executed within a shell. The docker image's ENTRYPOINT is used if this is not provided. Variable references $(VAR_NAME) are expanded using the container's environment. If a variable cannot be resolved, the reference in the input string will be unchanged. The $(VAR_NAME) syntax can be escaped with a double $$, ie: $$(VAR_NAME). Escaped references will never be expanded, regardless of whether the variable exists or not. Cannot be updated. More info: https://kubernetes.io/docs/tasks/inject-data-application/define-command-argument-container/#running-a-command-in-a-shell
    - env::Vector{IoK8sApiCoreV1EnvVar} : List of environment variables to set in the container. Cannot be updated.
    - envFrom::Vector{IoK8sApiCoreV1EnvFromSource} : List of sources to populate environment variables in the container. The keys defined within a source must be a C_IDENTIFIER. All invalid keys will be reported as an event when the container is starting. When a key exists in multiple sources, the value associated with the last source will take precedence. Values defined by an Env with a duplicate key will take precedence. Cannot be updated.
    - image::String : Docker image name. More info: https://kubernetes.io/docs/concepts/containers/images
    - imagePullPolicy::String : Image pull policy. One of Always, Never, IfNotPresent. Defaults to Always if :latest tag is specified, or IfNotPresent otherwise. Cannot be updated. More info: https://kubernetes.io/docs/concepts/containers/images#updating-images
    - lifecycle::IoK8sApiCoreV1Lifecycle : Lifecycle is not allowed for ephemeral containers.
    - livenessProbe::IoK8sApiCoreV1Probe : Probes are not allowed for ephemeral containers.
    - name::String : Name of the ephemeral container specified as a DNS_LABEL. This name must be unique among all containers, init containers and ephemeral containers.
    - ports::Vector{IoK8sApiCoreV1ContainerPort} : Ports are not allowed for ephemeral containers.
    - readinessProbe::IoK8sApiCoreV1Probe : Probes are not allowed for ephemeral containers.
    - resources::IoK8sApiCoreV1ResourceRequirements : Resources are not allowed for ephemeral containers. Ephemeral containers use spare resources already allocated to the pod.
    - securityContext::IoK8sApiCoreV1SecurityContext : SecurityContext is not allowed for ephemeral containers.
    - startupProbe::IoK8sApiCoreV1Probe : Probes are not allowed for ephemeral containers.
    - stdin::Bool : Whether this container should allocate a buffer for stdin in the container runtime. If this is not set, reads from stdin in the container will always result in EOF. Default is false.
    - stdinOnce::Bool : Whether the container runtime should close the stdin channel after it has been opened by a single attach. When stdin is true the stdin stream will remain open across multiple attach sessions. If stdinOnce is set to true, stdin is opened on container start, is empty until the first client attaches to stdin, and then remains open and accepts data until the client disconnects, at which time stdin is closed and remains closed until the container is restarted. If this flag is false, a container processes that reads from stdin will never receive an EOF. Default is false
    - targetContainerName::String : If set, the name of the container from PodSpec that this ephemeral container targets. The ephemeral container will be run in the namespaces (IPC, PID, etc) of this container. If not set then the ephemeral container is run in whatever namespaces are shared for the pod. Note that the container runtime must support this feature.
    - terminationMessagePath::String : Optional: Path at which the file to which the container's termination message will be written is mounted into the container's filesystem. Message written is intended to be brief final status, such as an assertion failure message. Will be truncated by the node if greater than 4096 bytes. The total message length across all containers will be limited to 12kb. Defaults to /dev/termination-log. Cannot be updated.
    - terminationMessagePolicy::String : Indicate how the termination message should be populated. File will use the contents of terminationMessagePath to populate the container status message on both success and failure. FallbackToLogsOnError will use the last chunk of container log output if the termination message file is empty and the container exited with an error. The log output is limited to 2048 bytes or 80 lines, whichever is smaller. Defaults to File. Cannot be updated.
    - tty::Bool : Whether this container should allocate a TTY for itself, also requires 'stdin' to be true. Default is false.
    - volumeDevices::Vector{IoK8sApiCoreV1VolumeDevice} : volumeDevices is the list of block devices to be used by the container. This is a beta feature.
    - volumeMounts::Vector{IoK8sApiCoreV1VolumeMount} : Pod volumes to mount into the container's filesystem. Cannot be updated.
    - workingDir::String : Container's working directory. If not specified, the container runtime's default will be used, which might be configured in the container image. Cannot be updated.
"""
mutable struct IoK8sApiCoreV1EphemeralContainer <: SwaggerModel
    args::Any # spec type: Union{ Nothing, Vector{String} } # spec name: args
    command::Any # spec type: Union{ Nothing, Vector{String} } # spec name: command
    env::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1EnvVar} } # spec name: env
    envFrom::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1EnvFromSource} } # spec name: envFrom
    image::Any # spec type: Union{ Nothing, String } # spec name: image
    imagePullPolicy::Any # spec type: Union{ Nothing, String } # spec name: imagePullPolicy
    lifecycle::Any # spec type: Union{ Nothing, IoK8sApiCoreV1Lifecycle } # spec name: lifecycle
    livenessProbe::Any # spec type: Union{ Nothing, IoK8sApiCoreV1Probe } # spec name: livenessProbe
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    ports::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1ContainerPort} } # spec name: ports
    readinessProbe::Any # spec type: Union{ Nothing, IoK8sApiCoreV1Probe } # spec name: readinessProbe
    resources::Any # spec type: Union{ Nothing, IoK8sApiCoreV1ResourceRequirements } # spec name: resources
    securityContext::Any # spec type: Union{ Nothing, IoK8sApiCoreV1SecurityContext } # spec name: securityContext
    startupProbe::Any # spec type: Union{ Nothing, IoK8sApiCoreV1Probe } # spec name: startupProbe
    stdin::Any # spec type: Union{ Nothing, Bool } # spec name: stdin
    stdinOnce::Any # spec type: Union{ Nothing, Bool } # spec name: stdinOnce
    targetContainerName::Any # spec type: Union{ Nothing, String } # spec name: targetContainerName
    terminationMessagePath::Any # spec type: Union{ Nothing, String } # spec name: terminationMessagePath
    terminationMessagePolicy::Any # spec type: Union{ Nothing, String } # spec name: terminationMessagePolicy
    tty::Any # spec type: Union{ Nothing, Bool } # spec name: tty
    volumeDevices::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1VolumeDevice} } # spec name: volumeDevices
    volumeMounts::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1VolumeMount} } # spec name: volumeMounts
    workingDir::Any # spec type: Union{ Nothing, String } # spec name: workingDir

    function IoK8sApiCoreV1EphemeralContainer(;args=nothing, command=nothing, env=nothing, envFrom=nothing, image=nothing, imagePullPolicy=nothing, lifecycle=nothing, livenessProbe=nothing, name=nothing, ports=nothing, readinessProbe=nothing, resources=nothing, securityContext=nothing, startupProbe=nothing, stdin=nothing, stdinOnce=nothing, targetContainerName=nothing, terminationMessagePath=nothing, terminationMessagePolicy=nothing, tty=nothing, volumeDevices=nothing, volumeMounts=nothing, workingDir=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1EphemeralContainer, Symbol("args"), args)
        setfield!(o, Symbol("args"), args)
        validate_property(IoK8sApiCoreV1EphemeralContainer, Symbol("command"), command)
        setfield!(o, Symbol("command"), command)
        validate_property(IoK8sApiCoreV1EphemeralContainer, Symbol("env"), env)
        setfield!(o, Symbol("env"), env)
        validate_property(IoK8sApiCoreV1EphemeralContainer, Symbol("envFrom"), envFrom)
        setfield!(o, Symbol("envFrom"), envFrom)
        validate_property(IoK8sApiCoreV1EphemeralContainer, Symbol("image"), image)
        setfield!(o, Symbol("image"), image)
        validate_property(IoK8sApiCoreV1EphemeralContainer, Symbol("imagePullPolicy"), imagePullPolicy)
        setfield!(o, Symbol("imagePullPolicy"), imagePullPolicy)
        validate_property(IoK8sApiCoreV1EphemeralContainer, Symbol("lifecycle"), lifecycle)
        setfield!(o, Symbol("lifecycle"), lifecycle)
        validate_property(IoK8sApiCoreV1EphemeralContainer, Symbol("livenessProbe"), livenessProbe)
        setfield!(o, Symbol("livenessProbe"), livenessProbe)
        validate_property(IoK8sApiCoreV1EphemeralContainer, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(IoK8sApiCoreV1EphemeralContainer, Symbol("ports"), ports)
        setfield!(o, Symbol("ports"), ports)
        validate_property(IoK8sApiCoreV1EphemeralContainer, Symbol("readinessProbe"), readinessProbe)
        setfield!(o, Symbol("readinessProbe"), readinessProbe)
        validate_property(IoK8sApiCoreV1EphemeralContainer, Symbol("resources"), resources)
        setfield!(o, Symbol("resources"), resources)
        validate_property(IoK8sApiCoreV1EphemeralContainer, Symbol("securityContext"), securityContext)
        setfield!(o, Symbol("securityContext"), securityContext)
        validate_property(IoK8sApiCoreV1EphemeralContainer, Symbol("startupProbe"), startupProbe)
        setfield!(o, Symbol("startupProbe"), startupProbe)
        validate_property(IoK8sApiCoreV1EphemeralContainer, Symbol("stdin"), stdin)
        setfield!(o, Symbol("stdin"), stdin)
        validate_property(IoK8sApiCoreV1EphemeralContainer, Symbol("stdinOnce"), stdinOnce)
        setfield!(o, Symbol("stdinOnce"), stdinOnce)
        validate_property(IoK8sApiCoreV1EphemeralContainer, Symbol("targetContainerName"), targetContainerName)
        setfield!(o, Symbol("targetContainerName"), targetContainerName)
        validate_property(IoK8sApiCoreV1EphemeralContainer, Symbol("terminationMessagePath"), terminationMessagePath)
        setfield!(o, Symbol("terminationMessagePath"), terminationMessagePath)
        validate_property(IoK8sApiCoreV1EphemeralContainer, Symbol("terminationMessagePolicy"), terminationMessagePolicy)
        setfield!(o, Symbol("terminationMessagePolicy"), terminationMessagePolicy)
        validate_property(IoK8sApiCoreV1EphemeralContainer, Symbol("tty"), tty)
        setfield!(o, Symbol("tty"), tty)
        validate_property(IoK8sApiCoreV1EphemeralContainer, Symbol("volumeDevices"), volumeDevices)
        setfield!(o, Symbol("volumeDevices"), volumeDevices)
        validate_property(IoK8sApiCoreV1EphemeralContainer, Symbol("volumeMounts"), volumeMounts)
        setfield!(o, Symbol("volumeMounts"), volumeMounts)
        validate_property(IoK8sApiCoreV1EphemeralContainer, Symbol("workingDir"), workingDir)
        setfield!(o, Symbol("workingDir"), workingDir)
        o
    end
end # type IoK8sApiCoreV1EphemeralContainer

const _property_map_IoK8sApiCoreV1EphemeralContainer = Dict{Symbol,Symbol}(Symbol("args")=>Symbol("args"), Symbol("command")=>Symbol("command"), Symbol("env")=>Symbol("env"), Symbol("envFrom")=>Symbol("envFrom"), Symbol("image")=>Symbol("image"), Symbol("imagePullPolicy")=>Symbol("imagePullPolicy"), Symbol("lifecycle")=>Symbol("lifecycle"), Symbol("livenessProbe")=>Symbol("livenessProbe"), Symbol("name")=>Symbol("name"), Symbol("ports")=>Symbol("ports"), Symbol("readinessProbe")=>Symbol("readinessProbe"), Symbol("resources")=>Symbol("resources"), Symbol("securityContext")=>Symbol("securityContext"), Symbol("startupProbe")=>Symbol("startupProbe"), Symbol("stdin")=>Symbol("stdin"), Symbol("stdinOnce")=>Symbol("stdinOnce"), Symbol("targetContainerName")=>Symbol("targetContainerName"), Symbol("terminationMessagePath")=>Symbol("terminationMessagePath"), Symbol("terminationMessagePolicy")=>Symbol("terminationMessagePolicy"), Symbol("tty")=>Symbol("tty"), Symbol("volumeDevices")=>Symbol("volumeDevices"), Symbol("volumeMounts")=>Symbol("volumeMounts"), Symbol("workingDir")=>Symbol("workingDir"))
const _property_types_IoK8sApiCoreV1EphemeralContainer = Dict{Symbol,String}(Symbol("args")=>"Vector{String}", Symbol("command")=>"Vector{String}", Symbol("env")=>"Vector{IoK8sApiCoreV1EnvVar}", Symbol("envFrom")=>"Vector{IoK8sApiCoreV1EnvFromSource}", Symbol("image")=>"String", Symbol("imagePullPolicy")=>"String", Symbol("lifecycle")=>"IoK8sApiCoreV1Lifecycle", Symbol("livenessProbe")=>"IoK8sApiCoreV1Probe", Symbol("name")=>"String", Symbol("ports")=>"Vector{IoK8sApiCoreV1ContainerPort}", Symbol("readinessProbe")=>"IoK8sApiCoreV1Probe", Symbol("resources")=>"IoK8sApiCoreV1ResourceRequirements", Symbol("securityContext")=>"IoK8sApiCoreV1SecurityContext", Symbol("startupProbe")=>"IoK8sApiCoreV1Probe", Symbol("stdin")=>"Bool", Symbol("stdinOnce")=>"Bool", Symbol("targetContainerName")=>"String", Symbol("terminationMessagePath")=>"String", Symbol("terminationMessagePolicy")=>"String", Symbol("tty")=>"Bool", Symbol("volumeDevices")=>"Vector{IoK8sApiCoreV1VolumeDevice}", Symbol("volumeMounts")=>"Vector{IoK8sApiCoreV1VolumeMount}", Symbol("workingDir")=>"String")
Base.propertynames(::Type{ IoK8sApiCoreV1EphemeralContainer }) = collect(keys(_property_map_IoK8sApiCoreV1EphemeralContainer))
Swagger.property_type(::Type{ IoK8sApiCoreV1EphemeralContainer }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1EphemeralContainer[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1EphemeralContainer }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1EphemeralContainer[property_name]

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

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


text:
"""
    DisjoinSet(size::Int)

Disjoint Set Data Structure.
"""
struct DisjointSet
    size::Int
    parent::Array{Int}
    min::Array{Int}
    max::Array{Int}
    DisjointSet(size) = new(size, fill(-1, size), collect(1:size), collect(1:size))
end

"""
    function setUnion!(disjointSet::DisjointSet, representative1::Int, representative2::Int)

Unite the sets of representative 1 and representative 2.
"""
function setUnion!(disjointSet::DisjointSet, representative1::Int, representative2::Int)
    #assert disjointSet.parent[representative1] < 0 && disjointSet.parent[representative1] < 0
    if (disjointSet.parent[representative1] < disjointSet.parent[representative2])
        disjointSet.parent[representative2] = representative1
        disjointSet.max[representative1] = max(disjointSet.max[representative1], disjointSet.max[representative2]);
        disjointSet.min[representative1] = min(disjointSet.min[representative1], disjointSet.min[representative2]);
    elseif (disjointSet.parent[representative2] < disjointSet.parent[representative1])
        disjointSet.parent[representative1] = representative2
        disjointSet.max[representative2] = max(disjointSet.max[representative1], disjointSet.max[representative2]);
        disjointSet.min[representative2] = min(disjointSet.min[representative1], disjointSet.min[representative2]);
    else
        disjointSet.parent[representative1] = representative2
        disjointSet.parent[representative2] = disjointSet.parent[representative2] - 1
        disjointSet.max[representative2] = max(disjointSet.max[representative1], disjointSet.max[representative2]);
        disjointSet.min[representative2] = min(disjointSet.min[representative1], disjointSet.min[representative2]);
    end

end

"""
    function findRepresentative!(disjointSet::DisjointSet, element::Int)

Find the representative of the element in the set.
"""
function findRepresentative!(disjointSet::DisjointSet, element::Int)
    root = element
    while (disjointSet.parent[root] >= 0)
        root = disjointSet.parent[root]
    end
    while (element != root)
        temporary_element = disjointSet.parent[element]
        disjointSet.parent[element] = root
        element = temporary_element
    end
    return root
end

"""
    function greatest!(disjointSet::DisjointSet, element::Int)

Find the gretest element in the set of element.
"""
function greatest!(disjointSet::DisjointSet, element::Int)
    representative = findRepresentative!(disjointSet, element)
    return disjointSet.max[representative]
end

"""
    function smallest!(disjointSet::DisjointSet, element::Int)

Find the smallest element in the set of element.
"""
function smallest!(disjointSet::DisjointSet, element::Int)
    representative = findRepresentative!(disjointSet, element)
    return disjointSet.min[representative]
end


text:
function observables(R, xroot, lroot, ϵroot, κ, norm, coef, energ, j, mass)

    @info "Extracting observables..."

    sinϕ = @. mass/√(energ^2+mass^2)
    cosϕ = @. √(1-sinϕ^2)
    vnjl² = @. (1-energ/√(energ^2+mass^2))/2

    # Average mass computed with vnjl² and sinϕ
    Mvavg = sum(@. mass*j*vnjl²)/sum(@. j*vnjl²)
    Msavg = sum(@. mass*j*sinϕ) / sum(@. j*sinϕ)

    # Average mass computed as an integral
    Mcoef(r) = sum(@. 1/2*(2*j+1)/(4*π)*sinϕ*norm^2*(sphericalbessel(xroot/R*r, lroot)^2+sphericalbessel(xroot/R*r, lroot+ϵroot)^2))
    M(r) = 4*Nc*Nf*G*ħ^3*Mcoef(r)
    Mint(r) = 4*π*r^2*M(r)
    Mavg = quadgk(Mint, 0, R)[1]
    Mavg = Mavg / (4*π*R^3/3)

    # Chiral condensate extracted from the radial mass
    # For dubious reasons, for certain values it returns a complex value with very small imaginary part...
    ψ̅ψ(r) = - complex(1/(2*G) * M(r))^(1/3)
    # ψ̅ψ(r) = - (1/(2*G) * M(r))^(1/3)
    ψ̅ψint(r) = 4*π*r^2*ψ̅ψ(r)
    ψ̅ψavg = quadgk(ψ̅ψint, 0, R)[1]
    ψ̅ψavg = ψ̅ψavg / (4*π*R^3/3)

    # Trying to brute force get rid of the imaginary part...
    if imag(ψ̅ψavg) != 0
        ψ̅ψavg = NaN
    else
        ψ̅ψavg = real(ψ̅ψavg)
    end

    # Radial dependence of mass
    Mr = r -> M(r)
    # Radial dependence of chiral condensate
    qq̅ = r -> ψ̅ψ(r)

    # Energy density of the chiral condensate
    ε₁ = - sum(@. j*energ*cosϕ)
    ε₂ = - sum([(coef[k₁, k₂]*sinϕ[k₁]*sinϕ[k₂]) for k₁ in 1:κ for k₂ in 1:κ])
    ε₂ = ε₂ *2*Nc*Nf*G*ħ^3
    ε = (ε₁ + ε₂) / (4*π*R^3/3)

    @info "Done!"
    @debug "Average mass: " Mavg
    @debug "Average value of chiral condensate: " ψ̅ψavg
    @debug "Energy density of a chiral pair: " ε

    return Mvavg, Msavg, Mr, Mavg, qq̅, ψ̅ψavg, ε

end

max_stars_repo_name: avramescudana/QuarkDroplets

max_issues_repo_name: avramescudana/QuarkDroplets

max_forks_repo_name: avramescudana/QuarkDroplets


text:
# This file is part of the IntervalArithmetic.jl package; MIT licensed

function sinh(a::Interval)
    isempty(a) && return a

    return @round(sinh(a.lo), sinh(a.hi))
end

function cosh(a::Interval)
    isempty(a) && return a

    return @round(cosh(mig(a)), cosh(mag(a)))
end

# function tanh(a::Interval{Float64})
#     isempty(a) && return a
#
#     float(tanh(big53(a)))
# end

function tanh(a::Interval{BigFloat})
    isempty(a) && return a

    return @round(tanh(a.lo), tanh(a.hi))
end

function asinh(a::Interval{BigFloat})
    isempty(a) && return a

    return @round(asinh(a.lo), asinh(a.hi))
end

function acosh(a::Interval{BigFloat})
    domain = Interval(one(eltype(a)), Inf)
    a = a ∩ domain
    isempty(a) && return a

    return @round(acosh(a.lo), acosh(a.hi))
end

function atanh(a::Interval{BigFloat})
    domain = Interval(-one(eltype(a)), one(eltype(a)))
    a = a ∩ domain

    isempty(a) && return a

    res_lo = @round_down(atanh(a.lo))
    res_hi = @round_up(atanh(a.hi))

    # The IEEE Std 1788-2015 does not allow intervals like of the
    # form Interval(∞,∞) and Interval(-∞,-∞)
    (res_lo == res_hi == Inf || res_lo == res_hi == -Inf) && return emptyinterval(a)  # IEEE 1788-2015 does not allow intervals consisting only of ∞, i.e. Interval(∞,∞) and Interval(-∞,-∞)

    return Interval(res_lo, res_hi)
end

# Float64 versions of functions missing from CRlibm:
for f in (:tanh, :asinh, :acosh, :atanh)

    @eval function ($f)(a::Interval{Float64})
        isempty(a) && return a

        atomic(Interval{Float64}, ($f)(big53(a)) )
    end
end


text:
using PackageCompiler

compile_package("Poisson")


text:
include("network.jl")

struct Discriminator
    initial
    model
end

@functor Discriminator
function Discriminator(
    in_channels::Int = 3,
    features::Any = [64, 128, 256, 512]
)
    layers = []
    channel = features[1]
    for index in range(2, length(features))
            if features[index] != last(features)
                push!(layers,
                    Block(channel, features[index], stride=2)
                )
            else
                push!(layers,
                    Block(channel, features[index], stride=1)
                )
            end
            channel = features[index]
    end
    push!(layers,
        Conv((4,4), channel => 1; stride=1, pad=1)
    )
return Discriminator(
    Chain(
        Conv((4,4), in_channels=>features[1]; stride=2, pad=1),
        x -> leakyrelu.(x, 0.2)
    ),
    Chain(layers...)
) |> gpu
end

function (net::Discriminator)(x)
    input = net.initial(x)
    return sigmoid.(net.model(input))
end

using Random
function test()
    img_channels = 3
    img_size = 100
    ## need to explicity type to avoid Slow fallback implementation
    ## https://discourse.julialang.org/t/flux-con-warning/49456
    x = randn(Float32, (img_size, img_size, img_channels, 5))
    preds = Discriminator()
    println(size(preds(x)))
end


text:
# This file is a part of SimilaritySearch.jl

using Random

"""
    BeamSearch(bsize::Integer=16, Δ::Float32)

BeamSearch is an iteratively improving local search algorithm that explores the graph using blocks of `bsize` elements and neighborhoods at the time.

- `bsize`: The size of the beam.
- `Δ`: Soft margin for accepting elements into the beam
- `maxvisits`: MAximum visits while searching, useful for early stopping without convergence
"""
@with_kw mutable struct BeamSearch <: LocalSearchAlgorithm
    bsize::Int32 = 8  # size of the search beam
    Δ::Float32 = 1.0  # soft-margin for accepting an element into the beam
    maxvisits::Int64 = 1000_000 # maximum visits by search, useful for early stopping without convergence, very high by default
end

Base.copy(bsearch::BeamSearch; bsize=bsearch.bsize, Δ=bsearch.Δ, maxvisits=bsearch.maxvisits) =
    BeamSearch(; bsize, Δ, maxvisits)

const GlobalBeamKnnResult = [KnnResultShift(32)]  # see __init__ function

function __init__beamsearch()
    for i in 2:Threads.nthreads()
        push!(GlobalBeamKnnResult, KnnResultShift(32))
    end
end

### local search algorithm

function beamsearch_queue(index::SearchGraph, q, res, objID, vstate)
    visited_ = 0
    @inbounds if !visited(vstate, objID)
        visit!(vstate, objID)
        visited_ += 1
        d = evaluate(index.dist, q, index[objID])
        push!(res, objID, d)
    end

    visited_
end

function beamsearch_init(bs::BeamSearch, index::SearchGraph, q, res, hints, vstate, bsize)
    visited_ = 0

    for objID in hints
        visited_ += beamsearch_queue(index, q, res, objID, vstate)
    end

    if length(res) == 0
        _range = 1:length(index)
        for i in 1:bsize
           objID = rand(_range)
           visited_ += beamsearch_queue(index, q, res, objID, vstate)
       end
    end

    visited_
end

function beamsearch_inner(bs::BeamSearch, index::SearchGraph, q, res, vstate, beam, Δ, maxvisits, visited_)
    beam_st = initialstate(beam)
    beam_st = push!(beam, beam_st, argmin(res), minimum(res))
    # @show res.id length(res, st) length(index) res st beam_st
    while length(beam, beam_st) > 0
        p, beam_st = popfirst!(beam, beam_st)
        prev_id = p.first

        @inbounds for childID in index.links[prev_id]
            if !visited(vstate, childID)
                visit!(vstate, childID)
                d = evaluate(index.dist, q, index[childID])
                push!(res, childID, d)
                visited_ += 1
                visited_ > maxvisits && @goto finish_search
                if d <= Δ * maximum(res)
                    #if length(index.links[childID]) > 1
                    beam_st = push!(beam, beam_st, childID, d)
                    # length(beam) == maxlength(beam) && continue
                    # sat_should_push(keys(beam), index, q, childID, d) && push!(beam, childID, d)
                    #end
                end
            end
        end
    end

    @label finish_search
    (res=res, cost=visited_)
end

"""
    search(bs::BeamSearch, index::SearchGraph, q, res, hints, pools; bsize=bs.bsize, Δ=bs.Δ, maxvisits=bs.maxvisits)

Tries to reach the set of nearest neighbors specified in `res` for `q`.
- `bs`: the parameters of `BeamSearch`
- `index`: the local search index
- `q`: the query
- `res`: The result object, it stores the results and also specifies the kind of query
- `hints`: Starting points for searching, randomly selected when it is an empty collection
- `pools`: A SearchGraphPools object with preallocated pools

Optional arguments (defaults to values in `bs`)
- `bsize`: Beam size
- `Δ`: exploration expansion factor
- `maxvisits`: Maximum number of nodes to visit (distance evaluations)

"""
function search(bs::BeamSearch, index::SearchGraph, q, res, hints, pools::SearchGraphPools; bsize=bs.bsize, Δ=bs.Δ, maxvisits=bs.maxvisits, vstate=getvstate(length(index), pools))
    # k is the number of neighbors in res
    visited_ = beamsearch_init(bs, index, q, res, hints, vstate, bsize)
    beam = getbeam(bsize, pools)
    beamsearch_inner(bs, index, q, res, vstate, beam, Δ, maxvisits, visited_)
end

max_stars_repo_name: sadit/NNS

max_issues_repo_name: sadit/NNS

max_forks_repo_name: sadit/NNS


text:
# This file is a part of GeDetROOT.jl, licensed under the MIT License (MIT).

__precompile__(false)

module GeDetROOT

include("sis3316_to_root.jl")

end # module


text:
function secant(f::Function, xi::Number, dx::Number, ea::Number, maxIterations::Integer)
  # Estimate the root of functions using the secant method.
  #
  # Keyword arguments:
  #  f -- Anonymous function to find the root of
  #  xi -- Initial estimate for the root
  #  dx -- Pertubation to use for estimating df/dx
  #  ea -- Minimum relative error for a solution
  #  maxIterations -- Maximum number of iterations to carry out
  #
  # Returns:
  #  An estimate of the root of the function

  xo = xi; # Current estimate of root

  for i = 1:maxIterations
    xoLast = xo;
    xo = xo - (f(xo) * dx) / (f(xo+dx) - f(xo));
    relativeError = abs((xo - xoLast) / xo) * 100;

    (relativeError < ea || isnan(relativeError)) && break;
  end

  return xo;
end

max_stars_repo_path: Julia/Root Estimation/secant.jl

max_stars_repo_name: alexjohnj/numerical-methods

max_issues_repo_path: Julia/Root Estimation/secant.jl

max_issues_repo_name: alexjohnj/numerical-methods

max_forks_repo_path: Julia/Root Estimation/secant.jl

max_forks_repo_name: alexjohnj/numerical-methods


text:
#  Note:  evalVolume is tested under test_linearFFD.jl since it is needed
#         everytime a control point is perturbed to compute the new location of
#         the geometric point

@testset "--- Checking B-spline Formulation ---" begin

  U = [0,0,0,0,2/5,3/5,3/5,1,1,1,1]
  order = 4
  p = order - 1
  u = 0.42
  nctl = length(U) - order

  @test ( nctl )== 7

  @testset "Checking knot span index evaluation" begin
    span = FFD.findSpan(u, U, order, nctl)
    @test ( span )== 5
  end

  @testset "Checking basis function evaluations" begin
    N = zeros(order)
    span = FFD.findSpan(u, U, order, nctl)
    FFD.basisFunctions(U, order, u, span, N)
    @test isapprox( N[1], 0.081) atol= 1e-15
    @test isapprox( N[2], 0.405) atol= 1e-15
    @test isapprox( N[3], 0.5136666666666666) atol= 1e-15
    @test isapprox( N[4], 0.00033333333333333153) atol= 1e-15
  end

  @testset "Checking curve point evaluation" begin
    P = 0:1/(nctl-1):1
    C = [0.0]
    FFD.evalCurve([u], U, order, P, C)
    @test isapprox( C[1], 0.40555555555555556) atol= 1e-16
  end

end  # End facts("--- Check Basis Function Calculations ---")

@testset "--- Checking B-spline Derivatives ---" begin

  @testset "Checking Basis Function Derivatives" begin

    U = [0,0,0,0,2/5,3/5,3/5,1,1,1,1]
    order = 4
    p = order - 1
    u = 0.42
    nctl = length(U) - order
    N = zeros(AbstractFloat, order)
    Nderiv = zeros(AbstractFloat, order)
    span = FFD.findSpan(u, U, order, nctl)
    FFD.derivBasisFunctions(u, U, order, span, N, Nderiv)

    @test isapprox( N[1], 0.081) atol= 1e-15
    @test isapprox( N[2], 0.405) atol= 1e-15
    @test isapprox( N[3], 0.5136666666666666) atol= 1e-15
    @test isapprox( N[4], 0.00033333333333333153) atol= 1e-15
    @test isapprox( Nderiv[1], -1.3500000000000003) atol= 1e-15
    @test isapprox( Nderiv[2], -2.25) atol= 1e-15
    @test isapprox( Nderiv[3], 3.55) atol= 1e-15
    @test isapprox( Nderiv[4], 0.04999999999999982) atol= 1e-15

  end  # End context("Checking Basis Function Derivatives")

end  # End facts("--- Checking B-spline Derivatives ---")


text:
# another tree test

using IncrementalInference
using Test

@testset "second tree formation test" begin

fg = initfg()

addVariable!(fg, :x1, ContinuousScalar)
addVariable!(fg, :x2, ContinuousScalar)
addVariable!(fg, :x3, ContinuousScalar)
addVariable!(fg, :x4, ContinuousScalar)
addVariable!(fg, :x5, ContinuousScalar)
addVariable!(fg, :l1, ContinuousScalar)
addVariable!(fg, :l2, ContinuousScalar)
addVariable!(fg, :l3, ContinuousScalar)

addFactor!(fg, [:x1;:l1], LinearConditional(Normal()))
addFactor!(fg, [:x1;:x2], LinearConditional(Normal()))
addFactor!(fg, [:x2;:l1], LinearConditional(Normal()))
addFactor!(fg, [:x2;:x3], LinearConditional(Normal()))
addFactor!(fg, [:x3;:x4], LinearConditional(Normal()))
addFactor!(fg, [:x4;:l2], LinearConditional(Normal()))
addFactor!(fg, [:x4;:x5], LinearConditional(Normal()))
addFactor!(fg, [:l2;:x5], LinearConditional(Normal()))
addFactor!(fg, [:x4;:l3], LinearConditional(Normal()))
addFactor!(fg, [:x5;:l3], LinearConditional(Normal()))

writeGraphPdf(fg, show=true)

eo = [:x1; :l3; :l1; :x5; :x2; :l2; :x4; :x3]

tree = buildTreeFromOrdering!(fg,eo)

# using Gadfly, Fontconfig, Cairo
# drawTree(tree, show=true, imgs=false)
# drawTree(tree, show=true, imgs=true)
#

end

## Comparible reference available from gtsam/testSymbolicEliminationTree.cpp

#        l1                  l2
#    /    |                /  |
# x1 --- x2 --- x3 --- x4 --- x5
#                          \  |
#                            l3

# (X(1), L(1));
# (X(1), X(2));
# (X(2), L(1));
# (X(2), X(3));
# (X(3), X(4));
# (X(4), L(2));
# (X(4), X(5));
# (L(2), X(5));
# (X(4), L(3));
# (X(5), L(3));

# EliminationTreeTester::MakeTree(list_of
#     (MakeNode(X(3), SymbolicFactorGraph(), list_of
#       (MakeNode(X(2), list_of(SymbolicFactor(X(2), X(3))), list_of
#         (MakeNode(L(1), list_of(SymbolicFactor(X(2), L(1))), list_of
#           (MakeNode(X(1), list_of(SymbolicFactor(X(1), L(1))) (SymbolicFactor(X(1), X(2)))))))))
#       (MakeNode(X(4), list_of(SymbolicFactor(X(3), X(4))), list_of
#         (MakeNode(L(2), list_of(SymbolicFactor(X(4), L(2))), list_of
#           (MakeNode(X(5), list_of(SymbolicFactor(X(4), X(5))) (SymbolicFactor(L(2), X(5))), list_of
#             (MakeNode(L(3), list_of(SymbolicFactor(X(4), L(3))) (SymbolicFactor(X(5), L(3))))))))))) )));

#


text:
module ImageDraw

# package code goes here
using ImageCore, Distances

include("core.jl")
include("line2d.jl")
include("ellipse2d.jl")
include("circle2d.jl")
include("paths.jl")
include("cross.jl")

#export methods
export
	draw,
	draw!,
	bresenham,
	xiaolin_wu

#export types
export
	#Drawable
	Drawable,

	#Point
	Point,

	#Lines
	LineNormal,
	LineTwoPoints,

	#LineSegment
	LineSegment,

	#Ellipse
	Ellipse,

	#Circles
	CirclePointRadius,
	CircleThreePoints,

	#Path
	Path,

	#Polygons
	Polygon,
	RegularPolygon,

	#Cross
	Cross

end # module


text:
module TestStaticArrays

using Mutabilities
using Mutabilities: isfrozen, isfrozenindex, isfrozenvalue
using StaticArrays
using Test

@testset "isfrozen" begin
    @test isfrozen(SVector(1))
    @test isfrozenindex(SVector(1))
    @test isfrozenvalue(SVector(1))

    @test !isfrozen(MVector(1))
    @test isfrozenindex(MVector(1))
    @test !isfrozenvalue(MVector(1))
end

@testset "melt" begin
    @test typeof(melt(SVector(1))) <: Vector{Int}
    @test typeof(melt(MVector(1))) <: Vector{Int}
    @test typeof(melt(SMatrix{1,1}(1))) <: MMatrix
    @test typeof(melt(MMatrix{1,1}(1))) <: MMatrix
end

@testset "meltvalue" begin
    @test typeof(meltvalue(SVector(1))) <: MVector
    @test typeof(meltvalue(MVector(1))) <: MVector
    @test typeof(meltvalue(SMatrix{1,1}(1))) <: MMatrix
    @test typeof(meltvalue(MMatrix{1,1}(1))) <: MMatrix
end

@testset "meltindex" begin
    @test typeof(meltindex(SVector(1))) <: Vector{Int}
    @test typeof(meltindex(MVector(1))) <: Vector{Int}
    @test typeof(meltindex(SMatrix{1,1}(1))) <: MMatrix
    @test typeof(meltindex(MMatrix{1,1}(1))) <: MMatrix
end

end  # module


text:
# Julia wrapper for header: /usr/include/scip/scip_solve.h
# Automatically generated using Clang.jl wrap_c

function SCIPtransformProb(scip)
    ccall((:SCIPtransformProb, libscip), SCIP_RETCODE, (Ptr{SCIP_},), scip)
end

function SCIPpresolve(scip)
    ccall((:SCIPpresolve, libscip), SCIP_RETCODE, (Ptr{SCIP_},), scip)
end

function SCIPsolve(scip)
    ccall((:SCIPsolve, libscip), SCIP_RETCODE, (Ptr{SCIP_},), scip)
end

function SCIPsolveParallel(scip)
    ccall((:SCIPsolveParallel, libscip), SCIP_RETCODE, (Ptr{SCIP_},), scip)
end

function SCIPsolveConcurrent(scip)
    ccall((:SCIPsolveConcurrent, libscip), SCIP_RETCODE, (Ptr{SCIP_},), scip)
end

function SCIPfreeSolve(scip, restart)
    ccall((:SCIPfreeSolve, libscip), SCIP_RETCODE, (Ptr{SCIP_}, UInt32), scip, restart)
end

function SCIPfreeReoptSolve(scip)
    ccall((:SCIPfreeReoptSolve, libscip), SCIP_RETCODE, (Ptr{SCIP_},), scip)
end

function SCIPfreeTransform(scip)
    ccall((:SCIPfreeTransform, libscip), SCIP_RETCODE, (Ptr{SCIP_},), scip)
end

function SCIPinterruptSolve(scip)
    ccall((:SCIPinterruptSolve, libscip), SCIP_RETCODE, (Ptr{SCIP_},), scip)
end

function SCIPrestartSolve(scip)
    ccall((:SCIPrestartSolve, libscip), SCIP_RETCODE, (Ptr{SCIP_},), scip)
end

function SCIPenableReoptimization(scip, enable)
    ccall((:SCIPenableReoptimization, libscip), SCIP_RETCODE, (Ptr{SCIP_}, UInt32), scip, enable)
end

function SCIPisReoptEnabled(scip)
    ccall((:SCIPisReoptEnabled, libscip), UInt32, (Ptr{SCIP_},), scip)
end

function SCIPgetReoptSolsRun(scip, run, sols, allocmem, nsols)
    ccall((:SCIPgetReoptSolsRun, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cint, Ptr{Ptr{SCIP_SOL}}, Cint, Ptr{Cint}), scip, run, sols, allocmem, nsols)
end

function SCIPresetReoptSolMarks(scip)
    ccall((:SCIPresetReoptSolMarks, libscip), Cvoid, (Ptr{SCIP_},), scip)
end

function SCIPcheckReoptRestart(scip, node, restart)
    ccall((:SCIPcheckReoptRestart, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_NODE}, Ptr{UInt32}), scip, node, restart)
end

function SCIPaddReoptDualBndchg(scip, node, var, newbound, oldbound)
    ccall((:SCIPaddReoptDualBndchg, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_NODE}, Ptr{SCIP_VAR}, Cdouble, Cdouble), scip, node, var, newbound, oldbound)
end

function SCIPgetReoptLastOptSol(scip)
    ccall((:SCIPgetReoptLastOptSol, libscip), Ptr{SCIP_SOL}, (Ptr{SCIP_},), scip)
end

function SCIPgetReoptOldObjCoef(scip, var, run, objcoef)
    ccall((:SCIPgetReoptOldObjCoef, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Cint, Ptr{Cdouble}), scip, var, run, objcoef)
end

function SCIPisInRestart(scip)
    ccall((:SCIPisInRestart, libscip), UInt32, (Ptr{SCIP_},), scip)
end


text:
@testset "Fast cis" begin
    @test @inferred(GNSSSignals.sin_fast(π / 2)) ≈ sin(π / 2)
    @test @inferred(GNSSSignals.sin_fast(-π / 2)) ≈ sin(-π / 2)
    @test @inferred(GNSSSignals.sin_fast(0.0)) ≈ sin(0.0)
    @test @inferred(GNSSSignals.sin_fast(1π)) ≈ sin(π) atol = eps(Float64)
    @test @inferred(GNSSSignals.sin_fast(-1π)) ≈ sin(-π) atol = eps(Float64)

    @test @inferred(GNSSSignals.cos_fast(π / 2)) ≈ cos(π / 2) atol = eps(Float64)
    @test @inferred(GNSSSignals.cos_fast(-π / 2)) ≈ cos(-π / 2) atol = eps(Float64)
    @test @inferred(GNSSSignals.cos_fast(0.0)) ≈ cos(0.0)
    @test @inferred(GNSSSignals.cos_fast(1π)) ≈ cos(π)
    @test @inferred(GNSSSignals.cos_fast(-1π)) ≈ cos(-π)

    x = -π:0.01:π
    @test sqrt(sum(abs2.(GNSSSignals.sin_vfast.(x) .- sin.(x))) / length(x)) < 0.0359
    @test sqrt(sum(abs2.(GNSSSignals.sin_fast.(x) .- sin.(x))) / length(x)) < 0.000597

    @test sqrt(sum(abs2.(GNSSSignals.cos_vfast.(x) .- cos.(x))) / length(x)) < 0.03589
    @test sqrt(sum(abs2.(GNSSSignals.cos_fast.(x) .- cos.(x))) / length(x)) < 0.000597

    @test sqrt(sum(abs2.(GNSSSignals.cis_vfast.(x) .- cis.(x))) / length(x)) < 0.0507
    @test sqrt(sum(abs2.(GNSSSignals.cis_fast.(x) .- cis.(x))) / length(x)) < 0.00085

    @test get_carrier_fast_unsafe.(x) == GNSSSignals.cis_fast.(x)
    @test get_carrier_vfast_unsafe.(x) == GNSSSignals.cis_vfast.(x)
end


text:
include("OSMPBF.jl")

using ProtoBuf, CodecZlib

const BLOCK = Dict(
    "OSMHeader" => OSMPBF.HeaderBlock(),
    "OSMData" => OSMPBF.PrimitiveBlock(),
)

function parseblob(io::IO)
    len = ntoh(read(io, Int32))
    h = read(io, len)
    bh = readproto(IOBuffer(h), OSMPBF.BlobHeader())
    blob_buf = IOBuffer(read(io, bh.datasize))
    b = readproto(blob_buf, OSMPBF.Blob())
    ok = transcode(ZlibDecompressor, b.zlib_data)
    if length(ok) != b.raw_size
        @warn("Uncompressed size $(length(ok)) bytes does not match $(b.raw_size).")
    end
    readproto(IOBuffer(ok), BLOCK[bh._type])
end

function process_block(osm::OSMData, block::OSMPBF.HeaderBlock)
    # We don't currently use that information in `OSMData`
    box = block.bbox
    osm.bounds = Bounds(box.bottom / 1e9, box.top / 1e9, box.left / 1e9, box.right / 1e9)
    return
end

function geo(offset, granularity, int)
    return 0.000000001 * (offset + (granularity * int))
end

function tag(osm, element, table, key_id, value_id)
    tag(osm, element, table[key_id + 1], table[value_id + 1])
end

function process_elements(osm::OSMData, nodes::OSMPBF.DenseNodes, table, lat_offset, lon_offset, granularity)
    ids = nodes.id
    cumsum!(ids, ids)
    lat = nodes.lat
    cumsum!(lat, lat)
    lon = nodes.lon
    cumsum!(lon, lon)
    for i in eachindex(ids)
        @assert !haskey(osm.nodes, ids[i])
        osm.nodes[ids[i]] = LLA(
            geo(lat_offset, granularity, lat[i]),
            geo(lon_offset, granularity, lon[i]),
        )
    end
    i = 1
    j = 1
    keys_vals = nodes.keys_vals
    while j <= length(keys_vals)
        if iszero(keys_vals[j])
            i += 1
            j += 1
        else
            tag(osm, ids[i], table, keys_vals[j], keys_vals[j+1])
            j += 2
        end
    end
end

function process_element(osm, pbf_node::OSMPBF.Node, table, lat_offset, lon_offset, granularity)
    id = pbf_node.id
    osm.nodes[id] = LLA(
        geo(lat_offset, granularity, pbf_node.lat),
        geo(lon_offset, granularity, pbf_node.lon),
    )
    keys = pbf_node.keys
    vals = pbf_node.vals
    for i in eachindex(keys)
        tag(osm, id, table, keys[i], vals[i])
    end
end

function process_element(osm, pbf_way::OSMPBF.Way, table, lat_offset, lon_offset, granularity)
    way = Way(pbf_way.id)
    # `refs` is Δ coded
    refs = pbf_way.refs
    cumsum!(refs, refs)
    way.nodes = refs
    if isdefined(pbf_way, :keys)
        keys = pbf_way.keys
        vals = pbf_way.vals
        for i in eachindex(keys)
            tag(osm, way, table, keys[i], vals[i])
        end
    end
    push!(osm.ways, way)
    return
end

function process_element(osm, pbf_relation::OSMPBF.Relation, table, lat_offset, lon_offset, granularity)
    relation = Relation(pbf_relation.id)
    memids = pbf_relation.memids
    cumsum!(memids, memids)
    types = pbf_relation.types
    for i in eachindex(types)
        push!(relation.members, Dict(
            "type" => OSMPBF.Relation_MemberType[types[i] + 1],
            "ref" => string(memids[i]),
        ))
    end
    if isdefined(pbf_relation, :keys)
        keys = pbf_relation.keys
        vals = pbf_relation.vals
        for i in eachindex(keys)
            tag(osm, relation, table, keys[i], vals[i])
        end
    end
    push!(osm.relations, relation)
    return
end

function process_elements(osm, elements::Vector, table, lat_offset, lon_offset, granularity)
    for element in elements
        process_element(osm, element, table, lat_offset, lon_offset, granularity)
    end
end

function process_block(osm, block::OSMPBF.PrimitiveBlock)
    table = String.(block.stringtable.s)
    granularity = block.granularity
    lat_offset = block.lat_offset
    lon_offset = block.lon_offset
    for group in block.primitivegroup
        for key in fieldnames(OSMPBF.PrimitiveGroup)
            if isdefined(group, key)
                process_elements(osm, getproperty(group, key), table, lat_offset, lon_offset, granularity)
            end
        end
    end
end

function parsePBF(filename::AbstractString)
    osm = OSMData()
    open(filename) do io
        while !eof(io)
            block = parseblob(io)
            process_block(osm, block)
        end
    end
    return osm
end


text:
using .Clustering
using .Plots
using .Colors
using .LinearAlgebra

function clustercovariance(Σ²)
    issymmetric(Σ²) || (Σ² = cov(Σ²'))
    Dr = 1.0.-abs.(Σ²)
    if !issymmetric(Dr)
        @warn "Correlation distance matrix is not symmetric, so not clustering"
    end
    Clustering.hclust(Dr; linkage=:average, branchorder=:optimal)
end

"""
    covarianceimage(f, F; [palette=[:cornflowerblue, :crimson, :forestgreen], colormode=:top, colorbargrad=:binary, donames=true, kwargs...])
Plot the covariance matrix of the columns of an `Array`, coloring elements by their contribution to each of the top 3 principal components.
Either provide as positional arguments a vector `f` of N row names and an N×_ matrix `F`, or just a `Catch22.FeatureArray`.

# Keyword Arguments
- `palette`: a vector containing a color for each principal component.
- `colorbargrad`: the color gradient for the colorbar.
- `donames`: override `f` and display no tick labels.
- `colormode`: how to color the covariance matrix. `:raw` gives no coloring by principal components, `:top` is a combination of the top three PC colors and `:all` is a combination of all PC colors, where PCN = :black if N > length(palette).
- `docluster`: whether to perform clustering on the covariance matrix.
- `kwargs...`: other `Plots` attributes.
"""
covarianceimage
@userplot CovarianceImage
Plots.@recipe function f(g::CovarianceImage; palette=[:cornflowerblue, :crimson, :forestgreen], colormode=:top, colorbargrad=:binary, donames=true, docluster=true, verbose=true, dendrogram=false)
    if g.args[1] isa AbstractFeatureArray || g.args[1] isa AbstractDimArray
        f, Σ² = string.(getdim(g.args[1], 1)), g.args[1] |> Array
    elseif length(g.args) == 2 && g.args[2] isa AbstractMatrix
        f, Σ² = g.args[1], g.args[2] |> Array
    else
        @error "Incorrect arguments; give row names and a matrix or an annotated DimArray"
    end

    issymmetric(Σ²) || (Σ² = cov(Σ²'))

    any(diag(Σ²) .== 0) && @warn "Covariance matrix is not positive definite, which may cause an error"
    linecolor --> nothing
    if docluster == true
        idxs = clustercovariance(Σ²).order
    elseif docluster isa Union{AbstractVector, Tuple}
        idxs = docluster # Precomputed indices
    else
        idxs = 1:size(Σ², 1)
    end
    Σ̂² = Σ²[idxs, idxs]
    A = abs.(Σ̂²)./max(abs.(Σ̂²)...)
    f̂ = f[idxs]
    N = min(length(palette), size(Σ̂², 1))
    if colormode == :raw # * Don't color by PC's
        H = abs.(Σ̂²)
        colorbar --> true
    else
        λ = (eigvals∘Symmetric∘Array)(Σ̂²)
        λi = sortperm(abs.(λ), rev=true)
        λ = λ[λi]
        P = (eigvecs∘Symmetric∘Array)(Σ̂²)[:, λi] # Now sorted by decreasing eigenvalue norm
        vidxs = sortperm(abs.(P[:, 1]), rev=true)
        verbose && isnothing(printstyled("Feature weights:\n", color=:red, bold=true)) && display(vcat(hcat("Feature", ["PC$i" for i ∈ 1:N]...) , hcat(f̂[vidxs], round.(P[vidxs, 1:N], sigdigits=3))))
        P = abs.(P)
        if colormode == :top # * Color by the number of PC's given by the length of the color palette
            P = P[:, 1:N]
            P̂ = P.^2.0./sum(P.^2.0, dims=2)
            # Square the loadings, since they are added in quadrature. Maybe not a completely faithful representation of the PC proportions, but should get the job done.
            𝑓′ = parse.(Colors.XYZ, palette[1:N]);
        elseif colormode == :all # * Color by all PC's. This can end up very brown
            Σ̂′² = Diagonal(abs.(λ))
            P̂ = P.^2.0./sum(P.^2.0, dims=2)
            p = fill(:black, size(P, 2))
            p[1:N] = palette[1:N]
            𝑓′ = parse.(Colors.XYZ, p);
            [𝑓′[i] = Σ̂′²[i, i]*𝑓′[i] for i ∈ 1:length(𝑓′)]
        end
        𝑓 = Vector{eltype(𝑓′)}(undef, size(P̂, 1))
        try # Load colors by PC weights
            𝑓 = P̂*𝑓′
        catch
            # Equivalent but slower
            @info "Iterating to load covariances"
            for ii ∈ 1:length(𝑓)
                𝑓[ii] = sum([P̂[ii, jj]*𝑓′[jj] for jj ∈ 1:length(𝑓′)])
            end
        end

        H = Array{Colors.XYZA}(undef, size(Σ̂²))
        for (i, j) ∈ Tuple.(CartesianIndices(H)) # Apply the correlations as transparencies
            J = (𝑓[i] + 𝑓[j])/2
            H[i, j] = Colors.XYZA(J.x, J.y, J.z, A[i, j])
        end
        H = convert.((Colors.RGBA,), H)
    end
    @series begin
        seriestype := :heatmap
        if colormode == :raw
            seriescolor --> colorbargrad
        end
        if backend() == Plots.GRBackend() # For some reason, GR does heatmaps differently
            xs = 0.5:1:size(H, 1)+0.5
        else
            xs = 1:size(H, 1)
        end
        (xs, xs, H)
    end
    # Plot dummy data and set attributes
    @series begin
        if dendrogram
            colorbar := false
        else
            colorbar --> true
        end
        seriestype := :scatter
        markersize := 0.0
        label := nothing
        legend := :none
        marker_z := [0, max(abs.(Σ²)...)]
        markercolor := colorbargrad
        (zeros(2), zeros(2))
    end
    for i ∈ 1:N
        @series begin
            seriestype := :shape
            if colormode != :raw
                label := "PC$i"
                legend --> :topright
                line_width := 20
            else
                label := nothing
                legend := nothing
            end
            colorbar_title --> "\$|\\Sigma^2|\$"
            colorbar_titlefontsize := 14
            xticks := :none
            size --> (800, 400)
            yflip --> true
            lims := (0.5, size(H, 1)+0.5)
            aspect_ratio := :equal
            legendfontsize := 8
            if donames
                yticks := (1:size(H, 1), replace.(string.(f̂), (r"_" => s"\\_",)))
            else
                yticks := nothing
            end
            grid := :none
            framestyle := :box
            seriescolor := palette[i]
            (Shape([0.0;], [0.0;]))
        end
    end
    if dendrogram
        @series begin
            inset_subplots := (1, bbox(0.895, 0.0, 0.1, 1.0))
            subplot := 2
            orientation := :horizontal
            yticks := nothing
            showaxis := false
            colorbar := nothing
            axis := nothing
            useheight := false
            yflip := true
            x := clustercovariance(Σ²)
        end
    end
end


text:
##
reload("MetabolicFluxAnalysis")
using MetabolicFluxAnalysis

# load file
#fluxi_1_f = open("datasets/antoniewicz_toy_network/toy_network_01.fluxi")
fluxi_1_f = open("testdata/toy_network_01.fluxi")
fluxi_1_s = readstring(fluxi_1_f)
fluxi1 = MetabolicFluxAnalysis.parseFluxi( fluxi_1_s )

#
str_json_01 = MetabolicFluxAnalysis.export_json(fluxi1)

print(str_json_01)

# try import..
fluxi_imported = MetabolicFluxAnalysis.import_json(str_json_01)

fluxi_1_reconstructed = MetabolicFluxAnalysis.import_json_to_fluxi(str_json_01)


text:
# ----------------------  Worker

function make_test_worker(s :: Integer)
    uFct = make_test_worker_utility();
    retireAge = TimeInt(45);
    retireDuration = TimeInt(15);
    wage = 3.5 + Int(s);
    retireIncome = 0.5 * wage;
    R = 1.04;
    xp = collect(range(0.0, 0.6, length = retireAge));
    kMin = -20.0;
    kMax = 75.0;
    hMin = 1.0;
    hMax = 9.5;
    return Worker(uFct, xp, wage, R, retireAge, retireDuration, retireIncome,
        kMin, kMax, hMin, hMax)
end

function validate_worker(wk :: Worker)
    isValid = (wk.wage > 0.0)  &&  (wk.R > 0.9)  &&  (wk.retireDuration > 5);
    if !(wk.xp[1] ≈ 0.0)
        @warn "Experience at t = 1 should be 0"
        isValid = false;
    end
    return isValid
end

## ---------------  Show

function settings_table(wk :: Worker)
    maxExper, maxAge = findmax(wk.xp);
    maxExper = round(maxExper, digits = 2);
    return [
        "Worker"  " ";
        "Retirement age"  "$(wk.retireAge)";
        "Retirement duration"  "$(wk.retireDuration)";
        "Max log experience"  "$maxExper at age $maxAge"
    ]
end

StructLH.describe(wk :: Worker) = settings_table(wk);

function Base.show(io :: IO,  wk :: Worker)
    wageStr = string(round(wk.wage, digits = 2));
    maxExper, maxAge = findmax(wk.xp);
    maxExper = round(maxExper, digits = 2);
    print(io, "Worker:")
    println(io, "  Wage: $wageStr    retireAge: $(wk.retireAge)");
    println(io, "  Max log experience $maxExper at exper $maxAge")
    return nothing
end

## -----------  Access properties

k_range(wk :: Worker) = (wk.kMin, wk.kMax);
h_range(wk :: Worker) = (wk.hMin, wk.hMax);
k_min(wk :: Worker) = wk.kMin;
k_max(wk :: Worker) = wk.kMax;
h_min(wk :: Worker) = wk.hMin;
h_max(wk :: Worker) = wk.hMax;

# Life-span of a worker with work start 1
max_life_span(wk :: Worker) = wk.retireAge + wk.retireDuration - 1;
retire_periods(wk :: Worker, workStartAge) = wk.retireDuration;

# These ages are relative to `workStartAge`
life_span(wk :: Worker, workStartAge) =
    work_periods(wk, workStartAge) .+ retire_periods(wk, workStartAge);
cons_periods(wk :: Worker, workStartAge) = life_span(wk, workStartAge);

last_work_age(wk :: Worker, workStartAge) =
    wk.retireAge .- workStartAge;
work_periods(wk :: Worker, workStartAge) = last_work_age(wk, workStartAge);
work_ages(wk :: Worker, workStartAge) = 1 : last_work_age(wk, workStartAge);
# Retirement ages (relative to workStartAge)
retire_ages(wk :: Worker, workStartAge) =
    last_work_age(wk, workStartAge) .+ (1 : retire_periods(wk, workStartAge));

# Input: ages relative to workStartAge
# is_retired(wk :: Worker, workStartAge :: Integer, tV) =

# function retire_income(wk :: Worker{F1}, workStartAge :: Integer, tV) where F1
#     trV = retire_ages(wk, workStartAge);
#     incomeV = zeros(F1, size(tV)...);
#     incomeV[is_retired(wk, workStartAge, tV)] .= wk.retireIncome;
#     return incomeV
# end

# Present value of retirement income, discounted to workStartAge
function pv_retire_income(wk :: Worker, workStartAge)
    tV = retire_ages(wk, workStartAge);
    pv = retire_income(wk) * pv_factor(wk.R, length(tV));
    pv *= ((1 / wk.R) ^ (tV[1] - 1));
    return pv
end

retire_income(wk :: Worker) = wk.retireIncome;

"""
	$(SIGNATURES)

Log experience productivity profile for given workStartAge until retirement.
Normalized to 1 at work start.

# Arguments
- `experV`: experience levels; work start age is 1.
"""
function log_exper_profile(wk :: Worker, workStartAge :: Integer;
    experV = nothing)
    T = work_periods(wk, workStartAge);
    @assert T > 0
    if isnothing(experV)
        return wk.xp[1 : T]
    else
        return wk.xp[experV]
    end
end

"""
	$(SIGNATURES)

Experience-wage profile. Normalized to 1 at work start.
"""
function exper_profile(wk :: Worker{F1}, workStartAge :: Integer; experV = nothing) where F1
    return exp.(log_exper_profile(wk, workStartAge; experV = experV))
end

# Earnings profile by experience
# Does not include retirement income
function earn_profile(wk :: Worker{F1}, workStartAge :: Integer, h :: F1;
    experV = nothing) where F1

    earnV = h .* wk.wage .* exper_profile(wk, workStartAge; experV = experV);
    return earnV
end

"""
    $(SIGNATURES)

Lifetime earnings for a vector of h. Discounted to `workStartAge`.
Does not include retirement income.
"""
function lifetime_earnings(w :: Worker{F1}, workStartAge :: Integer,
    h :: AbstractArray{F1}) where F1

    # T = work_periods(w, workStartAge);
    # This way it works for array `h` inputs.
    lty = h .* present_value(earn_profile(w, workStartAge, 1.0), w.R);
    @assert all(lty .> 0.0)  "Negative lifetime earnings"
    return lty
end

function lifetime_earnings(w :: Worker{F1}, workStartAge :: Integer, h :: F1) where F1
    lty = only(lifetime_earnings(w, workStartAge, [h]));
    @assert lty > 0.0
    return lty
end

# ---------------