text:
export initsoln, easy_sdp

struct LoadingProblem
    ptr::Ptr{Cvoid}
end

function allocate_loading_prob(pC::Ref{blockmatrix}, block_dims::Vector{CSDP_INT}, num_constraints::Integer, num_entries::Matrix{CSDP_INT}, printlevel::Integer)
    ptr = allocate_loading_prob(pC, fptr(block_dims), num_constraints, pointer(num_entries), printlevel)
    return LoadingProblem(ptr)
end
free_loading_prob(problem::LoadingProblem) = free_loading_prob(problem.ptr)
function free_loaded_prob(problem::LoadingProblem, X::blockmatrix, y::Vector{Cdouble}, Z::blockmatrix)
    free_loaded_prob(problem.ptr, X, fptr(y), Z)
end

function setconstant(problem::LoadingProblem, mat::Integer, ent::Cdouble)
    setconstant(problem.ptr, mat, ent)
end

function addentry(problem::LoadingProblem, mat::Integer, blk::Integer, indexi::Integer, indexj::Integer, ent::Cdouble, allow_duplicates::Bool)
    ret = addentry(problem.ptr, mat, blk, indexi, indexj, ent, allow_duplicates)
    return !iszero(ret)
end

function loaded_initsoln(problem::LoadingProblem, num_constraints::Integer, X::Ref{blockmatrix}, Z::Ref{blockmatrix})
    y = loaded_initsoln(problem.ptr, X, Z)
    return mywrap(y, num_constraints)
end
function initsoln(C::BlockMatrix, b::Vector{Cdouble}, As::Vector{constraintmatrix})
    m = length(As)
    X, y, Z = initsoln(CSDP_INT(size(C, 1)), CSDP_INT(m), C.csdp, fptr(b), fptr(As))
    mywrap(X), mywrap(y, m), mywrap(Z)
end
function initsoln(C::BlockMatrix, b::Vector{Cdouble}, As::Vector{ConstraintMatrix})
    initsoln(C, b, [A.csdp for A in As])
end

function setupAs!(As::Vector{ConstraintMatrix}, C::BlockMatrix)
    # See lib/easysdp.c
    # Fills info in each ConstraintMatrix that couldn't be determined before
    # since it needs info about block category (which is in C) and other constraints matrices

    nblocks = length(C.blocks)

    byblocks = fill(Ptr{sparseblock}(C_NULL), nblocks)

    for constr in length(As):-1:1
        A = As[constr]
        @assert nblocks == length(A.jblocks)
        for blk in 1:nblocks
            jblock = A.jblocks[blk]
            if jblock.csdp.numentries > 0.25 * jblock.csdp.blocksize && jblock.csdp.numentries > 15
                jblock.csdp.issparse = 0
            else
                jblock.csdp.issparse = 1
            end

            if C.blocks[jblock.csdp.blocknum].blockcategory == DIAG
                jblock.csdp.issparse = 1
            end

            jblock.csdp.nextbyblock = byblocks[blk]
            byblocks[blk] = pointer_from_objref(jblock.csdp)
        end
    end
    byblocks
end

sdp(C::BlockMatrix, b::Vector{Cdouble}, As::Vector{ConstraintMatrix}, X::BlockMatrix, y::Vector{Cdouble}, Z::BlockMatrix, params::Dict{Symbol}) = sdp(C, b, As, 0.0, X, y, Z, options(params)...)
sdp(C::BlockMatrix, b::Vector{Cdouble}, As::Vector{ConstraintMatrix}, X::BlockMatrix, y::Vector{Cdouble}, Z::BlockMatrix, printlevel, params) = sdp(C, b, As, 0.0, X, y, Z, printlevel, params)
sdp(C::BlockMatrix, b::Vector{Cdouble}, As::Vector{ConstraintMatrix}, co::Cdouble, X::BlockMatrix, y::Vector{Cdouble}, Z::BlockMatrix, params::Dict{Symbol}) = sdp(C, b, As, co, X, y, Z, options(params)...)
function sdp(C::BlockMatrix, b::Vector{Cdouble}, As::Vector{ConstraintMatrix}, co::Cdouble, X::BlockMatrix, y::Vector{Cdouble}, Z::BlockMatrix, printlevel, params)
    # I pass pointers pX, py and pZ to X, y and Z but only *pX, *py and *pZ are use in the code
    # so no need to worry, they won't change :)
    Xcsdp = X.csdp
    ycsdp = fptr(y)
    Zcsdp = Z.csdp
    Ascsdp = map(A->A.csdp, As)

    n = CSDP_INT(size(C, 1))
    k = CSDP_INT(length(As))
    nd = sizeof(Cdouble) * (n+1)
    kd = sizeof(Cdouble) * (k+1)
    md = max(nd, kd)
    ldam = (iseven(k) ? k+1 : k)
    ld = sizeof(Cdouble) * ldam * ldam

    byblocks = setupAs!(As, C)

    status, pobj, dobj = sdp(n,                                    # n::CSDP_INT
                             k,                                    # k::CSDP_INT
                             C.csdp,                               # C::blockmatrix
                             fptr(b),                              # a::Ptr{Cdouble}
                             co,                                   # constant_offset::Cdouble
                             fptr(Ascsdp),                         # constraints::Ptr{constraintmatrix}
                             fptr(byblocks),
                             Xcsdp,                                # X::blockmatrix
                             ycsdp,                                # y::Ptr{Cdouble}
                             Zcsdp,                                # Z::blockmatrix
                             pointer(Vector{Cdouble}(undef, md)),  # workvec1::Ptr{Cdouble}
                             pointer(Vector{Cdouble}(undef, md)),  # workvec2::Ptr{Cdouble}
                             pointer(Vector{Cdouble}(undef, md)),  # workvec3::Ptr{Cdouble}
                             pointer(Vector{Cdouble}(undef, md)),  # workvec4::Ptr{Cdouble}
                             pointer(Vector{Cdouble}(undef, md)),  # workvec5::Ptr{Cdouble}
                             pointer(Vector{Cdouble}(undef, md)),  # workvec6::Ptr{Cdouble}
                             pointer(Vector{Cdouble}(undef, md)),  # workvec7::Ptr{Cdouble}
                             pointer(Vector{Cdouble}(undef, md)),  # workvec8::Ptr{Cdouble}
                             pointer(Vector{Cdouble}(undef, md)),  # diagO::Ptr{Cdouble}
                             pointer(Vector{Cdouble}(undef, kd)),  # besty::Ptr{Cdouble}
                             pointer(Vector{Cdouble}(undef, ld)),  # O::Ptr{Cdouble}
                             pointer(Vector{Cdouble}(undef, kd)),  # rhs::Ptr{Cdouble}
                             pointer(Vector{Cdouble}(undef, kd)),  # dy::Ptr{Cdouble}
                             pointer(Vector{Cdouble}(undef, kd)),  # dy1::Ptr{Cdouble}
                             pointer(Vector{Cdouble}(undef, kd)),  # Fp::Ptr{Cdouble}
                             CSDP_INT(printlevel),                     # printlevel::CSDP_INT
                             params)                               # parameters::paramstruc
    # I can even assert that they won't change
    @assert Xcsdp == X.csdp
    @assert ycsdp == fptr(y)
    @assert Zcsdp == Z.csdp
    return status, pobj, dobj
end

loaded_sdp(problem::LoadingProblem, X::Ref{blockmatrix}, y::Vector{Cdouble}, Z::Ref{blockmatrix}, params::Dict{Symbol}) = loaded_sdp(problem, 0.0, X, y, Z, options(params)...)
loaded_sdp(problem::LoadingProblem, X::Ref{blockmatrix}, y::Vector{Cdouble}, Z::Ref{blockmatrix}, printlevel, params) = loaded_sdp(problem, 0.0, X, y, Z, printlevel, params)
loaded_sdp(problem::LoadingProblem, co::Cdouble, X::Ref{blockmatrix}, y::Vector{Cdouble}, Z::Ref{blockmatrix}, params::Dict{Symbol}) = loaded_sdp(problem, co, X, y, Z, options(params)...)
function loaded_sdp(problem::LoadingProblem, co::Cdouble, X::Ref{blockmatrix}, y::Vector{Cdouble}, Z::Ref{blockmatrix}, printlevel, params)
    # I pass pointers py to X, y and Z but only *pX, *py and *pZ are
    # used in the code so no need to worry, they won't change :)
    ycsdp = Ref{Ptr{Cdouble}}(fptr(y))

    status, pobj, dobj = loaded_sdp(
        problem.ptr,      # problem::Ptr{Cvoid}
        co,               # constant_offset::Cdouble
        X,                # pX::Ptr{blockmatrix}
        ycsdp,            # py::Ptr{Cdouble}
        Z,                # pZ::Ptr{blockmatrix}
        CSDP_INT(printlevel), # printlevel::CSDP_INT
        params)           # parameters::paramstruc
    # I can even assert that they won't change
    @assert ycsdp[] == fptr(y)
    return status, pobj, dobj
end

parametrized_sdp(C::BlockMatrix, b::Vector{Cdouble}, As::Vector{ConstraintMatrix}, X::BlockMatrix, y::Vector{Cdouble}, Z::BlockMatrix, params::Dict{Symbol}) = parametrized_sdp(C, b, As, 0.0, X, y, Z, options(params)...)
parametrized_sdp(C::BlockMatrix, b::Vector{Cdouble}, As::Vector{ConstraintMatrix}, X::BlockMatrix, y::Vector{Cdouble}, Z::BlockMatrix, printlevel, params) = parametrized_sdp(C, b, As, 0.0, X, y, Z, printlevel, params)
parametrized_sdp(C::BlockMatrix, b::Vector{Cdouble}, As::Vector{ConstraintMatrix}, co::Cdouble, X::BlockMatrix, y::Vector{Cdouble}, Z::BlockMatrix, params::Dict{Symbol}) = parametrized_sdp(C, b, As, co, X, y, Z, options(params)...)
function parametrized_sdp(C::BlockMatrix, b::Vector{Cdouble}, As::Vector{ConstraintMatrix}, co::Cdouble, X::BlockMatrix, y::Vector{Cdouble}, Z::BlockMatrix, printlevel, params)
    # I pass pointers pX, py and pZ to X, y and Z but only *pX, *py and *pZ are
    # used in the code so no need to worry, they won't change :)
    Xcsdp = X.csdp
    ycsdp = Ref{Ptr{Cdouble}}(fptr(y))
    Zcsdp = Z.csdp
    Ascsdp = map(A->A.csdp, As)

    status, pobj, dobj = parametrized_sdp(
        CSDP_INT(size(C, 1)), # n::CSDP_INT
        CSDP_INT(length(As)), # k::CSDP_INT
        C.csdp,           # C::blockmatrix
        fptr(b),          # a::Ptr{Cdouble}
        fptr(Ascsdp),     # constraints::Ptr{constraintmatrix}
        co,               # constant_offset::Cdouble
        ptr(Xcsdp),       # pX::Ptr{blockmatrix}
        ycsdp,            # py::Ptr{Cdouble}
        ptr(Zcsdp),       # pZ::Ptr{blockmatrix}
        CSDP_INT(printlevel), # printlevel::CSDP_INT
        params)           # parameters::paramstruc
    # I can even assert that they won't change
    @assert Xcsdp == X.csdp
    @assert ycsdp[] == fptr(y)
    @assert Zcsdp == Z.csdp
    return status, pobj, dobj
end

easy_sdp(C::BlockMatrix, b::Vector{Cdouble}, As::Vector{constraintmatrix}, X::BlockMatrix, y::Vector{Cdouble}, Z::BlockMatrix) = easy_sdp(C, b, As, 0.0, X, y, Z)
function easy_sdp(C::BlockMatrix, b::Vector{Cdouble}, As::Vector{constraintmatrix}, co::Cdouble, X::BlockMatrix, y::Vector{Cdouble}, Z::BlockMatrix)
    # I pass pointers pX, py and pZ to X, y and Z but only *pX, *py and *pZ are
    # used in the code so no need to worry, they won't change :)
    Xcsdp = X.csdp
    ycsdp = Ref{Ptr{Cdouble}}(fptr(y))
    Zcsdp = Z.csdp

    status, pobj, dobj = easy_sdp(CSDP_INT(size(C, 1)), # n::CSDP_INT
                                  CSDP_INT(length(As)), # k::CSDP_INT
                                  C.csdp,           # C::blockmatrix
                                  fptr(b),          # a::Ptr{Cdouble}
                                  fptr(As),         # constraints::Ptr{constraintmatrix}
                                  co,               # constant_offset::Cdouble
                                  ptr(Xcsdp),       # pX::Ptr{blockmatrix}
                                  ycsdp,            # py::Ptr{Cdouble}
                                  ptr(Zcsdp))       # pZ::Ptr{blockmatrix}
    # I can even assert that they won't change
    @assert Xcsdp == X.csdp
    @assert ycsdp[] == fptr(y)
    @assert Zcsdp == Z.csdp
    return status, pobj, dobj
end

function write_prob(fname::String, C::BlockMatrix, b::Vector{Cdouble}, As::Vector{constraintmatrix})
    write_prob(fname, CSDP_INT(size(C, 1)), CSDP_INT(length(As)), C.csdp, fptr(b), fptr(As))
end

function write_sol(fname::String, X::BlockMatrix, y::Vector{Cdouble}, Z::BlockMatrix)
    write_sol(fname, CSDP_INT(size(X, 1)), CSDP_INT(length(y)), X.csdp, fptr(y), Z.csdp)
end


text:
function absdiff(x::T,y::T) where {T<:Unsigned}
    d = max(x,y) - min(x,y)
    d, d
end
function absdiff(x::T,y::T) where {T<:Signed}
    d = x - y
    abs(d), d
end
# binary GCD (aka Stein's) algorithm
# about 1.7x (2.1x) faster for random Int64s (Int128s)
# Unfortunately, we need to manually annotate this as `@assume_effects :terminates_locally` to work around #41694.
# Since this is used in the Rational constructor, constant folding is something we do care about here.
function _gcd(ain::T, bin::T) where T<:Integer
    zb = trailing_zeros(bin)
    za = trailing_zeros(ain)
    a = abs(ain)
    b = abs(bin >> zb)
    k = min(za, zb)
    while a != 0
        a >>= za
        absd, diff = absdiff(a, b)
        za = trailing_zeros(diff)
        b = min(a, b)
        a = absd
    end
    r = b << k
    return r % T
end

ext: jl


text:
# Helper functions

using Dates

function makehomesavepath(foldername)
    rightnow = Dates.format(Dates.now(), "yyyy-mm-dd_HH-MM-SS")

    savepath = joinpath(homedir(), foldername * "_" * rightnow * "/")

    while isdir(savepath)
        sleep(1.)
        rightnow = Dates.format(Dates.now(), "yyyy-mm-dd_HH-MM-SS")

        savepath = joinpath(homedir(), foldername * "_" * rightnow * "/")
    end
    savepath
end

function checkifpathexists(result_path_base)
    if ispath(result_path_base)
        if "del" in ARGS
            println("rm($result_path_base)")
            rm(result_path_base; recursive=true)
        elseif "append" in ARGS
        # nothing to do, append new content to the result folder
        else
            error("Result folder already exists: $result_path_base")
        end
    end
end

nanmean(x) = mean(filter(!isnan,x))
nanmean(x,dims) = mapslices(nanmean, x, dims=dims)
nanstd(x) = std(filter(!isnan,x))
nanstd(x,dims) = mapslices(nanstd, x, dims=dims)
nanlength(x) = length(filter(!isnan,x))
nanlength(x,dims) = mapslices(nanlength, x, dims=dims)
function nanserror(x, dims)
    std = nanstd(x, dims)
    n = nanlength(x, dims)
    std ./ sqrt.(n)
end
function calc_stats_nan(dataarray) # dataarray : timepoints x instances

    m = vec(nanmean(dataarray, 2))
    s = vec(nanstd(dataarray, 2))
    serror = vec(nanserror(dataarray, 2))
    m, s, serror
end

function checkInf(variable)
    if isinf(variable) #|| ratioterm >= 1e16
        #println("I'm Inf")
        variable = 1e16
    end
    variable
end

searchdirendings(path, key) = filter(x->endswith(x, key), readdir(path))

getsubfolders(path) = [x for x in readdir(path) if isdir(joinpath(path,x))]

# Dataframe manipulations
function dftranspose(;
                    filepath =  "../R2adj/",
                    filename = "medianR2adj.csv",
                    titleRowsSymbols = [:REINFORCE, :AC, :SurpAC, :HybridAC, :SurpREINFORCE],
                    savepath = "./",
                    )

    df = CSV.read(joinpath(filepath, filename), header=false)

    dftrans = dfjusttranspose(df)

    dftrans = dfplacetitlesinfirstcolumn(dftrans, titleRowsSymbols=titleRowsSymbols)

    prefix = split(filename, ".")[end-1]
    CSV.write(joinpath(savepath, prefix * "_trans.csv"), dftrans, delim = " ")
    dftrans
end

function dfjusttranspose(df)
    dftrans = DataFrame([[names(df)]; collect.(eachrow(df))], [:column; Symbol.(axes(df, 1))])
end

function dfstatspercolumn(;
                filepath =  "../R2adj_fitruns/",
                filename = "medianR2adj.csv",
                titleRowsSymbols = [:run1, :run2, :run3, :run4, :run5, :run6],
                savepath = "./",
                )

    df = CSV.read(joinpath(filepath, filename), header=false)

    dfstats = dfjustdostats(df, titleRowsSymbols)

    prefix = split(filename, ".")[end-1]
    CSV.write(joinpath(savepath, prefix * "_statistics.csv"), dfstats, delim = " ")
    dfstats
end

function dfjustdostats(df, titleRowsSymbols)

    dfstats = DataFrame(
                column = titleRowsSymbols,
                _mean =  mean.(eachcol(df)),
                _std = std.(eachcol(df)),
                _serror = std.(eachcol(df))./sqrt.(length.(eachcol(df)))
                )
    dfstats
end

function dfplacetitlesinfirstcolumn(df;
                    titleRowsSymbols = [:REINFORCE, :AC, :SurpAC, :HybridAC]
                    )

    if ~isempty(titleRowsSymbols)
        if length(titleRowsSymbols) == size(df,1)
            df[:,1] = titleRowsSymbols
        else
            println("Length of row titles does not match!")
        end
    else
        println("No row titles")
    end
    df
end

function dfsort(df)

    dfstats = DataFrame(meanmedR2adj =  mean.(eachcol(df)))

    p = sortperm(dfstats[:,1])

     dsorted = df[:,p]
 end

function dfcenterwrtlastcolumn(df)

    dfcentered = df .- df[:,end]

end

# deleterows!(dftrans,1)
# rename!(dftrans, titleRowsSymbols)


text:
include("SA.jl")
# using ProfileView
# Function to be run from the command line
# Profile.clear()
# @profile SA.optimize(ARGS)
SA.optimize("") #garbage run to get things into memory
SA.optimize(ARGS)
# ProfileView.view()


text:
using Documenter, ArgParse2

makedocs(;
    modules=[ArgParse2],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/kmsquire/ArgParse2.jl/blob/{commit}{path}#L{line}",
    sitename="ArgParse2.jl",
    authors="Kevin Squire",
    assets=String[],
)

deploydocs(;
    repo="github.com/kmsquire/ArgParse2.jl",
)


text:
using Test
using GPX

@testset "Sample" begin
    @testset "write_gpx" begin
		include("../sample/write_gpx.jl")
		@test 1 == 1
    end

    @testset "read_gpx" begin
		include("../sample/read_gpx.jl")
		@test 1 == 1
    end

end

@testset "parse gpx" begin
    s = """<gpx xmlns="http://www.topografix.com/GPX/1/1" version="1.1" creator="GPX.jl">
	<metadata>
		<name>07/11/2019 LFBI (09:32) LFBI (11:34)</name>
		<author>
			<name>GPX.jl</name>
		</author>
		<time>2019-11-09T09:14:17.473Z</time>
	</metadata>
	<trk>
		<trkseg>
			<trkpt lat="46.58634900" lon="0.30904200">
				<ele>63.6000</ele>
				<time>2019-11-07T08:44:06Z</time>
				<desc>Vitesse : 0 kts
Altitude : 209 ft
Latitude : 46.586349
Longitude : 0.309042
Cap : 314°
</desc>
			</trkpt>
			<trkpt lat="46.58653100" lon="0.30869000">
				<ele>64.8000</ele>
				<time>2019-11-07T08:44:55Z</time>
				<desc>Vitesse : 11 kts
Altitude : 213 ft
Latitude : 46.586531
Longitude : 0.308690
Cap : 307°
</desc>
            </trkpt>
        </trkseg>
    </trk>
</gpx>"""
    gpx = parse_gpx_string(s)
	println(gpx)

	# Metadata
	# ToFix @test gpx.metadata.name == "07/11/2019 LFBI (09:32) LFBI (11:34)"
	# ToFix @test gpx.metadata.author.name == "GPX.jl"
	# ToFix @test gpx.metadata.time == ZonedDateTime("2019-11-09T09:14:17.473+00:00")

	# Tracks / Track segments / Track points
	@test length(gpx.tracks) == 1
	@test length(gpx.tracks[1].segments) == 1
	@test length(gpx.tracks[1].segments[1].points) == 2

	pt = gpx.tracks[1].segments[1].points[1]
	@test pt.lat == 46.58634900
	@test pt.lon == 0.30904200
	@test pt.ele == 63.6000
	@test pt.time == ZonedDateTime("2019-11-07T08:44:06.000+00:00")

	pt = gpx.tracks[1].segments[1].points[2]
	@test pt.lat == 46.58653100
	@test pt.lon == 0.30869000
	@test pt.ele == 64.8000
	@test pt.time == ZonedDateTime("2019-11-07T08:44:55.000+00:00")
end


text:
# This file was generated, do not modify it. # hide
r = report(m)
r.best_result


text:
# chemically equivalent nuclei are ones that have the same chemical shift.

### routines related to checking whether the J-coupling value between
#   chemically equivalent nuclei are the same.

#### misc.

"""
mapQtoHID(Q_g::Vector{Vector{Int}}, H_IDs::Vector{Int})

Convert entries in `Q_g` from taking on integer values between 1 to length(`H_IDs`) to taking on values in `H_IDs`.
"""
function mapQtoHID(Q_g::Vector{Vector{Int}}, H_IDs::Vector{Int})

    dict_g_to_H_IDs = Dict( collect(1:length(H_IDs)) .=> H_IDs)

    Q = Vector{Vector{Int}}(undef, length(Q_g))
    for i = 1:length(Q_g)

        Q[i] = Vector{Int}(undef, length(Q_g[i]))
        for k = 1:length(Q_g[i])
            Q[i][k] = dict_g_to_H_IDs[Q_g[i][k]]
        end
    end

    return Q
end


text:
# Speech waveform synthesis filters

poledf_delay_length(order) = order
poledf_delay(order) = zeros(poledf_delay_length(order))

function poledf(x::Cdouble, a::StridedVector{Cdouble},
                delay::StridedVector{Cdouble})
    order = length(a) - 1
    if length(delay) != poledf_delay_length(order)
        throw(DimensionMismatch("inconsistent delay length"))
    end
    ccall((:poledf, libSPTK), Cdouble,
          (Cdouble, Ptr{Cdouble}, Cint, Ptr{Cdouble}),
          x, a, length(a)-1, delay)
end

lmadf_delay_length(order, pd) = 2pd*(order + 1)
lmadf_delay(order, pd) = zeros(lmadf_delay_length(order, pd))

function lmadf(x::Cdouble, b::StridedVector{Cdouble}, pd,
               delay::StridedVector{Cdouble})
    assert_pade(pd)
    order = length(b) - 1
    if length(delay) != lmadf_delay_length(order, pd)
        throw(DimensionMismatch("inconsistent delay length"))
    end
    ccall((:lmadf, libSPTK), Cdouble,
          (Cdouble, Ptr{Cdouble}, Cint, Cint, Ptr{Cdouble}),
          x, b, length(b)-1, pd, delay)
end

lspdf_delay_length(order) = 2order + 1
lspdf_delay(order) = zeros(lspdf_delay_length(order))

function lspdf(x::Cdouble, f::StridedVector{Cdouble},
               delay::StridedVector{Cdouble})
    order = length(f) - 1
    if iseven(order)
        lspdf_even(x, f, delay)
    else
        lspdf_odd(x, f, delay)
    end
end

function lspdf_even(x::Cdouble, f::StridedVector{Cdouble},
                    delay::StridedVector{Cdouble})
    order = length(f) - 1
    if length(delay) != lspdf_delay_length(order)
        throw(DimensionMismatch("inconsistent delay length"))
    end
    ccall((:lspdf_even, libSPTK), Cdouble,
          (Cdouble, Ptr{Cdouble}, Cint, Ptr{Cdouble}),
          x, f, length(f) - 1, delay)
end

function lspdf_odd(x::Cdouble, f::StridedVector{Cdouble},
                   delay::StridedVector{Cdouble})
    order = length(f) - 1
    if length(delay) != lspdf_delay_length(order)
        throw(DimensionMismatch("inconsistent delay length"))
    end
    ccall((:lspdf_odd, libSPTK), Cdouble,
          (Cdouble, Ptr{Cdouble}, Cint, Ptr{Cdouble}),
          x, f, length(f) - 1, delay)
end

ltcdf_delay_length(order) = order + 1
ltcdf_delay(order) = zeros(ltcdf_delay_length(order))

function ltcdf(x::Cdouble, k::StridedVector{Cdouble},
               delay::StridedVector{Cdouble})
    order = length(k) - 1
    if length(delay) != ltcdf_delay_length(order)
        throw(DimensionMismatch("invalid delay length"))
    end
    ccall((:ltcdf, libSPTK), Cdouble,
          (Cdouble, Ptr{Cdouble}, Cint, Ptr{Cdouble}),
          x, k, length(k) - 1, delay)
end

glsadf_delay_length(order, stage) = order * (stage + 1) + 1
glsadf_delay(order, stage) = zeros(glsadf_delay_length(order, stage))

# NOTE: stage = -1/γ
function glsadf(x::Cdouble, c::StridedVector{Cdouble}, stage,
                delay::StridedVector{Cdouble})
    stage >= 1 || throw(ArgumentError("stage >= 1 (-1 <= γ < 0)"))
    order = length(c) - 1
    if length(delay) != glsadf_delay_length(order, stage)
        throw(DimensionMismatch("invalid delay length"))
    end
    ccall((:glsadf, libSPTK), Cdouble,
          (Cdouble, Ptr{Cdouble}, Cint, Cint, Ptr{Cdouble}),
          x, c, length(c) - 1, stage, delay)
end

# see mlsadf.c in original SPTK for this magic allocation
mlsadf_delay_length(order, pd) = 3*(pd+1) + pd*(order+2)
mlsadf_delay(order, pd) = zeros(mlsadf_delay_length(order, pd))

function mlsadf(x::Cdouble, b::StridedVector{Cdouble}, α, pd,
                delay::StridedVector{Cdouble})
    assert_pade(pd)
    order = length(b) - 1
    if length(delay) != mlsadf_delay_length(order, pd)
        throw(DimensionMismatch("invalid delay length"))
    end
    ccall((:mlsadf, libSPTK), Cdouble,
          (Cdouble, Ptr{Cdouble}, Cint, Cdouble, Cint, Ptr{Cdouble}),
          x, b, length(b)-1, α, pd, delay)
end

# see mglsadf.c in original SPTK for this magic allocation
mglsadf_delay_length(order, stage) = (order+1)*stage
mglsadf_delay(order, stage) = zeros(mglsadf_delay_length(order, stage))

function mglsadf(x::Cdouble, b::StridedVector{Cdouble}, α, stage,
                 delay::StridedVector{Cdouble})
    stage >= 1 || throw(ArgumentError("stage >= 1 (-1 <= γ < 0)"))
    order = length(b) - 1
    if length(delay) != mglsadf_delay_length(order, stage)
        throw(DimensionMismatch("invalid delay length"))
    end
    ccall((:mglsadf, libSPTK), Cdouble,
          (Cdouble, Ptr{Cdouble}, Cint, Cdouble, Cint, Ptr{Cdouble}),
          x, b, length(b)-1, α, stage, delay)
end


text:
dropmissing(x) = collect(skipmissing(x))

"""
    smallest(da::AbstractArray, n::Int) or smallest(df::DataFrame, varname::Symbol, n::Int)

List the `n` smallest values in `da` in a descending order. The default `n` is 5.
"""
function smallest(da::AbstractArray; n::Int = 5)
	return sort(collect(skipmissing(da)))[1:n]
end
smallest(df::DataFrame,varname::Symbol; n = 5) = smallest(df[!,varname], n = n)

"""
    largest(da::AbstractArray, n::Int) or largest(df::DataFrame, varname::Symbol, n::Int)

List the `n` largest values in `da` in an ascending order. The default `n` is 5.
"""
function largest(da::AbstractArray; n::Int = 5)
	return sort(collect(skipmissing(da)))[end-n+1:end]
end
largest(df::DataFrame,varname::Symbol; n = 5) = largest(df[!,varname], n = n)

"""
    univariate(da::AbstractArray) or univariate(df::DataFrame,varname::Symbol)

Produce univariate statistics from the `da` and return a DataFrame with two columns,
:Statistic and :Value. Computed statistics include `N_Total` (number of rows
in the input DataArray), `N_Missing` (number of rows with missing values), `N_Used` (number of non-missing rows),
`Sum` (sum), `Mean` (mean), `SD` (standard deviation), `Var` (variance), `Min` (minimum),
`P25` (25th percentile), `Median` (median), `P75` (75th percentile), `Max` (maximum),
`Skewness` (skewness), and `Kurtosis` (kurtosis).
"""
function univariate(da::AbstractVector)

    da2 = dropmissing(da)

	return DataFrame(
		Statistic = [:N_Total,
			:N_Missing,
			:N_Used,
			:Sum,
			:Mean,
			:SD,
			:Var,
			:Min,
			:P25,
			:Median,
			:P75,
			:Max,
			:Skewness,
			:Kurtosis],
		Value = [size(da,1),
			size(da,1) - size(da2,1),
			size(da2,1),
			sum(da2),
			mean(da2),
			sqrt(var(da2)),
			var(da2),
			minimum(da2),
			quantile(da2,0.25),
			quantile(da2,0.5),
			quantile(da2,0.75),
			maximum(da2),
			skewness(da2),
			kurtosis(da2)]
	)
end
univariate(df::DataFrame,var::Symbol) = univariate(df[!,var])

function strval(val::AbstractFloat)
  return @sprintf("%.2f",val)
end
function strval(val::AbstractFloat, decimals::Int)
    if decimals == 1
        return @sprintf("%.1f",val)
    elseif decimals == 2
        return @sprintf("%.2f",val)
    elseif decimals == 3
        return @sprintf("%.3f",val)
    elseif decimals == 4
        return @sprintf("%.4f",val)
    elseif decimals == 5
        return @sprintf("%.5f",val)
    elseif decimals == 6
        return @sprintf("%.6f",val)
    else
        return @sprintf("%f",val)
    end
end
function strval(val::Integer)
  return @sprintf("%.0d",val)
end

function getdictval(dt::Dict,val)
    return haskey(dt,val) ? dt[val] : val
end

"""
    tabstat(df::DataFrame,varname::Symbol,groupvar::Symbol;s::Vector{Function} = [N,mean,sd,minimum,p25,median,p75,maximum ],wt::Union{Nothing,Symbol}=nothing)

Produce a DataFrame that contains summary statistics for each `groupvar` subgroup
of the `varname` column in the `df`. The following are computed: `n` (total non-missing rows),
`mean` (mean), `sd` (standard deviation), `min` (minimum), `p25` (25th percentile),
`median` (median), `p75` (75th percentile), and `max` (maximum).
"""
function tabstat(indf::DataFrame, var1::Symbol, groupvar::Symbol; s::Vector{Function} = [N,mean,sd,minimum,p25,median,p75,maximum ],wt::Union{Nothing,Symbol}=nothing)

    if length(s) == 0
        error("No statistic functions were specified.")
    end

    # strip prepended "Stella." from the function list
    namevec = [Symbol(replace(string(x),r"(Stella|Statistics)\." => "")) for x in s]

    # number of levels in the groupvar
    lev = DataFrames.levels(indf[!,groupvar])
    gvnum = length(lev)

    # output DataFrame
    outdf = DataFrame()

    # groupvar
    outdf[!,groupvar] = Vector{Union{Missing,eltype(groupvar)}}(undef, gvnum)

    # stats as variables
    for j = 1:length(namevec)

        # for N, create a vector of integers
        if namevec[j] == :N
            outdf[!,namevec[j]] = Vector{Union{Missing,Int}}(undef, gvnum)
        # for all other stats, create float64 vectors
        else
            outdf[!,namevec[j]] = Vector{Union{Missing,Float64}}(undef, gvnum)
        end
    end

    # get stat for each level
    for subdf in groupby(indf, groupvar)

        # level index to be used as the row number in outdf
        gidx = findfirst(x->x==subdf[1,groupvar],lev)

        # non-missing values only
        da = collect(skipmissing(subdf[var1]))

        # if array is empty, there is nothing to process
        if length(da) == 0
            continue
        end

        # groupvar value
        outdf[gidx,groupvar] = subdf[1,groupvar]

        # obtain stats
        for j = 1:length(namevec)
            outdf[gidx,namevec[j]] = s[j](da)
        end
    end

    sort!(outdf,groupvar)

    return outdf
end

#----------------------------------------------------------------------------
# stats by subdataframe - begin
#----------------------------------------------------------------------------

"""
    p5(x::AbstractVector)
    p10(x::AbstractVector)
    p25(x::AbstractVector)
    p50(x::AbstractVector)
    p75(x::AbstractVector)
    p90(x::AbstractVector)
    p95(x::AbstractVector)

Produce percentile values (e.g., p5() produces 5% percentile, p10 10% percentile, etc)
from `x` AbstractVector. If `x` is a DataArray, it must be stripped of NA's before used in these functions.
Other percentile functions not shown here can easily be created by, for example,
`p15(x) = quantile(x,.15)`.
"""
p5(x) = quantile(x,.5)
p10(x) = quantile(x,.1)
p25(x) = quantile(x,.25)
p50(x) = quantile(x,.5)
p75(x) = quantile(x,.75)
p90(x) = quantile(x,.9)
p95(x) = quantile(x,.95)

"""
    cv(x::AbstractVector)

Is an alias of `variation()`.
"""
cv(x) = variation(x)

"""
    se(x::AbstractVector)

Is an alias of `sem(x)`.
"""
se(x) = sem(x)

"""
    sd(x::AbstractVector)

Is an alias of `std(x)`.
"""
sd(x) = std(x)

"""
    N(x::AbstractVector)

Is an alias of `size(x,1)`.
"""
N(x) = size(x,1)
# substat - pass functions to produce the stats by subframe
# functions: mean, var, std, sum, count, kurtosis, skewness, median, percentiles, maximum, minimum
#            cov, se

"""
    substat(df::DataFrame, varname::Symbol, groupvars::Array{Symbol,1}, func::Function)

Produces a vector of the same length as the `df` that contains summary statistics of
the `varname` column computed by `func` for each subgroup stratified by `groupvars`.
Multiple columns can be used as `groupvar`'s if provided in an array of Symbols.
Any function that can take one numeric vector as input can be used. Percentile functions
such as those that take two arguments (e.g., `quantile(x,.25)` for a 25th percentile of x)
should first be converted to a function that takes one argument (e.g, `p25(x) = quantile(x,.25)`).
The following percentile functions are already defined: `p5`, `p10`, `p25`, 'p50', `p75`, and `p90`.
This function emulates the Stata's `egen functions` such as `egen testmean = mean(test), by(groupvar)`.

```
julia> df[:testmean] = substat(df,:test, :class, mean)
```
"""
function substat(df::DataFrame, varname::Symbol, groupvars::Vector{Symbol}, func::Function)

    if (Missings.T(eltype(df[varname])) <: Number) == false
        error("Only numeric variables are allowed.")
    end

    df2 = by(df,groupvars) do subdf
        da = Vector(collect(skipmissing(subdf[varname])))
        if length(da) == 0
            DataFrame(x1 = missing)
        else
            DataFrame(x1 = func(da))
        end
    end

    df4 = DataFrame()
    df4[groupvars] = df[groupvars]
    df4[:___obs___] = collect(1:size(df4,1))

    df5 = join(df4, df2, on = groupvars, kind = :left)
    sort!(df5,cols=[:___obs___])
    return df5[:x1]
end
substat(df::DataFrame, varname::Symbol, groupvar::Symbol, func::Function) = substat(df,varname,[groupvar],func)

#----------------------------------------------------------------------------
# eform
#----------------------------------------------------------------------------
function eform(glmout::StatsModels.RegressionModel)

    # family and link function
    if isa(glmout.model,GeneralizedLinearModel)
        distrib = glmout.model.rr.d
        linkfun = GLM.Link(glmout.model.rr)
    else
        error("GLM model is required as the argument")
    end

	coeftable2 = coeftable(glmout)

	# estimates
	coeftable2.cols[1] = exp.(coeftable2.cols[1])

	# standard errors
	coeftable2.cols[2] = coeftable2.cols[1] .* coeftable2.cols[2]

	# rename column1 to OR
    coeftable2.colnms[1] = coeflab(distrib,linkfun)

	return coeftable2
end
function eform(glmout::StatsModels.RegressionModel, labels::Label)

    # family and link function
    if isa(glmout.model,GeneralizedLinearModel)
        distrib = glmout.model.rr.d
        linkfun = GLM.Link(glmout.model.rr)
    else
        error("GLM model is required as the argument")
    end

	coeftable2 = coeftable(glmout)

	# estimates
	coeftable2.cols[1] = exp.(coeftable2.cols[1])

	# standard errors
	coeftable2.cols[2] = coeftable2.cols[1] .* coeftable2.cols[2]

	# rename column1 to OR
	coeftable2.colnms[1] = coeflab(distrib,linkfun)

	# parse the row names and change variable names and values
	for i in 2:length(coeftable2.rownms)
		# parse row name and split into a tuple (varname, value)
		if occursin(": ",coeftable2.rownms[i])
			(varname,value) = split(coeftable2.rownms[i],": ")
		else
			(varname,value) = (coeftable2.rownms[i],"")
		end

        # get variable label from the label dictionary
		varlabel = varlab(labels,Symbol(varname))
        if varlabel == ""
            varlabel = varname
        end
        if value != ""
		    vlabel = vallab(labels,Symbol(varname),parse(Int,value))
        end

        # If value is 1 and value label is Yes, it is a binary variable
		# do not print
		if value == 1 && ismatch(r"^ *yes *$"i,vlabel)
			coeftable2.rownms[i] = varlabel
		elseif value != ""
			coeftable2.rownms[i] = string(varlabel, ": ", vlabel)
        else
            coeftable2.rownms[i] = varlabel
		end
	end

	return coeftable2
end

function coeflab(d::UnivariateDistribution,l::Link)
    if (isa(d,Bernoulli) || isa(d,Binomial)) && isa(l,LogitLink)
        return "OR"
    elseif (isa(d,Binomial) || isa(d,NegativeBinomial)) && isa(l,LogLink)
        return "RR"
    elseif isa(d,Poisson) && isa(l,LogLink)
        return "IRR"
    else
        return "exp(Est)"
    end
end

import GLM.coeftable
function GLM.coeftable(r::StatsModels.RegressionModel,labels::Label)
    coeftable2 = coeftable(r)

    # parse the row names and change variable names and values
	for i in 2:length(coeftable2.rownms)
		# parse row name and split into a tuple (varname, value)
		if occursin(": ",coeftable2.rownms[i])
			(varname,value) = split(coeftable2.rownms[i],": ")
		else
			(varname,value) = (coeftable2.rownms[i],"")
		end

        # get variable label from the label dictionary
		varlabel = varlab(labels,Symbol(varname))
        if varlabel == ""
            varlabel = varname
        end
        if value != ""
		    vlabel = vallab(labels,Symbol(varname),parse(Int,value))
        end

        # If value is 1 and value label is Yes, it is a binary variable
		# do not print
		if value == 1 && ismatch(r"^ *yes *$"i,vlabel)
			coeftable2.rownms[i] = varlabel
		elseif value != ""
			coeftable2.rownms[i] = string(varlabel, ": ", vlabel)
        else
            coeftable2.rownms[i] = varlabel
		end
	end

	return coeftable2
end

#--------------------------------------------------------------------------
# pairwise correlations
#--------------------------------------------------------------------------
struct pwcorr_return
    r::AbstractArray
    pval::AbstractArray
    N::AbstractArray
    colnms::Vector
end

"""
    pwcorr(a::AbstractArray)
    pwcorr(a::AbstractVector...)
    pwcorr(df::DataFrame, args::Symbol...)

Produces a n x n matrix that contain correlation coefficients between all `n` columns.
When the option `p = true` is specified, it produces a second n x n matrix containing
p-values which is calculated using the Fisher transformation. When the option `out = true`
is set, pwcorr prints the correlation coeffients and p-values (if specified) to the console
in addition to returning them in arrays.
"""
function pwcorr(a::AbstractArray)

    cols = size(a,2)
    N = zeros(Int64,cols,cols)
    r = zeros(Float64,cols,cols)
    pval = zeros(Float64,cols,cols)
    if typeof(a) == DataFrame
        colnames = names(a)
    else
        colnames = collect(1:size(a,2))
    end

    for j = 1:cols
        for i = j:cols
            if i != j
                if typeof(a) == DataFrame
                    x = Array(a[completecases(a[[i,j]]),[i,j]])
                else
                    x = a[:,[i,j]]
                end
                r[i,j] = r[j,i] = cor(x)[1,2]
                rows = N[i,j] = N[j,i] = size(x,1)
                z = (sqrt(rows - 3) / 2)*log((1+r[i,j])/(1-r[i,j]))
                pval[i,j] = pval[j,i] = Distributions.ccdf(Distributions.Normal(),z) / 2.0
            end
        end
    end
    return pwcorr_return(r, pval, N, colnames)
end
pwcorr(a::AbstractArray...) = pwcorr(hcat(a...))
pwcorr(a::DataFrame, args::Vector{Symbol}; out=true) = pwcorr(df[args], out = out)

#--------------------------------------------------------------------------
# ranksum(), signrank(), signtest()
#--------------------------------------------------------------------------
#
# struct ranksum_return
#     df::DataFrame
#     varname::Symbol
#     t::Float64
#     U::Float64
#     s²::Float64
#     var_t::Float64
#     z::Float64
#     pvalue::Float64
#     porder::Float64
# end
#
# function ranksum(df::DataFrame,varname::Symbol; group::Symbol = nothing)
#     if group == nothing
#         error("`group` variable is required")
#     end
#
#     df2 = df[completecases(df[[varname,group]]),[varname,group]]
#     df2[:__ranking] = tiedrank(df[varname])
#
#     df3 = by(df2,group) do subdf
#         DataFrame(
#             obs = size(subdf,1),
#             ranksum = sum(subdf[:__ranking]),
#             expected = size(subdf,1)*(size(df2,1)+1) / 2
#         )
#     end
#
#     # test statistic
#     t = df3[1,:ranksum]
#     n1 = df3[1,:obs]
#     U = t - n1*(n1+1) / 2
#     meanrank = mean(df2[:__ranking])
#     s² = var(df2[:__ranking])
#     vart = df3[1,:obs]*df3[2,:obs]*s² / size(df2,1)
#     z = (t - df3[1,:expected]) / sqrt(vart)
#     pval = ccdf()
#     porder = U / (df3[1,:obs]*df3[2,:obs])
#     pval = 2*Distributions.cdf(Distributions.Normal(), z)
#     return ranksum_return(df3,varname,t,U,s²,vart,z,pval,porder)
# end
# #
# # function print(r::ranksum_return)
# #     group = names(r.df)[1]
# #     print(r.df,"\n\n")
# #     @printf("         z = %.3f\n",r.z)
# #     @printf("Prob > |z| = %.3f\n",r.pvalue)
# #     @printf("P{%s(%s == %s)} >  P{%s(%s == %s)} = %.3f\n",
# #         string(r.varname),
# #         string(group),
# #         string(r.df[1,1]),
# #         string(r.varname),
# #         string(group),
# #         string(r.df[2,1]),
# #         r.porder
# #         )
# # end
#
# struct signrank_return
#     df::DataFrame
#     t::Float64
#     var_t::Float64
#     z::Float64
#     p::Float64
# end
#
# function signrank(df::DataFrame,var1::Symbol,var2::Symbol)
#
#     # complete cases only
#     df2 = df[completecases(df[[var1,var2]]),[var1,var2]]
#
#     # difference
#     d = df[var1] .- df[var2]
#
#     # sign
#     s = sign.(d)
#
#     # rank the absolute difference
#     r = tiedrank(abs(d)).*s
#
#     # test statistic
#     t = sum(r)
#
#     # dataframe output
#     n = size(df2,1)
#     n0 = sum(s .== 0)
#     et = (n*(n + 1)/2 - sum(s .== 0))/2
#     df = DataFrame(
#         sign = ["positive","negative","zero"],
#         obs = [
#             sum(s .== 1),
#             sum(s .== -1),
#             n0
#         ],
#         sum_ranks = [
#             sum((s .== 1) .* r),
#             abs(sum((s .== -1) .* r)),
#             n0
#         ],
#         expected = [et,et,n0]
#     )
#
#     varadjt = sum(r.^2) / 4
#     varunadjt = (n*(n + 1)*(2n + 1))/24
#     Δvarzeroadj = -1 * n0*(n0 + 1)*(2*n0 + 1) / 24
#     Δvartiesadj = varadjt - varunadjt - Δvarzeroadj
#
#     # z-value
#     z = (df[1,:sum_ranks] - et) / sqrt(varadjt)
#
#     # p-value
#     pval = 2*Distributions.cdf(Distributions.Normal(), z)
#
#     return signrank_return(df, et, varadjt, z, pval)
# end
# #
# # function print(sr::signrank_return)
# #     print(sr.df,"\n\n")
# #     @printf("Adj variance = %.3f\n",sr.var_t)
# #     @printf("           z = %.3f\n",sr.z)
# #     @printf("  Prob > |z| = %.3f\n",sr.p)
# # end
#
#
# struct signtest_return
#     df::DataFrame
#     p_left::Float64
#     p_right::Float64
#     p_both::Float64
# end
#
# function signtest(df::DataFrame,var1::Symbol,var2::Symbol)
#
#     # complete cases only
#     df2 = df[completecases(df[[var1,var2]]),[var1,var2]]
#
#     # difference
#     d = df[var1] .- df[var2]
#
#     # sign
#     s = sign.(d)
#
#     # test statistic
#     nplus = sum(s .> 0)
#     n0 = sum(s .== 0)
#     n_nonzero = size(df2,1) - n0
#     et = n_nonzero / 2
#
#     # output dataframe
#     df3 = DataFrame(
#         sign = ["positive","negative","zero"],
#         observed = [
#             sum(s .== 1),
#             sum(s .== -1),
#             n0
#         ],
#         expected = [et,et,n0]
#     )
#
#     # p1 - Ha: median of var1 - var2 > 0
#     p1 = Distributions.cdf(Distributions.Binomial(n_nonzero,.5),df3[2,:observed])
#
#     # p2 - Ha: median of var1 - var2 < 0
#     p2 = Distributions.cdf(Distributions.Binomial(n_nonzero,.5),df3[1,:observed])
#
#     # p2 - Ha: median of var1 - var2 != 0
#     p3 = min(1.0, 2.0*Distributions.cdf(Distributions.Binomial(n_nonzero,.5),df3[1,:observed]))
#
#     return signtest_return(df3, p1, p2, p3)
# end
# #
# # function print(sr::signtest_return)
# #     print(sr.df,"\n\n")
# #     @printf("Pr(#positive ≥ %d) = %.3f\n",sr.df[1,:observed],sr.p_left)
# #     @printf("Pr(#negative ≥ %d) = %.3f\n",sr.df[2,:observed],sr.p_right)
# #     @printf("Pr(#positive ≥ %d or #negative ≥ %d) = %.3f\n",sr.df[1,:observed],sr.df[2,:observed],sr.p_both)
# # end
# #

function dir(str::String)
  if occursin("*",str) || occursin("?",str)
    printdir(glob(str))
  elseif str == "." || str == ".."
    printdir(glob(string(str,"/*")))
  elseif str == ""
    printdir(readdir())
  else
    printdir(readdir(str))
  end
end
dir() = dir("")

function printdir(vstr::Vector{String})

  for i=1:length(vstr)
    vstr[i] = isdir(vstr[i]) ? string(vstr[i],"/") : vstr[i]
  end
  maxlen = maximum(length.(vstr))

  for i=1:length(vstr)
    println(rpad(vstr[i],maxlen),"  ", humanReadableByteCountBin(stat(vstr[i]).size))
  end
end

units   = ["K",  "M",  "G",  "T",  "P"]
bin_div = [0x1p10, 0x1p20, 0x1p30, 0x1p40, 0x1p40]

function humanReadableByteCountBin(bytes::Int64)
    b = bytes == typemin(Int64) ? typemax(Int64) : abs(bytes)
    if b < 1024
        return @sprintf("%d B", b)
    end
    sf = 40
    for i in 1:5
        if b < 0xfffcccccccccccc >> sf
            if i == 5
                b >>= 10
            end
            return @sprintf("%.1f %sB", b / bin_div[i], units[i])
        end
        sf -= 10
    end
    b = (b >> 20) / 0x1p40
    return @sprintf("%.1f EB", b)
end


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl
=#

matrixregression_insts = [
    vcat((4, 3), [(n, m) for n in ns]) # includes compile run
    for (m, ns) in ((15, vcat(50:50:150, 250, 500:500:3000)),
        (30, vcat(50:50:150, 250:250:1500)))
    ]

insts = OrderedDict()
insts["nat"] = (nothing, matrixregression_insts)
insts["ext"] = (:SOCExpPSD, matrixregression_insts)
return (MatrixRegressionJuMP, insts)


text:
import Base.show

struct ParameterSet
    arg_names::Vector{Symbol}
    arg_defaults::Vector
    arg_ranges::Vector
    arg_types::Vector{<:Type}
    n_args::Int
    #TODO: parameter constraints (e.g. ensure one parameter always greater than another)
    #TODO: refactor out the arg_ prefix (its redundant if they all start with it)
    function ParameterSet(arg_names::Vector{Symbol},
                          arg_defaults::Vector,
                          arg_ranges::Vector=[x:x for x in arg_defaults])
        @assert length(arg_names) == length(arg_defaults) == length(arg_ranges)
        @assert eltype.(arg_defaults) == eltype.(arg_ranges)
        arg_types::Vector{<:Type} = eltype.(arg_defaults)
        return new(arg_names, arg_defaults, arg_ranges, arg_types, length(arg_names))
    end
end

function count_runs(ps::ParameterSet)::Int
    n_runs = 1
    @inbounds for i in 1:ps.n_args
        n_runs *= length(ps.arg_ranges[i])
    end
    return n_runs
end

function generate_dict(ps::ParameterSet, arg_values::Vector)::Dict{Symbol,Any}
    out_dict = Dict{Symbol,Any}()
    for j in 1:ps.n_args
        out_dict[ps.arg_names[j]] = arg_values[j]
    end
    return out_dict
end

function show(io::IO, ps::ParameterSet)::Nothing
    println(io, "# Parameters:")
    @inbounds for i in 1:ps.n_args
        println(io, TAB, "($i) $(ps.arg_names[i])  →  $(ps.arg_defaults[i])  ∈  {$(string(ps.arg_ranges[i]))} :: $(ps.arg_types[i])")
    end
end

function generate_combinations(ps::ParameterSet)
    A = collect(Iterators.product(ntuple(i->ps.arg_ranges[i], length(ps.arg_ranges))...))
    B = A[:]
    T = Tuple(Array{arg_type}(undef, size(B,1)) for arg_type in ps.arg_types)
    for j in 1:length(ps.arg_ranges), i in 1:size(B,1)
        T[j][i] = B[i][j]
    end
    return hcat(T...)
end

max_stars_repo_name: lithiferous/Strathreads

max_issues_repo_name: lithiferous/Strathreads

max_forks_repo_name: lithiferous/Strathreads


text:
#####################################
#DiracVector#########################
#####################################

type DiracVector{S<:Single, T<:Union(Number, ScalarExpr)} <: Dirac
	coeffs::SparseMatrixCSC{T, Int}
	basis::AbstractBasis{S}
	function DiracVector(coeffs, basis)
		if S<:Ket
			@assert size(coeffs)==(length(basis),1) "Dimensions of coefficient array do not match basis"
			return new(coeffs, basis)
		else
			@assert size(coeffs)==(1,length(basis)) "Dimensions of coefficient array do not match basis"
			return new(coeffs, basis)
		end
	end
end

DiracVector{S,T<:Union(Number, ScalarExpr)}(coeffs::SparseMatrixCSC{T}, basis::AbstractBasis{S}) = DiracVector{S,T}(coeffs, basis)
DiracVector{S}(coeffs::SparseMatrixCSC, basis::AbstractBasis{S}) = DiracVector{S,ScalarExpr}(convert(SparseMatrixCSC{ScalarExpr},coeffs), basis)
DiracVector{S}(coeffs::SparseMatrixCSC{Any}, basis::AbstractBasis{S}) = DiracVector{S,ScalarExpr}(convert(SparseMatrixCSC{ScalarExpr},coeffs), basis)

DiracVector(coeffs::Array, basis::AbstractBasis) = DiracVector(qsparse(coeffs), basis)

DiracVector(coeffs, s::State...) = DiracVector(coeffs, basis(collect(s)))
dvec = DiracVector

convert{S,T}(::Type{DiracVector{S,T}}, dv::DiracVector{S}) = dvec(convert(SparseMatrixCSC{T}, dv.coeffs), dv.basis)

#####################################
#Access Functions####################
#####################################
bsym(d::DiracVector) = bsym(d.basis)

#####################################
#Boolean Functions###################
#####################################
isequal(a::DiracVector, b::DiracVector) = isequal(a.coeffs, b.coeffs) && a.basis==b.basis
==(a::DiracVector, b::DiracVector) = a.coeffs==b.coeffs && a.basis==b.basis
isdual{K<:Ket,B<:Bra}(a::DiracVector{K}, b::DiracVector{B}) = isdual(a.basis, b.basis) && a.coeffs==b.coeffs'
isdual{K<:Ket,B<:Bra}(a::DiracVector{B}, b::DiracVector{K}) = isdual(b,a)
isdual(a::DiracVector, b::DiracVector) = false

#####################################
#Array-like Functions################
#####################################
copy(d::DiracVector) = dvec(copy(d.coeffs), copy(d.basis))

ctranspose(d::DiracVector) = dvec(d.coeffs', d.basis')
size(d::DiracVector, args...) = size(d.coeffs, args...)
getindex{K<:Ket}(d::DiracVector{K}, x) = d.coeffs[x,1]
getindex{B<:Bra}(d::DiracVector{B}, x) = d.coeffs[1,x]
getindex(d::DiracVector, x, y) = d.coeffs[x,y]
setindex!{K<:Ket}(d::DiracVector{K}, y, x) = setindex!(d.coeffs, y, x, 1)
setindex!{B<:Bra}(d::DiracVector{B}, y, x) = setindex!(d.coeffs, y, 1, x)
setindex!(d::DiracVector, y, x...) = setindex!(d.coeffs, y, x...)

for op=(:endof, :ndims, :eltype, :length, :find, :findn, :findnz, :nnz, :ndims, :countnz)
	@eval ($op)(d::DiracVector) = ($op)(d.coeffs)
end

consnz_dvec{K<:Ket}(d::DiracVector{K}, nz_info) = dvec(nz_info[3], basis(d.basis[nz_info[1]]))
consnz_dvec{B<:Bra}(d::DiracVector{B}, nz_info) = dvec(nz_info[3], basis(d.basis[nz_info[2]]))

filternz(d::DiracVector) = consnz_dvec(d, findnz(d))
in(s::State, d::DiracVector) = in(s,d.basis)
in(n, dm::DiracVector) = in(n, dm.coeffs)

#####################################
#Dict-like Functions#################
#####################################
getpos(d::DiracVector, s::State) = get(d.basis, s)
get(d::DiracVector, s::State, notfound) = in(s,d.basis) ? d[getpos(d, s)] : notfound
get(d::DiracVector, s::State) = d[getpos(d, s)]
in(s::State, d::DiracVector) = in(s, d.basis)

#####################################
#Show Functions######################
#####################################

summary{K<:Ket}(d::DiracVector{K}) = "$(size(d,1))x1 $(typeof(d))"
summary{B<:Bra}(d::DiracVector{B}) = "1x$(size(d,2)) $(typeof(d))"

function showcompact(io::IO, d::DiracVector)
	if length(d)==0
		print(io, "$(typeof(d))[]")
	else
		tempio = IOBuffer()
		print(tempio, [" + $(d.coeffs[i])$(d.basis[i])" for i in find(d)]...)
		print(io, takebuf_string(tempio)[3:end])
	end
end
function show{S}(io::IO, d::DiracVector{S})
	if length(d)==0
		print(io, "$(typeof(d))[]")
	else
		println(io, summary(d))
		table = cell(length(d), 2)
		if length(d)>=52
			for i=1:25
				table[i,1]= d.coeffs[i]
				table[i,2]= d.basis[i]
			end
			table[26:(length(d)-25),:] = 0 # prevents access to undefined reference
			for i=(length(d)-25):length(d)
				table[i,1]= d.coeffs[i]
				table[i,2]= d.basis[i]
			end
		else
			for i=1:length(d)
				table[i,1]= d.coeffs[i]
				table[i,2]= d.basis[i]
			end
		end
		temp_io = IOBuffer()
		if S<:Ket
			show(temp_io, table)
		else
			show(temp_io, [transpose(table[:,2]), transpose(table[:,1])])
		end
		io_str = takebuf_string(temp_io)
		io_str = io_str[searchindex(io_str, "\n")+1:end]
		print(io, io_str)
	end
end

function showsp{S<:Union(Bra, Ket),T<:Union(Float64,Int)}(io::IO, dv::DiracVector{S,T})
	numpad = length("$(round(maximum(dv.coeffs.nzval), 4))")
	showsp_dv(io, dv, print_el = ((io,el)->print(io, rpad("$(round(el, 4))", numpad))))
end

showsp(io::IO, dv::DiracVector) = showsp_dmat(io, dv)

function showsp_dv(io::IO, dv::DiracVector; print_el=((io,el)->showcompact(io, el)))
 	print(io, "$(size(dv,1))x$(size(dv,2)) DiracVector with $(nnz(dv)) $(eltype(dv)) entries: ")
	rows = Base.tty_rows()
	S = dv.coeffs
    half_screen_rows = div(rows - 8, 2)
    pad = ndigits(max(S.m,S.n))
    k = 0
    sep = "\n\t"
    for col = 1:S.n, k = S.colptr[col] : (S.colptr[col+1]-1)
        if k < half_screen_rows || k > nnz(S)-half_screen_rows
            print(io, sep, '[', rpad(S.rowval[k], pad), ", ", lpad(col, pad), "]  =  ")
            print_el(io, S.nzval[k])
            print(io, "  $(dv.basis[S.rowval[k]])")
        elseif k == half_screen_rows
            print(io, sep, '\u22ee')
        end
        k += 1
    end
end

#####################################
#Function-passing Functions##########
#####################################

find(f::Function, d::DiracVector) = find(f, d.coeffs)
findstates(f::Function, d::DiracVector) = find(f, d.basis)

map(f::Function, d::DiracVector)=dvec(map(f, d.coeffs), d.basis)

function loadcoeffs!(coeffs, d, newbasis)
	for i=1:length(newbasis)
		coeffs[i] = get(d, newbasis[i])
	end
	return coeffs
end

function filtercons{S,T}(newbasis::AbstractBasis, d::DiracVector{S,T})
	return dvec(loadcoeffs!(S<:Ket ? Array(T, length(newbasis)) : Array(T, 1, length(newbasis)), d, newbasis), newbasis)
end

filterstates{S,T}(f::Function, d::DiracVector{S,T}) = filtercons(filter(f, d.basis), d)
filtercoeffs(f::Function, d::DiracVector) = filtercons(basis(d.basis[find(map(f, d.coeffs))]), d)

qeval(f::Function, d::DiracVector) = dvec(vcat([qeval(f,x) for x in d.coeffs]...), d.basis) #vcat forces correct typing, but this is hacky

#####################################
#Arithmetic Operations###############
#####################################

for op=(:.*,:.-,:.+,:./,:.^)
	@eval ($op)(a::DiracVector, b::DiracVector) = dvec(($op)(a.coeffs,b.coeffs), a.basis)
	@eval ($op)(n::DiracCoeff, d::DiracVector) = dvec(($op)(n,d.coeffs), d.basis)
	@eval ($op)(d::DiracVector, n::DiracCoeff) = dvec(($op)(d.coeffs,n), d.basis)
end

/(dv::DiracVector, c::DiracCoeff) = dvec(dv.coeffs/c, dv.basis)
*(dv::DiracVector, c::DiracCoeff) = dvec(dv.coeffs*c, dv.basis)
*(c::DiracCoeff, dv::DiracVector) = dvec(c*dv.coeffs, dv.basis)

kron(c::DiracCoeff, d::DiracVector) = c*d
kron(d::DiracVector, c::DiracCoeff) = d*c

-(d::DiracVector) = -1*d

ireduce(s::State, d::DiracVector) = reduce(+,[d[i]*inner(s,d.basis[i]) for i=1:length(d)])
ireduce(d::DiracVector, s::State) = reduce(+,[d[i]*inner(d.basis[i],s) for i=1:length(d)])
ireduce(a::DiracVector, b::DiracVector) = reduce(+,[a[i]*b[j]*inner(a.basis[i],b.basis[j]) for i=1:length(a), j=1:length(b)])
ireduce(a::DiracVector, b::DiracVector, eqbasis::Bool) = reduce(+,[a[i]*b[i]*inner(a.basis[i],b.basis[i]) for i=1:length(a)])

inner{b,T1,T2}(s::State{Bra{b,T1}}, d::DiracVector{Ket{b,T2}}) = get(d, s', 0)
inner{b,T1,T2}(d::DiracVector{Bra{b,T1}},s::State{Ket{b,T2}}) = get(d, s', 0)

inner{B<:Bra, K<:Ket}(s::State{B}, d::DiracVector{K}) = ireduce(s,d)
inner{B<:Bra, K<:Ket}(d::DiracVector{B}, s::State{K}) = ireduce(d,s)
inner{B<:Bra, K<:Ket, N1<:Number, N2<:Number}(a::DiracVector{B, N1}, b::DiracVector{K, N2}) = isdual(a.basis, b.basis) ? (a.coeffs*b.coeffs)[1] : ireduce(a,b)
inner{B<:Bra, K<:Ket}(a::DiracVector{B}, b::DiracVector{K}) = isdual(a.basis, b.basis) ? ireduce(a,b,true) : ireduce(a,b)

for t=(:Bra, :Ket)
	@eval begin
	tensor{A<:($t), B<:($t)}(a::DiracVector{A}, b::DiracVector{B}) = dvec(kron(a.coeffs, b.coeffs), tensor(a.basis, b.basis))
	tensor{A<:($t), B<:($t)}(a::DiracVector{A}, b::State{B}) = dvec(a.coeffs, tensor(a.basis, b))
	tensor{A<:($t), B<:($t)}(a::State{B}, b::DiracVector{A}) = dvec(b.coeffs, tensor(a, b.basis))
	kron{A<:($t), B<:($t)}(a::DiracVector{A}, b::DiracVector{B}) = tensor(a,b)
	kron{A<:($t), B<:($t)}(a::DiracVector{A}, b::State{B}) = tensor(a,b)
	kron{A<:($t), B<:($t)}(a::State{B}, b::DiracVector{A}) = tensor(a,b)
	end
end

kron{K<:Ket, B<:Bra}(a::DiracVector{K}, b::State{B}) = dmat(a.coeffs, a.basis, basis(b))
kron{K<:Ket, B<:Bra}(a::State{K}, b::DiracVector{B}) = dmat(b.coeffs, basis(a), b.basis)
kron{K<:Ket, B<:Bra}(a::DiracVector{K}, b::DiracVector{B}) = dmat(a.coeffs*b.coeffs, a.basis, b.basis)
kron{K<:Ket, B<:Bra}(a::DiracVector{B}, b::DiracVector{K}) = kron(b,a)

kron(c::DiracCoeff, s::State) = dvec([c], basis(s))
kron(s::State, c::DiracCoeff) = kron(c,s)

*{K<:Ket, B<:Bra}(a::DiracVector{K}, b::DiracVector{B}) = kron(a,b)
*{K<:Ket, B<:Bra}(a::State{K}, b::DiracVector{B}) = kron(a,b)
*{K<:Ket, B<:Bra}(a::DiracVector{K}, b::State{B}) = kron(a,b)
*(c::DiracCoeff, s::State) = kron(c,s)
*(s::State, c::DiracCoeff) = kron(c,s)
-(s::State) = kron(-1,s)

*{B<:Bra, K<:Ket}(a::DiracVector{B}, b::DiracVector{K}) = inner(a,b)
*{B<:Bra, K<:Ket}(a::DiracVector{B}, b::State{K}) = inner(a,b)
*{B<:Bra, K<:Ket}(a::State{B}, b::DiracVector{K}) = inner(a,b)

function addstate!(d,s)
	d[getpos(d,s)] = 1+get(d, s)
	return d
end

function +{K<:Ket}(d::DiracVector{K}, s::State{K})
	if in(s, d.basis)
		return addstate!(copy(d),s)
	else
		return dvec(vcat(d.coeffs, speye(1)), bcat(d.basis,s))
	end
end

function +{B<:Bra}(d::DiracVector{B}, s::State{B})
	if in(s, d.basis)
		return addstate!(copy(d),s)
	else
		return dvec(hcat(d.coeffs, speye(1)), bcat(d.basis,s))
	end
end

function +{K<:Ket}(s::State{K}, d::DiracVector{K})
	if in(s, d.basis)
		return addstate!(copy(d),s)
	else
		return dvec(vcat(speye(1), d.coeffs), bcat(s,d.basis))
	end
end

function +{B<:Bra}(s::State{B}, d::DiracVector{B})
	if in(s, d.basis)
		return addstate!(copy(d),s)
	else
		return dvec(hcat(speye(1),d.coeffs), bcat(s,d.basis))
	end
end

function add_dvec!(d::DiracVector, y::DiracVector)
	for i=1:length(y)
		d = d+y.basis[i]
		d[getpos(d, y.basis[i])] = get(d, y.basis[i]) + y[i] - 1
	end
	return d
end

for t=(:Bra,:Ket)
	@eval begin
	function +{S<:($t)}(x::DiracVector{S}, y::DiracVector{S})
		if x.basis==y.basis
			return dvec(x.coeffs+y.coeffs, x.basis)
		else
			return add_dvec!(copy(x), y)
		end
	end
	-{S1<:($t),S2<:($t)}(d::DiracVector{S1}, s::State{S2}) = d+(-s)
	-{S1<:($t),S2<:($t)}(s::State{S1}, d::DiracVector{S2}) = s+(-d)
	-{S1<:($t),S2<:($t)}(a::DiracVector{S1}, b::DiracVector{S2}) = a+(-b)
	end
end

+{K1<:Ket,K2<:Ket}(a::State{K1}, b::State{K2}) = a==b ? 2*a : dvec([1, 1], basis([a,b]))
+{B1<:Bra,B2<:Bra}(a::State{B1}, b::State{B2}) = a==b ? 2*a : dvec([1 1], basis([a,b]))
-{K1<:Ket,K2<:Ket}(a::State{K1}, b::State{K2}) = a==b ? dvec(spzero(1),basis(a)) : dvec([1, -1], basis([a,b]))
-{B1<:Bra,B2<:Bra}(a::State{B1}, b::State{B2}) = a==b ? dvec(spzero(1),basis(a)) : dvec([1 -1], basis([a,b]))

norm(d::DiracVector, p::Int=2) = reduce(+,map(i->abs(i)^p, d.coeffs))^(1/p)
norm{S<:Single,N<:Number}(d::DiracVector{S,N}, p::Int=2) = norm(d.coeffs)
normalize(d::DiracVector) = dvec((1/norm(d))*d.coeffs, d.basis)


text:
module QDatesCLI

using Comonicon
using QDates
using Dates

abstract type ARGType end

struct ShowQDate <: ARGType
    src::String
    qdt::QDate
end

struct ShowQDateAndDayname <: ARGType
    src::String
    qdt::QDate
end

struct ShowDate <: ARGType
    src::String
    dt::Date
end

struct ShowError <: ARGType
    src::String
    msg::String
end

function showresult(res::ShowQDate)
    println("$(res.src)は$(res.qdt)です。")
end

function showresult(res::ShowQDateAndDayname)
    println("$(res.src)は$(res.qdt)、$(QDates.dayname(res.qdt))です。")
end

function showresult(res::ShowDate)
    println("$(res.src)は西暦$(res.dt)です。")
end

function showresult(res::ShowError)
    println("$(res.src): $(res.msg)")
end

function parseargs(args::String...)
    results = ARGType[]
    if length(args) == 0
        push!(results, ShowQDateAndDayname("今日（$(Dates.today())）", QDates.today()))
    else
        for arg in args
            m = match(r"^(\d{3,4})-(\d{1,2})-(\d{1,2})$", arg)
            if !isnothing(m)
                try
                    dt = Date(parse(Int, m[1]), parse(Int, m[2]), parse(Int, m[3]))
                    push!(results, ShowQDate(arg, QDate(dt)))
                catch ex
                    if ex isa ArgumentError
                        push!(results, ShowError(arg, "不正な日付です。"))
                    else
                        push!(results, ShowError(arg, sprint(showerror, ex)))
                        rethrow(ex)
                    end
                end
                continue
            end
            m = match(r"^(\d{3,4})/(\d{1,2})/(\d{1,2})$", arg)
            if !isnothing(m)
                try
                    dt = Date(parse(Int, m[1]), parse(Int, m[2]), parse(Int, m[3]))
                    push!(results, ShowQDate(arg, QDate(dt)))
                catch ex
                    if ex isa ArgumentError
                        push!(results, ShowError(arg, "不正な日付、または対応範囲外の日付です。"))
                    else
                        push!(results, ShowError(arg, sprint(showerror, ex)))
                        rethrow(ex)
                    end
                end
                continue
            end
            m = match(r"^旧(\d{3,4})年(閏)?(\d{1,2})月(\d{1,2})日$", arg)
            if !isnothing(m)
                try
                    qdt = QDate(parse(Int, m[1]), parse(Int, m[3]), !isnothing(m[2]), parse(Int, m[4]))
                    push!(results, ShowDate(arg, Date(qdt)))
                catch ex
                    if ex isa ArgumentError
                        push!(results, ShowError(arg, "不正な日付、または対応範囲外の日付です。"))
                    else
                        push!(results, ShowError(arg, sprint(showerror, ex)))
                        rethrow(ex)
                    end
                end
                continue
            end
            m = match(r"^(一昨[昨々]日)|(一昨日)|(昨日)|(今日)|(明日)|(明後日)|(明[明々]後日)$", arg)
            if !isnothing(m)
                try
                    segment = Day(findfirst(!isnothing, m.captures) - 4)
                    push!(results, ShowQDate(arg, QDates.today() + segment))
                catch ex
                    if ex isa ArgumentError
                        push!(results, ShowError(arg, "不正な日付、または対応範囲外の日付です。"))
                    else
                        push!(results, ShowError(arg, sprint(showerror, ex)))
                        rethrow(ex)
                    end
                end
                continue
            end
            push!(results, ShowError(arg, "有効な日付として認識できません。"))
        end
    end
    results
end

"""
Show Kyūreki (旧暦) or Gregorian date for the corresponding date.
Multiple `<date_like>` arguments may be specified, and displays all the results.
When no arguments are specified, show Kyūreki information of `今日` (today).

- `<date_like>`: one of:
    + Gregorian date (`yyyy-mm-dd`|`yyyy/mm/dd`)
    + Kyūreki date (`旧yyyy年mm月dd日`|`旧yyyy年閏mm月dd日`)
    + one of date_keywords ("今日"|"明日"|"明後日"|"明々後日"|"昨日"|"一昨日"|"一昨々日")

# Usage

qdate [<date_like> [<date_like> ...]]

# Arguments

- `args`: Some `<date_like>`s.
"""
@main function qdate(args::String...)
    for res in parseargs(args...)
        showresult(res)
    end
end
# @main function main(arg)
#     println("arg: $arg")
# end

end # module


text:
using Distributions
using Base.Test

N = 10^5

d = fit(DiscreteUniform, rand(DiscreteUniform(10, 15), N))
@test isa(d, DiscreteUniform)
@test min(d) == 10
@test max(d) == 15

d = fit(Bernoulli, rand(Bernoulli(0.7), N))
@test isa(d, Bernoulli)
@test_approx_eq_eps mean(d) 0.7 0.01

d = fit(Beta, rand(Beta(1.3, 3.7), N))
@test isa(d, Beta)
@test_approx_eq_eps d.alpha 1.3 0.1
@test_approx_eq_eps d.beta  3.7 0.1

d = fit(Binomial, 100, rand(Binomial(100, 0.3), N))
@test isa(d, Binomial)
@test d.size == 100
@test_approx_eq_eps d.prob 0.3 0.01

d = fit(Exponential, rand(Exponential(0.5), N))
@test isa(d, Exponential)
@test_approx_eq_eps mean(d) 0.5 0.01

x = rand(Normal(11.3, 5.7), N)
d = fit(Normal, x)
@test isa(d, Normal)
@test_approx_eq_eps mean(d) 11.3 0.1
@test_approx_eq_eps std(d) 5.7 0.1

d = fit_mle(Normal, suffstats(Normal, x))
@test isa(d, Normal)
@test_approx_eq_eps mean(d) 11.3 0.1
@test_approx_eq_eps std(d) 5.7 0.1

d = fit(Uniform, rand(Uniform(1.2, 10.3), N))
@test isa(d, Uniform)
@test 1.2 <= min(d) <= max(d) <= 10.3
@test_approx_eq_eps min(d) 1.2 0.02
@test_approx_eq_eps max(d) 10.3 0.02

d = fit(Gamma, rand(Gamma(3.9, 2.1), N))
@test isa(d, Gamma)
@test_approx_eq_eps d.shape 3.9 0.1
@test_approx_eq_eps d.scale 2.1 0.2

d = fit(Geometric, rand(Geometric(0.3), N))
@test isa(d, Geometric)
@test_approx_eq_eps d.prob 0.3 0.01

d = fit(Laplace, rand(Laplace(5.0, 3.0), N))
@test isa(d, Laplace)
@test_approx_eq_eps d.location 5.0 0.1
@test_approx_eq_eps d.scale 3.0 0.2

d = fit(Poisson, rand(Poisson(8.2), N))
@test isa(d, Poisson)
@test_approx_eq_eps mean(d) 8.2 0.2


text:
# # Landcover data

# In this example, we will look at landcover data in Corsica

using SimpleSDMLayers
using Plots

# To avoid loading the (rather large) dataset of land cover at once, we will
# rely on a bounding box:

boundaries = (left=8.25, right=10.0, bottom=41.2, top=43.2)

# The numbers corresponding to the layers can be found in the documentation for
# `SimpleSDMPredictor` - urban area is 9.

urban = SimpleSDMPredictor(EarthEnv, LandCover, 9; boundaries...)

# This dataset is returning data as `UInt8` (as it represents a proportion of
# the pixel occupied by the type), but this is not something that can be plotted
# efficiently (because we rely on `NaN` to indicate no data, and `UInt8` has no
# `NaN`). So in the next step, we will manipulate this object a little bit to
# have something more workable.

urban = convert(Float16, urban)

# Why `Float16`? It is the smallest floating point type with a `NaN` value
# (`NaN16`). We will replace the values of 0 by `nothing`, to only see the
# pixels with *at least* some urban cover:

replace!(urban, zero(eltype(urban)) => nothing)

# With this done, we can plot the results:

plot(urban; c=:berlin, clim=(0, 100))

# Unsuprisingly to anyone who had the chance to visit Corsica, it is not a very
# densely urbanized island. This is a good time to question whether we can look
# at (i) which landcover type dominates within each pixel, and (ii) how
# heterogeneous the land use within each pixel is.

# First, we will download all values for the landcover layers, including open
# water (12).

landcover =
    convert.(
        Float16, SimpleSDMPredictor(EarthEnv, LandCover, 1:12; full=false, boundaries...)
    )

plot(
    plot.(landcover, leg=false, c=:oleron, clim=(0, 100))...;
    leg=false,
    grid=:none,
    frame=:none,
)

# To perform the actual analysis, we will define a `shannon` function, which
# will return the entropy of the land use categories:

function shannon(x)
    v = filter(n -> n > zero(eltype(x)), x)
    length(v) == 0 && return NaN
    v = v ./ sum(v)
    return -sum(v .* log2.(v))
end

entropy = mosaic(shannon, landcover)

# We can also get the index of the most common layer within the pixel:

consensus = mosaic(x -> last(findmax(x)), landcover)

# We may not be *that* interested in fully open water, so let's define a mask to
# remove it:

openwater = broadcast(!isequal(12), consensus)

# All that's left to do is to plot after applying this mask, and we now get a
# map of most common land cover type (left), and land cover heterogeneity
# (right)

p1 = plot(mask(openwater, consensus); c=:terrain, frame=:none)
p2 = plot(mask(openwater, entropy); c=:lapaz, frame=:none)
plot(p1, p2)


text:
using Luxor

Drawing(625,625,"logo.svg")

origin()
background("transparent")
Luxor.scale(1.7)
setfont("Lithos Pro Black", 180)
N,F,C = juliacircles()
settext("F",F+Point(0,12),halign="center",valign="center")
settext("N",N+Point(0,12),halign="center",valign="center")
settext("C",C+Point(0,12),halign="center",valign="center")

finish()

Drawing(625,625,"logo-dark.svg")

origin()
background("transparent")
Luxor.scale(1.7)
setfont("Lithos Pro Black", 180)
N,F,C = juliacircles()
setcolor("white")
settext("F",F+Point(0,12),halign="center",valign="center")
settext("N",N+Point(0,12),halign="center",valign="center")
settext("C",C+Point(0,12),halign="center",valign="center")

finish()
preview()


text:
using DSGE, Test, FileIO, Random, ModelConstructors
m = AnSchorfheide()
@testset "Test moments" begin
    @test DSGE.moments(m.parameters[findfirst(x -> x.key==:τ, m.parameters)]) == (2.0, 0.5)
    @test DSGE.moments(m.parameters[findfirst(x -> x.key==:e_y, m.parameters)]) == (0.1159846, 0.0)
    @test DSGE.moments(m.parameters[findfirst(x -> x.key==:σ_R, m.parameters)]) == (0.4, 4.0)
end

# This script currently just tests sample_λ, compute_Eλ

fp = dirname(@__FILE__)
save_λmat_noparallel = load("$(fp)/../reference/moments_poolmodel_outputs.jld2", "lammat_noparallel")
save_λmat_parallel = load("$(fp)/../reference/moments_poolmodel_outputs.jld2", "lammat_parallel")
Random.seed!(1793)
m = PoolModel("ss1")
θs = load("$(fp)/../reference/moments_poolmodel_inputs.jld2", "thetas")[1:2,:]
pred_dens = load("$(fp)/../reference/moments_poolmodel_inputs.jld2", "pred_dens")
λmat_noparallel = sample_λ(m, pred_dens, θs, 1)
λmat_parallel   = sample_λ(m, pred_dens, θs, 1; parallel = true)

@testset "Check sample_λ works correctly" begin
    @test @test_matrix_approx_eq save_λmat_noparallel λmat_noparallel
    @test @test_matrix_approx_eq save_λmat_parallel λmat_parallel
end

λvec = 0.5 * ones(2)
θchange = [.8 0. 1.] .* ones(2)
λhat_plush, λhat_t = compute_Eλ(m, 4, λvec, θchange)
λhat_plush_parallel, λhat_t_parallel = compute_Eλ(m, 4, λvec, θchange; parallel = true)
@testset "Check compute_Eλ works correctly" begin
    @test λhat_plush == 0.5
    @test λhat_plush_parallel == 0.5
    @test λhat_t == 0.5
    @test λhat_t_parallel == 0.5
end

Random.seed!(1793)
sm = PoolModel("ss1"; weight_type = :static)
sm <= Setting(:saveroot, "$(fp)/../reference/")
sm <= Setting(:hessian_path, "$(fp)/../reference/mh_hessian_poolmodel.h5")
sm <= Setting(:n_mh_simulations, 1)
save_sλmat_noparallel = load("$(fp)/../reference/moments_poolmodel_outputs.jld2", "slammat_noparallel_one")
save_sλmat_parallel = load("$(fp)/../reference/moments_poolmodel_outputs.jld2", "slammat_parallel_one")
sλmat_noparallel = sample_λ(sm, pred_dens, 1)
sλmat_parallel   = sample_λ(sm, pred_dens, 1; parallel = true)
rm(joinpath(saveroot(sm), "output_data/poolmodel"); recursive = true)

@testset "Check sample_λ works for a static pool" begin
    @test sλmat_noparallel == save_sλmat_noparallel
    @test sλmat_parallel == save_sλmat_parallel
end


text:
########################################
## File Name: belief_transition_test.jl
## Author: Haruki Nishimura (hnishimura@stanford.edu)
## Date Created: 2020/05/13
## Description: Test code for src/belief_transition_models.jl
########################################

using LinearAlgebra
using Random
using Distributions

@testset "Belief Transition Test" begin

# covObservePos
p = randn(2);
qi = randn(2);
C = SACBP.covObservePos(p, qi)
@test all(isapprox(C, diagm((0.001 + norm(p - qi)*0.01)*ones(2))))

# UKFPosRange with BelMvNormal
ukfModel = UKFPosRange();
bPos = BelMvNormal(0.,zeros(2),100*Matrix(1.0I, 2, 2));
bVecPos = [bPos for ii = 1:100];
p_robot = PhysPos(0.,[5.,7.]);
s = AugState(p_robot,bPos);
sV = AugState(p_robot,bVecPos);

b_new = trans(ukfModel,bPos,p_robot,0.01,2.,Matrix(1.0I, 2, 2))
@test all(isapprox(b_new.params, [3.24374, 4.62673, 94.7848, -7.45299, -7.45299, 89.3794], atol=1e-4))

Jp = trans_jacobi_auto_p(ukfModel,bPos,p_robot,0.01,2.,Matrix(1.0I, 2, 2))
@test all(isapprox(Jp, [0.851283   0.243418
                        0.251906   1.05732
                       -2.23007   -0.0795886
                       -1.63071   -1.25742
                       -1.63071   -1.25742
                       -0.114878  -3.42515  ], atol=1e-4))

Jb = trans_jacobi_auto_b(ukfModel,bPos,p_robot,0.01,2.,Matrix(1.0I, 2, 2))
@test all(isapprox(Jb, [0.148717  -0.243418    0.00303772  0.0   0.0437606  -0.0139536
                       -0.251906  -0.0573231  -0.0206755   0.0   0.027708    0.00430951
                        2.23007    0.0795886   0.975254    0.0  -0.147259    0.0310205
                        1.63071    1.25742     0.00498744  0.5   0.345869    0.00524383
                        1.63071    1.25742     0.00498744  0.5   0.345869    0.00524383
                        0.114878   3.42515     0.0645743   0.0  -0.140094    0.951848  ], atol=1e-4))

# UKFPosRange with VecBelMVNormal
yvec = [2. for ii = 1:100];
b_new_array = trans(ukfModel,bVecPos,p_robot,0.01,yvec,Matrix(1.0I, 2, 2));
@test all([b_new_array[ii].params == b_new.params for ii = 1:100])

Jp_array = trans_jacobi_auto_p(ukfModel,bVecPos,p_robot,0.01,yvec,Matrix(1.0I, 2, 2));
@test all([isapprox(Jp_array[ii], Jp) for ii = 1:100])

Jb_array = trans_jacobi_auto_b(ukfModel,bVecPos,p_robot,0.01,yvec,Matrix(1.0I, 2, 2));
@test all([isapprox(Jb_array[ii], Jb) for ii = 1:100])

# UKFPosRange with AugState and BelMvNormal
s_new = trans(ukfModel,s,0.01,2.,Matrix(1.0I, 2, 2))
@test s_new.p.pos == [5., 7.]
@test s_new.b.params == b_new.params
@test all(isapprox.(trans_jacobi_auto_p(ukfModel,s,0.01,2.,Matrix(1.0I,2,2)), Jp))
@test all(isapprox.(trans_jacobi_auto_b(ukfModel,s,0.01,2.,Matrix(1.0I,2,2)), Jb))

# UKFPosRange with AugState and VecBelMvNormal
sV_new = trans(ukfModel,sV,0.01,yvec,Matrix(1.0I,2,2))
@test sV_new.p.pos == [5., 7.]
@test all([all(isapprox.(b.params, b_new.params)) for b in sV_new.b])
@test all([all(isapprox.(J, Jp)) for J in trans_jacobi_auto_p(ukfModel,sV,0.01,yvec,Matrix(1.0I,2,2))])
@test all([all(isapprox.(J, Jb)) for J in trans_jacobi_auto_b(ukfModel,sV,0.01,yvec,Matrix(1.0I,2,2))])

trans(ukfModel,bVecPos,p_robot,0.01,Matrix(1.0I,2,2),MersenneTwister(1234))
trans(ukfModel,sV,0.01,Matrix(1.0I,2,2),MersenneTwister(1234))

# CPredictManipulate2D
xMP = PhysManipulate2D(0.,[0.1,0.2],0.5,[0.3,0.4],0.6,3.,5.,[10.,11.],0.7);
uMP = MControl2D(0.,6.,7.,8.);

predictEKFMP = CPredictManipulate2D();
bMP = BelMvNormal(xMP.t,vec(xMP),10*Matrix(1.0I, 11, 11));

bMP_new = trans(predictEKFMP,bMP,uMP,0.01,Matrix(1.0I, 11, 11))
@test bMP_new.t == 0.01
d = Distributions.MvNormal(bMP_new)
@test all(isapprox.(d.μ, [0.103, 0.204, 0.506, 0.3193, 0.4224, 0.491658, 3.0, 5.0, 10.0, 11.0, 0.7], atol=1e-4))

J1 = trans_jacobi_auto(predictEKFMP,bMP,uMP,Matrix(1.0I,11,11))
J2 = trans_jacobi(predictEKFMP,bMP,uMP,Matrix(1.0I,11,11))
@test all(isapprox.(J1, J2))

Hu = trans_u_coeff(predictEKFMP,bMP,Matrix(1.0I,11,11))

# DUpdateManipulate2D
updateEKFMP = DUpdateManipulate2D();
bMP_new2 = trans(updateEKFMP, bMP, uMP, ones(9), Matrix(1.0I, 9, 9))
@test bMP_new2.t == 0.0
d2 = Distributions.MvNormal(bMP_new2)
@test all(isapprox.(d2.μ, [-2.60194, -11.0021, 0.478504, 0.552211, 0.657014, 0.0276191, 2.31804, 8.64641, 9.45285, 9.87449, 0.595059], atol=1e-4))

trans_jacobi_auto(updateEKFMP,bMP,uMP,ones(9),Matrix(1.0I, 9, 9))

# CDEKFManipulate2D
cdEKFMP = CDEKFManipulate2D();
b_new_2 = trans(cdEKFMP,bMP,uMP,0.01,0.2,ones(9),Matrix(1.0I, 11, 11),Matrix(1.0I, 9, 9))
@test b_new_2.t == 0.2
d3 = Distributions.MvNormal(b_new_2);
@test all(isapprox.(d3.μ, [-4.91794, -9.74958, 0.293202, -1.55987, 4.92263, -0.437129, 3.12938, 6.42346, 9.27099, 9.46519, 0.814831], atol=1e-4))

end


text:
push!(LOAD_PATH, ".")

module GAInitialization

using Population

abstract type InitRandom end

function initializePopulation!(pop::Population.AbstractPopulation, pop_size::Integer, ::Type{InitRandom})
  Population.populateRandom!(pop, pop_size)
  Population.evalFitness!(pop)
end

end

max_stars_repo_name: Markhyz/Julia-Genetic-Algorithm

max_issues_repo_name: Markhyz/Julia-Genetic-Algorithm

max_forks_repo_name: Markhyz/Julia-Genetic-Algorithm


text:
module Binance

using PyCall, DataFrames

export
    set_auth_token,
    parsekeyfile,
    Client,
    get_all_tickers,
    get_order_book,
    get_symbol_info,
    get_exchange_info,
    get_ticker,
    get_symbol_ticker,
    get_symbols,
    base_symbols,
    split_symbols,
    get_symbol_map,
    get_order_book

@pyimport binance.client as bn

const KEYFILE = joinpath(@__DIR__, "auth_token.txt")

const BASE_SYMBOLS = ["BTC", "ETH", "BNB", "USDT"]

Client() = Client(parsekeyfile()...)
Client(api_key::AbstractString, api_secret::AbstractString) = bn.Client(api_key, api_secret)

function set_auth_token(api_key::AbstractString, api_secret::AbstractString, keyfile::AbstractString=KEYFILE)
    f = open(keyfile, "w")
    println(f, api_key)
    println(f, api_secret)
    close(f)
end

function parsekeyfile(keyfile::AbstractString=KEYFILE)
    f = open(keyfile)
    api_key = readline(f)
    api_secret = readline(f)
    close(f)
    (api_key, api_secret)
end

get_all_tickers(client) = client[:get_all_tickers]()
get_order_book(client, sym::Union{Symbol,String}) = client[:get_order_book](symbol=sym)
get_symbol_info(client, sym::Union{Symbol,String}) = client[:get_symbol_info](sym)
get_exchange_info(client) = client[:get_exchange_info]()
get_ticker(client,sym::Union{Symbol,String}) = client[:get_ticker](symbol=sym)
get_symbol_ticker(client,sym::Union{Symbol,String}) = client[:get_symbol_ticker](symbol=sym)

function get_symbols(client)
    resp = get_all_tickers(client)
    [d[] for d in resp]
end

base_symbols() = BASE_SYMBOLS

function get_all_tickers(::Type{DataFrame}, client; verbose::Bool=false)
    df = DataFrame([String, Float64],[:symbol, :price], 0)
    for x in get_all_tickers(client)
        verbose && @show x["symbol"]
        push!(df, [x["symbol"], parse(Float64,x["price"])])
    end
    df
end
function get_all_tickers(::Type{Dict}, client,
                         df::DataFrame=get_all_tickers(DataFrame, client);
                         verbose::Bool=false)
    d = Dict{String,Float64}()
    for x in eachrow(df)
        verbose && @show x[:symbol]
        d[x[:symbol]] = x[:price]
    end
    d
end

function get_symbol_map(::Type{DataFrame}, client; verbose::Bool=false)
    df = DataFrame([String, String, String],[:symbol, :from_symbol, :to_symbol], 0)
    for x in get_all_tickers(client)
        verbose && @show x["symbol"]
        from_symbol, to_symbol = split_symbol(client, x["symbol"])
        push!(df, [x["symbol"], from_symbol, to_symbol])
    end
    df
end
function get_symbol_map(::Type{Dict}, client,
                        df::DataFrame=get_symbol_map(DataFrame, client);
                        verbose::Bool=false)
    d = Dict{String,Tuple{String,String}}()
    for x in eachrow(df)
        d[x[:symbol]] = (x[:from_symbol],x[:to_symbol])
    end
    d
end
function split_symbol(client, s::String)
    d = get_symbol_info(client, s)
    (d["quoteAsset"], d["baseAsset"]) #from_symbol, to_symbol
end

end # module


text:
using ImplicitEquations
using Test

f(x,y) = y-x
g(x,y) = y+x

ImplicitEquations.GRAPH(f ⩵ 0, -5, 5, -5, 5, 2^4, 2^4)
ImplicitEquations.GRAPH(f ≦ 0, -5, 5, -5, 5, 2^4, 2^4)
ImplicitEquations.GRAPH(f ≪  g, -5, 5, -5, 5, 2^4, 2^4)

ImplicitEquations.GRAPH(f ≶  g, -5, 5, -5, 5, 2^4, 2^4)
ImplicitEquations.GRAPH(f ≷  g, -5, 5, -5, 5, 2^4, 2^4)
ImplicitEquations.GRAPH(Neq(f, g), -5, 5, -5, 5, 2^4, 2^4)

ImplicitEquations.GRAPH(Eq(f, g), -5, 5, -5, 5, 2^4, 2^4)
ImplicitEquations.GRAPH(f ⩵ g, -5, 5, -5, 5, 2^4, 2^4)

# Issue #30: relax assumptions in predicates
struct F end
(::F)(x,y) = x^2 + y^2
ImplicitEquations.GRAPH(f ⩵ 0, -5, 5, -5, 5, 2^4, 2^4)


text:
# ---
# title: 127. Word Ladder
# id: problem127
# author: AquaIndigo
# date: 2020-11-07
# difficulty: Medium
# categories: Breadth-first Search
# link: <https://leetcode.com/problems/word-ladder/description/>
# hidden: true
# ---
#
# Given two words ( _beginWord_ and _endWord_ ), and a dictionary's word list,
# find the length of shortest transformation sequence from _beginWord_ to
# _endWord_ , such that:
#
#   1. Only one letter can be changed at a time.
#   2. Each transformed word must exist in the word list.
#
# **Note:**
#
#   * Return 0 if there is no such transformation sequence.
#   * All words have the same length.
#   * All words contain only lowercase alphabetic characters.
#   * You may assume no duplicates in the word list.
#   * You may assume _beginWord_ and _endWord_ are non-empty and are not the same.
#
# **Example 1:**
#
#
#
#     Input:
#     beginWord = "hit",
#     endWord = "cog",
#     wordList = ["hot","dot","dog","lot","log","cog"]
#
#     Output: 5
#
#     Explanation: As one shortest transformation is "hit" -> "hot" -> "dot" -> "dog" -> "cog",
#     return its length 5.
#
#
# **Example 2:**
#
#
#
#     Input:
#     beginWord = "hit"
#     endWord = "cog"
#     wordList = ["hot","dot","dog","lot","log"]
#
#     Output:  0
#
#     Explanation:  The endWord "cog" is not in wordList, therefore no possible ** ** transformation.
#
#
#
## @lc code=start
using LeetCode
using DataStructures
function ladder_length(begin_word::String, end_word::String, word_list::Vector{String})::Int
    function isadj(s1, s2)
        flg = false
        for i in 1:length(s1)
            if s1[i] != s2[i]
                if flg
                    return false
                else
                    flg = true
                end
            end
        end
        return true
    end
    if !(end_word in word_list)
        return 0
    end
    !(begin_word in word_list) && push!(word_list, begin_word)
    s, t = findall(x -> x == begin_word, word_list)[1],
    findall(x -> x == end_word, word_list)[1]
    ## println(s, ", ", t)
    ### wl = collect.(word_list)
    len = length(word_list)
    edges = [Set{Int}() for i in 1:len]
    ### construct adj list
    for i in 1:len
        for j in (i + 1):len
            if isadj(word_list[i], word_list[j])
                push!(edges[i], j)
                push!(edges[j], i)
            end
        end
    end
    ## BFS
    qs, qt = Queue{Int}(), Queue{Int}()
    dists = [0 for i in 1:len]
    visited = [0 for i in 1:len]
    enqueue!(qs, s)
    enqueue!(qt, t)
    visited[s] = 1
    visited[t] = 2
    dists[s] = dists[t] = 1
    while !isempty(qs) && !isempty(qt)
        rt1, rt2 = dequeue!(qs), dequeue!(qt)
        dis1, dis2 = dists[rt1] + 1, dists[rt2] + 1
        ## println(rt1, " ", dis1)
        for neib in edges[rt1]
            if visited[neib] == 2
                return dis1 + dists[neib] - 1
            elseif visited[neib] == 0
                visited[neib] = 1
                enqueue!(qs, neib)
                dists[neib] = dis1
            end
        end
        for neib in edges[rt2]
            if visited[neib] == 1
                return dis2 + dists[neib] - 1
            elseif visited[neib] == 0
                visited[neib] = 2
                enqueue!(qt, neib)
                dists[neib] = dis2
            end
        end
    end
    return 0
end
## @lc code=end


text:
export
    Dynamics,
    step,
    substep,
    build_pid_policy

#=
This file defines the dynamics type. It is used to
step the state of the system forward in time.
=#

#=
Dynamics type contains constant values used for progressing the
state of the system forward in time.
=#
type Dynamics
    dt::Int64
    dt_step::Float64
    pid_omega::Float64
    pid_b_std::Float64
    pid_b_std_coc::Float64
    use_noise::Bool
    rng::MersenneTwister
    function Dynamics(dt::Int64, dt_step::Float64, pid_omega::Float64,
                pid_b_std::Float64, pid_b_std_coc::Float64, use_noise::Bool,
                rng::MersenneTwister)
        return new(dt, dt_step, pid_omega, pid_b_std, pid_b_std_coc, use_noise, rng)
    end
end

#=
Description:
This function steps the state of the system forward using the actions of the uavs

Parameters:
- dynamics: Dynamics type containing necessary fields
- state: State type containing all variable information necessary for progressing
        the system except for the agent actions
- actions: The actions of the uavs in the system

Return Values:
- next_state: the next state of the system after stepping forward in time
- min_horizontal_dist: minimum horizontal distance experienced by the uavs
=#
function Base.step(dynamics::Dynamics, state::State, actions::Vector{Float64})
    @assert state.num_uavs == length(actions) "Each uav must have one action. \
        state: $state\tactions: $actions"

    # generate pid controllers for the different uavs using their actions
    policies = Array(Function, state.num_uavs)
    for idx in 1:state.num_uavs
        action = actions[idx]

        # correct COC to 0 degrees nominal turn
        coc_flag = false
        if action == COC_ACTION
            coc_flag = true
            action = 0.0
        end

        # add noise to bank angle commanded
        if dynamics.use_noise
            if coc_flag
                bank_angle = action + dynamics.pid_b_std_coc * randn(dynamics.rng)
            else
                bank_angle = action + dynamics.pid_b_std * randn(dynamics.rng)
            end
        else
            bank_angle = action
        end

        policies[idx] = build_pid_policy(dynamics, bank_angle)
    end

    # step forward each uav independently
    # collect two things:
    # first, the final result of the dynamics applied to each uav
    # second, all the intermediate states also
    # use the intermediate states to determine the min_horizontal_dist
    next_uav_states = Array(UAVState, state.num_uavs)
    intermediate_uav_states = Array(Array{UAVState}, state.num_uavs)

    for idx in 1:state.num_uavs
        next_states = step(state.uav_states[idx], policies[idx], dynamics.dt, dynamics.dt_step)
        next_uav_states[idx] = next_states[end]
        intermediate_uav_states[idx] = next_states
    end

    # get the min_horizontal_dist to return
    min_horizontal_dist = get_min_horizontal_distance(intermediate_uav_states)

    # create and return a new State object using the list of UAVStates
    return State(next_uav_states), min_horizontal_dist
end

#=
Description:
This function steps forward the state of a single uav using a
policy (i.e., pid controller)

Parameters:
- uav_state: a UAVState object for the uav to be stepped forward
- policy: the pid controller for this uav
- dt: time steps to step forward
- dt_step: number of dynamics updates per timestep

Return Value:
- next_states: UAVState objects containing the intermediate states of the uav
=#
function Base.step(uav_state::UAVState, policy::Function, dt::Int64, dt_step::Float64)
    # convert to vector form, call substep, and return
    uav_state_vector = convert(Vector{Float64}, uav_state)
    next_states = substep(uav_state_vector, policy, dt, dt_step)
    return next_states
end

#=
Description:
Steps forward the state by the smallest unit of time

Parameters:
- uav_state_vector: the state of the uav as a vector instead of an object
- policy: the pid controller for this uav
- dt: time steps to step forward
- dt_step: number of dynamics updates per timestep

Return Values:
- next states: the next states in UAVState format
=#
function substep(uav_state_vector::Vector{Float64}, policy::Function, dt::Int64, dt_step::Float64)
    # use the actual speed of the uav rather than the speed constant
    speed = get_speed(uav_state_vector)
    G = 9.8

    # dynamics equation to solve
    function dxdt(t::Float64, state::Vector{Float64})
        dx = zeros(length(state))
        dx[1] = speed * cos(state[5])                   # xdot
        dx[3] = speed * sin(state[5])                   # ydot
        dx[5] = G * tan(state[7]) / speed               # psidot
        dx[7] = state[8]                                # phidot
        dx[8] = policy(state)                           # phiddot
        dx[2] = -speed * sin(state[5]) * dx[7]          # xddot
        dx[4] = speed * cos(state[5]) * dx[7]           # yddot
        dx[6] = G / speed * sec(state[8])^2 * dx[8]     # psiddot
        return dx
    end

    # create a list of the times to step through
    times = 0:dt_step:dt

    # step through them solving for next state each timestep
    # initial state is the incoming state vector
    state = uav_state_vector

    # collect the states as row vectors in an array
    states = Array(UAVState, length(times) - 1)
    for tidx = 1:length(times) - 1
        tspan = [times[tidx], times[tidx + 1]]
        _, xout = ODE.ode45(dxdt, state, tspan)

        # collect states to return
        states[tidx] = convert(UAVState, xout[end])

        # set the state for use in the next loop
        state = xout[end]
    end

    # return all intermediate states for min horizontal computation
    return states
end

#=
Description:
Builds a pid controller given an action and is an effective policy for the uav

Parameters:
- dynamics: the dynamics of the environment in this case provide constant values
- action: the bank angle the uav has decided to take (in radians)

Return Value:
- pid_policy: a function that given a state returns an action
=#
function build_pid_policy(dynamics::Dynamics, action::Float64)
    function pid_policy(uav_state_vector::Vector{Float64})
        policy =  2 * dynamics.pid_omega * -uav_state_vector[8] +
                    dynamics.pid_omega^2 * (action - uav_state_vector[7])
        return policy
    end
    return pid_policy
end


text:
# ---
# title: 395. Longest Substring with At Least K Repeating Characters
# id: problem395
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Divide and Conquer, Recursion, Sliding Window
# link: <https://leetcode.com/problems/longest-substring-with-at-least-k-repeating-characters/description/>
# hidden: true
# ---
#
# Given a string `s` and an integer `k`, return _the length of the longest
# substring of_ `s` _such that the frequency of each character in this substring
# is greater than or equal to_ `k`.
#
#
#
# **Example 1:**
#
#
#
#     Input: s = "aaabb", k = 3
#     Output: 3
#     Explanation: The longest substring is "aaa", as 'a' is repeated 3 times.
#
#
# **Example 2:**
#
#
#
#     Input: s = "ababbc", k = 2
#     Output: 5
#     Explanation: The longest substring is "ababb", as 'a' is repeated 2 times and 'b' is repeated 3 times.
#
#
#
#
# **Constraints:**
#
#   * `1 <= s.length <= 104`
#   * `s` consists of only lowercase English letters.
#   * `1 <= k <= 105`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end