text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing

#    @assert precompile(HTTP.request, (String, String,))
#    @assert precompile(HTTP.request, (String, URI, Headers, Vector{UInt8},))
#=
    @assert precompile(HTTP.URIs.parseurlchar, (UInt8, Char, Bool,))
    @assert precompile(HTTP.status, (HTTP.Response,))
    @assert precompile(HTTP.Cookies.pathmatch, (HTTP.Cookies.Cookie, String,))
    @assert precompile(HTTP.onheaderfield, (HTTP.Parser, Array{UInt8, 1}, Int64, Int64,))
    @assert precompile(HTTP.isjson, (Array{UInt8, 1}, UInt64, Int64,))
    @assert precompile(HTTP.onheadervalue, (HTTP.Parser, HTTP.Response, Array{UInt8, 1}, Int64, Int64, Bool, String,))
    @assert precompile(HTTP.isjson, (Array{UInt8, 1}, Int64, Int64,))
    @assert precompile(HTTP.onurlbytes, (HTTP.Parser, Array{UInt8, 1}, Int64, Int64,))
    @assert precompile(HTTP.onurl, (HTTP.Parser,))
    @assert precompile(HTTP.Response, (Int64, String,))
    @assert precompile(HTTP.URIs.getindex, (Array{UInt8, 1}, HTTP.URIs.Offset,))
    @assert precompile(HTTP.iscompressed, (Array{UInt8, 1},))
    @assert precompile(HTTP.Cookies.readcookies, (Base.Dict{String, String}, String,))
    @assert precompile(HTTP.tocameldash!, (String,))
    @assert precompile(HTTP.canonicalizeheaders, (Dict{String,String},))
    @assert precompile(HTTP.URIs.http_parse_host_char, (HTTP.URIs.http_host_state, Char,))
    @assert precompile(HTTP.Form, (Base.Dict{String, Any},))
    @assert precompile(HTTP.Cookies.hasdotsuffix, (String, String,))
    @assert precompile(HTTP.onheadervalue, (HTTP.Parser, Array{UInt8, 1}, Int64, Int64,))
    @assert precompile(HTTP.Cookies.parsecookievalue, (String, Bool,))
    @assert precompile(HTTP.processresponse!, (HTTP.Client, HTTP.Connection{Base.TCPSocket}, HTTP.Response, String, HTTP.Method, Task, Bool, Float64, Bool))
    @assert precompile(HTTP.Request, (HTTP.Method, HTTP.URIs.URI, Base.Dict{String, String}, HTTP.FIFOBuffers.FIFOBuffer,))
    @assert precompile(HTTP.Cookies.isCookieDomainName, (String,))
    @assert precompile(HTTP.getbytes, (Base.TCPSocket, Float64))
#    @assert precompile(HTTP.FIFOBuffers.write, (HTTP.FIFOBuffers.FIFOBuffer, Array{UInt8, 1}, Int64, Int64,))
    @assert precompile(HTTP.URIs.port, (HTTP.URIs.URI,))
    @assert precompile(HTTP.read, (HTTP.Form,))
    @assert precompile(HTTP.get, (HTTP.Nitrogen.ServerOptions, Symbol, Int64,))
    @assert precompile(HTTP.ignorewhitespace, (String, Int64, Int64,))
    @assert precompile(HTTP.Response, ())
    @assert precompile(HTTP.Cookies.string, (String, Array{HTTP.Cookies.Cookie, 1}, Bool,))
    @assert precompile(HTTP.FIFOBuffers.read, (HTTP.FIFOBuffers.FIFOBuffer, Int64,))
    @assert precompile(HTTP.processresponse!, (HTTP.Client, HTTP.Connection{MbedTLS.SSLContext}, HTTP.Response, String, HTTP.Method, Task, Bool, Float64, Bool))
    @assert precompile(HTTP.Form, (Base.Dict{String, String},))
    #@assert precompile(HTTP.FIFOBuffers.String, (HTTP.FIFOBuffers.FIFOBuffer,))
    @assert precompile(HTTP.update!, (HTTP.RequestOptions, HTTP.RequestOptions,))
    @assert precompile(HTTP.restofstring, (String, Int64, Int64,))
    @assert precompile(HTTP.ismatch, (Type{HTTP.MP4Sig}, Array{UInt8, 1}, Int64,))
    @assert precompile(HTTP.connectandsend, (HTTP.Client, Type{HTTP.http}, String, String, HTTP.Request, HTTP.RequestOptions, Bool,))
    @assert precompile(HTTP.URIs.unescape, (String,))
    @assert precompile(HTTP.Request, ())
    @assert precompile(HTTP.isjson, (String, Int64, Int64,))
    @assert precompile(HTTP.URIs.splitpath, (String,))
    @assert precompile(HTTP.ignorewhitespace, (Array{UInt8, 1}, Int64, Int64,))
    @assert precompile(HTTP.Cookies.String, (HTTP.Cookies.Cookie, Bool,))
    @assert precompile(HTTP.Cookies.shouldsend, (HTTP.Cookies.Cookie, Bool, String, String,))
    @assert precompile(HTTP.Client, (Nothing, HTTP.RequestOptions,))
    @assert precompile(HTTP.Cookies.parsecookievalue, (Base.SubString{String}, Bool,))
    @assert precompile(HTTP.ignorewhitespace, (Array{UInt8, 1}, UInt64, Int64,))
    @assert precompile(HTTP.Cookies.hash, (HTTP.Cookies.Cookie, UInt,))
    @assert precompile(HTTP.ismatch, (HTTP.HTMLSig, Array{UInt8, 1}, Int64,))
    @assert precompile(HTTP.redirect, (HTTP.Response, HTTP.Client, HTTP.Request, HTTP.RequestOptions, Bool, Array{HTTP.Response, 1}, Int64, Bool,))
    @assert precompile(HTTP.get, (HTTP.RequestOptions, Symbol, Int64,))
    @assert precompile(HTTP.Cookies.isequal, (HTTP.Cookies.Cookie, HTTP.Cookies.Cookie,))
    @assert precompile(HTTP.sniff, (String,))
    @assert precompile(HTTP.restofstring, (Array{UInt8, 1}, UInt64, Int64,))
    @assert precompile(HTTP.stalebytes!, (Base.TCPSocket,))
#    @assert precompile(HTTP.FIFOBuffers.write, (HTTP.FIFOBuffers.FIFOBuffer, UInt8,))
    @assert precompile(HTTP.eof, (HTTP.Form,))
    @assert precompile(HTTP.Cookies.readsetcookie, (String, String,))
    @assert precompile(HTTP.Response, (Int64, HTTP.Request,))
    @assert precompile(HTTP.URIs.http_parser_parse_url, (Array{UInt8, 1}, Int64, Int64, Bool,))
    @assert precompile(HTTP.Response, (String,))
#    @assert precompile(HTTP.FIFOBuffers.read, (HTTP.FIFOBuffers.FIFOBuffer, Type{Tuple{UInt8, Bool}},))
    @assert precompile(HTTP.Response, (Int64, Base.Dict{String, String}, String,))
    @assert precompile(HTTP.ismatch, (HTTP.Masked, Array{UInt8, 1}, Int64,))
    @assert precompile(HTTP.Cookies.sanitizeCookieValue, (String,))
    @assert precompile(HTTP.URIs.escape, (Int64, Int64,))
    @assert precompile(HTTP.URIs.http_parse_host, (Array{UInt8, 1}, HTTP.URIs.Offset, Bool,))
    @assert precompile(HTTP.read, (HTTP.Form, Int64,))
    @assert precompile(HTTP.nb_available, (HTTP.Multipart{Base.IOStream},))
    @assert precompile(HTTP.URIs.host, (HTTP.URIs.URI,))
    @assert precompile(HTTP.Request, (HTTP.Method, HTTP.URIs.URI, Base.Dict{String, String}, HTTP.Form,))
    @assert precompile(HTTP.getbytes, (MbedTLS.SSLContext, Float64))
    @assert precompile(HTTP.ismatch, (HTTP.Exact, Array{UInt8, 1}, Int64,))
    @assert precompile(HTTP.getconnections, (Type{HTTP.https}, HTTP.Client, String,))
    @assert precompile(HTTP.ismatch, (Type{HTTP.TextSig}, Array{UInt8, 1}, Int64,))
    @assert precompile(HTTP.restofstring, (Array{UInt8, 1}, Int64, Int64,))
    @assert precompile(HTTP.dead!, (HTTP.Connection{Base.TCPSocket},))
    @assert precompile(HTTP.addcookies!, (HTTP.Client, String, HTTP.Request, Bool,))
#    @assert precompile(HTTP.FIFOBuffers.length, (HTTP.FIFOBuffers.FIFOBuffer,))
    @assert precompile(HTTP.Cookies.domainandtype, (String, String,))
    @assert precompile(HTTP.sniff, (Array{UInt8, 1},))
    @assert precompile(HTTP.mark, (HTTP.Multipart{Base.IOStream},))
    @assert precompile(HTTP.seek, (HTTP.Form, Int64,))
    @assert precompile(HTTP.onheadervalue, (HTTP.Parser, HTTP.Request, Array{UInt8, 1}, Int64, Int64, Bool, String))
    @assert precompile(HTTP.FIFOBuffers.write, (HTTP.FIFOBuffers.FIFOBuffer, String,))
    @assert precompile(HTTP.URIs.escape, (String, String,))
    @assert precompile(HTTP.dead!, (HTTP.Connection{MbedTLS.SSLContext},))
    @assert precompile(HTTP.URIs.isvalid, (HTTP.URIs.URI,))
#    @assert precompile(HTTP.FIFOBuffers.read, (HTTP.FIFOBuffers.FIFOBuffer, Type{UInt8},))
    @assert precompile(HTTP.getconnections, (Type{HTTP.http}, HTTP.Client, String,))
    @assert precompile(HTTP.Cookies.validCookieDomain, (String,))
    @assert precompile(HTTP.headers, (HTTP.Response,))
    @assert precompile(HTTP.request, (HTTP.Response,))
    @assert precompile(HTTP.URIs.escape, (String, Array{String, 1},))
    @assert precompile(HTTP.connectandsend, (HTTP.Client, Type{HTTP.https}, String, String, HTTP.Request, HTTP.RequestOptions, Bool,))
    @assert precompile(HTTP.body, (HTTP.Response,))
    @assert precompile(HTTP.http_should_keep_alive, (HTTP.Parser, HTTP.Request,))
    @assert precompile(HTTP.length, (HTTP.Form,))
#    @assert precompile(HTTP.FIFOBuffers.position, (HTTP.FIFOBuffers.FIFOBuffer,))
    @assert precompile(HTTP.URIs.escape, (String,))
    @assert precompile(HTTP.busy!, (HTTP.Connection{Base.TCPSocket},))
    @assert precompile(HTTP.connect, (HTTP.Client, HTTP.http, String, String, HTTP.RequestOptions, Bool,))
    @assert precompile(HTTP.string, (HTTP.Request,))
    @assert precompile(HTTP.parse!, (HTTP.Request, HTTP.Parser, Array{UInt8, 1},))
    @assert precompile(HTTP.parse!, (HTTP.Request, HTTP.Parser, Array{UInt8, 1}, Int64, String, HTTP.Method, Task))
    @assert precompile(HTTP.parse!, (HTTP.Response, HTTP.Parser, Array{UInt8, 1}, Int64, String, HTTP.Method, Task))
    @assert precompile(HTTP.onbody, (HTTP.Request, Task, Array{UInt8, 1}, Int64, Int64,))
    @assert precompile(HTTP.take!, (HTTP.Response,))
    @assert precompile(HTTP.Response, (String,))
    @assert precompile(HTTP.Cookies.Cookie, (Base.SubString{String}, Base.SubString{String},))
    @assert precompile(HTTP.initTLS!, (Type{HTTP.http}, String, HTTP.RequestOptions, Base.TCPSocket,))
    @assert precompile(HTTP.cookies, (HTTP.Response,))
    @assert precompile(HTTP.history, (HTTP.Response,))
    @assert precompile(HTTP.setconnection!, (Type{HTTP.https}, HTTP.Client, String, HTTP.Connection{MbedTLS.SSLContext},))
    @assert precompile(HTTP.setconnection!, (Type{HTTP.http}, HTTP.Client, String, HTTP.Connection{Base.TCPSocket},))
    @assert precompile(HTTP.eof, (HTTP.Multipart{Base.IOStream},))
    @assert precompile(HTTP.position, (HTTP.Form,))
    @assert precompile(HTTP.haskey, (Type{HTTP.http}, HTTP.Client, String,))
    @assert precompile(HTTP.parse!, (HTTP.Response, HTTP.Parser, Array{UInt8, 1},))
    @assert precompile(HTTP.reset, (HTTP.Multipart{Base.IOStream},))
#    @assert precompile(HTTP.FIFOBuffers.seek, (HTTP.FIFOBuffers.FIFOBuffer, Tuple{Int64, Int64, Int64},))
#    @assert precompile(HTTP.FIFOBuffers.wait, (HTTP.FIFOBuffers.FIFOBuffer,))
    @assert precompile(HTTP.initTLS!, (Type{HTTP.https}, String, HTTP.RequestOptions, Base.TCPSocket,))
    @assert precompile(HTTP.stalebytes!, (MbedTLS.SSLContext,))
    @assert precompile(HTTP.sniff, (Base.IOStream,))
    @assert precompile(HTTP.request, (HTTP.Client, HTTP.Request, HTTP.RequestOptions, Bool, Array{HTTP.Response, 1}, Int, Bool,))
    @assert precompile(HTTP.read, (HTTP.Multipart{Base.IOStream}, Int64,))
    @assert precompile(HTTP.isjson, (Array{UInt8, 1},))
#    @assert precompile(HTTP.FIFOBuffers.close, (HTTP.FIFOBuffers.FIFOBuffer,))
    @assert precompile(HTTP.headers, (HTTP.Request,))
    @assert precompile(HTTP.busy!, (HTTP.Connection{MbedTLS.SSLContext},))
    @assert precompile(HTTP.seek, (HTTP.Form, Tuple{Int64, Int64, Int64},))
#    @assert precompile(HTTP.FIFOBuffers.eof, (HTTP.FIFOBuffers.FIFOBuffer,))
    @assert precompile(HTTP.contenttype, (HTTP.Masked,))
    @assert precompile(HTTP.get, (HTTP.RequestOptions, Symbol, MbedTLS.SSLConfig,))
    @assert precompile(HTTP.contenttype, (HTTP.Exact,))
#    @assert precompile(HTTP.FIFOBuffers.FIFOBuffer, (HTTP.FIFOBuffers.FIFOBuffer,))
    @assert precompile(HTTP.haskey, (Type{HTTP.https}, HTTP.Client, String,))
#    @assert precompile(HTTP.FIFOBuffers.readavailable, (HTTP.FIFOBuffers.FIFOBuffer,))
    @assert precompile(HTTP.string, (HTTP.Response, HTTP.Nitrogen.ServerOptions,))
    @assert precompile(HTTP.string, (HTTP.Response, HTTP.RequestOptions,))
    @assert precompile(HTTP.string, (HTTP.Request, HTTP.RequestOptions,))
    @assert precompile(HTTP.Request, (String,))
    @assert precompile(HTTP.ismatch, (Type{HTTP.JSONSig}, Array{UInt8, 1}, Int64,))
    @assert precompile(HTTP.hasmessagebody, (HTTP.Request,))
#    @assert precompile(HTTP.FIFOBuffers.write, (HTTP.FIFOBuffers.FIFOBuffer, Array{UInt8, 1},))
    @assert precompile(HTTP.readavailable, (HTTP.Form,))
    @assert precompile(HTTP.get, (String,))
    @assert precompile(HTTP.URL, (String,))
    @assert precompile(HTTP.request, (HTTP.Client, HTTP.Method, HTTP.URI,))
    @assert precompile(HTTP.RequestOptions, ())
#    @assert precompile(HTTP.Request, (HTTP.Method, HTTP.URI, Dict{String, String}, HTTP.FIFOBuffer))
    @assert precompile(HTTP.request, (HTTP.Request,))
    @assert precompile(HTTP.request, (HTTP.Client, HTTP.Request))
    @assert precompile(HTTP.request, (HTTP.Client, HTTP.Request, HTTP.RequestOptions, Bool, Vector{HTTP.Response}, Int, Bool))
    @static if VERSION < v"0.7-DEV"
        @assert precompile(HTTP.Client, (Base.AbstractIOBuffer{Array{UInt8, 1}}, HTTP.RequestOptions,))
        @assert precompile(HTTP.URIs.printuri, (Base.AbstractIOBuffer{Array{UInt8, 1}}, String, String, String, String, String, String, String,))
#        @assert precompile(HTTP.FIFOBuffers.FIFOBuffer, (Base.AbstractIOBuffer{Array{UInt8, 1}},))
        @assert precompile(HTTP.startline, (Base.AbstractIOBuffer{Array{UInt8, 1}}, HTTP.Response,))
        @assert precompile(HTTP.writemultipartheader, (Base.AbstractIOBuffer{Array{UInt8, 1}}, HTTP.Multipart{Base.IOStream},))
        @assert precompile(HTTP.print, (Base.AbstractIOBuffer{Array{UInt8, 1}}, HTTP.Method,))
        @assert precompile(HTTP.sniff, (Base.AbstractIOBuffer{Array{UInt8, 1}},))
        @assert precompile(HTTP.headers, (Base.AbstractIOBuffer{Array{UInt8, 1}}, HTTP.Response,))
        @assert precompile(HTTP.headers, (Base.AbstractIOBuffer{Array{UInt8, 1}}, HTTP.Request,))
        @assert precompile(HTTP.body, (Base.AbstractIOBuffer{Array{UInt8, 1}}, HTTP.Response, HTTP.Nitrogen.ServerOptions,))
        @assert precompile(HTTP.body, (Base.AbstractIOBuffer{Array{UInt8, 1}}, HTTP.Response, HTTP.RequestOptions,))
        @assert precompile(HTTP.URIs.print, (Base.AbstractIOBuffer{Array{UInt8, 1}}, HTTP.URIs.ParsingStateCode,))
        @assert precompile(HTTP.URIs.print, (Base.AbstractIOBuffer{Array{UInt8, 1}}, HTTP.URIs.http_host_state,))
        @assert precompile(HTTP.showcompact, (Base.AbstractIOBuffer{Array{UInt8, 1}}, HTTP.Request,))
        @assert precompile(HTTP.URIs.show, (Base.AbstractIOBuffer{Array{UInt8, 1}}, HTTP.URIs.URI,))
        @assert precompile(HTTP.writemultipartheader, (Base.AbstractIOBuffer{Array{UInt8, 1}}, Base.IOStream,))
        @assert precompile(HTTP.URIs.print, (Base.AbstractIOBuffer{Array{UInt8, 1}}, HTTP.URIs.URI,))
        @assert precompile(HTTP.body, (Base.AbstractIOBuffer{Array{UInt8, 1}}, HTTP.Request, HTTP.RequestOptions,))
        @assert precompile(HTTP.startline, (Base.AbstractIOBuffer{Array{UInt8, 1}}, HTTP.Request,))
    else
        @assert precompile(HTTP.Client, (Base.GenericIOBuffer{Array{UInt8, 1}}, HTTP.RequestOptions,))
        @assert precompile(HTTP.URIs.printuri, (Base.GenericIOBuffer{Array{UInt8, 1}}, String, String, String, String, String, String, String,))
#        @assert precompile(HTTP.FIFOBuffers.FIFOBuffer, (Base.GenericIOBuffer{Array{UInt8, 1}},))
        @assert precompile(HTTP.startline, (Base.GenericIOBuffer{Array{UInt8, 1}}, HTTP.Response,))
        @assert precompile(HTTP.writemultipartheader, (Base.GenericIOBuffer{Array{UInt8, 1}}, HTTP.Multipart{Base.IOStream},))
        @assert precompile(HTTP.print, (Base.GenericIOBuffer{Array{UInt8, 1}}, HTTP.Method,))
        @assert precompile(HTTP.sniff, (Base.GenericIOBuffer{Array{UInt8, 1}},))
        @assert precompile(HTTP.headers, (Base.GenericIOBuffer{Array{UInt8, 1}}, HTTP.Response,))
        @assert precompile(HTTP.headers, (Base.GenericIOBuffer{Array{UInt8, 1}}, HTTP.Request,))
        @assert precompile(HTTP.body, (Base.GenericIOBuffer{Array{UInt8, 1}}, HTTP.Response, HTTP.Nitrogen.ServerOptions,))
        @assert precompile(HTTP.body, (Base.GenericIOBuffer{Array{UInt8, 1}}, HTTP.Response, HTTP.RequestOptions,))
        @assert precompile(HTTP.URIs.print, (Base.GenericIOBuffer{Array{UInt8, 1}}, HTTP.URIs.ParsingStateCode,))
        @assert precompile(HTTP.URIs.print, (Base.GenericIOBuffer{Array{UInt8, 1}}, HTTP.URIs.http_host_state,))
        @assert precompile(HTTP.showcompact, (Base.GenericIOBuffer{Array{UInt8, 1}}, HTTP.Request,))
        @assert precompile(HTTP.URIs.show, (Base.GenericIOBuffer{Array{UInt8, 1}}, HTTP.URIs.URI,))
        @assert precompile(HTTP.writemultipartheader, (Base.GenericIOBuffer{Array{UInt8, 1}}, Base.IOStream,))
        @assert precompile(HTTP.URIs.print, (Base.GenericIOBuffer{Array{UInt8, 1}}, HTTP.URIs.URI,))
        @assert precompile(HTTP.body, (Base.GenericIOBuffer{Array{UInt8, 1}}, HTTP.Request, HTTP.RequestOptions,))
        @assert precompile(HTTP.startline, (Base.GenericIOBuffer{Array{UInt8, 1}}, HTTP.Request,))
    end
=#
end
_precompile_()


text:
mutable struct MemoryBuffer{L,T} <: DenseVector{T}
    data::NTuple{L,T}
    @inline function MemoryBuffer{L,T}(::UndefInitializer) where {L,T}
        @assert isbitstype(T) "Memory buffers must point to bits types, but `isbitstype($T) == false`."
        new{L,T}()
    end
end

struct BCache{T<:Union{UInt,Nothing}}
    p::Ptr{Cvoid}
    i::T
end

ext: jl


text:
# Demonstrate/test functions in cornerfeatures

using ImageProjectiveGeometry
using PyPlot
using PerceptualColourMaps

img = PyPlot.imread("lena.tif");
figure(10), PyPlot.imshow(img);
PyPlot.set_cmap(PyPlot.ColorMap("gray"))
PyPlot.title("Input image")

println("Harris corner detector")
(cimg, r, c) = ImageProjectiveGeometry.harris(img, 1,k= 0.04, radius=2, N=100, subpixel=false, img=img)
cimg[cimg .< 0] = 0;   # eliminate -ve values
figure(2), PyPlot.imshow(cimg.^0.25);
PyPlot.set_cmap(PyPlot.ColorMap("gray"))
ImageProjectiveGeometry.keypause()

println("Noble corner detector")
(cimg, r, c) = ImageProjectiveGeometry.noble(img, 1, radius=2, N=100, subpixel=false, img=img);
figure(2), PyPlot.imshow(cimg.^0.25);
PyPlot.set_cmap(PyPlot.ColorMap("gray"))
ImageProjectiveGeometry.keypause()

println("Shi-Tomasi corner detector")
(cimg, r, c) = ImageProjectiveGeometry.shi_tomasi(img, 1, radius=2, N=100, subpixel=false, img=img);
figure(2), PyPlot.imshow(cimg.^0.25);
PyPlot.set_cmap(PyPlot.ColorMap("gray"))
ImageProjectiveGeometry.keypause()

println("-ve Hessian features")
hdet = ImageProjectiveGeometry.hessianfeatures(img, 1)
(r, c) = ImageProjectiveGeometry.nonmaxsuppts(-hdet, N=100, img=img)
figure(2), PyPlot.imshow(-hdet);
PyPlot.set_cmap(PyPlot.ColorMap("gray"))
ImageProjectiveGeometry.keypause()

println("+ve Hessian features")
hdet = ImageProjectiveGeometry.hessianfeatures(img, 1)
(r, c) = ImageProjectiveGeometry.nonmaxsuppts(+hdet, N=100, img=img)
figure(2), PyPlot.imshow(hdet);
PyPlot.set_cmap(PyPlot.ColorMap("gray"))
ImageProjectiveGeometry.keypause()

println("Coherence")
cimg = ImageProjectiveGeometry.coherence(img, 2)
(r, c) = ImageProjectiveGeometry.nonmaxsuppts(cimg, N=10, img=img)
figure(2), PyPlot.imshow(cimg);
PyPlot.set_cmap(PyPlot.ColorMap("gray"))
ImageProjectiveGeometry.keypause()

println("Fast radial")
(S, So) = ImageProjectiveGeometry.fastradial(img, [1, 3, 5, 7])
(r, c) = ImageProjectiveGeometry.nonmaxsuppts(-S, N=10, img=img)
figure(2), PyPlot.imshow(-S);
PyPlot.set_cmap(PyPlot.ColorMap("gray"))

figure(3), PyPlot.imshow(-So);
PyPlot.set_cmap(PyPlot.ColorMap("gray"))


text:
@with_kw mutable struct DLevelSetData{T <: Real}
    Q::T = 1e10
    Q̃::T = 1e10
    θ::T = -1e10
    timestamp::Int = 1
    iterations::Int = 0
end

@with_kw mutable struct DLevelSetParameters{T <: Real}
    κ::T = 0.3
    τ::T = 1e-6
    λ::T = 0.5
    log::Bool = true
end

struct DLevelSet{T <: Real, A <: AbstractVector, M <: LQSolver, S <: LQSolver} <: AbstractLShapedSolver{T,A,M,S}
    structuredmodel::JuMP.Model
    solverdata::DLevelSetData{T}

    # Master
    mastersolver::M
    projectionsolver::M
    mastervector::A
    c::A
    x::A
    Q_history::A

    # Subproblems
    nscenarios::Int
    subobjectives::Vector{A}
    finished::Vector{Int}

    # Workers
    subworkers::Vector{SubWorker{T,A,S}}
    work::Vector{Work}
    decisions::Decisions{A}
    cutqueue::CutQueue{T}

    # Regularizer
    ξ::A
    Q̃_history::A
    levels::A

    # Cuts
    θs::A
    cuts::Vector{SparseHyperPlane{T}}
    θ_history::A

    # Params
    parameters::DLevelSetParameters{T}
    progress::ProgressThresh{T}

    @implement_trait DLevelSet HasLevels
    @implement_trait DLevelSet IsParallel

    function (::Type{DLevelSet})(model::JuMP.Model,ξ₀::AbstractVector,mastersolver::AbstractMathProgSolver,subsolver::AbstractMathProgSolver; kw...)
        if nworkers() == 1
            warn("There are no worker processes, defaulting to serial version of algorithm")
            return LevelSet(model,ξ₀,mastersolver,subsolver; kw...)
        end
        length(ξ₀) != model.numCols && error("Incorrect length of starting guess, has ",length(ξ₀)," should be ",model.numCols)
        !haskey(model.ext,:SP) && error("The provided model is not structured")

        T = promote_type(eltype(ξ₀),Float32)
        c_ = convert(AbstractVector{T},JuMP.prepAffObjective(model))
        c_ *= model.objSense == :Min ? 1 : -1
        mastervector = convert(AbstractVector{T},copy(ξ₀))
        x₀_ = convert(AbstractVector{T},copy(ξ₀))
        ξ₀_ = convert(AbstractVector{T},copy(ξ₀))
        A = typeof(x₀_)

        msolver = LQSolver(model,mastersolver)
        psolver = LQSolver(model,mastersolver)
        M = typeof(msolver)
        S = LQSolver{typeof(LinearQuadraticModel(subsolver)),typeof(subsolver)}
        n = StochasticPrograms.nscenarios(model)

        lshaped = new{T,A,M,S}(model,
                               DLevelSetData{T}(),
                               msolver,
                               psolver,
                               mastervector,
                               c_,
                               x₀_,
                               A(),
                               n,
                               Vector{A}(),
                               Vector{Int}(),
                               Vector{SubWorker{T,A,S}}(nworkers()),
                               Vector{Work}(nworkers()),
                               RemoteChannel(() -> DecisionChannel(Dict{Int,A}())),
                               RemoteChannel(() -> Channel{QCut{T}}(4*nworkers()*n)),
                               ξ₀_,
                               A(),
                               A(),
                               A(fill(-Inf,n)),
                               Vector{SparseHyperPlane{T}}(),
                               A(),
                               DLevelSetParameters{T}(;kw...),
                               ProgressThresh(1.0, "Distributed Leveled L-Shaped Gap "))
        lshaped.progress.thresh = lshaped.parameters.τ
        push!(lshaped.subobjectives,zeros(n))
        push!(lshaped.finished,0)
        push!(lshaped.Q_history,Inf)
        push!(lshaped.Q̃_history,Inf)
        push!(lshaped.θ_history,-Inf)

        init!(lshaped,subsolver)

        return lshaped
    end
end
DLevelSet(model::JuMP.Model,mastersolver::AbstractMathProgSolver,subsolver::AbstractMathProgSolver; kw...) = DLevelSet(model,rand(model.numCols),mastersolver,subsolver; kw...)

function (lshaped::DLevelSet)()
    # Reset timer
    lshaped.progress.tfirst = lshaped.progress.tlast = time()
    # Start workers
    active_workers = init_workers!(lshaped)
    # Start procedure
    while true
        status = iterate!(lshaped)
        if status != :Valid
            close_workers!(lshaped,active_workers)
            return status
        end
    end
end


text:
#!/bin/bash
# -*- mode: julia -*-
#=
exec ${JULIA:-julia} --color=yes --startup-file=no \
    -e "include(popfirst!(ARGS))" "${BASH_SOURCE[0]}" \
    "$@"
=#

csv_path, = ARGS

using SnoopCompile

SnoopCompile.@snoopc ["-i"] csv_path begin
    using IPython
    IPython.test_ipython_jl(inprocess=true)
    IPython.start_ipython()
end


text:
# TODO:
# - [ ] add documentation.

default_monitors() = [[:p, :beta, :gammaC, :gammaT, :etaC, :etaT, :mu, :nu,
                       :omega, :psi, :tau, :sigma, :phi]]

complete_monitors() = [default_monitors()[1],
                      [:lambdaC, :lambdaT, :vC, :vT, :zetaC, :zetaT]]

default_flags() = [UpdateBetaWithSkewT(), UpdateLambdaWithSkewT()]

function default_callback_fn(state::State, data::Data, prior::Prior,
                             iter::Int, pbar::MCMC.ProgressBars.ProgressBar)
  ll = round(loglike_latent_var(state, data), digits=3)
  MCMC.ProgressBars.set_postfix(pbar, loglike=ll)
  return Dict(:loglike => ll)
end

"""
- `init::State`: Initial state.
- `data::Data`: Data.
- `prior::Prior`: Priors.
- `tuners::Tuners`: Tuners.
- `nsamps::Vector{Int}`: Number of posterior samples for each monitor. Needs to
   be in descending order. Defaults to `[1000]`.
- `nburn::Int`: Number of burn-in iterations. Defaults to `1000`.
- `thin::Int`: Thinning factor for first monitor. Defaults to `1`.
- `monitors::Vector{Vector{Symbol}}`: Parameters to keep samples of. Defaults
   to `default_monitors()`.
- `fix::Vector{Symbol}`: Parameters to hold constant. Defaults to `Symbol[]`.
- `flags::Vector{Flag}`: Available options include `:update_beta_with_skewt`
   and `:update_lambda_with_skewt`.
   Defaults to `default_flags()`.
- `verbose::Int`: How much info to print. Defaults to 1.
- `callback_fn`: A function with signature `(state::State, data::Data,
   prior::Prior, iter::Int, pbar::ProgressBar)` that returns a dicionary of
   summary statistics. Enables one to add metrics to progress bar.
"""
function fit(init::State, data::Data, prior::Prior, tuners::Tuners;
             nsamps::Vector{Int}=[1000], nburn::Int=1000, thin::Int=1,
             monitors::Vector{Vector{Symbol}}=default_monitors(),
             callback_fn::Function=default_callback_fn,
             fix::Vector{Symbol}=Symbol[],
             rep_beta_flipped::Int=0,
             rep_aux::Int=1, seed=nothing,
             flags::Vector{Flag}=Flag[], verbose::Int=1)

  # Set random seed.
  seed == nothing || Random.seed!(seed)

  # `fix` and `flags` should be unique.
  fix = unique(fix)
  flags = unique(flags)

  # Print settings for sanity check.
  if verbose > 0
    println("seed: ", seed)
    println("flags: ", flags)
    println("rep_beta_flipped: ", rep_beta_flipped)
    println("rep_aux: ", rep_aux)
    println("fix: ", fix)
    println("monitors: ", monitors)
  end

  isfixed(sym::Symbol) = sym in fix

  function update!(state)
    prev_beta = state.beta
    update_state!(state, data, prior, tuners, fix=fix, flags=flags)
    if state.beta != prev_beta
      for _ in 1:rep_beta_flipped
        update_state!(state, data, prior, tuners, flags=flags,
                      fix=[fix; [:p, :beta]])
      end
    else
      for _ in 1:rep_aux
        update_state!(state, data, prior, tuners, flags=flags,
                      fix=[fix; [:p, :beta, :lambda, :gamma, :eta, :mu]])
      end
    end
  end

  function _callback_fn(state::State, iter::Int,
                       pbar::MCMC.ProgressBars.ProgressBar)
    return callback_fn(state, data, prior, iter, pbar)
  end

  chain, laststate, summarystats = MCMC.gibbs(
      init, update!, nsamps=nsamps, nburn=nburn, thin=thin, monitors=monitors,
      callback_fn=_callback_fn)

  return chain, laststate, summarystats
end

max_stars_repo_name: luiarthur/CytofDensityEstimation

max_issues_repo_name: luiarthur/CytofDensityEstimation

max_forks_repo_name: luiarthur/CytofDensityEstimation


text:
@testset "Newton polygon" begin

    @polyvar x y
    f = 4 + 4x + 4y + 4x*y + 2x^2 + 2y^2

    newtonpoly = newtonpolygon(f)
    @test lattices(newtonpoly) == [[2, 0], [1, 1], [0, 2], [1, 0], [0, 1], [0, 0]]
    @test vertices(newtonpoly) == [[0, 2], [0, 0], [2, 0]]
    @test vertexindices(newtonpoly) == [3, 6, 1]
    @test facets(newtonpoly) == [[[0, 2], [0, 0]],[[0, 0], [2, 0]], [[2, 0], [0, 2]]]
    @test newtonpoly.subdivision == [[2, 5, 6, 4], [3, 5, 2], [2, 4, 1]]

    @test_throws AssertionError newtonpolygon(1 + x)

    @test_throws AssertionError vertices(newtonpolygon(x + y))

    @polyvar z w
    f = 1 + z * cis(0.8π) + z^3*w*cis(1.6π) + z*w^2*cis(0.4π) + w*cis(1.2π)
    n = newtonpolygon(archimedean_tropical_polynomial(f))
    @test length(n.subdivision) == 1
    @test length(n.subdivision[1]) == 5

    @polyvar x y
    f = x^2*y + y^2 + 3x^2*y^3 + y^4 + x^4*y^4
    newt = newtonpolygon(f)
    @test vertices(newt) == [[0, 4], [0, 2], [2, 1], [4, 4]]
    @test PolynomialAmoebas.on_boundary(newt, [2, 3]) == false
    @test PolynomialAmoebas.on_boundary(newt, [1, 1.5])
    @test PolynomialAmoebas.iscontained(newt, [2, 3])
end


text:
module SpatialDependence

    """
    SpatialDependence
    A Julia package for spatial dependence (spatial autocorrelation), spatial weights matrices, and exploratory spatial data analysis (ESDA).
    [SpatialDependence repository](https://github.com/javierbarbero/SpatialDependence.jl).
    """

    using GeoInterface: AbstractGeometry, AbstractPoint, AbstractPolygon, AbstractMultiPolygon, Point, Polygon, MultiPolygon, coordinates, geotype, shapecoords
    using NearestNeighbors: KDTree, knn, inrange
    using PlotUtils: palette
    using Random: shuffle, AbstractRNG, default_rng
    using RecipesBase
    using Tables: istable

    import Base: length
    import SparseArrays: SparseMatrixCSC, sparse
    import StatsBase: counts, levels, maximum, median, minimum, nobs, percentile, sample, score, standardize, weights, zscore, ZScoreTransform
    import Statistics: mean, median, std, quantile

    export
        # Types
        SpatialWeights,

        # Spatial Weights Creation functions
        dnearneigh,
        knearneigh,
        polyneigh,

        # Spatial Weights functions
        cardinalities,
        neighbors,
        weights,
        nislands,
        islands,
        wtransform,
        wtransform!,
        wtransformation,
        slag,
        sparse,
        reggeomlattice,
        meancenter,
        centroid,

        # Spatial autocorrelation
        AbstractGlobalSpatialAutocorrelation,
        GlobalMoran,
        moran,
        GlobalGeary,
        geary,
        score,
        scoreperms,
        pvalue,
        AbstractLocalSpatialAutocorrelation,
        LocalMoran,
        localmoran,
        LocalGeary,
        localgeary,
        GetisOrd,
        getisord,
        issignificant,

        # Re-export from Statistics
        nobs,
        mean,
        median,
        std,
        score,
        zscore,

        # Choropleth Maps and Classification<
        AbstractMapClassificator,
        AbstractGraduatedMapClassificator,
        EqualIntervals,
        Quantiles,
        CustomBreaks,
        GraduatedMapClassification,

        AbstractStatisticalMapClassificator,
        NaturalBreaks,
        BoxPlot,
        StdMean,
        Percentiles,

        AbstractUniqueMapClassificator,
        Unique,
        UniqueMapClassification,

        AbstractCoLocationMapClassificator,
        CoLocation,
        CoLocationMapClassification,

        mapclassify,
        maplabels,

        assignments,
        counts,
        bounds,
        levels

    include("sweights/sweights.jl")
    include("sweights/dnearneigh.jl")
    include("sweights/knearneigh.jl")
    include("sweights/polyneigh.jl")
    include("sweights/lattice.jl")
    include("sweights/centroid.jl")

    include("slag.jl")

    include("scor/global.jl")
    include("scor/moran.jl")
    include("scor/geary.jl")
    include("scor/local.jl")
    include("scor/localmoran.jl")
    include("scor/localgeary.jl")
    include("scor/getisord.jl")

    include("maps/classification.jl")
    include("maps/graduated.jl")
    include("maps/statistical.jl")
    include("maps/unique.jl")
    include("maps/colocation.jl")

    include("plotrecipes.jl")
    include("plotchoropleth.jl")

end


text:
doc"""
    TDist(ν)

The *Students T distribution* with `ν` degrees of freedom has probability density function

$f(x; d) = \frac{1}{\sqrt{d} B(1/2, d/2)}
\left( 1 + \frac{x^2}{d} \right)^{-\frac{d + 1}{2}}$

```julia
TDist(d)      # t-distribution with d degrees of freedom

params(d)     # Get the parameters, i.e. (d,)
dof(d)        # Get the degrees of freedom, i.e. d
```

External links

[Student's T distribution on Wikipedia](https://en.wikipedia.org/wiki/Student%27s_t-distribution)

"""
immutable TDist <: ContinuousUnivariateDistribution
    ν::Float64

    TDist(ν::Real) = (@check_args(TDist, ν > zero(ν)); new(ν))
end

@distr_support TDist -Inf Inf

#### Parameters

dof(d::TDist) = d.ν
params(d::TDist) = (d.ν,)

#### Statistics

mean(d::TDist) = d.ν > 1.0 ? 0.0 : NaN
median(d::TDist) = 0.0
mode(d::TDist) = 0.0

function var(d::TDist)
    ν = d.ν
    ν > 2.0 ? ν / (ν - 2.0) :
    ν > 1.0 ? Inf : NaN
end

skewness(d::TDist) = d.ν > 3.0 ? 0.0 : NaN

function kurtosis(d::TDist)
    ν = d.ν
    ν > 4.0 ? 6.0 / (ν - 4.0) :
    ν > 2.0 ? Inf : NaN
end

function entropy(d::TDist)
    h = 0.5 * d.ν
    h1 = h + 0.5
    h1 * (digamma(h1) - digamma(h)) + 0.5 * log(d.ν) + lbeta(h, 0.5)
end

#### Evaluation & Sampling

@_delegate_statsfuns TDist tdist ν

rand(d::TDist) = StatsFuns.Rmath.tdistrand(d.ν)

function cf(d::TDist, t::Real)
    t == 0 && return complex(1.0)
    h = d.ν * 0.5
    q = d.ν * 0.25
    t2 = t*t
    complex(2*(q*t2)^q*besselk(h,sqrt(d.ν)*abs(t))/gamma(h))
end

gradlogpdf(d::TDist, x::Float64) = -((d.ν + 1.0) * x) / (x^2 + d.ν)


text:
# replay a game step by step
export GameReplayer

function GameReplayer()
    board = initboard()
    i = 1
    return function()
        display(h[i])
        display(policy(h[i]))
        i+=1
    end
end


text:
# include("6a_comparison.jl")

include("_common.jl")
include("_calib_funcs.jl")

using Streamfall: load_calibration

using StatsPlots, DataStructures
using Infiltrator
using HypothesisTests

ensemble_figs = joinpath(FIG_PATH, "ensemble_results", "cmd") * "/"
mkpath(ensemble_figs)

ensemble_res_path = joinpath(DATA_PATH, "ensemble_results") * "/"
mkpath(ensemble_res_path)

baseline_calib = joinpath(DATA_PATH, "baselines") * "/"
calib_param_path = joinpath(DATA_PATH, "calib_params") * "/"

function read_params(fn)
    calibrated_params = open(fn, "r") do f
        p_ins = Float64[]
        for p in eachline(f)
            p_ins = collect(map(x -> parse(Float64, x), split(p, ",")))
        end

        return p_ins
    end

    return calibrated_params
end

function separate_state_periods(activity)
    low = findall(activity .== 1)
    mid = findall(activity .== 2)
    high = findall(activity .== 3)

    return low, mid, high
end

function separate_state_periods(activity, t_min, t_max=nothing)
    if !isnothing(t_max)
        active = activity[t_min:t_max]
    else
        active = activity[t_min:end]
    end

    return separate_state_periods(active)
end

function build_ensemble_flow(Qb, Qm, state_factors, activity)
    Qe = deepcopy(Qb)

    # Weight towards whichever has the better score in calibration period
    low_factor, mid_factor, high_factor = state_factors

    # Create ensemble outflow for entire simulation period
    low_store, mid_store, high_store = separate_state_periods(activity)

    Qe[high_store] = (Qm[high_store] .* high_factor) .+ (Qb[high_store] .* (1.0 - high_factor))
    Qe[mid_store] = (Qm[mid_store] .* mid_factor) .+ (Qb[mid_store] .* (1.0 - mid_factor))
    Qe[low_store] = (Qm[low_store] .* low_factor) .+ (Qb[low_store] .* (1.0 - low_factor))

    return Qe
end

function combine_ensemble(x, Qb, Qm, active_params)
    if all(isapprox.(x, 1.0)) | all(isapprox.(x, 0.0; atol=1e-6))
        return -9999.9, -9999.9, -9999.9
    end

    # Create ensemble outflow for entire simulation period
    Qe = build_ensemble_flow(Qb, Qm, x, active_params)

    # Subset to calibration period only
    calib_low_store, calib_mid_store, calib_high_store = separate_state_periods(active_params, calib_ts+1, CALIB_LN)
    Qe_calib = Qe[calib_ts+1:CALIB_LN]
    # overall_mKGE = Streamfall.mKGE(Qo_calib, Qe_calib)
    # overall_A2k = KSampleADTest(Qo_calib, Qe_calib).A²k
    low_score = Streamfall.mKGE(Qo_calib[calib_low_store], Qe_calib[calib_low_store])
    mid_score = Streamfall.mKGE(Qo_calib[calib_mid_store], Qe_calib[calib_mid_store])
    high_score = Streamfall.mKGE(Qo_calib[calib_high_store], Qe_calib[calib_high_store])

    # low_A2k = KSampleADTest(Qo_calib[calib_low_store], Qe_calib[calib_low_store]).A²k
    # mid_A2k = KSampleADTest(Qo_calib[calib_mid_store], Qe_calib[calib_mid_store]).A²k
    # high_A2k = KSampleADTest(Qo_calib[calib_high_store], Qe_calib[calib_high_store]).A²k
    # mean_A2k = mean([low_A2k, mid_A2k, high_A2k])

    # ensemble_valid = Qe[CALIB_LN+1:end]
    # valid_low_store, valid_mid_store, valid_high_store = separate_state_periods(active_params, CALIB_LN)
    # high_score = Streamfall.mKGE(Qo_valid[valid_high_store], ensemble_valid[valid_high_store])
    # if length(valid_high_store) == 0
    #     return (-9999.9, -9999.9, -9999.9, -9999.9)
    # end

    return (low_score, mid_score, high_score)
end

"""Create base-base ensemble instance"""
function base_base_score(x; model1="split_mean_NmKGE", model2="NmKGE")
    local active_params = state_active_params[model1]
    local Qb = base_Q[model1]
    local Qm = base_Q[model2]

    return combine_ensemble(repeat(x, 3), Qb, Qm, active_params)
end

"""Create base-state ensemble instance"""
function base_state_score(x; model1="split_mean_NmKGE", model2="NmKGE")
    local active_params = state_active_params[model2]
    local Qb = base_Q[model1]
    local Qm = state_Qe[model2]

    return combine_ensemble(x, Qb, Qm, active_params)
end

"""Create state-state ensemble instance"""
function state_state_score(x; model1="split_mean_NmKGE", model2="NmKGE")
    local active_params = state_active_params[model1]
    local Qb = state_Qe[model1]
    local Qm = state_Qe[model2]

    return combine_ensemble(x, Qb, Qm, active_params)
end

function insert_NmKGE_scores(df, approach, Qo, Qm)
    local Qm_calib = Qm[calib_ts+1:CALIB_LN]
    local Qm_valid = Qm[CALIB_LN+1:end]

    if length(valid_low) == 0
        low_score = missing
    else
        low_score = Streamfall.NmKGE(Qo_valid[valid_low], Qm_valid[valid_low])
    end

    if length(valid_high) == 0
        high_score = missing
    else
        high_score = Streamfall.NmKGE(Qo_valid[valid_high], Qm_valid[valid_high])
    end

    push!(df, [approach,
          Streamfall.NmKGE(Qo, Qm),

          Streamfall.NmKGE(Qo_calib, Qm_calib),
          Streamfall.NmKGE(Qo_calib[calib_low], Qm_calib[calib_low]),
          Streamfall.NmKGE(Qo_calib[calib_mid], Qm_calib[calib_mid]),
          Streamfall.NmKGE(Qo_calib[calib_high], Qm_calib[calib_high]),

          Streamfall.NmKGE(Qo_valid, Qm_valid),
          low_score,
          Streamfall.NmKGE(Qo_valid[valid_mid], Qm_valid[valid_mid]),
          high_score,
    ])
end

function insert_A2k_scores(df, approach, Qo, Qm)
    local Qm_calib = Qm[calib_ts+1:CALIB_LN]
    local Qm_valid = Qm[CALIB_LN+1:end]

    AD_low = KSampleADTest(Qo_calib[calib_low], Qm_calib[calib_low]).A²k
    AD_mid = KSampleADTest(Qo_calib[calib_mid], Qm_calib[calib_mid]).A²k
    AD_high = KSampleADTest(Qo_calib[calib_high], Qm_calib[calib_high]).A²k

    if length(valid_low) == 0
        low_score = missing
    else
        low_score = KSampleADTest(Qo_valid[valid_low], Qm_valid[valid_low]).A²k
    end

    if length(valid_high) == 0
        high_score = missing
    else
        high_score = KSampleADTest(Qo_valid[valid_high], Qm_valid[valid_high]).A²k
    end

    push!(df, [approach,
          KSampleADTest(Qo, Qm).A²k,  # Overall A2k score
          mean([AD_low, AD_mid, AD_high]),  # Mean of states in calibration period

          KSampleADTest(Qo_calib, Qm_calib).A²k,
          AD_low,
          AD_mid,
          AD_high,

          KSampleADTest(Qo_valid, Qm_valid).A²k,
          low_score,
          KSampleADTest(Qo_valid[valid_mid], Qm_valid[valid_mid]).A²k,
          high_score,
    ])
end

calib_ts = 3651
state = :storage
state_str = String(state)

target_idx = [1,2,4,5,6,7,8]
quantiles = [0.0, 0.1, 0.9]
# quantiles = [0.0, 0.16, 0.84]

Qo = FULL_DATASET[:, "410730_Q"]
Qo_calib = CALIB_OBS[calib_ts+1:CALIB_LN]
Qo_valid = VALID_OBS

# Recreate results if they do not exist
if (!@isdefined baseline_NmKGE_results) | !isfile(joinpath(ensemble_res_path, "overall_baseline_NmKGE_results.csv"))
    # Record of ensemble flow
    base_Q = Dict()
    state_Qe = Dict()

    # catchment state indicated by each state instance
    state_active_params = Dict()

    # Records of CMD
    base_cmd = Dict()
    state_cmd = Dict()

    for app in APPROACHES
        # baseline
        base_sn, base_id = setup_network("baselines/cotter_baseline_IHACRES_$(app).yml")
        base_node = base_sn[base_id]
        Streamfall.run_basin!(base_sn, FULL_CLIMATE)

        base_Q[app] = base_node.outflow
        base_cmd[app] = base_node.storage[2:end]

        # State-based
        # Read in parameters
        calibrated_state_params = nothing
        try
            fn = "$(calib_param_path)log_$(state_str)_$(app)_10_year.txt"
            calibrated_state_params = read_params(fn)
        catch err
            if err isa SystemError
                @info "Skipping $(app) as file not available"
                continue
            end
        end

        sn, n_id = setup_network("baselines/cotter_baseline_IHACRES_$app.yml")
        node = sn[n_id]

        _, x0, __ = param_info(node; with_level=false)
        mod_params = update_partial(x0, target_idx, calibrated_state_params)

        active_params, param_idxs = online_burn_state_node!(sn, n_id, FULL_CLIMATE, state, mod_params, quantiles; burn_in=calib_ts, log_transform=true, releases=nothing)

        state_Qe[app] = node.outflow
        state_cmd[app] = node.storage
        state_active_params[app] = active_params
    end

    base_base_Q = Dict()
    base_state_Qe = Dict()
    state_state_Qe = Dict()

    base_base_mix = Dict()
    base_state_mix = Dict()
    state_state_mix = Dict()

    # Optimize ensembles and build comparison DFs
    for app in APPROACHES
        for state_approach in APPROACHES
            # Compare baseline-state_based combinations
            local active_params = state_active_params[app]
            local Qb = base_Q[app]
            local Qm = state_Qe[state_approach]

            ensemble_app = (app, state_approach)
            @info "Mixing " ensemble_app

            opt_func = (x) -> base_state_score(x; model1=app, model2=state_approach)
            opt = bbsetup(opt_func; SearchRange=[(0.0, 1.0), (0.0, 1.0), (0.0, 1.0)],
                        MaxTime=900, TraceInterval=30,
                        MaxStepsWithoutProgress=50_000,
                        Method=:borg_moea,
                        ϵ=0.05,
                        FitnessScheme=ParetoFitnessScheme{3}(is_minimizing=false))

            res = bboptimize(opt)
            bc_res = best_candidate(res)
            base_state_Qe[ensemble_app] = build_ensemble_flow(Qb, Qm, bc_res, active_params)
            base_state_mix[ensemble_app] = bc_res

            # Compare pairs built with state-based models only
            # Only compare unique combinations
            if app == state_approach
                continue
            end

            # Create base-base ensemble
            local Qb = base_Q[app]
            local Qm = base_Q[state_approach]
            opt_func = (x) -> base_base_score(x; model1=app, model2=state_approach)
            opt = bbsetup(opt_func; SearchRange=[(0.0, 1.0)],
                        MaxTime=900, TraceInterval=30,
                        MaxStepsWithoutProgress=50_000,
                        Method=:borg_moea,
                        ϵ=0.05,
                        FitnessScheme=ParetoFitnessScheme{3}(is_minimizing=false))
            res = bboptimize(opt)
            mixer = best_candidate(res)[1]
            bc_res = [mixer, mixer, mixer]
            base_base_Q[ensemble_app] = build_ensemble_flow(Qb, Qm, bc_res, active_params)
            base_base_mix[ensemble_app] = bc_res

            # Build state-based pairs
            local Qb = state_Qe[app]
            local Qm = state_Qe[state_approach]
            opt_func = (x) -> state_state_score(x; model1=app, model2=state_approach)
            opt = bbsetup(opt_func; SearchRange=[(0.0, 1.0), (0.0, 1.0), (0.0, 1.0)],
                        MaxTime=900, TraceInterval=30,
                        MaxStepsWithoutProgress=50_000,
                        Method=:borg_moea,
                        ϵ=0.05,
                        FitnessScheme=ParetoFitnessScheme{3}(is_minimizing=false))

            res = bboptimize(opt)
            bc_res = best_candidate(res)
            state_state_Qe[ensemble_app] = build_ensemble_flow(Qb, Qm, bc_res, active_params)
            state_state_mix[ensemble_app] = bc_res
        end
    end
end

if !isfile(joinpath(ensemble_res_path, "overall_baseline_NmKGE_results.csv"))
    # Collate results using a single identified state series (e.g., by state-based NmKGE)
    baseline_NmKGE_results = DataFrame(Approach=[], Overall_NmKGE=[],
                                       Overall_calib=[], Wet_calib=[], Mid_calib=[], Dry_calib=[],
                                       Overall_valid=[], Wet_valid=[], Mid_valid=[], Dry_valid=[])
    base_base_NmKGE_results = deepcopy(baseline_NmKGE_results)
    base_state_NmKGE_results = deepcopy(baseline_NmKGE_results)
    state_NmKGE_results = deepcopy(baseline_NmKGE_results)
    state_state_NmKGE_results = deepcopy(baseline_NmKGE_results)

    baseline_A2k_results = DataFrame(Approach=[], Overall_A2k=[], Mean_A2k=[],
                                     Overall_calib=[], Wet_calib=[], Mid_calib=[], Dry_calib=[],
                                     Overall_valid=[], Wet_valid=[], Mid_valid=[], Dry_valid=[])
    base_base_A2k_results = deepcopy(baseline_A2k_results)
    base_state_A2k_results = deepcopy(baseline_A2k_results)
    state_A2k_results = deepcopy(baseline_A2k_results)
    state_state_A2k_results = deepcopy(baseline_A2k_results)

    # active_params = state_active_params["NmKGE"]
    # calib_low, calib_mid, calib_high = separate_state_periods(active_params, calib_ts+1, CALIB_LN)
    # valid_low, valid_mid, valid_high = separate_state_periods(active_params, CALIB_LN+1)
    common_cmd = log.(base_cmd["NmKGE"])
    common_cmd_calib = common_cmd[1:calib_ts]
    common_thresholds = quantile(common_cmd_calib, quantiles[2:end])
    calib_low = findall(common_cmd .< common_thresholds[1])
    calib_mid = findall(common_thresholds[1] .<= common_cmd .< common_thresholds[2])
    calib_high = findall(common_thresholds[2] .<= common_cmd)

    active_params = copy(common_cmd)
    active_params[calib_low] .= 1
    active_params[calib_mid] .= 2
    active_params[calib_high] .= 3

    calib_low, calib_mid, calib_high = separate_state_periods(active_params, calib_ts+1, CALIB_LN)
    valid_low, valid_mid, valid_high = separate_state_periods(active_params, CALIB_LN+1)

    for app in APPROACHES
        # Assess baseline instances
        Qs = base_Q[app]
        insert_NmKGE_scores(baseline_NmKGE_results, app, Qo, Qs)
        insert_A2k_scores(baseline_A2k_results, app, Qo, Qs)

        # State-based instances
        Qs = state_Qe[app]
        insert_NmKGE_scores(state_NmKGE_results, app, Qo, Qs)
        insert_A2k_scores(state_A2k_results, app, Qo, Qs)

        for s_app in APPROACHES
            approach = (app, s_app)

            # base-state instances
            Qs = base_state_Qe[approach]
            insert_NmKGE_scores(base_state_NmKGE_results, approach, Qo, Qs)
            insert_A2k_scores(base_state_A2k_results, approach, Qo, Qs)

            if app == s_app
                # skip non-unique combinations
                continue
            end

            # base-base instances
            Qs = base_base_Q[approach]
            insert_NmKGE_scores(base_base_NmKGE_results, approach, Qo, Qs)
            insert_A2k_scores(base_base_A2k_results, approach, Qo, Qs)

            # State-State ensembles
            Qs = state_state_Qe[approach]
            insert_NmKGE_scores(state_state_NmKGE_results, approach, Qo, Qs)
            insert_A2k_scores(state_state_A2k_results, approach, Qo, Qs)
        end
    end

    # Export to file
    CSV.write(joinpath(ensemble_res_path, "overall_baseline_NmKGE_results.csv"), baseline_NmKGE_results)
    CSV.write(joinpath(ensemble_res_path, "overall_state_NmKGE_results.csv"), state_NmKGE_results)
    CSV.write(joinpath(ensemble_res_path, "overall_base-base_NmKGE_results.csv"), base_base_NmKGE_results)
    CSV.write(joinpath(ensemble_res_path, "overall_base-state_NmKGE_results.csv"), base_state_NmKGE_results)
    CSV.write(joinpath(ensemble_res_path, "overall_state-state_NmKGE_results.csv"), state_state_NmKGE_results)

    CSV.write(joinpath(ensemble_res_path, "overall_baseline_A2k_results.csv"), baseline_A2k_results)
    CSV.write(joinpath(ensemble_res_path, "overall_state_A2k_results.csv"), state_A2k_results)
    CSV.write(joinpath(ensemble_res_path, "overall_base-base_A2k_results.csv"), base_base_A2k_results)
    CSV.write(joinpath(ensemble_res_path, "overall_base-state_A2k_results.csv"), base_state_A2k_results)
    CSV.write(joinpath(ensemble_res_path, "overall_state-state_A2k_results.csv"), state_state_A2k_results)

    # Save outflows
    CSV.write(joinpath(ensemble_res_path, "baseline_outflows.csv"), DataFrame(base_Q))
    CSV.write(joinpath(ensemble_res_path, "state_outflows.csv"), DataFrame(state_Qe))
    CSV.write(joinpath(ensemble_res_path, "base_base_outflows.csv"),
                DataFrame([k=>v for (k,v) in zip(["$(b) - $(m)"
                for (b, m) in keys(base_base_Q)], values(base_base_Q))]))
    CSV.write(joinpath(ensemble_res_path, "base_state_outflows.csv"),
                DataFrame([k=>v for (k,v) in zip(["$(b) - $(m)"
                for (b, m) in keys(base_state_Qe)], values(base_state_Qe))]))
    CSV.write(joinpath(ensemble_res_path, "state_state_outflows.csv"),
                DataFrame([k=>v for (k,v) in zip(["$(b) - $(m)"
                for (b, m) in keys(state_state_Qe)], values(state_state_Qe))]))

    # Save instance weightings
    CSV.write(joinpath(ensemble_res_path, "base_base_mix.csv"),
                DataFrame([k=>v for (k,v) in zip(["$(b) - $(m)"
                for (b, m) in keys(base_base_mix)], values(base_base_mix))]))
    CSV.write(joinpath(ensemble_res_path, "base_state_mix.csv"),
                DataFrame([k=>v for (k,v) in zip(["$(b) - $(m)"
                for (b, m) in keys(base_state_mix)], values(base_state_mix))]))
    CSV.write(joinpath(ensemble_res_path, "state_state_mix.csv"),
                DataFrame([k=>v for (k,v) in zip(["$(b) - $(m)"
                for (b, m) in keys(state_state_mix)], values(state_state_mix))]))

    CSV.write(joinpath(ensemble_res_path, "state_active_params.csv"),
                DataFrame(state_active_params))

    CSV.write(joinpath(ensemble_res_path, "base_NmKGE_state.csv"), DataFrame("NmKGE"=>active_params))
end


text:
module DitherPunk

using ImageCore
using ImageCore: NumberLike, Pixel, GenericImage, GenericGrayImage, MappedArrays
using ImageCore.Colors: DifferenceMetric
using Random
using IndirectArrays
using OffsetArrays
using Requires

abstract type AbstractDither end

include("compat.jl")
include("utils.jl")
include("api/binary.jl")
include("api/color.jl")
include("threshold.jl")
include("ordered.jl")
include("ordered_imagemagick.jl")
include("error_diffusion.jl")
include("closest_color.jl")
include("eval.jl")

export dither, dither!
# Threshold dithering
export ConstantThreshold, WhiteNoiseThreshold
# Ordered dithering
export OrderedDither
export Bayer, ClusteredDots, CentralWhitePoint, BalancedCenteredPoint, Rhombus
export IM_checks, IM_h4x4a, IM_h6x6a, IM_h8x8a, IM_h4x4o, IM_h6x6o, IM_h8x8o
export IM_c5x5, IM_c6x6, IM_c7x7
# Error diffusion
export ErrorDiffusion
export SimpleErrorDiffusion, FloydSteinberg, JarvisJudice, Stucki, Burkes
export Sierra, TwoRowSierra, SierraLite, Atkinson, Fan93, ShiauFan, ShiauFan2
# Closest color
export ClosestColor
# Other utilities
export upscale

# Conditional dependencies using Requires.jl
function __init__()
    @require ColorSchemes = "35d6a980-a343-548e-a6ea-1d62b119f2f4" begin
        include("colorschemes.jl")
    end
    @require Clustering = "aaaa29a8-35af-508c-8bc3-b662a17a0fe5" begin
        include("clustering.jl")
    end
    @require UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228" begin
        include("braille.jl")
        export braille, brailleprint
    end
end

end # module


text:
@testset "Certification" begin
    @testset "Simple 1: Input = Vector{Expression}" begin
        @var x y
        # define the polynomials
        f₁ = (x^4 + y^4 - 1) * (x^2 + y^2 - 2) + x^5 * y
        f₂ = x^2 + 2x * y^2 - 2 * y^2 - 1 / 2
        F = [f₁, f₂]
        result = solve(F)

        # Input: System(F)
        cert = certify(System(F), result)
        @test nnonsingular(result) == 18
        @test ncertified(cert) == 18
        @test ndistinct_certified(cert) == 18
        @test nreal_certified(cert) == 4
        @test ncomplex_certified(cert) == 14
        @test ndistinct_real_certified(cert) == 4
        @test ndistinct_complex_certified(cert) == 14
        save("tmp_cert.txt", cert)
        @test !isempty(read("tmp_cert.txt", String))

        # Input: F
        cert = certify(F, result)
        @test nnonsingular(result) == 18
        @test ncertified(cert) == 18
        @test ndistinct_certified(cert) == 18
        @test nreal_certified(cert) == 4
        @test ndistinct_real_certified(cert) == 4

        # Control Display
        cert = certify(F, result, show_progress = false)

        # Double solutions
        S = solutions(result)
        cert = certify(F, [S; S])
        @test ncertified(cert) == 36
        @test ndistinct_certified(cert) == 18
        @test nreal_certified(cert) == 8
        @test ndistinct_real_certified(cert) == 4
    end

    @testset "Simple 2: Input = Vector{MP.Polynomial}" begin
        @polyvar x y
        # define the polynomials
        f₁ = (x^4 + y^4 - 1) * (x^2 + y^2 - 2) + x^5 * y
        f₂ = x^2 + 2x * y^2 - 2 * y^2 - 1 / 2
        F = [f₁, f₂]
        result = solve(F)

        # Input: System(F)
        cert = certify(System(F), result)
        @test nnonsingular(result) == 18
        @test ncertified(cert) == 18
        @test ndistinct_certified(cert) == 18
        @test nreal_certified(cert) == 4
        @test ndistinct_real_certified(cert) == 4

        # Input: F
        cert = certify(F, result)
        @test nnonsingular(result) == 18
        @test ncertified(cert) == 18
        @test ndistinct_certified(cert) == 18
        @test nreal_certified(cert) == 4
        @test ndistinct_real_certified(cert) == 4

        # Control Display
        cert = certify(F, result, show_progress = false)

        # Double solutions
        S = solutions(result)
        cert = certify(System(F), [S; S])
        @test ncertified(cert) == 36
        @test ndistinct_certified(cert) == 18
        @test nreal_certified(cert) == 8
        @test ndistinct_real_certified(cert) == 4
    end

    @testset "Parameters: Input = Vector{Expression}" begin
        @var x y
        f = (x^4 + y^4 - 1) * (x^2 + y^2 - 2) + x^5 * y
        # define new variables u₁, u₂ and λ₁
        @var λ[1:1] u[1:2]
        # define the jacobian of F
        J = differentiate([f], [x, y])
        # J' defines the transpose of J
        F = [[x, y] - u - J' * λ; f]
        u₀ = [-0.32, -0.1]
        res = solve(F, parameters = u, target_parameters = u₀)

        # Input: System(F)
        C = System(F, parameters = u)
        cert = certify(C, res; target_parameters = u₀)
        @test ncertified(cert) == 36
        @test ndistinct_certified(cert) == 36
        @test nreal_certified(cert) == 8
        @test ndistinct_real_certified(cert) == 8

        # Input: F
        cert = certify(F, res, u₀; parameters = u)
        @test ncertified(cert) == 36
        @test ndistinct_certified(cert) == 36
        @test nreal_certified(cert) == 8
        @test ndistinct_real_certified(cert) == 8

        # Invalid solutions
        cert = certify(C, [100 .* randn(ComplexF64, 3) for i = 1:10], u₀)
        @test ncertified(cert) < 10
    end

    @testset "Parameters: Input = Vector{MP.Polynomial}" begin
        @polyvar x y
        f = (x^4 + y^4 - 1) * (x^2 + y^2 - 2) + x^5 * y
        # define new variables u₁, u₂ and λ₁
        @polyvar λ[1:1] u[1:2]
        # define the jacobian of F
        J = differentiate([f], [x, y])
        # J' defines the transpose of J
        F = [[x, y] - u - J' * λ; f]
        u₀ = [-0.32, -0.1]
        res = solve(F, parameters = u, target_parameters = u₀)

        # Input: System(F)
        C = System(F, parameters = u)
        cert = certify(C, res; target_parameters = u₀)
        @test ncertified(cert) == 36
        @test ndistinct_certified(cert) == 36
        @test nreal_certified(cert) == 8
        @test ndistinct_real_certified(cert) == 8

        # Input: F
        cert = certify(F, res; parameters = u, target_parameters = u₀)
        @test ncertified(cert) == 36
        @test ndistinct_certified(cert) == 36
        @test nreal_certified(cert) == 8
        @test ndistinct_real_certified(cert) == 8
    end

    @testset "positive" begin
        @var x y
        f = System([x^2 + y^2 - 1, x - y])
        res = solve(f; compile = false, start_system = :total_degree)
        cert = certify(f, res, compile = false)
        @test count(is_positive, certificates(cert)) == 1
        @test count(s -> is_positive(s, 1), certificates(cert)) == 1
        @test count(is_real, certificates(cert)) == 2
    end

    @testset "3264" begin
        F = steiner()
        real_conics = [
            10124547 // 662488724,
            8554609 // 755781377,
            5860508 // 2798943247,
            -251402893 // 1016797750,
            -25443962 // 277938473,
            1 // 1,
            520811 // 1788018449,
            2183697 // 542440933,
            9030222 // 652429049,
            -12680955 // 370629407,
            -24872323 // 105706890,
            1 // 1,
            6537193 // 241535591,
            -7424602 // 363844915,
            6264373 // 1630169777,
            13097677 // 39806827,
            -29825861 // 240478169,
            1 // 1,
            13173269 // 2284890206,
            4510030 // 483147459,
            2224435 // 588965799,
            33318719 // 219393000,
            92891037 // 755709662,
            1 // 1,
            8275097 // 452566634,
            -19174153 // 408565940,
            5184916 // 172253855,
            -23713234 // 87670601,
            28246737 // 81404569,
            1 // 1,
        ]
        real_sols = read_solutions(joinpath(@__DIR__, "data/3264_real_sols.txt"))
        cert = certify(F, real_sols, real_conics; compile = true)
        @test ndistinct_real_certified(cert) == 3264
        cert = certify(F, real_sols, real_conics; compile = false)
        @test ndistinct_real_certified(cert) == 3264
        @test_throws ArgumentError certify(F, real_sols; compile = false)
        @test_throws ArgumentError certify(F, real_sols; compile = true)
    end

    @testset "certify uses complex inversion" begin
        @var x y
        F = System([y / x + x^2 - 3], parameters = [y])
        monres = monodromy_solve(F)
        cert = certify(F, monres)
        @test ndistinct_certified(cert) >= 1
    end
end


text:
#
#   SVD decomposition of a matrix truncated to a rank
#
using Manopt, Manifolds
import LinearAlgebra: norm, svd, Diagonal

"""
    truncated_svd(A, p)

return a singular value decomposition of a real valued matrix A truncated to
rank p.

# Input
* `A` – a real-valued matrix A of size mxn
* `p` – an integer p in min { m, n }

# Output
* `U` – an orthonormal matrix of size mxp
* `V` – an orthonormal matrix of size nxp
* `S` – a diagonal matrix of size pxp with nonnegative and decreasing diagonal
        entries
"""
function truncated_svd(A::Array{Float64,2} = randn(42, 60), p::Int64 = 5)
    (m, n) = size(A)

    if p > min(m,n)
        throw( ErrorException("The Rank p=$p must be smaller than the smallest dimension of A = $min(m, n).") )
    end

    M = ProductManifold(Grassmann(m,p), Grassmann(n,p))

    function cost(X::ProductRepr)
        return cost([submanifold_components(X)...])
    end
    function cost(X::Array{Matrix{Float64},1})
        return -0.5 * norm(transpose(X[1]) * A * X[2])^2
    end

    function egrad(X::Array{Matrix{Float64},1})
        U = X[1]
        V = X[2]
        AV = A*V
        AtU = transpose(A)*U
        return [ -AV*(transpose(AV)*U), -AtU*(transpose(AtU)*V) ];
    end

    function rgrad(M::ProductManifold, X::ProductRepr)
        eG = egrad([submanifold_components(M,X)...])
        x = [submanifold_components(M,X)...]
        return Manifolds.ProductRepr(project.(M.manifolds, x, eG)...)
    end

    function e2rHess(M::Grassmann, x, ξ, eGrad::Matrix{T},Hess::Matrix{T}) where T<:Union{U, Complex{U}} where U<:AbstractFloat
        pxHess = project(M,x,Hess)
        xtGrad = x'*eGrad
        ξxtGrad = ξ*xtGrad
        return pxHess - ξxtGrad
    end

    function eHess(X::Array{Matrix{Float64},1}, H::Array{Matrix{Float64},1})
        U = X[1]
        V = X[2]
        Udot = H[1]
        Vdot = H[2]
        AV = A*V
        AtU = transpose(A)*U
        AVdot = A*Vdot
        AtUdot = transpose(A)*Udot
        return [ -(AVdot*transpose(AV)*U + AV*transpose(AVdot)*U + AV*transpose(AV)*Udot),
                 -(AtUdot*transpose(AtU)*V + AtU*transpose(AtUdot)*V + AtU*transpose(AtU)*Vdot)
            ]
    end

    function rhess(M::ProductManifold, X::ProductRepr, H::ProductRepr)
        x = [submanifold_components(M,X)...]
        h = [submanifold_components(M,H)...]
        eG = egrad(x)
        eH = eHess(x,h)
        return Manifolds.ProductRepr(e2rHess.(M.manifolds, x, h, eG, eH)...)
    end

    x = random_point(M)
    X = trust_regions(M, cost, rgrad, x, rhess;
        Δ_bar=4*sqrt(2*p),
        debug = [:Iteration, " ", :Cost, " | ", DebugEntry(:Δ), "\n", 1, :Stop]
    )

    U = get_component(M,X,1)
    V = get_component(M,X,2)

    Spp = transpose(U)*A*V
    SVD = svd(Spp)
    U = U*SVD.U
    S = SVD.S
    V = V*SVD.V

    return [U, S, V]
end

A=[1. 2. 3. 4.; 5. 6. 7. 8.; 9. 10. 11. 12.; 13. 14. 15. 16.]

truncated_svd(A,2)


text:
using PooksoftOptionsKit
using PyPlot
using Dates

# setup components -
callOptionContract = PSCallOptionContract("AAPL", Date(2020,6,25), 100.0, 1.20, 1; sense=:sell)
equityObject = PSEquityAsset("AAPL", 80.0, 100, Date(2020,6,4))
putOptionContract = PSPutOptionContract("AAPL", Date(2020,6,25), 60.0, 5.00, 1; sense=:buy)

# create set -
assetSet = Set{PSAbstractAsset}()
#push!(assetSet, callOptionContract)
push!(assetSet, equityObject)
push!(assetSet, putOptionContract)

# calculate the PL -
result = compute_complex_trade_profit_and_loss_at_expiration(assetSet, 0.0, 140.0; number_of_price_steps=1000)

# plot -
profit_array = result.value
plot(profit_array[:,1], profit_array[:,2], "r--")     # buyer - P/L
plot(profit_array[:,1], zeros(1000))


text:
mutable struct Edge
    rev::Int
    from::Int
    to::Int
    cap::Int
end

function push!(G::Vector{Vector{Edge}}, from, to, cap)
    fromrev = length(G[from]) + 1
    torev = length(G[to]) + 1
    Base.push!(G[from], Edge(torev, from, to, cap))
    Base.push!(G[to], Edge(fromrev, to, from, 0))
end

function run_flow!(G::Vector{Vector{Edge}}, e::Edge, f)
    e.cap -= f
    G[e.to][e.rev].cap += f
    return G
end

function fodfs(G, v, t, f, seens)
    v == t && return f

    seens[v] = true
    for e ∈ G[v]
        if seens[e.to] || e.cap == 0
            # Skip
        else
            flow = fodfs(G, e.to, t, min(f, e.cap), seens)
            if flow == 0
                # Skip
            else
                run_flow!(G, e, flow)
                return flow
            end
        end
    end

    return 0
end

function FordFulkerson(G, s, t)
    INF = 10^9
    res = 0

    while true
        flow = fodfs(G, s, t, INF, falses(length(G)))
        if flow == 0
            return res
        else
            res += flow
        end
    end

    return 0
end

function main()
    N = parse(Int, readline())
    as = parse.(Int, split(readline()))
    INF = 10^9

    offset = 0
    G = [Vector{Edge}() for _ = 1:N+2]
    for (i, a) = enumerate(as)
        if a ≥ 0
            push!(G, N+1, i, 0)
            push!(G, i, N+2, a)
            offset += a
        else
            push!(G, N+1, i, -a)
            push!(G, i, N+2, 0)
        end
    end

    for i = 1:N
        for j = i+1:N
            if j % i == 0
                push!(G, i, j, INF)
            end
        end
    end

    res = FordFulkerson(G, N+1, N+2)
    println(offset - res)
end

main()


text:
using Documenter

# manual/unicode-input.md
download("http://www.unicode.org/Public/9.0.0/ucd/UnicodeData.txt", joinpath(@__DIR__, "UnicodeData.txt"))

makedocs(
    clean = true,
    doctest = false, # TODO: Fix doctest
    sitename="Julia Taiwan 中文文件",
    pages = [
        "首頁" => "index.md",
        "手冊" => [
            "manual/getting-started.md",
            "manual/variables.md",
            "manual/integers-and-floating-point-numbers.md",
            "manual/mathematical-operations.md",
            "manual/complex-and-rational-numbers.md",
            "manual/strings.md",
            "manual/functions.md",
            "manual/control-flow.md",
            "manual/variables-and-scoping.md",
            "manual/types.md",
            "manual/methods.md",
            "manual/constructors.md",
            "manual/conversion-and-promotion.md",
            "manual/interfaces.md",
            "manual/modules.md",
            "manual/documentation.md",
            "manual/metaprogramming.md",
            "manual/arrays.md",
            "manual/missing.md",
            "manual/networking-and-streams.md",
            "manual/parallel-computing.md",
            "manual/running-external-programs.md",
            "manual/calling-c-and-fortran-code.md",
            "manual/handling-operating-system-variation.md",
            "manual/environment-variables.md",
            "manual/embedding.md",
            "manual/code-loading.md",
            "manual/profile.md",
            "manual/stacktraces.md",
            "manual/performance-tips.md",
            "manual/workflow-tips.md",
            "manual/style-guide.md",
            "manual/faq.md",
            "manual/noteworthy-differences.md",
            "manual/unicode-input.md",
        ],
        "Base" => [
        ],
        "標準程式庫" => [
        ],
        "開發者手冊" => [
        ],
    ],
    Documenter.HTML(
        prettyurls = !("local" in ARGS),
        canonical = "https://julialangtw.github.io/JuliaZHTW.jl/latest/",
    ),
    assets = ["assets/favicon.ico"],
)

if "deploy" in ARGS
deploydocs(
    repo = "github.com/JuliaLangTW/JuliaZHTW.jl.git",
)
end


text:
export psth, segment, NoGrp, Trialwise, Binwise

# enumeration for grouping dimention
@enum GrpDim NoGrp Trialwise Binwise

# ============================================================================ #
function psth(ts::Vector{Float64}, evt::AbstractArray{<:Real},
    kbin::AbstractArray{Int,1}, bin_size::Float64, grp_dim::GrpDim = NoGrp)

    pad = kbin[1] < 0 ? -kbin[1] : 0

    # make sure all time stamps are sufficiently > 0, otherwise we end up
    # trying to use 0 as an index...
    if (ts[1] < bin_size) || (evt[1] < bin_size)
        ts = ts .+ bin_size
        evt = evt .+ bin_size
    end

    # convert timestamps to indicies based on bin size
    # pad indicies to make sure the <kbin> index vector can be subtracted from
    # event indicies without producing indicies < 1
    kspk = round.(Int, ts ./ bin_size) .+ pad

    # if ts and evt refer to the same memory address (i.e. we're doing an
    # autocorrelation) then just have kevt and kspk refer to the same array
    if ts === evt
        kevt = kspk
    else
        kevt = round.(Int, evt ./ bin_size) .+ pad
    end

    # maximum time that we'll be dealing with, gives us the length of the
    # "binary spike vector"
    max_time = kevt[end] + kbin[end]
    kspk = kspk[kspk .<= max_time]

    # Int8 is acceptable so long as we are assiging 1 to the <kspk> elements
    # (see warning below)
    bsv = zeros(Int, max_time)

    @inbounds @fastmath for k in eachindex(kspk)
        bsv[kspk[k]] += 1
    end

    data = zeros(length(kbin), length(kevt))
    @inbounds @fastmath for k in eachindex(evt)
        data[:,k] .= bsv[kevt[k] .+ kbin]
    end

    if grp_dim != NoGrp
        # -------------------------------------------------------------------- #
        # NOTE: this is a crappy, repetitive, inelegant solution
        # -------------------------------------------------------------------- #
        tmp = zeros(Int, max_time)
        tmp[kspk] = 1:length(kspk)

        tmp_data = zeros(Int, length(kbin), length(kevt))

        @inbounds @fastmath for k in eachindex(kevt)
            tmp_data[:,k] .= tmp[kevt[k] .+ kbin]
        end

        # collapse across the non-requested dimention
        if grp_dim == Trialwise
            grp_length = length(evt)
            kdim = 2
        else
            grp_length = length(kbin)
            kdim = 1
        end

        grp = Vector{Vector{Int}}(undef, size(tmp_data, kdim))
        @inbounds for k in eachindex(grp)
            grp[k] = rmzeros(selectdim(tmp_data, kdim, k))
        end
    else
        grp = Vector{Vector{Int}}()
    end

    return data, grp
end
# ============================================================================ #
rmzeros(x::AbstractArray{<:Real}) = x[x .> 0]
# ============================================================================ #
function segment(ts::Vector{Float64}, evt::Vector{Float64}, label::Vector{L},
    dur::Float64, bin_size::Float64) where L
    return segment(ts, evt, label, 0.0, dur, bin_size)
end
# ---------------------------------------------------------------------------- #
function segment(ts::Vector{Float64}, evt::Vector{Float64}, label::Vector{L},
    pre::Float64, post::Float64, bin_size::Float64) where L

    if length(label) != length(evt)
        error("Length of event label and event timestamp vectors must be equal")
    end

    kbins = round.(Int, -abs(pre)/bin_size) : round(Int, post/bin_size)
    data = psth(ts, evt, kbins, bin_size)[1]

    unique_label = sort(unique(label))

    nlabel = length(unique_label)
    grouped_data = Vector{Matrix{Float64}}(undef, nlabel)

    @inbounds for k = 1:nlabel
        group = findall(isequal(unique_label[k]), label)
        grouped_data[k] = data[:, group]
    end

    return grouped_data, unique_label
end
# ============================================================================ #


text:
using Distributions

import Tables

abstract type AbstractMatch <: Tables.AbstractRow end

Base.getproperty(m::AbstractMatch, nm::Symbol) = getfield(m, nm)
Base.propertynames(m::AbstractMatch) = fieldnames(typeof(m))
Tables.schema(m::AbstractMatch) = Tables.Schema(propertynames(m), Tuple(typeof.(m)))

abstract type AbstractMatches <: Tables.AbstractColumns end

Base.getproperty(m::AbstractMatches, nm::Symbol) = getfield(m, nm)
Base.propertynames(m::AbstractMatches) = fieldnames(typeof(m))
Tables.schema(m::AbstractMatches) = Tables.Schema(propertynames(m), Tuple(typeof(i[1]) for i in m))


text:
using NumDisplay

println("Starting...")

d = DisplayDirect([27, 22, 10, 9], (7, 8, 25, 24, 23, 18, 15, 14), inverted_sectors = true)
shutdown_mode_off(d)
decode_mode(d, 0x00)

###

println("High level method write_number()...")

data = [0,1,2,3]
write_number(d, data)
println("$data\n")
sleep(3)

data = [4,5,6,7]
write_number(d, data)
println("$data\n")
sleep(3)

data = [8,9,10,11]
write_number(d, data)
println("$data\n")
sleep(3)

data = [6,6,6,nothing]
write_number(d, data)
println("$data\n")
sleep(3)

write_number(d, data, 1)
println("$data with dot at 1\n")
sleep(3)

write_number(d, data, 2)
println("$data with dot at 2\n")
sleep(3)

write_number(d, data, 3)
println("$data with dot at 3\n")
sleep(3)

write_number(d, data, 4)
println("$data with dot at 4\n")
sleep(3)

data = 3141
write_number(d, data, 4)
println("$data with dot at 4\n")
sleep(3)

println("High level method write_symbols()...")

data = [0b01001001, 0b01001001, 0b01001001, 0b01001001]
write_symbols(d, data)
println("$data\n")
sleep(3)

data = ['I', 'L', 'U', 'J']
write_symbols(d, data)
println("$data\n")
sleep(3)

data = "Hell"
write_symbols(d, data)
println("$data\n")
sleep(3)

###

println("Low level methods")

decode_mode(d, 0b00000001)
println("decode_mode for digit 1 only\n")
sleep(3)

decode_mode(d)
println("decode_mode for all digits\n")
sleep(3)

decode_mode(d, 0b0)
println("decode_mode for none\n")
sleep(3)

println("change intensity")
for i in 1:16
    set_intensity(d, i)
    println("$i/16")
    sleep(0.5)
end
for i in collect(16:(-1):1)
    set_intensity(d, i)
    println("$i/16")
    sleep(0.5)
end
sleep(1)

println("change limit")
for i in collect(4:(-1):1)
    set_limit(d, i)
    println("$i/4")
    sleep(0.5)
end

for i in 1:4
    set_limit(d, i)
    println("$i/4")
    sleep(0.5)
end
sleep(1)

println("shutdown_mode_on()")
shutdown_mode_on(d)
sleep(3)

println("shutdown_mode_off()")
shutdown_mode_off(d)
sleep(3)

println("test_mode_on()")
test_mode_on(d)
sleep(3)

println("test_mode_off()")
test_mode_off(d)
sleep(3)

println("switch to decode and set 3,3,3,3 at positions 4,3,2,1")
decode_mode(d)
write_digit(d, UInt8(3), 4)
sleep(1)
write_digit(d, UInt8(3), 3)
sleep(1)
write_digit(d, UInt8(3), 2)
sleep(1)
write_digit(d, UInt8(3), 1)
sleep(1)

println("Stop.")


text:
module test_loltoolsdatadragon_datadragon

using Test
using LOLToolsDataDragon.DataDragon

@test DataDragon.version >= v"9.24.1"

champ = DataDragon.en_US.Champions[114]
@test champ.name == "Fiora"
spell = DataDragon.en_US.SummonerSpells[12]
@test spell.name == "Teleport"
item = DataDragon.en_US.Items[3111]
@test item.name == "Mercury's Treads"
@test item.gold == (base = 350, purchasable = true, total = 1100, sell = 770)
@test item.stats == (FlatMovementSpeedMod = 45, FlatSpellBlockMod = 25)

champ = DataDragon.ko_KR.Champions[114]
@test champ.name == "피오라"
spell = DataDragon.ko_KR.SummonerSpells[12]
@test spell.name == "순간이동"
item = DataDragon.ko_KR.Items[3111]
@test item.name == "헤르메스의 발걸음"
@test item.gold == (base = 350, purchasable = true, total = 1100, sell = 770)
@test item.stats == (FlatMovementSpeedMod = 45, FlatSpellBlockMod = 25)

end # module test_loltoolsdatadragon_datadragon


text:
# asd

a=1

# %% [markdown]
# asd
# + asd
# asd

a=1

# [markdown]
# # title
# content

c=1

# MD # title
# MD content

d=1

"""
    bar(x[, y])

julia standard markdown doc
"""
function bar(x, y)
    x
end

#=

a=1

# %% [markdown]
# asd
# + asd
# asd

a=1

# [markdown]
# # title
# content

c=1

# MD # title
# MD content

d=1

"""
    bar(x[, y])

julia standard markdown doc
"""
function bar(x, y)
    x
end

=#

max_stars_repo_name: zhaouv/vscode-markdown-everywhere

max_issues_repo_name: zhaouv/vscode-markdown-everywhere

max_forks_repo_name: zhaouv/vscode-markdown-everywhere


text:
## Look at musical artists
using DelimitedFiles
using SparseArrays
names = readlines("artistsim.names")
data = readdlm("artistsim.smat")
##
A = sparse(Int.(data[2:end,1]).+1,
           Int.(data[2:end,2]).+1,
           Int.(data[2:end,3]),
           Int(data[1,1]),Int(data[1,1]))
##
A = max.(A,A')
##
using PyCall
using Conda
pyimport_conda("igraph","python-igraph","conda-forge")
##
using PyCall
igraph = pyimport("igraph")

function igraph_layout(A::SparseMatrixCSC{T}, layoutname::AbstractString="drl") where T
    ei,ej,ew = findnz(A)
    edgelist = [(ei[i]-1,ej[i]-1) for i = 1:length(ei)]
    nverts = size(A)
    G = igraph.Graph(nverts, edges=edgelist, directed=false)
    xy = G.layout(layoutname)
    xy = [Float64(xy[i][j]) for i in 1:length(xy),  j in 1:length(xy[1])]
end

xy = igraph_layout(A, "drl")
##
using GraphRecipes
using LinearAlgebra
using Plots
ei, ej = findnz(triu(A,1))
graphplot(ei, ej, x =xy[:,1], y=xy[:,2],
  markercolor=:black, markerstrokecolor=:white,
  size=(1200,1200),dpi=300,
  markersize=2, linecolor=1, linealpha=0.01, linewidth=0.5,
  markeralpha=0.2,
  axis_buffer=0.02, background=nothing)
savefig("artistsim-plot.png")


text:
include("register.jl")
include("instruct.jl")
include("blocks.jl")
include("patch.jl")


text:
include(joinpath(@__DIR__, "..", "src", "ADMM.jl"))
include(joinpath(@__DIR__, "..", "src", "ClassificationInverseValueProblem.jl"))
using .ADMM
using .ClassificatinoInerseValueProblem

using LinearAlgebra
using JuMP, OSQP
using Test

function softthreshold(x::T, λ::T)::T where {T <: AbstractFloat}
    @assert λ >= 0 "Argument λ must exceed 0"
    x > λ && return (x - λ)
    x < -λ && return (x + λ)
    return zero(T)
end

"Solve LASSO problem for small matrix A ∈ R^{N×d} and sparse vector x."
function test_admm()
    N, d = 4, 5
    A = randn(N, d)
    x = zeros(d)
    x[1] = 1.0
    y = A * x
    λ = 0.001  # Lasso regularization weight
    ρ = 1.0  # ADMM proximal weight
    ϵ = 1e-5  # ADMM tolerance

    prox_f = x -> softthreshold.(x, λ / ρ)
    prox_g = x -> (A'A + ρ * I) \ (A'y + ρ * x)

    xr = ADMM.admm(zero(x), zero(x), prox_f, prox_g, ϵ)
    @test x ≈ xr atol = 1e-2
end

test_admm()

"""Minimizes f(δ, π) = δ + 1/(2λ) ‖ π - πk ‖^2, s.t., π'b - z ≥ δ, δ ≥ 0."""
function test_min_fpidelta()
    n = 20
    πk = randn(n)
    bi = randn(n)
    zi = rand()
    λ = 1.0

    model = Model(OSQP.Optimizer)
    @variable(model, δi >= 0)
    @variable(model, πi[1:n])
    @objective(model, Min, δi + 1 / (2 * λ) * sum((πi - πk) .^ 2))
    @constraint(model, con, πi' * bi - zi - δi >= 0.0)
    optimize!(model)

    πi_sol1 = value.(πi)
    πi_sol2 = ClassificatinoInerseValueProblem.min_fpidelta(πk, bi, zi, λ)
    @test πi_sol1 ≈ πi_sol2 atol = 1e-4 * n
end

test_min_fpidelta()

"""Minimizes f(ϵ, c; ck, cα, λ, xi) = ϵ'1 + 1/(2λ) ‖c - ck‖^2 + α/2 ‖c - cα‖^2,
s.t., c'xi - zi ≤ ϵi, ϵi ≥ 0, for all i.

"""
function test_fcepsilon()
    n = 20  # dimensions
    k = 10  # number of measurements
    xi = rand(n, k)  # feature vector
    zi = 2 * rand(k)  # measurements
    ck = rand(n)
    cα = rand(n)
    λ = 1.0
    α = 1.0

    model = Model(OSQP.Optimizer)
    @variable(model, ϵi[1:k] >= 0)
    @variable(model, c[1:n])
    @objective(
        model,
        Min,
        sum(ϵi) + 1.0 / 2.0 / λ * sum((c - ck) .^ 2) + α / 2.0 * sum((c - cα) .^ 2)
    )
    @constraint(model, con[i=1:k], xi[:, i]'c - zi[i] - ϵi[i] <= 0.0)
    optimize!(model)

    c_sol1 = value.(c)
    c_sol2 = ClassificatinoInerseValueProblem.min_fcepsilon(xi, zi, ck, cα, λ, α)
    @test c_sol1 ≈ c_sol2 rtol = 1e-3 * n
end

test_fcepsilon()

max_stars_repo_name: jzazo/InverseValueProblem

max_issues_repo_name: jzazo/InverseValueProblem

max_forks_repo_name: jzazo/InverseValueProblem


text:
using SodShockTube
using CairoMakie
using Test
using Documenter

doctest(SodShockTube)

@testset "SodShockTube.jl" begin
    problem = ShockTubeProblem(
        geometry = (0.0, 1.0, 0.5),
        left_state = (ρ = 1.0, u = 0.0, p = 1.0),
        right_state = (ρ = 0.125, u = 0.0, p = 0.1),
        t = 0.2,
        γ = 1.4
    );
    xs = LinRange(0.0, 1.0, 500);
    positions, regions, values = solve(problem, xs)
    @test positions["Shock"] ≈ 0.8504311464060357
    @test positions["Contact Discontinuity"] ≈ 0.6854905240097902
    @test positions["Head of rarefaction"] ≈ 0.26335680867601535
    @test positions["Foot of rarefaction"] ≈ 0.4859454374877634

    @test all(regions["Region 1"] .≈ (1.0, 1.0, 0.0))
    @test regions["Region 2"] == "RAREFACTION"
    @test all(regions["Region 3"] .≈ (0.30313017804679177, 0.42631942817462254, 0.9274526200369384))
    @test all(regions["Region 4"] .≈ (0.30313017805064707, 0.26557371170530725, 0.92745262004895057))
    @test all(regions["Region 5"] .≈ (0.1, 0.125, 0.0))

    f = Figure(resolution = (1000, 1000))
    ax_ρ = Axis(f[1,1], xlabel = "x", ylabel = "ρ", title = "Density")
    ax_u = Axis(f[2,1], xlabel = "x", ylabel = "u", title = "Velocity")
    ax_p = Axis(f[1,2], xlabel = "x", ylabel = "p", title = "Pressure")
    ax_E = Axis(f[2,2], xlabel = "x", ylabel = "E", title = "Stagnation Energy")

    opts = (;linewidth = 4)

    lines!(ax_ρ, values.x, values.ρ; opts...)
    lines!(ax_u, values.x, values.u; opts...)
    lines!(ax_p, values.x, values.p; opts...)
    lines!(ax_E, values.x, values.e; opts...)
end

max_stars_repo_name: archermarx/SodShockTube

max_issues_repo_name: archermarx/SodShockTube

max_forks_repo_name: archermarx/SodShockTube


text:
abstract type PandocFormat <: WeaveFormat end

function render_figures(docformat::PandocFormat, chunk)
    fignames = chunk.figures
    length(fignames) > 0 || (return "")

    caption = chunk.options[:fig_cap]
    label = get(chunk.options, :label, nothing)
    result = ""
    figstring = ""
    attribs = ""
    width = chunk.options[:out_width]
    height = chunk.options[:out_height]

    # Build figure attibutes
    attribs = String[]
    isnothing(width) || push!(attribs, "width=$width")
    isnothing(height) || push!(attribs, "height=$height")
    isnothing(label) || push!(attribs, "#fig:$label")
    attribs = isempty(attribs) ? "" : "{" * join(attribs, " ") * "}"

    if !isnothing(caption)
        result *= "![$caption]($(fignames[1]))$attribs\n"
        for fig in fignames[2:end]
            result *= "![]($fig)$attribs\n"
            println("Warning, only the first figure gets a caption\n")
        end
    else
        for fig in fignames
            result *= "![]($fig)$attribs\\ \n\n"
        end
    end
    return result
end

Base.@kwdef mutable struct Pandoc <: PandocFormat
    description = "Pandoc Markdown"
    extension = "md"
    codestart = "~~~~{.julia}"
    codeend = "~~~~~~~~~~~~~\n\n"
    termstart = codestart
    termend = codeend
    outputstart = "~~~~"
    outputend = "~~~~\n\n"
    # Prefer png figures for markdown conversion, svg doesn't work with latex
    mimetypes = ["image/png", "image/jpg", "image/svg+xml", "text/markdown", "text/plain"]
    fig_ext = ".png"
    out_width = nothing
    out_height = nothing
    fig_pos = nothing
    fig_env = nothing
    # specials
    preserve_header = true
end
register_format!("pandoc", Pandoc())

const DEFAULT_PANDOC_OPTIONS = String[]

Base.@kwdef mutable struct Pandoc2PDF <: PandocFormat
    description = "PDF via intermediate Pandoc Markdown"
    extension = "md"
    codestart = "~~~~{.julia}"
    codeend = "~~~~~~~~~~~~~\n\n"
    termstart = codestart
    termend = codeend
    outputstart = "~~~~"
    outputend = "~~~~\n\n"
    # Prefer png figures for markdown conversion, svg doesn't work with latex
    mimetypes = ["image/png", "image/jpg", "image/svg+xml", "text/markdown", "text/plain"]
    fig_ext = ".png"
    out_width = nothing
    out_height = nothing
    fig_pos = nothing
    fig_env = nothing
    # specials
    preserve_header = true
    header_template = normpath(TEMPLATE_DIR, "pandoc2pdf_header.txt")
    pandoc_options = DEFAULT_PANDOC_OPTIONS
end
register_format!("pandoc2pdf", Pandoc2PDF())

function set_format_options!(docformat::Pandoc2PDF; pandoc_options = DEFAULT_PANDOC_OPTIONS, _kwargs...)
    docformat.pandoc_options = pandoc_options
end


text:
module ViZDoom

using CxxWrap
@wrapmodule(joinpath(@__DIR__, "..", "deps", "usr", "ViZDoom-1.1.6", "bin", "libvizdoomjl.so"), :ViZDoom)

function __init__()
    @initcxx
end

include("util.jl")
export get_scenario_path, set_game

include("games/games.jl")
export basic_game

end # module


text:
#Example command line to run this program on multiple processors:
# julia -p 4 "run_evolve_goals.jl"

@everywhere include("EvolveGoals.jl")

# Example runs.  Uncomment one of these lines
#run_evolve_goals("prun",2,3,10,4)
#run_evolve_goals("rrun",4,15,200,3,4)
#run_evolve_goals("erun",3,14,50,6)
run_evolve_goals("erun",3,14,50,6,4)
#run_evolve_goals("erun",3,14,50,20)

# In run_evolve_goals("frun",3,14,50,6),
#   numinputs = 3
#   random seed = 14
#   genome size = 50
#   number of repetitions per goal = 6
# Evolves all goals of length 3 with 6 repetitions per goal

# In run_evolve_goals("frun",4,14,200,3,4),
#   numinputs = 4
#   random seed = 14
#   genome size = 200
#   number of repetitions per goal = 3
#   number of random goals = 4

# An alternative that does not involve running this file
# julia -p 2  -L EvolveGoals.jl -e "run_evolve_goals(\"prun\",3,16,20,3,4)"


text:
function DefaultIterativeLinearSolverEvent(arg0::Object, arg1::jint, arg2::RealVector, arg3::RealVector, arg4::RealVector, arg5::jdouble)
    return DefaultIterativeLinearSolverEvent((Object, jint, RealVector, RealVector, RealVector, jdouble), arg0, arg1, arg2, arg3, arg4, arg5)
end

function DefaultIterativeLinearSolverEvent(arg0::Object, arg1::jint, arg2::RealVector, arg3::RealVector, arg4::jdouble)
    return DefaultIterativeLinearSolverEvent((Object, jint, RealVector, RealVector, jdouble), arg0, arg1, arg2, arg3, arg4)
end

function get_iterations(obj::IterationEvent)
    return jcall(obj, "getIterations", jint, ())
end

function get_norm_of_residual(obj::DefaultIterativeLinearSolverEvent)
    return jcall(obj, "getNormOfResidual", jdouble, ())
end

function get_residual(obj::DefaultIterativeLinearSolverEvent)
    return jcall(obj, "getResidual", RealVector, ())
end

function get_right_hand_side_vector(obj::DefaultIterativeLinearSolverEvent)
    return jcall(obj, "getRightHandSideVector", RealVector, ())
end

function get_solution(obj::DefaultIterativeLinearSolverEvent)
    return jcall(obj, "getSolution", RealVector, ())
end

function get_source(obj::EventObject)
    return jcall(obj, "getSource", Object, ())
end

function provides_residual(obj::DefaultIterativeLinearSolverEvent)
    return jcall(obj, "providesResidual", jboolean, ())
end

function to_string(obj::EventObject)
    return jcall(obj, "toString", JString, ())
end


text:
module GMC_NS
    using  Distributions, Serialization, UnicodePlots
    import ProgressMeter: AbstractProgress, Progress, @showprogress, next!, move_cursor_up_while_clearing_lines, printover, durationstring
    import Printf: @sprintf
    import Random: rand
    import Statistics: unscaled_covzm, det, eigen
    import Serialization: serialize, deserialize
    import LinearAlgebra: dot, normalize
    import BioBackgroundModels: lps
    import StatsFuns: logaddexp, logsumexp, logsubexp
    import Base: show
    import Measurements: measurement
    import ConjugatePriors: NormalGamma, NormalInverseGamma
    import NGRefTools: MarginalTDist
    import KernelDensity: kde, AbstractKDE

    #GMC settings vector fields: [GMC_Nmin::Int64, GMC_τ_death::Float64, GMC_init_τ::Float64, GMC_tune_μ::Int64, GMC_tune_α::Float64, GMC_tune_PID::NTuple{3,Float64}, GMC_timestep_η::Float64, GMC_reflect_η::Float64, GMC_exhaust_σ::Float64, GMC_chain_κ::Int64]

    #should GMC_reflect_η be smaller?
    GMC_DEFAULTS=Vector{Any}([50,1e-5,1., 4,.25,(.3,0.,.0), .25, .005, 10., typemax(Int64)])
    export GMC_DEFAULTS

    include("ensemble/GMC_NS_Ensemble.jl")
    export GMC_NS_Ensemble
    include("model/GMC_NS_Model.jl")
    export GMC_NS_Model, GMC_NS_Model_Record
    include("model/model_utilities.jl")
    include("GMC/galilean_trajectory.jl")
    include("utilities/t_Tuner.jl")
    include("utilities/ellipsoid.jl")
    include("nested_sampler/search_patterns.jl")
    include("nested_sampler/nested_step.jl")
    include("nested_sampler/converge_ensemble.jl")
    export converge_ensemble!
    include("utilities/coordinate_utils.jl")
    export box_bound!, to_prior, to_unit_ball, box_reflect!
    include("ensemble/ensemble_utilities.jl")
    export ensemble_history, clean_ensemble_dir, measure_evidence, reset_ensemble!, move_ensemble!, copy_ensemble, rewind_ensemble, show_models, show_models_e, get_model, rectify_ensemble!, posterior_kde
    include("utilities/ns_progressmeter.jl")
    include("utilities/progress_displays.jl")
    export tuning_display,evidence_display,convergence_display,info_display,lh_display,liwi_display,ensemble_display,model_display,model_obs_display
    include("utilities/stats.jl")
    include("model/normal/Normal_Model.jl")
    include("model/normal/Normal_Ensemble.jl")
    export Normal_Ensemble
    include("model/normal/LogNormal_Model.jl")
    include("model/normal/LogNormal_Ensemble.jl")
    export LogNormal_Ensemble
    include("model/eggbox/Eggbox_Model.jl")
    include("model/eggbox/Eggbox_Ensemble.jl")
    export Eggbox_Ensemble
end # module


text:
function lbinom(n::T, s::T) where T <: Integer
    if n < 60
        log(binomial(n,s))
    else
        logabsbinomial(n,s)[1]
    end
end

# A simple implementation of a DFT to avoid introducing a dependency
# on an external FFT package just for this one distribution
# adapted from Distributions.jl
function _dft!(y::Vector{T}, x::Vector{T}) where T
    n = length(x)
    y .= zero(T)
    @inbounds for j = 0:n-1, k = 0:n-1
        y[k+1] += x[j+1] * cis(-π * T(2 * mod(j * k, n)) / n)
    end
end

#adopted from DSP.jl
function _filt!(out::AbstractArray{T,1}, b::AbstractArray{T,1}, x::AbstractArray{T,1}, si::AbstractArray{T,1}) where T <: Real
    silen = length(si)
    @inbounds for i in 1:length(x)
        xi = x[i]
        val = muladd(xi, b[1], si[1])
        for j in 1:(silen-1)
            si[j] = muladd(xi, b[j+1], si[j+1])
        end
        si[silen] = b[silen+1]*xi
        out[i] = val
    end
end

function _filt_reg!(out::AbstractArray{T,1}, b::AbstractArray{T,1}, x::AbstractArray{T,1}, si::AbstractArray{T,1}, k1::D, k2::D) where {T <: Real, D <: Integer}
    silen = length(si)
    @inbounds for j in 1:(silen-1)
        adj1 = exp(lbinom(k1,j))
        si[j] = x[1]*b[j+1]*adj1
    end
    adj1 = exp(lbinom(k1,silen))
    si[silen] = b[silen+1]*x[1]*adj1
    out[1] = one(T)

    @inbounds for i in 2:length(x)
        xi = x[i]
        adj1 = exp(lbinom(k2,i-1)-lbinom(k1+k2,i-1))
        adj2 = (i-1) / (k1+k2-i+2)
        val = muladd(xi, b[1]*adj1, si[1]*adj2)
        for j in 1:(silen-1)
            adj1 = exp(lbinom(k1,j)+lbinom(k2,i-1)-lbinom(k1+k2,i-1))
            si[j] = muladd(xi, b[j+1]*adj1, si[j+1]*adj2)
        end
        adj1 = exp(lbinom(k1,silen)+lbinom(k2,i-1)-lbinom(k1+k2,i-1))
        si[silen] = b[silen+1]*xi*adj1
        out[i] = val
    end
end


text:
@with_kw struct LocalProjection
    df                  ::DataFrame
    horizons            ::Array{Int,1} = [1]
    models              ::Array{FixedEffectModel,1} = Array{FixedEffectModel,1}(undef, length(collect(horizons)))
    transformfunc       ::Symbol = :identity
    grwthoffset         ::Int = 0
    formula             ::FormulaTerm
    regstyle            ::Symbol = :rates
    ynames              ::Array{Symbol,1} = Array{Symbol,1}(undef, length(collect(horizons)))
end

#######################################################################
# LocalProjection functions
#######################################################################
function makedepvars(lp::LocalProjection)
    if lp.regstyle == :rates
        FΔ!(lp.df, Symbol(lp.formula.lhs), lp.horizons, getfield(Main, lp.transformfunc), lp.grwthoffset)
        lp.ynames .= [Symbol("FΔ", i, lp.formula.lhs) for i in lp.horizons]
    elseif lp.regstyle == :levels
        F!(lp.df, Symbol(lp.formula.lhs), lp.horizons, getfield(Main, lp.transformfunc))
        lp.ynames .= [Symbol("F", i, lp.formula.lhs) for i in lp.horizons]
    end
end

function runregs(lp::LocalProjection)
    for (index, value) in enumerate(lp.horizons)
        form = Term(lp.ynames[index]) ~ lp.formula.rhs
        lp.models[index] = reg(lp.df, form)
    end
end

function coef(lp::LocalProjection, xname::Symbol)
    idx = findfirst(Symbol.(coefnames(lp.models[1])) .== xname)
    β = [m.coef[idx] for m in lp.models]
    return β
end

function se(lp::LocalProjection, xname::Symbol)
    idx = findfirst(Symbol.(coefnames(lp.models[1])) .== xname)
    se = [sqrt(vcov(m)[idx, idx]) for m in lp.models]
    return se
end

function coef_se(lp::LocalProjection, xname::Symbol)
    return DataFrame(horizons = lp.horizons, β = coef(lp, xname), se = se(lp, xname))
end


text:
using Test
using WorldBankData
using DataFrames
using JSON

# http://api.worldbank.org/v2/countries/US/indicators/NY.GNP.PCAP.CD?date=1990:1997&format=json
testjsonstr = "[{\"page\":1,\"pages\":1,\"per_page\":50,\"total\":8,\"sourceid\":\"2\",\"lastupdated\":\"2020-07-01\"},[{\"indicator\":{\"id\":\"NY.GNP.PCAP.CD\",\"value\":\"GNI per capita, Atlas method (current US\$)\"},\"country\":{\"id\":\"US\",\"value\":\"United States\"},\"countryiso3code\":\"USA\",\"date\":\"1997\",\"value\":31270,\"unit\":\"\",\"obs_status\":\"\",\"decimal\":0},{\"indicator\":{\"id\":\"NY.GNP.PCAP.CD\",\"value\":\"GNI per capita, Atlas method (current US\$)\"},\"country\":{\"id\":\"US\",\"value\":\"United States\"},\"countryiso3code\":\"USA\",\"date\":\"1996\",\"value\":30270,\"unit\":\"\",\"obs_status\":\"\",\"decimal\":0},{\"indicator\":{\"id\":\"NY.GNP.PCAP.CD\",\"value\":\"GNI per capita, Atlas method (current US\$)\"},\"country\":{\"id\":\"US\",\"value\":\"United States\"},\"countryiso3code\":\"USA\",\"date\":\"1995\",\"value\":29040,\"unit\":\"\",\"obs_status\":\"\",\"decimal\":0},{\"indicator\":{\"id\":\"NY.GNP.PCAP.CD\",\"value\":\"GNI per capita, Atlas method (current US\$)\"},\"country\":{\"id\":\"US\",\"value\":\"United States\"},\"countryiso3code\":\"USA\",\"date\":\"1994\",\"value\":27650,\"unit\":\"\",\"obs_status\":\"\",\"decimal\":0},{\"indicator\":{\"id\":\"NY.GNP.PCAP.CD\",\"value\":\"GNI per capita, Atlas method (current US\$)\"},\"country\":{\"id\":\"US\",\"value\":\"United States\"},\"countryiso3code\":\"USA\",\"date\":\"1993\",\"value\":26390,\"unit\":\"\",\"obs_status\":\"\",\"decimal\":0},{\"indicator\":{\"id\":\"NY.GNP.PCAP.CD\",\"value\":\"GNI per capita, Atlas method (current US\$)\"},\"country\":{\"id\":\"US\",\"value\":\"United States\"},\"countryiso3code\":\"USA\",\"date\":\"1992\",\"value\":25680,\"unit\":\"\",\"obs_status\":\"\",\"decimal\":0},{\"indicator\":{\"id\":\"NY.GNP.PCAP.CD\",\"value\":\"GNI per capita, Atlas method (current US\$)\"},\"country\":{\"id\":\"US\",\"value\":\"United States\"},\"countryiso3code\":\"USA\",\"date\":\"1991\",\"value\":24270,\"unit\":\"\",\"obs_status\":\"\",\"decimal\":0},{\"indicator\":{\"id\":\"NY.GNP.PCAP.CD\",\"value\":\"GNI per capita, Atlas method (current US\$)\"},\"country\":{\"id\":\"US\",\"value\":\"United States\"},\"countryiso3code\":\"USA\",\"date\":\"1990\",\"value\":24060,\"unit\":\"\",\"obs_status\":\"\",\"decimal\":0}]]"

dfref = DataFrame(iso2c=["US" for _ in range(1, length=8)],
    country=["United States" for _ in range(1, length=8)],
    year=[ 1997.0, 1996.0, 1995.0, 1994.0, 1993.0, 1992.0, 1991.0, 1990.0],
    NY_GNP_PCAP_CD=[31270.0, 30270.0, 29040.0, 27650.0, 26390.0, 25680.0, 24270.0, 24060.0])

testjson = JSON.parse(testjsonstr)[2]

df = WorldBankData.parse_wdi(["NY.GNP.PCAP.CD"], testjson, 1990, 1997, false)

@test df[!,:year] == dfref[!,:year]
@test sort(names(df)) == sort(names(dfref))
@test df == dfref