text:
# Handers for execute_request and related messages, which are
# the core of the Jupyter protocol: execution of Julia code and
# returning results.

import Base.Libc: flush_cstdio

const text_plain = MIME("text/plain")
const image_svg = MIME("image/svg+xml")
const image_png = MIME("image/png")
const image_jpeg = MIME("image/jpeg")
const text_markdown = MIME("text/markdown")
const text_html = MIME("text/html")
const text_latex = MIME("text/latex") # Jupyter expects this
const text_latex2 = MIME("application/x-latex") # but this is more standard?
const application_vnd_vega_v3 = MIME("application/vnd.vega.v3+json")
const application_vnd_vegalite_v2 = MIME("application/vnd.vegalite.v2+json")
const application_vnd_dataresource = MIME("application/vnd.dataresource+json")

include("magics.jl")

# return a String=>Any dictionary to attach as metadata
# in Jupyter display_data and pyout messages
metadata(x) = Dict()

# return a String=>String dictionary of mimetype=>data
# for passing to Jupyter display_data and execute_result messages.
function display_dict(x)
    data = Dict{String,Any}("text/plain" => limitstringmime(text_plain, x))
    if showable(application_vnd_vegalite_v2, x)
        data[string(application_vnd_vegalite_v2)] = JSON.JSONText(limitstringmime(application_vnd_vegalite_v2, x))
    elseif showable(application_vnd_vega_v3, x) # don't send vega if we have vega-lite
        data[string(application_vnd_vega_v3)] = JSON.JSONText(limitstringmime(application_vnd_vega_v3, x))
    end
    if showable(application_vnd_dataresource, x)
        data[string(application_vnd_dataresource)] = JSON.JSONText(limitstringmime(application_vnd_dataresource, x))
    end
    if showable(image_svg, x)
        data[string(image_svg)] = limitstringmime(image_svg, x)
    end
    if showable(image_png, x)
        data[string(image_png)] = limitstringmime(image_png, x)
    elseif showable(image_jpeg, x) # don't send jpeg if we have png
        data[string(image_jpeg)] = limitstringmime(image_jpeg, x)
    end
    if showable(text_markdown, x)
        data[string(text_markdown)] = limitstringmime(text_markdown, x)
    elseif showable(text_html, x)
        data[string(text_html)] = limitstringmime(text_html, x)
    elseif showable(text_latex, x)
        data[string(text_latex)] = limitstringmime(text_latex, x)
    elseif showable(text_latex2, x)
        data[string(text_latex)] = limitstringmime(text_latex2, x)
    end
    return data
end

# queue of objects to display at end of cell execution
const displayqueue = Any[]

# remove x from the display queue
function undisplay(x)
    i = Compat.findfirst(isequal(x), displayqueue)
    i !== nothing && splice!(displayqueue, i)
    return x
end

if VERSION < v"0.7.0-DEV.3498" # julia #25544
    import Base.REPL: ip_matches_func
else
    import Base: ip_matches_func
end

function show_bt(io::IO, top_func::Symbol, t, set)
    # follow PR #17570 code in removing top_func from backtrace
    eval_ind = Compat.findlast(addr->ip_matches_func(addr, top_func), t)
    eval_ind !== nothing && (t = t[1:eval_ind-1])
    Base.show_backtrace(io, t)
end

# wrapper for showerror(..., backtrace=false) since invokelatest
# doesn't support keyword arguments.
showerror_nobt(io, e, bt) = showerror(io, e, bt, backtrace=false)

@static if VERSION < v"0.7.0-DEV.4724"
		# https://github.com/JuliaLang/Compat.jl/pull/572
    # thanks to https://github.com/jipolanco/WriteVTK.jl/commit/ea046493c7787ae651b5f629455aaf7bdaa000b1
    rsplit(s::AbstractString; limit::Integer=0, keepempty::Bool=false) =
        Base.rsplit(s)
    rsplit(s::AbstractString, splitter; limit::Integer=0, keepempty::Bool=false) =
        Base.rsplit(s, splitter; limit=limit, keep=keepempty)
    split(s::AbstractString; limit::Integer=0, keepempty::Bool=false) =
        Base.split(s)
    split(s::AbstractString, splitter; limit::Integer=0, keepempty::Bool=false) =
        Base.split(s, splitter; limit=limit, keep=keepempty)
end

# return the content of a pyerr message for exception e
function error_content(e, bt=catch_backtrace(); backtrace_top::Symbol=:include_string, msg::AbstractString="")
    tb = map(x->String(x), split(sprint(show_bt,
                                        backtrace_top,
                                        bt, 1:typemax(Int)),
                                  "\n", keepempty=true))

    ename = string(typeof(e))
    evalue = try
        # Peel away one LoadError layer that comes from running include_string on the cell
        isa(e, LoadError) && (e = e.error)
        sprint((io, e, bt) -> invokelatest(showerror_nobt, io, e, bt), e, bt)
    catch
        "SYSTEM: show(lasterr) caused an error"
    end
    pushfirst!(tb, evalue) # fperez says this needs to be in traceback too
    if !isempty(msg)
        pushfirst!(tb, msg)
    end
    Dict("ename" => ename, "evalue" => evalue,
                 "traceback" => tb)
end

#######################################################################

# global variable so that display can be done in the correct Msg context
execute_msg = Msg(["julia"], Dict("username"=>"jlkernel", "session"=>uuid4()), Dict())
# global variable tracking the number of bytes written in the current execution
# request
const stdio_bytes = Ref(0)

if VERSION < v"0.7.0-DEV.3589" # julia #25738
    import Base.Docs: helpmode # added in 0.6, see julia #19858
else
    import REPL: helpmode
end

# use a global array to accumulate "payloads" for the execute_reply message
const execute_payloads = Dict[]

const stdout_name = isdefined(Base, :stdout) ? "stdout" : "STDOUT"

function execute_request(socket, msg)
    code = msg.content["code"]
    @vprintln("EXECUTING ", code)
    global execute_msg = msg
    global n, In, Out, ans
    stdio_bytes[] = 0
    silent = msg.content["silent"]
    store_history = get(msg.content, "store_history", !silent)
    empty!(execute_payloads)

    if !silent
        n += 1
        send_ipython(publish[],
                     msg_pub(msg, "execute_input",
                             Dict("execution_count" => n,
                                          "code" => code)))
    end

    silent = silent || REPL.ends_with_semicolon(code)
    if store_history
        In[n] = code
    end

    # "; ..." cells are interpreted as shell commands for run
    code = replace(code, r"^\s*;.*$" =>
                   m -> string(replace(m, r"^\s*;", "Base.repl_cmd(`"),
                               "`, ", stdout_name, ")"))

    # a cell beginning with "? ..." is interpreted as a help request
    hcode = replace(code, r"^\s*\?" => "")

    try
        for hook in preexecute_hooks
            invokelatest(hook)
        end

        ans = result = if hcode != code # help request
            Core.eval(Main, helpmode(hcode))
        else
            #run the code!
            occursin(magics_regex, code) ? magics_help(code) :
                include_string(current_module[], code, "In[$n]")
        end

        if silent
            result = nothing
        elseif result !== nothing
            if store_history
                Out[n] = result
            end
        end

        user_expressions = Dict()
        for (v,ex) in msg.content["user_expressions"]
            user_expressions[v] = invokelatest(parse, ex)
        end

        for hook in postexecute_hooks
            invokelatest(hook)
        end

        # flush pending stdio
        flush_all()

        undisplay(result) # dequeue if needed, since we display result in pyout
        invokelatest(display) # flush pending display requests

        if result !== nothing
            result_metadata = invokelatest(metadata, result)
            result_data = invokelatest(display_dict, result)
            send_ipython(publish[],
                         msg_pub(msg, "execute_result",
                                 Dict("execution_count" => n,
                                              "metadata" => result_metadata,
                                              "data" => result_data)))

        end
        send_ipython(requests[],
                     msg_reply(msg, "execute_reply",
                               Dict("status" => "ok",
                                            "payload" => execute_payloads,
                                            "execution_count" => n,
                                            "user_expressions" => user_expressions)))
        empty!(execute_payloads)
    catch e
        bt = catch_backtrace()
        try
            # flush pending stdio
            flush_all()
            for hook in posterror_hooks
                invokelatest(hook)
            end
        catch
        end
        empty!(displayqueue) # discard pending display requests on an error
        content = error_content(e,bt)
        send_ipython(publish[], msg_pub(msg, "error", content))
        content["status"] = "error"
        content["execution_count"] = n
        send_ipython(requests[], msg_reply(msg, "execute_reply", content))
    end
end


text:
## Fichier qui créé els graphiques grace aux sorties d'un ReactiveGRASP

#
#Version modifiée à partire de celle de Mr Gandibleux
#

function simulation2(matrix,cost, n, m, ite, coupe, alphaset,f)

    zinit = zeros(Int64, ite) # zero
    zls   = zeros(Int64, ite) # zero
    zbest = zeros(Int64, ite) # zero

    x     = zeros(Int64, nbDivisionRun)
    zmax  = Matrix{Int64}(undef,nbInstances , nbDivisionRun); zmax[:] .= typemin(Int64)  # -Inf entier
    zmoy  = zeros(Float64, nbInstances, nbDivisionRun) # zero
    zmin  = Matrix{Int64}(undef,nbInstances , nbDivisionRun) ; zmin[:] .= typemax(Int64)  # +Inf entier
    tmoy  = zeros(Float64, nbInstances)  # zero

    # calcule la valeur du pas pour les divisions

    println("Experimentation ReactiveGRASP-SPP avec :")
    println("  nbIterationGrasp  = ", ite)

    println(" ")
    cpt = 0

    # run non comptabilise (afin de produire le code compile)
    #zinit, zls, zbest = graspSPP(allfinstance[1], 0.5, 1)

    print("  ",f," : ")
    evol_p,liste_zavg,liste_zmax,liste_zmin,z_rouge,z_vert,ligne_verte = ReactiveGRASP(matrix,cost, n, m, ite, coupe, alphaset)
    gr = GRASP(cost,matrix,n,m,alpha)

    println(" ")

    #Pkg.add("PyPlot") # Mandatory before the first use of this package
    println(" ");println("  Graphiques de synthese")
#    using PyPlot
    instancenb = 1
    plotRunGrasp(f, zinit, zls, zbest)
    plotAnalyseGrasp(f, 1:ite, liste_zavg, liste_zmin, liste_zmax )
    plotCPUt(allfinstance, tmoy)
end

ext: jl


text:
export QdldlKKTSolver, CholmodKKTSolver
# -------------------------------------
# abstract type defs
# -------------------------------------
abstract type AbstractKKTSolver end

# NB: all concrete examples of this type should
# implement refactor! and solve! methods and should
# implement a constructor taking (P,A,sigma,rho) as
# arguments

# -------------------------------------
# some internal utility functions
# -------------------------------------

function _kktutils_check_dims(P, A, sigma, rho)

    n = size(P, 1)
    m = size(A, 1)

    size(A,2) == n || throw(DimensionMismatch())

    length(rho)   == m || length(rho)   == 1 || throw(DimensionMismatch())
    length(sigma) == n || length(sigma) == 1 || throw(DimensionMismatch())

    return m, n
end

function _kktutils_make_kkt(P, A, sigma, rho, shape::Symbol=:F)

    S = length(sigma) == 1 ? (sigma[1]) * I : Diagonal(sigma)
    n = size(P, 1)
    m = size(A, 1)

    if length(rho) == 1
        rho = rho .* ones(m)
    end

    if     shape == :F
        #compute the full KKT matrix
        K = [P + S SparseMatrixCSC(A'); A -I]

    elseif shape == :U
        #upper triangular
        K = [triu(P) + S  SparseMatrixCSC(A'); spzeros(eltype(A), m, n)  -I]

    elseif shape == :L
        #lower triangular
        K = [tril(P)+S  spzeros(eltype(A), n, m); A  -I]

    else
        error("Bad matrix shape description")
    end

    @inbounds @simd for i = (n + 1):(n + m)
        K[i, i] = -1.0 / rho[i - n]
    end

    return K

end

function update_kkt_matrix(K, n::Int, m::Int, rho)
    if length(rho) == m
        @inbounds @simd for i = (n + 1):(n + m)
            K[i, i] = -1.0 / rho[i - n]
        end
    elseif length(rho) == 1
        @inbounds @simd for i = (n + 1):(n + m)
            K[i, i] = -1.0 / rho[1]
        end
    else
        throw(DimensionMismatch)
    end
end

# -------------------------------------
# QDLDL solver
# -------------------------------------

struct QdldlKKTSolver{Tv, Ti} <: AbstractKKTSolver

    m::Integer
    n::Integer
    ldlfact::QDLDL.QDLDLFactorisation{Tv, Ti}

    function QdldlKKTSolver(P::SparseMatrixCSC{Tv, Ti}, A::SparseMatrixCSC{Tv, Ti}, sigma, rho) where {Tv, Ti}

        m,n = _kktutils_check_dims(P, A, sigma, rho)

        #NB: qdldl uses triu internally, but it reorders
        #with AMD first.  This way is memory inefficient
        #but saves having to permute the rhs/lhs each
        #time we solve.
        K   = _kktutils_make_kkt(P, A, sigma, rho, :F)
        ldlfact = qdldl(K)

        #check for exactly n positive eigenvalues
        positive_inertia(ldlfact) == n || error("Objective function is not convex.")

        new{Tv, Ti}(m, n, ldlfact)
    end
end

function solve!(s::QdldlKKTSolver, lhs, rhs)
    lhs .= rhs
    QDLDL.solve!(s.ldlfact, lhs)
end

function update_rho!(s::QdldlKKTSolver, rho)
    QDLDL.update_diagonal!(s.ldlfact, (s.n+1):(s.n+s.m),(-1 ./ rho))
end

# -------------------------------------
# Julia Native solver (CHOLMOD based)
# -------------------------------------
mutable struct CholmodKKTSolver{Tv, Ti} <: AbstractKKTSolver

    fact::SuiteSparse.CHOLMOD.Factor{Tv}
    K::SparseMatrixCSC{Tv, Ti}
    m::Integer
    n::Integer

    function CholmodKKTSolver(P::SparseMatrixCSC{Tv, Ti}, A::SparseMatrixCSC{Tv, Ti}, sigma, rho) where {Tv, Ti}

         m,n  = _kktutils_check_dims(P, A, sigma, rho)
         K    = _kktutils_make_kkt(P, A, sigma, rho, :F)
         fact = ldlt(K)

        return new{Tv, Ti}(fact, K, m, n)

    end
end

solve!(s::CholmodKKTSolver, lhs, rhs) = lhs .= s.fact \ rhs

function update_rho!(s::CholmodKKTSolver, rho)
    update_kkt_matrix(s.K, s.n, s.m, rho)
    #complete restart
    s.fact = ldlt(s.K)
end

free_memory!(s::AbstractKKTSolver) = nothing


text:
include("include2.jl")

foo(s)


text:
function sphere_streamline(linebuffer, ∇ˢf, pt, h, n)
    push!(linebuffer, pt)
    df = normalize(∇ˢf(pt[1], pt[2], pt[3]))
    push!(linebuffer, normalize(pt .+ h*df))
    for k=2:n
        cur_pt = last(linebuffer)
        push!(linebuffer, cur_pt)
        df = normalize(∇ˢf(cur_pt...))
        push!(linebuffer, normalize(cur_pt .+ h*df))
    end
    return
end

"""
    StreamLines

TODO add function signatures
TODO add descripton

## Attributes
$(ATTRIBUTES)
"""
@recipe(StreamLines, points, directions) do scene
    Attributes(
        h = 0.01f0,
        n = 5,
        color = :black,
        linewidth = 1
    )
end


text:
module HoJMTK

using Reexport
@reexport using ModelingToolkit
@reexport using DifferentialEquations

include("Electrical/Electrical.jl")
end


text:
module VectorsOfBlockArrayCooTests

using Test
using Gridap.Arrays
using BlockArrays
using FillArrays
using LinearAlgebra
using Gridap.Arrays: BlockArrayCooMap

l = 10
b11 = [  i*[1 2; 3 4]  for i in 1:l ]
b21 = [  i*[5 6; 7 8; 9 10] for i in 1:l]
b32 = [  i*[5 6 7 8; 7 8 1 3] for i in 1:l]
blockids = [(1,1),(2,1),(3,2)]
ax = Fill((blockedrange([2,3,2]), blockedrange([2,4])),l)
al = lazy_map(BlockArrayCooMap((3,2),blockids),ax,b11,b21,b32)

#using BenchmarkTools
#c = array_cache(al)
#@btime getindex!($c,$al,1)

b1 = [  i*[3,4] for i in 1:l]
b2 = [  i*[3,4,5,6] for i in 1:l]
blockids = [(1,),(2,)]
ax = Fill((blockedrange([2,4]),),l)
bl = lazy_map(BlockArrayCooMap((2,),blockids),ax,b1,b2)

b2 = [  i*[3,4,5,6] for i in 1:l]
blockids = [(2,)]
ax = Fill((blockedrange([2,4]),),l)
bl = lazy_map(BlockArrayCooMap((2,),blockids),ax,b2)

c1 = [  i*[7,8] for i in 1:l]
c2 = [  i*[3,4,5] for i in 1:l]
blockids = [(1,),(2,)]
ax = Fill((blockedrange([2,3,2]),),l)
cl = lazy_map(BlockArrayCooMap((3,),blockids),ax,c1,c2)

dl = lazy_map(*,al,bl)
test_array(dl,[ a*b for (a,b) in zip(al,bl) ])

#using BenchmarkTools
#cache = array_cache(dl)
#@btime getindex!($cache,$dl,2)

dl = lazy_map(MulAddMap(2,3),al,bl,cl)
test_array(dl,[ 2*a*b + 3*c for (a,b,c) in zip(al,bl,cl) ])

#using BenchmarkTools
#cache = array_cache(dl)
#@btime getindex!($cache,$dl,2)

dl = lazy_map(transpose,al)
test_array(dl,transpose.(al))

dl = lazy_map(*,al,lazy_map(transpose,al))
test_array(dl,[ a*transpose(a) for a in al ])
#print_op_tree(dl)

#using BenchmarkTools
#cache = array_cache(dl)
#@btime getindex!($cache,$dl,2)

# in-homogeneous case

l1 = 4
l2 = 6
a11 = vcat([i*[1 2; 3 4] for i in 1:l1 ],[i*[1 2 3; 3 4 5; 6 7 8] for i in 1:l2])
a12 = vcat([10*i*[1 2; 3 4] for i in 1:l1 ],[10*i*[1 2 3; 3 4 5; 6 7 8] for i in 1:l2])
blockids = [(1,1),(1,2)]
ax1 = (blockedrange([2,2]),blockedrange([2,2]))
ax2 = (blockedrange([3,3]),blockedrange([3,3]))
ax = CompressedArray([ax1,ax2],vcat(fill(1,l1),fill(2,l2)))
al = lazy_map(BlockArrayCooMap((2,2),blockids),ax,a11,a12)

b2 = vcat([i*[1,2] for i in 1:l1 ],[i*[1,2,3] for i in 1:l2])
blockids = [(2,)]
ax1 = (blockedrange([2,2]),)
ax2 = (blockedrange([3,3]),)
ax = CompressedArray([ax1,ax2],vcat(fill(1,l1),fill(2,l2)))
bl = lazy_map(BlockArrayCooMap((2,),blockids),ax,b2)

dl = lazy_map(*,al,bl)
test_array(dl,[ a*b for (a,b) in zip(al,bl) ])

dl = lazy_map(MulAddMap(2,3),al,bl,bl)
test_array(dl,[ 2*a*b + 3*c for (a,b,c) in zip(al,bl,bl) ])

dl = lazy_map(transpose,al)
test_array(dl,transpose.(al))

dl = lazy_map(*,al,lazy_map(transpose,al))
test_array(dl,[ a*transpose(a) for a in al ])

# Blocks of blocks (in-homogeneous case)

blockids = [(1,2),]
_ax1 = (blockedrange([2,2]),blockedrange([2,2]))
_ax2 = (blockedrange([3,3]),blockedrange([3,3]))
ax1 = (append_ranges([_ax1[1],_ax1[1]]),append_ranges([_ax1[2],_ax1[2]]))
ax2 = (append_ranges([_ax2[1],_ax2[1]]),append_ranges([_ax2[2],_ax2[2]]))
ax = CompressedArray([ax1,ax2],vcat(fill(1,l1),fill(2,l2)))
aBl = lazy_map(BlockArrayCooMap((2,2),blockids),ax,al)

blocks = (bl,)
blockids = [(1,)]
_ax1 = (blockedrange([2,2]),)
_ax2 = (blockedrange([3,3]),)
ax1 = (append_ranges([_ax1[1],_ax1[1]]),)
ax2 = (append_ranges([_ax2[1],_ax2[1]]),)
ax = CompressedArray([ax1,ax2],vcat(fill(1,l1),fill(2,l2)))
bBl = lazy_map(BlockArrayCooMap((2,),blockids),ax,bl)

dl = lazy_map(*,aBl,bBl)
test_array(dl,[ a*b for (a,b) in zip(aBl,bBl) ])

dl = lazy_map(MulAddMap(2,3),aBl,bBl,bBl)
test_array(dl,[ 2*a*b + 3*c for (a,b,c) in zip(aBl,bBl,bBl) ])

#using BenchmarkTools
#cache = array_cache(dl)
#@btime getindex!($cache,$dl,2)

dl = lazy_map(*,aBl,lazy_map(transpose,aBl))
test_array(dl,[ a*transpose(a) for a in aBl ])

end # module


text:
"""
    edit_with(y::Dict{Any,Any}; kwargs...)
    edit_with(df::DataFrame, editor::T) where T<:Edit
    edit_with(df::DataFrame, lst::Array{T}) where T<:Edit
    edit_with(df::DataFrame, x::Describe, file::T) where T<:File
    edit_with(file::T, y::Dict{Any,Any}; kwargs...)
    edit_with(files::Array{T,N} where N, y::Dict{Any,Any}; kwargs...) where T<:File
This function edits the input DataFrame `df` and returns the resultant DataFrame.

# Arguments
- `df::DataFrame` on which to perform the edit.
- `editor::T where T<:Edit`: DataType containing information about which edit to perform.
    The following edit options are available and detailed below. If given a dictionary of
    edits, they will be made in this order:
    - [`SLiDE.Drop`](@ref): Remove information from the DataFrame -- either an entire column
        or rows containing specified values.
    - [`SLiDE.Rename`](@ref): Change column name `from` -> `to`.
    - [`SLiDE.Group`](@ref): Use to edit files containing data in successive dataframes with
        an identifying header cell or row.
    - [`SLiDE.Match`](@ref): Extract values from the specified column into a column or
        columns based on the specified regular expression.
    - [`SLiDE.Melt`](@ref): Normalize the dataframe by 'melting' columns into rows,
        lengthening the dataframe by duplicating values in the column `on` into new rows and
        defining 2 new columns:
        1. `var` with header names from the original dataframe.
        2. `val` with column values from the original dataframe.
    - [`SLiDE.Add`](@ref): Add new column `col` filled with `val`.
    - [`SLiDE.Map`](@ref): Define an `output` column containing values based on those in an
        `input` column. The mapping columns `from` -> `to` are contained in a .csv `file` in
        the coremaps directory. The columns `input` and `from` should contain the same
        values, as should `output` and `to`.
    - [`SLiDE.Replace`](@ref): Replace values in `col` `from` -> `to`.
    - [`SLiDE.Operate`](@ref): Perform an arithmetic operation across multiple DataFrame columns or rows.
    - [`SLiDE.Describe`](@ref): This DataType is required when multiple DataFrames will be
        appended into one output file (say, if multiple sheets from an XLSX file are
        included). Before the DataFrames are appended, a column `col` will be added and
        filled with the value in the file descriptor.
    - [`SLiDE.Order`](@ref): Rearranges columns in the order specified by `cols` and sets
        them to the specified type.
- `file::T where T <: File`: Data file containing information to read.
- `files::Array{T} where T <: File`: List of data files.
- `y::Dict{Any,Any}`: Dictionary containing all editing structures among other values read
    from the yaml file. Dictionary keys must correspond EXACTLY with SLiDE.Edit DataType
    names, or the edits will not be made.

# Keywords
- `shorten::Bool = false` or `shorten::Int`: if an integer length is specified, the
    DataFrame will be shortened to the input value. This is meant to aid troubleshooting
    during development.

# Returns
- `df::DataFrame` including edit(s)
"""
function edit_with(df::DataFrame, x::Add)
    df = copy(df)
    # If adding the length of a string...
    if typeof(x.val) == String && occursin("length", x.val)
        m = match(r"(?<col>\S*) length", x.val)

        # If this is not indicating a column length to add, add the value and exit.
        if (m === nothing || !(Symbol(m[:col]) in propertynames(df)))
            df[!, x.col] .= x.val
            return df
        end
        # If possible, return the length of characters  in the string.
        col_len = Symbol(m[:col])
        df[!, x.col] .= [ismissing(val_len) ? missing : length(convert_type(String, val_len))
            for val_len in df[:,col_len]]
    else
        df[!, x.col] .= x.val
    end
    return df
end

function edit_with(df::DataFrame, x::Drop)
    df = copy(df)
    if x.val === "all" && x.operation == "occursin"
        df = edit_with(df, Drop.(propertynames(df)[occursin.(x.col, propertynames(df))], "all", "=="))
    end

    !(x.col in propertynames(df)) && (return df)
    if x.val === "all"  # Drop entire column to remove dead weight right away.
        df = df[:, setdiff(propertynames(df), [x.col])]
    else  # Drop rows using an operation or based on a value.
        if x.val === missing
            dropmissing!(df, x.col)
        # elseif x.val === "unique"
        #     unique!(df, x.col)
        else
            df[!,x.col] .= convert_type.(typeof(x.val), df[:,x.col])
            df = if x.operation == "occursin"
                df[.!broadcast(datatype(x.operation), x.val, df[:,x.col]), :]
            else
                df[.!broadcast(datatype(x.operation), df[:,x.col], x.val), :]
            end
        end
    end
    return df
end

function edit_with(df::DataFrame, x::Group)
    df = copy(df)
    # First, add a column to the original DataFrame indicating where the data set begins.
    cols = unique([propertynames(df); x.output])
    df[!,:start] = (1:size(df)[1]) .+ 1

    # # Next, create a DataFrame describing where to "split" the input DataFrame.
    # # Editing with a map will remove all rows that do not contain relevant information.
    # # Add a column indicating where each data set STOPS, assuming all completely blank rows
    # # were removed by read_file().
    df_split = edit_with(copy(df), convert_type(Map, x); kind = :inner)
    sort!(unique!(df_split), :start)
    df_split[!, :stop] .= vcat(df_split[2:end, :start] .- 2, [size(df)[1]])

    # Add a new, blank output column to store identifying information about the data block.
    # Then, fill this column based on the identifying row numbers in df_split.
    for out in x.output
        df[!,out] .= ""
        [df[row[:start]:row[:stop], out] .= row[out] for row in eachrow(df_split)]
    end

    # Finally, remove header rows (these will be blank in the output column),
    # as well as the column describing where the sub-DataFrames begin.
    df = edit_with(df, Drop.(x.output, "", "=="))
    return df[:, cols]
end

function edit_with(df::DataFrame, x::Map; kind = :left)
    df = copy(df)
    # Save all input column propertynames, read the map file, and isolate relevant columns.
    # # This prevents duplicate columns in the final DataFrame.
    cols = unique([propertynames(df); x.output])
    df_map = copy(read_file(x))
    df_map = unique(df_map[:,unique([x.from; x.to])])

    # If there are duplicate columns in from/to, differentiate between the two to save results.
    duplicates = intersect(x.from, x.to)
    if length(duplicates) > 0
        (ii_from, ii_to) = (occursin.(duplicates, x.from), occursin.(duplicates, x.to));
        x.from[ii_from] = Symbol.(x.from[ii_from], :_0)
        [df_map[!,Symbol(col, :_0)] .= df_map[:,col] for col in duplicates]
    end

    # Rename columns in the mapping DataFrame to temporary values in case any of these
    # columns were already present in the input DataFrame.
    temp_to = Symbol.(:to_, 1:length(x.to))
    temp_from = Symbol.(:from_, 1:length(x.from))
    df_map = edit_with(df_map, Rename.([x.to; x.from], [temp_to; temp_from]))

    # Ensure the input and mapping DataFrames are consistent in type. Types from the mapping
    # DataFrame are used since all values in each column should be of the same type.
    for (col, col_map) in zip(x.input, temp_from)
        try
            new_type = eltypes(dropmissing(df_map[:,[col_map]]))
            # new_type = eltype.(eachcol(dropmissing(df_map[:,[col_map]])))
            df[!,col] .= convert_type.(new_type, df[:,col])
        catch
            df_map[!,col_map] .= convert_type.(String, df_map[:,col_map])
        end
    end

    join_cols = Pair.(x.input, temp_from)

    x.kind == :inner && (df = innerjoin(df, df_map, on = join_cols; makeunique = true))
    x.kind == :outer && (df = outerjoin(df, df_map, on = join_cols; makeunique = true))
    x.kind == :left  && (df = leftjoin(df, df_map, on = join_cols; makeunique = true))
    x.kind == :right && (df = rightjoin(df, df_map, on = join_cols; makeunique = true))
    x.kind == :semi  && (df = semijoin(df, df_map, on = join_cols; makeunique = true))

    # Remove all output column propertynames that might already be in the DataFrame. These will be
    # overwritten by the columns from the mapping DataFrame. Finally, remane mapping "to"
    # columns from their temporary to output values.
    df = df[:, setdiff(propertynames(df), x.output)]
    df = edit_with(df, Rename.(temp_to, x.output))
    return df[:,cols]
end

function edit_with(df::DataFrame, x::Match)
    df = copy(df)
    if x.on == r"expand range"
        ROWS, COLS = size(df)
        cols = propertynames(df)
        df = [[DataFrame(Dict(cols[jj] =>
                cols[jj] == x.input ? _expand_range(df[ii,jj]) : df[ii,jj]
            for jj in 1:COLS)) for ii in 1:ROWS]...;]
    else
        # Ensure all row values are strings and can be matched with a Regex, and do so.
        # Temporarily remove missing values, just in case.
        df[:,x.input] .= convert_type.(String, df[:,x.input])
        col = edit_with(copy(df), Replace(x.input, missing, ""))[:,x.input]
        m = match.(x.on, col)

        # Add empty columns for all output columns not already in the DataFrame.
        # Where there is a match, fill empty cells. If values in the input column,
        # leave cells without a match unchanged.
        df = edit_with(df, Add.(setdiff(x.output, propertynames(df)), ""))
        [m[ii] != nothing && ([df[ii,out] = m[ii][out] for out in x.output])
            for ii in 1:length(m)]
    end
    return df
end

function edit_with(df::DataFrame, x::Melt)
    df = copy(df)
    on = intersect(x.on, propertynames(df))
    df = melt(df, on, variable_name = x.var, value_name = x.val)
    df[!, x.var] .= convert_type.(String, df[:, x.var])
    return df
end

function edit_with(df::DataFrame, x::Operate)
    df = copy(df)
    # If it is a ROW-WISE operation,
    if x.axis == :row
        df = by(df, x.input, x.output => datatype(x.operation))
        df = edit_with(df, Rename.(setdiff(propertynames(df), x.input), ensurearray(x.output)))
    end

    # If it is a COLUMN-WISE operation,
    if x.axis == :col
        cols = [setdiff(propertynames(df), unique([x.from; x.to; x.input; x.output])); x.output; x.from]

        # Isolate columns to be operated on.
        # Append original columns that might be replaced "_0" to preserve information.
        df_val = convert_type.(Float64, copy(df[:,x.input]))
        x.output in x.input && (df = edit_with(df, Rename(x.output, Symbol(x.output, :_0))))
        df[!,x.output] .= broadcast(datatype(x.operation), [col for col in eachcol(df_val)]...)

        # Adjust labeling columns: If both from/to descriptive columns are distinct and
        # in the DataFrame, Replace the column values from -> to.
        for (from, to) in zip(x.from, x.to)
            if length(intersect(propertynames(df), [from,to])) == 2
                df_comment = dropmissing(unique(df[:, [from; to]]))
                df[!, Symbol(from, :_0)] .= df[:,from]
                df = edit_with(df, Replace.(from, df_comment[:,from], df_comment[:,to]))
            end
        end
    end
    # !!!! How to handle floating point arithmetic? (ex: 1.1 + 0.1 = 1.2000000000000002)
    df[!,x.output] .= round.(df[:,x.output], digits=11)
    return df
end

function edit_with(df::DataFrame, x::Order)
    df = copy(df)
    # If not all columns are present, return the DataFrame as is. Such is the case when a
    # descriptor column must be added when appending multiple data sets in one DataFrame.
    if size(intersect(x.col, propertynames(df)))[1] < size(x.col)[1]
        return df
    # If all of the columns are present in the original DataFrame,
    # reorder the DataFrame columns and set them to the specified type.
    else
        df = df[!, x.col]  # reorder
        [df[!, c] .= convert_type.(t, df[!, c]) for (c, t) in zip(x.col, x.type)]  # convert
        return df
    end
end

function edit_with(df::DataFrame, x::Rename)
    df = copy(df)
    x.from in propertynames(df) && (df = rename(df, x.from => x.to))
    x.to == :upper && (df = edit_with(df, Rename.(propertynames(df), uppercase.(propertynames(df)))))
    x.to == :lower && (df = edit_with(df, Rename.(propertynames(df), lowercase.(propertynames(df)))))
    return df
end

function edit_with(df::DataFrame, x::Replace)
    df = copy(df)
    !(x.col in propertynames(df)) && (return df)

    if x.from === missing && Symbol(x.to) in propertynames(df)
        df[ismissing.(df[:,x.col]),x.col] .= df[ismissing.(df[:,x.col]), Symbol(x.to)]
        return df
    end

    if x.to === Not && eltype(df[:,x.col]) == Bool
        df[!,x.col] .= .!df[:,x.col]
    end

    df[!,x.col] .= if x.to === "lower"  lowercase.(df[:,x.col])
    elseif x.to === "upper"             uppercase.(df[:,x.col])
    elseif x.to === "uppercasefirst"    uppercasefirst.(lowercase.(df[:,x.col]))
    elseif x.to === "titlecase"         titlecase.(df[:,x.col])
    else
        replace(strip.(copy(df[:,x.col])), x.from => x.to)
    end
    return df
end

function edit_with(df::DataFrame, x::Stack)
    df = copy(df)
    df = [[edit_with(df[:, occursin.(indicator, propertynames(df))],
        [Rename.(propertynames(df)[occursin.(indicator, propertynames(df))], x.col);
            Add(x.var, replace(string(indicator), "_" => " "))]
    ) for indicator in x.on]...;]
    return dropmissing(df)
end

function edit_with(df::DataFrame, lst::Array{T}) where T<:Edit
    [df = edit_with(df, x) for x in lst]
    return df
end

function edit_with(df::DataFrame, x::Describe, file::T) where T<:File
    return edit_with(copy(df), Add(x.col, file.descriptor))
end

function edit_with(file::T, y::Dict{Any,Any}; shorten = false) where T<:File
    df = read_file(y["PathIn"], file; shorten = shorten)
    # Specify the order in which edits must occur. "Drop" is included twice, once at the
    # beginning and once at the end. First, drop entire columns. Last, drop specific values.
    EDITS = ["Rename", "Group", "Stack", "Match", "Melt", "Add", "Map", "Replace", "Drop", "Operate"]

    # Find which of thyese edits are represented in the yaml file of defined edits.
    KEYS = intersect(EDITS, collect(keys(y)))
    "Drop" in KEYS && pushfirst!(KEYS, "Drop")

    [df = edit_with(df, y[k]) for k in KEYS]

    # Add a descriptor to identify the data from the file that was just added.
    # Then, reorder the columns and set them to the correct types.
    # This ensures consistency when concattenating.
    "Describe" in keys(y) && (df = edit_with(df, y["Describe"], file))
    "Order" in keys(y)    && (df = edit_with(df, y["Order"]))
    return df
end

function edit_with(files::Array{T}, y::Dict{Any,Any}; shorten = false) where T<:File
    return [[edit_with(file, y; shorten = shorten) for file in files]...;]
end

function edit_with(y::Dict{Any,Any}; shorten = false)
    # Find all dictionary keys corresponding to file names and save these in a list.
    file = convert_type(Array, find_oftype(y, File))
    df = edit_with(file, y; shorten = shorten)
    # return _sort_datastream(df)
end

"""
    _sort_datastream(df::DataFrame)
Returns the edited DataFrame, stored in a nicely-sorted order. This is most helpful for
mapping and developing. Sorting isn't *necessary* and we could remove this function to save
some time for users.
"""
function _sort_datastream(df::DataFrame)
    colidx = 1:size(df,2)
    isvalue = istype(df, AbstractFloat)
    ii = colidx[.!isvalue]

    # If it's a mapping dataframe...s
    if length(ii) == length(setdiff(propertynames(df),[:factor]))
        :state_code in propertynames(df) && (ii = intersect(colidx[occursin.(:code, propertynames(df))], ii))
        ii = intersect(sortperm(length.(unique.(eachcol(df)))), ii)
        splice!(ii, 2:1, colidx[isvalue])
    end
    return sort(df, ii)
end

"""
    _expand_range()
"""
function _expand_range(x::T) where T <: AbstractString
    if occursin("-", x)
        if all(string(strip(x)) .!= ["31-33", "44-45", "48-49"])
            x = split(x, "-")
            x = ensurearray(convert_type(Int, x[1]):convert_type(Int, x[1][1:end-1] * x[end][end]))
        end
    else
        x = convert_type(Int, x)
    end
    return x
end

function _expand_range(x::String)
    if match(r"\D", x) !== nothing
        m = String.([m.match for m in collect(eachmatch(r"\D.*?", x))])
        # length(setdiff(m, [",","-"," "])) .== 0 && (x = [_expand_range.(split(x, ","))...;])
        x = length(setdiff(m, [",","-"," "])) .== 0 ? [_expand_range.(split(x, ","))...;] : missing
    else
        x = convert_type(Int, x)
    end
    return x
end

_expand_range(x::Missing) = x
_expand_range(x::Int) = x

"""
    fill_zero(keys_unique::NamedTuple; value_colnames)
    fill_zero(keys_unique::NamedTuple, df::DataFrame)
    fill_zero(df::DataFrame...)
    fill_zero(d::Dict...)
    fill_zero(keys_unique, d::Dict)

# Arguments
- `keys_unique::Tuple`: A list of arrays whose permutations should be included in the
    resultant dictionary.
- `keys_unique::NamedTuple`: A list of arrays whose permutations should be included in the
    resultant dictionary. The NamedTuple's keys correspond to the DataFrame columns where
    they will be stored.
- `d::Dict...`: The dictionary/ies to edit.
- `df::DataFrame...`: The DataFrame(s) to edit.

# Keyword Arguments
- `value_colnames::Any = :value`: "value" column labels to add and set to zero when creating
    a new DataFrame. Default is `:value`.

# Usage
This function can be used to fill zeros in either a dictionary or DataFrame.
- Options for DataFrame editing:
    - If only (a) dictionary/ies is/are input, the dictionaries will be edited such that
        they all contain all permutations of their key values. All dictionaries in a
        resultant list of dictionaries will be the same length.
    - If a dictionary is input with a list of keys, it will be edited to ensure that it
        includes all permutations.
    - If only a list of keys is input, a new dictionary will be created, containing all key
        permutations with values initialized to zero.
- Options for DataFrame editing:
    - If only (a) DataFrame(s) is/are input, the DataFrame(s) will be edited such that
        they all contain all permutations of their key values. All DataFrames in a
        resultant list of DataFrames will be the same length.
    - If a DataFrame is input with a NamedTuple, it will be edited to ensure that it
        includes all permutations of the NamedTuple's values.
    - If only a NamedTuple is input, a new DataFrame will be created, containing all key
        permutations with values initialized to zero.

# Returns
- `d::Dict...` if input included dictionaries and/or Tuples
- `df::DataFrame...` if input included DataFrames and/or NamedTuples
"""
function fill_zero(keys_fill::NamedTuple; value_colnames = :value)
    df_fill = DataFrame(permute(keys_fill))
    return edit_with(df_fill, Add.(convert_type.(Symbol, value_colnames), 0.))
end

function fill_zero(keys_fill::Tuple; permute_keys::Bool = true)
    permute_keys && (keys_fill = permute(keys_fill))
    return Dict(k => 0. for k in keys_fill)
end

function fill_zero(d::Vararg{Dict}; permute_keys::Bool = true)
    d = copy.(ensurearray(d))
    # Find all keys present in the input dictionary/ies and ensure all are present.
    keys_fill = unique([collect.(keys.(d))...;])
    d = [fill_zero(keys_fill, x; permute_keys = permute_keys) for x in d]
    return length(d) == 1 ? d[1] : Tuple(d)
end

function fill_zero(keys_fill::NamedTuple, df::DataFrame)
    df = copy(df)
    df_fill = fill_zero(keys_fill)
    df = fill_zero(df, df_fill)[1]
    return df
end

function fill_zero(df::Vararg{DataFrame}; permute_keys::Bool = true)
    df = copy.(ensurearray(df))
    # Save propertynames of columns containing values to fill zeros later.
    # Find descriptor columns to permute OR make consistent across input DataFrames.
    value_colnames = find_oftype.(df, AbstractFloat)
    cols = intersect(setdiff.(propertynames.(df), value_colnames)...)

    # Find a unique list of descriptor keys in the input DataFrame(s). Permute as desired.
    df_fill = sort(unique([[x[:,cols] for x in df]...;]))
    permute_keys && (df_fill = permute(df_fill))

    # For each DataFrame in the list, join the input DataFrame to DataFrame keys_all on the
    # descriptor columns shared by both DataFrames. Using a left join will add "missing"
    # where a descriptor was not already present, which will be replaced by zero.
    [df[ii] = edit_with(leftjoin(df_fill, df[ii], on = cols),
        Replace.(value_colnames[ii], missing, 0.0)) for ii in 1:length(df)]
    return length(df) == 1 ? df[1] : Tuple(df)
end

function fill_zero(keys_fill::Any, d::Dict; permute_keys::Bool = true)
    d = copy(d)
    # If permuting keys, find all possible permutations of keys that should be present
    # and determine which are missing. Then add missing keys to the dictionary and return.
    permute_keys && (keys_fill = permute(keys_fill))
    keys_missing = setdiff(keys_fill, collect(keys(d)))

    [push!(d, k => 0.) for k in keys_missing]
    return d
end

function fill_with(inp::Any, val::Any)
    df = fill_zero(inp);
    df = edit_with(df, Replace(:value, 0.0, val))
    return df
end

"""
    extrapolate_year(df::DataFrame, yr::Array{Int64,1}; kwargs...)
    extrapolate_year(df::DataFrame, set::Dict; kwargs...)

# Arguments:
- `df::DataFrame` that might be in need of extrapolation.
- `yr::Array{Int64,1}`: List of years overwhich extrapolation is possible (depending on the kwargs)
- `set::Dict` containing list of years, identified by the key `:yr`.

# Keyword Arguments:
- `backward::Bool = true`: Do we extrapolate backward in time?
- `forward::Bool = true`: Do we extrapolate forward in time?

# Returns:
- `df::DataFrame` extrapolated in time.
"""
function extrapolate_year(
    df::DataFrame,
    yr::Array{Int64,1};
    backward::Bool = true,
    forward::Bool = true
)
    df = copy(df)
    yr_diff = setdiff(yr, unique(df[:,:yr]))
    length(yr_diff) == 0 && (return df)

    cols = setdiff(propertynames(df), [:yr])
    cols_ans = propertynames(df)

    df_ext = []

    if backward
        yr_min = minimum(df[:,:yr])
        df_min = filter_with(df, (yr = yr_min,))[:,cols]

        yr_back = yr_diff[yr_diff .< yr_min]
        df_back = crossjoin(DataFrame(yr = yr_back), df_min)[:,cols_ans]

        push!(df_ext, df_back)
    end

    if forward
        yr_max = maximum(df[:,:yr])
        df_max = filter_with(df, (yr = yr_max,))[:,cols]

        yr_forward = yr_diff[yr_diff .> yr_max]
        df_forward = crossjoin(DataFrame(yr = yr_forward), df_max)[:,cols_ans]

        push!(df_ext, df_forward)
    end

    return sort([df_ext...; df])
end

function extrapolate_year(
    df::DataFrame,
    set::Dict;
    backward::Bool = true,
    forward::Bool = true
)
    extrapolate_year(df, set[:yr]; forward = forward, backward = backward)
end

"""
    extrapolate_region(df::DataFrame, yr::Array{Int64,1}; kwargs...)
    extrapolate_year(df::DataFrame, set::Dict; kwargs...)

# Arguments:
- `df::DataFrame` that might be in need of extrapolation.
- `r::Pair = "md" => "dc"`: `Pair` indicating a region (`r.first`) to extrapolate to another
    region (`r.second`). A suggested regional extrapolation: MD data will be used to
    approximate DC data in the event that it is missing.

# Keyword Argument:
- `overwrite::Bool = false`: If data in the target region `r.second` is already present,
    should it be overwritten?

# Returns:
- `df::DataFrame` extrapolated in time.
"""
function extrapolate_region(df::DataFrame, r::Pair = "md" => "dc"; overwrite = false)
    df = copy(df)
    if !overwrite
        r = r.first => setdiff(ensurearray(r.second), unique(df[:,:r]))
        length(r.second) == 0 && (return df)
    else
        df = edit_with(df, Drop.(:r, r.second, "=="))
    end

    cols = setdiff(propertynames(df), [:r])
    df_close = crossjoin(DataFrame(r = r.second), filter_with(df, (r = r.first,))[:,cols])

    return sort([df_close; df])
end

"""
    filter_with(df::DataFrame, set::Any; kwargs...)

# Arguments:
- `df::DataFrame` to filter.
- `set::Dict` or `set::NamedTuple`: Values to keep in the DataFrame.

# Keyword Arguments:
- `extrapolate::Bool = false`: Add missing regions/years to the DataFrame?
    If `extrapolate` is set to true, the following `kwargs` become relevant:
    - When extrapolating over years,
        - `backward::Bool = true`: Do we extrapolate backward in time?
        - `forward::Bool = true`: Do we extrapolate forward in time?
    - When extrapolating across regions,
        - `r::Pair = "md" => "dc"`: `Pair` indicating a region (`r.first`) to extrapolate to
            another region (`r.second`). A suggested regional extrapolation: MD data will be
            used to approximate DC data in the event that it is missing.
        - `overwrite::Bool = false`: If data in the target region `r.second` is already present,
            should it be overwritten?

# Returns:
- `df::DataFrame` with only the desired keys.
"""
function filter_with(
    df::DataFrame,
    set::Any;
    extrapolate::Bool = false,
    forward::Bool = true,
    backward::Bool = true,
    r::Pair = "md" => "dc",
    overwrite::Bool = false
)
    df_ans = copy(df)
    cols_ans = propertynames(df_ans)

    # Find keys that reference both column names in the input DataFrame df and
    # values in the set Dictionary. Then, created a DataFrame containing all permutations.
    cols = find_oftype(df, Not(AbstractFloat))
    cols_set = intersect(cols, collect(keys(set)))
    vals_set = [intersect(unique(df[:,k]), ensurearray(set[k])) for k in cols_set]
    # vals_set = [set[k] for k in cols_set]

    # Drop values that are not in the current set.
    df_set = DataFrame(permute(NamedTuple{Tuple(cols_set,)}(vals_set,)))
    df_ans = innerjoin(df_ans, df_set, on = cols_set)

    if extrapolate
        :yr in cols_set && (df_ans = extrapolate_year(df_ans, set; forward = forward, backward = backward))
        :r in cols_set  && (df_ans = extrapolate_region(df_ans, r; overwrite = overwrite))
    end

    return sort(df_ans[:,cols_ans])
end

max_stars_repo_name: matbesancon/SLiDE

max_issues_repo_name: matbesancon/SLiDE

max_forks_repo_name: matbesancon/SLiDE


text:
#
# Circular Restricted Three-body Problem Systems
#


text:
"""
Stuct that defines the pricing MDP problem
"""
struct PMDPProblem{nₚ, nᵣ, nₐ}
   P::SVector{nₚ, Product}      # Products array
   c₀::SVector{nᵣ, Int64}       # Initial capacity of resources
   T::Timestep                  # Number of timesteps
   D::DiscreteCountingProcess   # Demand for each product
   B::SVector{nₚ, Distribution} # Budget distributions for each product
   A::SVector{nₐ, Action}       # Action array (including reject action)
   info::NamedTuple
   objective::Symbol            # objective for optimization

   function PMDPProblem(non_empty_P::AbstractArray{<:Product}, c₀, D, B, real_actions, objective; info=(;))
      T = problem_selling_horizon(non_empty_P)
      A = SVector(real_actions..., REJECT_ACTION)
      @assert objective in [:revenue, :utilization]
      @assert are_unique(non_empty_P)
      new{length(non_empty_P), length(c₀), length(A)}(non_empty_P,c₀,T,D,B,A,info,objective)
   end
end

"""Read number of resources and products from type parameters"""
Base.size(pp::PMDPProblem{nₚ, nᵣ, nₐ}) where {nₚ, nᵣ, nₐ} = (nₚ, nᵣ, nₐ)
n_products(pp::PMDPProblem) = size(pp)[1]
n_resources(pp::PMDPProblem) = size(pp)[2]
n_actions(pp::PMDPProblem) = size(pp)[3]

function problem_selling_horizon(P::AbstractArray{<:Product})::Timestep
   maximum([p.selling_period_end for p in P])
end

selling_period_end(p::PMDPProblem) = p.T

empty_product(p::PMDPProblem) = Product(falses(n_resources(p)), selling_period_end(p))


text:
module LichessBot

using Chess
using Chess.Book
using HTTP
using JSON
using Random
using Dates
using DelimitedFiles
import Base.isready, Base.isempty
@debug "LichessBot: import finished"
baseUrl = "https://lichess.org/api"
token = get(ENV, "LICHESS_TOKEN", "")
defaultHeader = Dict("Authorization" => "Bearer $token")
myId =
    isempty(token) ? "" :
    try
        JSON.parse(String(HTTP.get(baseUrl * "/account", defaultHeader).body))["id"]
    catch c
        ""
    end
mutable struct engineGame{G<:Union{Game,SimpleGame}}
    g::G
    wtime::Millisecond
    btime::Millisecond
    winc::Millisecond
    binc::Millisecond
    #movestogo::Integer
    depth::Integer
    engineGame{G}(g::G) where {G<:Union{Game,SimpleGame}} =
        new(g, Minute(1), Minute(1), Millisecond(0), Millisecond(0), 1) #=0,=#
end # mutable struct
include("moveEvaluation.jl")
include("moveSelection.jl")
include("playGames.jl")
include("pst.jl")
include("utils.jl")
include("uci.jl")
export runBot, setDefaultToken!, upgradeToBot, negamax

end


text:
export triangulation, triangulation_indices

# Implementation of Cohen & Hickey algorithm for triangulation; see [CH79] and [BEK00].
#
# [CH79] Cohen, J., & Hickey, T. (1979). Two algorithms for determining volumes of convex polyhedra. Journal of the ACM (JACM), 26(3), 401-414.
# [BEK00] Büeler, B., Enge, A., & Fukuda, K. (2000). Exact volume computation for polytopes: a practical study. In Polytopes—combinatorics and computation (pp. 131-154). Birkhäuser, Basel.

# TODO parallelize: this is inherently parallelizable
function _triangulation(Δs, Δ, v_idx, h_idx, incident_idx, is_weak_adjacent)
    isempty(v_idx) && return
    v = first(v_idx)
    Δ = push!(copy(Δ), v)
    if isone(length(v_idx))
        push!(Δs, Δ)
        return
    end
    tail = true
    for h in h_idx
        if !(v in incident_idx[h])
            tail = false
            # The adjacency may be outside the current face under scrutiny but that's ok,
            # it will simply end up calling `_triangulation` with an empty `v_idx`.
            weak_adjacent = [hj for hj in h_idx if hj != h && is_weak_adjacent[(h, hj)]]
            active = [point for point in v_idx if point in incident_idx[h]]
            _triangulation(Δs, Δ, active, weak_adjacent, incident_idx, is_weak_adjacent)
        end
    end
end
function triangulation_indices(p::Polyhedron)
    hasrays(p) && error("Triangulation only supported for polytope.")
    v_idx = eachindex(points(p))
    h_idx = eachindex(halfspaces(p))
    Δ = eltype(v_idx)[]
    Δs = typeof(Δ)[]
    incident_idx = Dict(h => Set(incidentpointindices(p, h)) for h in h_idx)
    is_weak_adjacent = Dict((hi, hj) => !isempty(incident_idx[hi] ∩ incident_idx[hj]) for hi in h_idx for hj in h_idx)
    _triangulation(Δs, Δ, v_idx, h_idx, incident_idx, is_weak_adjacent)
    return Δs
end
function triangulation(p::Polyhedron)
    return map(Δ -> vrep(get.(p, Δ)), triangulation_indices(p))
end


text:
"""
    Interval(first, last, [inclusivity::Inclusivity]) -> Interval
    Interval(first, last, [closed_left::Bool, closed_right::Bool]) -> Interval

An `Interval` represents a non-iterable range or span of values (non-interable because,
unlike a `StepRange`, no step is defined).

An `Interval` can be closed (both `first` and `last` are included in the interval), open
(neither `first` nor `last` are included), or half-open. This openness is defined by an
`Inclusivity` value, which defaults to closed.

### Example

```julia
julia> i = Interval(0, 100, true, false)
Interval{Int64}(0, 100, Inclusivity(true, false))

julia> in(0, i)
true

julia> in(50, i)
true

julia> in(100, i)
false

julia> intersect(Interval(0, 25, false, false), Interval(20, 50, true, true)
Interval{Int64}(20, 25, Inclusivity(true, false))
```

### Infix Constructor: `..`

A closed `Interval` can be constructed with the `..` infix constructor:

```julia
julia> Dates.today() - Dates.Week(1) .. Dates.today()
Interval{Date}(2018-01-24, 2018-01-31, Inclusivity(true, true))
```

### Note on Ordering

The `Interval` constructor will compare `first` and `last`; if it finds that
`first > last`, they will be reversed to ensure that `first < last`. This simplifies
calls to `in` and `intersect`:

```julia
julia> i = Interval(Date(2016, 8, 11), Date(2013, 2, 13), false, true)
Interval{Date}(2013-02-13, 2016-08-11, Inclusivity(true, false))
```

Note that the `Inclusivity` value is also reversed in this case.

See also: [`AnchoredInterval`](@ref), [`Inclusivity`](@ref)
"""
struct Interval{T} <: AbstractInterval{T}
    first::T
    last::T
    inclusivity::Inclusivity

    function Interval{T}(f::T, l::T, inc::Inclusivity) where T
        # Ensure that `first` preceeds `last`.
        f, l, inc = if f ≤ l
            f, l, inc
        elseif l ≤ f
            l, f, Inclusivity(last(inc), first(inc))
        else
            throw(ArgumentError("Unable to determine an ordering between: $f and $l"))
        end

        return new(f, l, inc)
    end
end

Interval{T}(f, l, inc::Inclusivity) where T = Interval{T}(convert(T, f), convert(T, l), inc)
Interval{T}(f, l, x::Bool, y::Bool) where T = Interval{T}(f, l, Inclusivity(x, y))
Interval{T}(f, l) where T = Interval{T}(f, l, true, true)

Interval(f::T, l::T, inc...) where T = Interval{T}(f, l, inc...)
Interval(f, l, inc...) = Interval(promote(f, l)..., inc...)

(..)(first, last) = Interval(first, last)

# In Julia 0.7 constructors no longer automatically fall back to using `convert`
Interval(interval::AbstractInterval) = convert(Interval, interval)
Interval{T}(interval::AbstractInterval) where T = convert(Interval{T}, interval)

# Endpoint constructors
function Interval{T}(left::LeftEndpoint{T}, right::RightEndpoint{T}) where T
    Interval{T}(left.endpoint, right.endpoint, left.included, right.included)
end

Interval(left::LeftEndpoint{T}, right::RightEndpoint{T}) where T = Interval{T}(left, right)

# Empty Intervals
Interval{T}() where T = Interval{T}(zero(T), zero(T), Inclusivity(false, false))
Interval{T}() where T <: TimeType = Interval{T}(T(0), T(0), Inclusivity(false, false))

function Interval{T}() where T <: ZonedDateTime
    return Interval{T}(T(0, tz"UTC"), T(0, tz"UTC"), Inclusivity(false, false))
end

Base.copy(x::Interval{T}) where T = Interval{T}(x.first, x.last, x.inclusivity)

function Base.hash(interval::AbstractInterval, h::UInt)
    h = hash(LeftEndpoint(interval), h)
    h = hash(RightEndpoint(interval), h)
    return h
end

##### ACCESSORS #####

Base.first(interval::Interval) = interval.first
Base.last(interval::Interval) = interval.last
span(interval::Interval) = interval.last - interval.first
inclusivity(interval::AbstractInterval) = interval.inclusivity
isclosed(interval::AbstractInterval) = isclosed(inclusivity(interval))
Base.isopen(interval::AbstractInterval) = isopen(inclusivity(interval))

##### CONVERSION #####

function Base.convert(::Type{T}, i::Interval{T}) where T
    first(i) == last(i) && isclosed(i) && return first(i)
    throw(DomainError(i, "The interval is not closed with coinciding endpoints"))
end

# Date/DateTime attempt to convert to Int64 instead of falling back to convert(T, ...)
Compat.Dates.Date(interval::Interval{Date}) = convert(Date, interval)
Compat.Dates.DateTime(interval::Interval{DateTime}) = convert(DateTime, interval)

##### DISPLAY #####

function Base.show(io::IO, interval::T) where T <: Interval
    if get(io, :compact, false)
        print(io, interval)
    else
        print(io, "$T(")
        show(io, interval.first)
        print(io, ", ")
        show(io, interval.last)
        print(io, ", ")
        show(io, interval.inclusivity)
        print(io, ")")
    end
end

function Base.print(io::IO, interval::AbstractInterval)
    # Print to io in order to keep properties like :limit and :compact
    if get(io, :compact, false)
        io = IOContext(io, :limit=>true)
    end

    print(
        io,
        first(inclusivity(interval)) ? "[" : "(",
        first(interval),
        " .. ",
        last(interval),
        last(inclusivity(interval)) ? "]" : ")",
    )
end

##### ARITHMETIC #####

Base.:+(a::T, b) where {T <: Interval} = T(first(a) + b, last(a) + b, inclusivity(a))

Base.:+(a, b::Interval) = b + a
Base.:-(a::Interval, b) = a + -b
Base.:-(a, b::Interval) = a + -b

function Base.:-(a::Interval{T}) where T
    inc = inclusivity(a)
    Interval{T}(-last(a), -first(a), Inclusivity(last(inc), first(inc)))
end

##### EQUALITY #####

function Base.:(==)(a::AbstractInterval, b::AbstractInterval)
    return LeftEndpoint(a) == LeftEndpoint(b) && RightEndpoint(a) == RightEndpoint(b)
end

function Base.isequal(a::AbstractInterval, b::AbstractInterval)
    le = isequal(LeftEndpoint(a), LeftEndpoint(b))
    re = isequal(RightEndpoint(a), RightEndpoint(b))
    return le && re
end

# While it might be convincingly argued that this should define < instead of isless (see
# https://github.com/invenia/Intervals.jl/issues/14), this breaks sort.
Base.isless(a::AbstractInterval, b) = LeftEndpoint(a) < b
Base.isless(a, b::AbstractInterval) = a < LeftEndpoint(b)

less_than_disjoint(a::AbstractInterval, b) = RightEndpoint(a) < b
less_than_disjoint(a, b::AbstractInterval) = a < LeftEndpoint(b)

function Base.:isless(a::AbstractInterval, b::AbstractInterval)
    return LeftEndpoint(a) < LeftEndpoint(b)
end

function less_than_disjoint(a::AbstractInterval, b::AbstractInterval)
    return RightEndpoint(a) < LeftEndpoint(b)
end

greater_than_disjoint(a, b) = less_than_disjoint(b, a)

"""
    ≪(a::AbstractInterval, b::AbstractInterval) -> Bool
    less_than_disjoint(a::AbstractInterval, b::AbstractInterval) -> Bool

Less-than-and-disjoint comparison operator. Returns `true` if `a` is less than `b` and they
are disjoint (they do not overlap).

```
julia> 0..10 ≪ 10..20
false

julia> 0..10 ≪ 11..20
true
```
"""
≪(a, b) = less_than_disjoint(a, b)
# ≪̸(a, b) = !≪(a, b)

"""
    ≫(a::AbstractInterval, b::AbstractInterval) -> Bool
    greater_than_disjoint(a::AbstractInterval, b::AbstractInterval) -> Bool

Greater-than-and-disjoint comparison operator. Returns `true` if `a` is greater than `b` and
they are disjoint (they do not overlap).

```
julia> 10..20 ≫ 0..10
false

julia> 11..20 ≫ 0..10
true
```
"""
≫(a, b) = greater_than_disjoint(a, b)
# ≫̸(a, b) = !≫(a, b)

##### SET OPERATIONS #####

Base.isempty(i::AbstractInterval) = LeftEndpoint(i) > RightEndpoint(i)
Base.in(a, b::AbstractInterval) = !(a ≫ b || a ≪ b)

function Base.in(a::AbstractInterval, b::AbstractInterval)
    # Intervals should be compared with set operations
    throw(ArgumentError("Intervals can not be compared with `in`. Use `issubset` instead."))
end

function Base.issubset(a::AbstractInterval, b::AbstractInterval)
    return LeftEndpoint(a) ≥ LeftEndpoint(b) && RightEndpoint(a) ≤ RightEndpoint(b)
end

Base.:⊈(a::AbstractInterval, b::AbstractInterval) = !issubset(a, b)
Base.:⊉(a::AbstractInterval, b::AbstractInterval) = !issubset(b, a)

function overlaps(a::AbstractInterval, b::AbstractInterval)
    left = max(LeftEndpoint(a), LeftEndpoint(b))
    right = min(RightEndpoint(a), RightEndpoint(b))

    return left <= right
end

function contiguous(a::AbstractInterval, b::AbstractInterval)
    left = max(LeftEndpoint(a), LeftEndpoint(b))
    right = min(RightEndpoint(a), RightEndpoint(b))

    return right.endpoint == left.endpoint && left.included != right.included
end

function Base.intersect(a::AbstractInterval{T}, b::AbstractInterval{T}) where T
    !overlaps(a,b) && return Interval{T}()
    left = max(LeftEndpoint(a), LeftEndpoint(b))
    right = min(RightEndpoint(a), RightEndpoint(b))

    return Interval{T}(left, right)
end

# There is power in a union.
"""
    union(intervals::AbstractVector{<:AbstractInterval})

Flattens a vector of overlapping intervals into a new, smaller vector containing only
non-overlapping intervals.
"""
function Base.union(intervals::AbstractVector{<:AbstractInterval})
    return union!(convert(Vector{AbstractInterval}, intervals))
end

"""
    union!(intervals::AbstractVector{<:Union{Interval, AbstractInterval}})

Flattens a vector of overlapping intervals in-place to be a smaller vector containing only
non-overlapping intervals.
"""
function Base.union!(intervals::Union{AbstractVector{<:Interval}, AbstractVector{AbstractInterval}})
    sort!(intervals)

    i = 2
    n = length(intervals)
    while i <= n
        prev = intervals[i - 1]
        curr = intervals[i]

        # If the current and previous intervals don't meet then move along
        if !overlaps(prev, curr) && !contiguous(prev, curr)
            i = i + 1

        # If the two intervals meet then we absorb the current interval into
        # the previous one.
        else
            intervals[i - 1] = merge(prev, curr)
            deleteat!(intervals, i)
            n -= 1
        end
    end

    return intervals
end

"""
    superset(intervals::AbstractArray{<:AbstractInterval}) -> Interval

Create the smallest single interval which encompasses all of the provided intervals.
"""
function superset(intervals::AbstractArray{<:AbstractInterval})
    left = minimum(LeftEndpoint.(intervals))
    right = maximum(RightEndpoint.(intervals))

    return Interval(left, right)
end

function Base.merge(a::AbstractInterval, b::AbstractInterval)
    if !overlaps(a, b) && !contiguous(a, b)
        throw(ArgumentError("$a and $b are neither overlapping or contiguous."))
    end

    left = min(LeftEndpoint(a), LeftEndpoint(b))
    right = max(RightEndpoint(a), RightEndpoint(b))
    return Interval(left, right)
end

##### TIME ZONES #####

function TimeZones.astimezone(i::Interval{ZonedDateTime}, tz::TimeZone)
    return Interval(astimezone(first(i), tz), astimezone(last(i), tz), inclusivity(i))
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

# -------------------
# VIEWS WITH INDICES
# -------------------

Base.view(domain::Domain, inds) = DomainView(domain, inds)
Base.view(data::Data, inds) = DataView(data, inds)

# specialize view to avoid infinite loops
Base.view(v::DataView, inds::AbstractVector{Int}) =
  DataView(getfield(v, :data), getfield(v, :inds)[inds])

# ----------------------
# VIEWS WITH GEOMETRIES
# ----------------------

"""
    view(domain, geometry)

Return a view of the `domain` containing all elements that
are inside the `geometry`.
"""
Base.view(domain::Domain, geometry::Geometry) =
  view(domain, viewindices(domain, geometry))

function Base.view(data::Data, geometry::Geometry)
  D   = typeof(data)
  dom = domain(data)
  tab = values(data)

  # retrieve subdomain
  inds   = viewindices(dom, geometry)
  subdom = view(dom, inds)

  # retrieve subtable
  tinds  = _linear(dom, inds)
  subtab = viewtable(tab, tinds)

  # data table for elements
  vals = Dict(paramdim(dom) => subtab)

  constructor(D)(subdom, vals)
end

# convert from Cartesian to linear indices if needed
@traitfn _linear(domain::D, inds) where {D; IsGrid{D}} =
  vec(LinearIndices(size(domain))[inds])
@traitfn _linear(domain::D, inds) where {D; !IsGrid{D}} = inds

"""
    viewindices(domain, geometry)

Return the indices of the `domain` that are inside the `geometry`.
"""
@traitfn function viewindices(domain::D, geometry::Geometry) where {D; !IsGrid{D}}
  pred(i) = _isinside(domain[i], geometry)
  filter(pred, 1:nelements(domain))
end

_isinside(p::Point, geometry) = p ∈ geometry
_isinside(g::Geometry, geometry) = g ⊆ geometry

@traitfn function viewindices(domain::D, box::Box) where {D; IsGrid{D}}
  # grid properties
  or = coordinates(minimum(domain))
  sp = spacing(domain)
  sz = size(domain)

  # intersection of boxes
  □ = boundingbox(domain) ∩ box
  lo, up = coordinates.(extrema(□))

  # Cartesian indices of new corners
  ilo = @. max(ceil(Int,  (lo - or) / sp) + 1,  1)
  iup = @. min(floor(Int, (up - or) / sp)    , sz)

  CartesianIndex(Tuple(ilo)):CartesianIndex(Tuple(iup))
end

# ----------
# UTILITIES
# ----------

"""
    slice(object, xmin:xmax, ymin:ymax, ...)

Slice the `object` using real coordinate ranges `xmin:xmax`, `ymin:ymax`, ...

### Notes

This function is equivalent to `view(object, Box(first.(ranges), last.(ranges))`.

In Julia the range `0.5:10.0` is materialized as `[0.5, ..., 9.5]` so it won't
necessarily include the right value. This behavior is different than the more
intuitive behavior of `view(object, Box((0.5,0.5), (10.0,10.0))`.
"""
slice(object, ranges...) = view(object, Box(first.(ranges), last.(ranges)))


text:
# Divergence calculation

"""
	divergence_ndgrid!(hx, hy, hz, vec, div)

Generic divergence `div` of vector `vec` with step lengths `hx`, `hy`, `hz`.
"""
function divergence_ndgrid!(hx, hy, hz, vec, div)

   siz = size(vec)::NTuple{4,Int64}

   px = zeros(Float64, Base.front(siz))
   qy = zeros(Float64, Base.front(siz))
   rz = zeros(Float64, Base.front(siz))

   n = size(hx, 1)
   # Right now do nothing for the ghost cells maybe needed later!
   # Take central differences on interior points
   if n > 2
      compute_px(siz, px, n, hx, vec)
   end

   n = size(hy, 1)
   if n > 2
      compute_qy(siz, qy, n, hx, vec)
   end

   n = size(hz, 1)
   if n > 2
      compute_rz(siz, rz, n, hx, vec)
   end

   @. div = px + qy + rz
   return
end

function divergence_ndgrid(hx, hy, hz, vec)

   siz = size(vec)

   div = zeros(Float64, Base.front(siz))

   divergence_ndgrid!(hx, hy, hz, vec, div)

   return div
end

function compute_px(siz, px, n, hx, vec)
   @inbounds for i = 1:n-2, j = 1:siz[2], k = 1:siz[3]
      px[i+1,j,k] = (vec[i+2,j,k,1] - vec[i,j,k,1])/(hx[i+2] - hx[i])
   end
end

function compute_qy(siz, qy, n, hy, vec)
   @inbounds for i = 1:siz[1], j = 1:n-2, k = 1:siz[3]
      qy[i,j+1,k] = (vec[i,j+2,k,2] - vec[i,j,k,2])/(hy[i+2] - hy[i])
   end
end

function compute_rz(siz, rz, n, hz, vec)
   @inbounds for i = 1:siz[1], j = 1:siz[2], k = 1:n-2
      rz[i,j,k+1] = (vec[i,j,k+2,3] - vec[i,j,k,3])/(hz[i+2] - hz[i])
   end
end

"""
   divergence!(param, vec, div)

Calculate the divergence of vectors specialized to my grid size.
Always assume starting with i -> j -> k for 1/2/3D!
Right now do nothing for the ghost cells. Maybe needed later!
Take central differences on interior points.
"""
@inline function divergence!(param, vec, div)

   x, y, z = param.x, param.y, param.z
   nI, nJ, nK, nG = param.nI, param.nJ, param.nK, param.nG
   iMin, iMax, jMin, jMax, kMin, kMax =
   param.iMin, param.iMax, param.jMin, param.jMax, param.kMin, param.kMax

   div .= 0.0

   @inbounds for k = kMin:kMax, j = jMin:jMax, i = iMin:iMax
      div[i-nG,j-nG,k-nG] = (vec[i+1,j,k,1] - vec[i-1,j,k,1]) / (x[i+1] - x[i-1])
   end

   if nJ > 1
      @inbounds for k = kMin:kMax, j = jMin:jMax, i = iMin:iMax
         div[i-nG,j-nG,k-nG] += (vec[i,j+1,k,2] - vec[i,j-1,k,2]) / (y[i+1] - y[i-1])
      end
   end

   if nK > 1
      @inbounds for k = kMin:kMax, j = jMin:jMax, i = iMin:iMax
         div[i-nG,j-nG,k-nG] += (vec[i,j,k+1,3] - vec[i,j,k-1,3]) / (z[i+1] - z[i-1])
      end
   end

   return
end


text:
using LinearAlgebra
using Random

Random.seed!(843)

include("feast.jl")
using ..feastNonlinear,..feastUtil
using IterativeSolvers

using IterativeSolvers

#eigenvalue problem dimension:
n=500

#eigenvector subspaces:
R=rand(n,n)
(X1,R)=qr(R)
R=rand(n,n)
(X2,R)=qr(R)

#quadratic eigenvalues:
L1=zeros(n,n)
L2=zeros(n,n)
#eigenvalues are 1,2,3,...,2n
for i in 1:n
    L1[i,i]=i
    L2[i,i]=i+n
end

#generate quadratic factor matrices:
B1=X1*L1*X1'
B2=X2*L2*X2'

#quadratic eigenvcalue problem matrices: T(z)=z^2*M+z*D+K
M=Matrix(I,n,n)
D=-B1-B2
K=B1*B2

#FEAST parameters
m0=5 #subspace dimension
x0=rand(ComplexF64,n,m0) #eigenvector initial guess
nc=4 #number of integration quadrature points
emid=50 #center point of contour
ra=1.5 #contour radius 1
rb=1.5 #contour radius 2
eps=1e-8 #convergence residual tolerance
maxit=50 #maximum number of FEAST iterations

#use quadratic FEAST
println("Using quadratic feast")
(lest,xest)=nlfeast_quad(M,D,K,x0,nc,emid,ra,rb,eps,maxit)

#Using Beyn
println("Using Beyn")
Tf(z)=(z^2*M+z*D+K)
(lest,xest)=beyn(Tf,x0,nc,emid,ra,rb)
ncmax=100
res=M*xest*diagm(0 => lest.^2)+D*xest*diagm(0 => lest)+K*xest
maxres=getMaxResInside(lest,xest,res,emid,ra,rb)
println("res=$maxres")

#Using general nonlinear feast, i.e. feast with beyn
println("Using general nlfeast")
Tf(z)=(z^2*M+z*D+K)
(lest, xest)=nlfeast_beyn(Tf,x0,nc,emid,ra,rb,eps,maxit)

#using iterative general nonlinear FEAST with GMRES
println("Using general iterative nlfeast")
alpha=0.0 #target accuracy for linear system solutions
maxitls = 50 # maximum linear system iterations
(lest, xest)=inlfeast_beyn(Tf,x0,alpha,maxitls,nc,emid,ra,rb,eps,maxit)

max_stars_repo_name: rleegates/feastjl

max_issues_repo_name: rleegates/feastjl

max_forks_repo_name: rleegates/feastjl


text:
using AdventOfCode2018
using AdventOfCode2018.Day8

input = get_aoc_input(8)
numbers = parse_input(input)
tree = build_tree(numbers)

println("Solution part 1: ", part_1(tree))
println("Solution part 2: ", part_2(tree))


text:
@testset "Heuristics" begin
    @testset "farmer" begin
        NS = 3  # number of scenarios
        probability = ones(3) / 3
        CROPS = 1:3 # set of crops (wheat, corn and sugar beets, resp.)
        PURCH = 1:2 # set of crops to purchase (wheat and corn, resp.)
        SELL = 1:4  # set of crops to sell (wheat, corn, sugar beets under 6K and those over 6K)
        Cost = [150 230 260]    # cost of planting crops
        Budget = 500            # budget capacity
        Purchase = [238 210]    # purchase price
        Sell = [170 150 36 10]  # selling price
        Yield = [3.0 3.6 24.0; 2.5 3.0 20.0; 2.0 2.4 16.0]
        Minreq = [200 240 0]    # minimum crop requirement

        @testset "MILP" begin
            # This creates a Lagrange dual problem for each scenario s.
            function create_scenario_model(s::Int64)
                m = Model(GLPK.Optimizer)
                @variable(m, 0 <= x[i=CROPS] <= 500, Int)
                @variable(m, y[j=PURCH] >= 0)
                @variable(m, w[k=SELL] >= 0)

                @objective(m, Min,
                    probability[s] * sum(Cost[i] * x[i] for i=CROPS)
                    + probability[s] * sum(Purchase[j] * y[j] for j=PURCH)
                    - probability[s] * sum(Sell[k] * w[k] for k=SELL))

                @constraint(m, sum(x[i] for i=CROPS) <= Budget)
                @constraint(m, [j=PURCH], Yield[s,j] * x[j] + y[j] - w[j] >= Minreq[j])
                @constraint(m, Yield[s,3] * x[3] - w[3] - w[4] >= Minreq[3])
                @constraint(m, w[3] <= 6000)
                return m
            end

            models = Dict{Int,JuMP.Model}(s => create_scenario_model(s) for s in 1:NS)
            coupling_variables = Vector{DD.CouplingVariableRef}()
            for s in 1:NS
                model = models[s]
                xref = model[:x]
                for i in CROPS
                    push!(coupling_variables, DD.CouplingVariableRef(s, i, xref[i]))
                end
            end

            @testset "All block heuristic with trust region method" begin
                # Create DualDecomposition instance.
                algo = DD.LagrangeDual()

                # Add Lagrange dual problem for each scenario s.
                for s in 1:NS
                    DD.add_block_model!(algo, s, models[s])
                end

                # Set nonanticipativity variables as an array of symbols.
                DD.set_coupling_variables!(algo, coupling_variables)

                # Lagrange master method
                LM = DD.BundleMaster(BM.TrustRegionMethod, GLPK.Optimizer)

                # add heuristic
                DD.add!(DD.AllBlockHeuristic, algo)

                # Solve the problem with the solver; this solver is for the underlying bundle method.
                DD.run!(algo, LM)

                @show DD.primal_objective_value(algo)
                @show DD.primal_solution(algo)
                @test isapprox(DD.primal_objective_value(algo), -107701.0, rtol=1e-3)
            end

            @testset "Rounding heuristic with trust region method" begin
                # Create DualDecomposition instance.
                algo = DD.LagrangeDual()

                # Add Lagrange dual problem for each scenario s.
                for s in 1:NS
                    DD.add_block_model!(algo, s, models[s])
                end

                # Set nonanticipativity variables as an array of symbols.
                DD.set_coupling_variables!(algo, coupling_variables)

                # Lagrange master method
                LM = DD.BundleMaster(BM.TrustRegionMethod, GLPK.Optimizer)

                # add heuristic
                DD.add!(DD.RoundingHeuristic, algo)

                # Solve the problem with the solver; this solver is for the underlying bundle method.
                DD.run!(algo, LM)

                @show DD.primal_objective_value(algo)
                @show DD.primal_solution(algo)
                @test isapprox(DD.primal_objective_value(algo), -107342., rtol=1e-3)
            end
        end
    end
end


text:
#SpiceData: Base types & core functions
#-------------------------------------------------------------------------------

#==Aliases
===============================================================================#
const DataWord = UInt32

#==Constants
===============================================================================#
const SIGNAME_BUFSIZE = 256 #Maximum supported length

const WRITEBLOCK_SYNCWORD = DataWord(0x4)
const WRITEBLOCK_HEADERSIZE = 4*sizeof(DataWord)

#TODO: Not convinced this code is really "block type"... or if ids match their meaning.
const BLOCKTYPEID_HEADER = DataWord(0x70)
const BLOCKTYPEID_DATA = DataWord(0x80)

const DATATYPEID_VOLTAGE = 1
const DATATYPEID_CURRENT = 8

const SWEEPDATA_LASTPOINT = 1e30 #Can detect end of data with this.

#==Main Types
===============================================================================#
abstract type Endianness end #Of data file
struct BigEndian <: Endianness; end
struct LittleEndian <: Endianness; end
const NetworkEndianness = BigEndian #Informative only
#=Comment:
Apparently SPICE files used to use network endianness (big-endian), but are now
little-endian.
=#

abstract type SpiceFormat end
struct Format_Unknown <: SpiceFormat; end
struct Format_9601 <: SpiceFormat; end #x: 32-bit floats, y: 32-bit floats
struct Format_9602 <: SpiceFormat; end #x: 64-bit floats, y: 64-bit floats
struct Format_2001 <: SpiceFormat; end #x: 64-bit floats, y: 64-bit floats
struct Format_2013 <: SpiceFormat; end #x: 64-bit floats, y: 32-bit floats
#=Comment:
I believe 9602 is a non-standard format used by CppSim.
=#

struct BlockHeader
	typeid::DataWord #TODO: is this really type id?
	_size::DataWord #Number of bits in current block
end

mutable struct BlockReader{T<:Endianness}
	io::IO
	header::BlockHeader
	endpos::Int
end
#NOTE: Parameterized so we can specialize (dispatch) on endianness.

#Convenience:
Endianness(::BlockReader{E}) where E<:Endianness = E

#TODO: Not sure if the SpiceTags are named correctly:
mutable struct SpiceTags
	id::String
	date::String
	time::String
	comments::String
end
SpiceTags() = SpiceTags("", "", "", "")

#SPICE file reader: Main object
#-------------------------------------------------------------------------------
mutable struct DataReader
	io::IOStream
	filepath::String #Informative only
	format::SpiceFormat
	sweepname::String
	signalnames::Vector{String}
	sweep::Vector
	tags::SpiceTags
	datastart::Int
	rowsize::Int
	endianness::Endianness
end

#==Exceptions
===============================================================================#
function corruptword_exception(io::IO, w::DataWord, expected::DataWord)
	pos = position(io) - sizeof(DataWord)
	pos = hex(pos)
	w = hex(w)
	expected = hex(expected)
	return "Corrupt word 0x$w @ 0x$pos (expected 0x$expected)"
end

function stringboundary_exception(io::IO, )
	hpos = hex(position(io))
	return "Reading string across block boundary: 0x$hpos"
end

#==Helper Functions
===============================================================================#

printable(v::String) = isprint(v) ? v : ""

#Debug: show block header info
function _show(io::IO, hdr::BlockHeader, pos::Int)
	print(io, "Block: 0x", hex(WRITEBLOCK_SYNCWORD))
	print(io, " 0x", hex(hdr.typeid))
	print(io, " 0x", hex(WRITEBLOCK_SYNCWORD))
	print(io, " 0x", hex(hdr._size))
	print(io, " (start 0x", hex(pos), ")")
	println(io)
end

#-------------------------------------------------------------------------------
xtype(::Format_9601) = Float32
ytype(::Format_9601) = Float32

xtype(::Format_9602) = Float64
ytype(::Format_9602) = Float64

xtype(::Format_2001) = Float64
ytype(::Format_2001) = Float64

xtype(::Format_2013) = Float64
ytype(::Format_2013) = Float32

#-------------------------------------------------------------------------------
_reorder(v::T, ::BigEndian) where T = ntoh(v)
_reorder(v::T, ::LittleEndian) where T = ltoh(v)

#IO reads
#-------------------------------------------------------------------------------
_read(io::IO, ::Type{T}, endianness::Endianness) where T<:Real =
	_reorder(read(io, T), endianness)

#Read in a WRITEBLOCK_SYNCWORD & validate:
function readsyncword(io::IO, endianness::Endianness)
	w = _read(io, DataWord, endianness)
	if w != WRITEBLOCK_SYNCWORD
		throw(corruptword_exception(io, w, WRITEBLOCK_SYNCWORD))
	end
end

#Read in a block header:
function _read(io::IO, ::Type{BlockHeader}, endianness::Endianness)
	readsyncword(io, endianness)
	typeid = _read(io, DataWord, endianness)
	readsyncword(io, endianness)
	_size = _read(io, DataWord, endianness)
	return BlockHeader(typeid, _size)
end

#Block reader
#-------------------------------------------------------------------------------
bytesleft(r::BlockReader) = (r.endpos - position(r.io))
canread(r::BlockReader, nbytes::Int) = bytesleft(r) >= nbytes

function nextblock(r::BlockReader{E}) where E
	seek(r.io, r.endpos)
	sz = _read(r.io, DataWord, E())
	if sz != r.header._size
		hpos = hex(position(r.io) - 1)
		throw("Inconsistent block size @ 0x$hpos.")
	end
	r.header = _read(r.io, BlockHeader, E())
	r.endpos = position(r.io) + r.header._size
	return r
end

function _skip(r::BlockReader, offset::Int)
	while offset > 0
		rmg = offset - bytesleft(r)
		if rmg > 0
			nextblock(r)
			offset = rmg
		else
			return skip(r.io, offset)
		end
	end
end

function _read(r::BlockReader{E}, ::Type{T}) where {E, T<:Number}
	#NOTE: don't check if bytesleft<0... checked by reading BlockHeader
	if bytesleft(r) < 1
		nextblock(r)
	end

	if !canread(r, sizeof(T))
		hpos = hex(position(r.io))
		throw("Cannot read $T @ 0x$hpos")
	end

	return _read(r.io, T, E())
end

#Read in fixed-length string:
function _read(r::BlockReader, ::Type{String}, nchars::Int)
	if !canread(r, nchars)
		throw(stringboundary_exception(r.io))
	end
	buf = Array{UInt8}(undef, nchars)
	readbytes!(r.io, buf)
	return String(buf)
end

#Read in space-delimited string:
function readsigname(r::BlockReader)
	DELIM = UInt8(' ') #WANTCONST
	buf = Array{UInt8}(undef, SIGNAME_BUFSIZE)

	#TODO: improve test?
	isdelim(v::UInt8) = (v != DELIM)

	lastchar = DELIM
	while DELIM == lastchar
		lastchar = _read(r, UInt8)
	end

	if !isdelim(lastchar)
		hpos = hex(position(r.io)-1)
		throw("Invalid string @ 0x$hpos")
	end

	i = 1
	while isdelim(lastchar)
		buf[i] = lastchar
		lastchar = _read(r, UInt8)
		i+=1
		if i > SIGNAME_BUFSIZE
			throw("Insufficient buffer size: 'SIGNAME_BUFSIZE'")
		end
	end

	buf[i] = 0
	return unsafe_string(pointer(buf))
end

#==Constructors
===============================================================================#
#"Construct" a BlockReader, by reading a header from IO:
function BlockReader(io::IO, endianness::Endianness; start::Int=0)
	seek(io, start)
	hdr = _read(io, BlockHeader, endianness)
	endpos = position(io) + hdr._size
	return BlockReader{typeof(endianness)}(io, hdr, endpos)
end

#==Main functions
===============================================================================#
#Detect endianness from first word:
function _read(io::IO, ::Type{Endianness})
	w = read(io, DataWord)
	for endianness in [LittleEndian(), BigEndian()]
		if WRITEBLOCK_SYNCWORD == _reorder(w, endianness)
			return endianness
		end
	end
	throw(corruptword_exception(io, w, WRITEBLOCK_SYNCWORD))
end

#Read in SPICE data file format:
function _read(r::BlockReader, ::Type{SpiceFormat})
	versiontxt = strip(_read(r, String, 8))

	try
		version = parse(Int, versiontxt)
		if 9601 == version
			return Format_9601()
		elseif 9602 == version
			return Format_9602()
		elseif 2001 == version
			return Format_2001()
		elseif 2013 == version
			return Format_2013()
		end
	catch
	end

	versiontxt = printable(versiontxt)
	throw("SPICE data format not recognized: '$versiontxt'")
end

#Read in signal names:
function readnames(r::BlockReader, datacolumns::Int)
	for i in 1:datacolumns
		sigtype = _read(r, String, 8)
		try
			parse(Int, sigtype)
		catch
			hpos = hex(position(r.io) - 8)
			sigtype = printable(sigtype)
			throw("Non-numerical signal type '$sigtype' @ 0x$hpos")
		end
	end

	sweepname = readsigname(r)
	nsigs = datacolumns - 1

	signalnames = Array{String}(undef, nsigs)
	for i in 1:length(signalnames)
		signalnames[i] = readsigname(r)
	end

	return (sweepname, signalnames)
end

#Read in signal data to vector d.
function readsignal!(r::BlockReader, d::Vector{T}, offset::Int, rowsize::Int) where T
	rowskip = rowsize - sizeof(T)
	_skip(r, offset)
	lastrowcomplete = false
	npoints = 0
	lastpos = 0

	try
		while npoints <= length(d)
			val = _read(r, T)
			npoints += 1
			d[npoints] = val
			lastrowcomplete = false
#lastpos=position(r.io)
			_skip(r, rowskip)
			lastrowcomplete = true
		end
	catch
	end

	#When reading main sweep (offset == 0):
	#Get rid of last value if last row is not completely written:
	if 0 == offset && !lastrowcomplete
#hpos = hex(lastpos); chpos = hex(position(r.io))
#throw("INCOMPLETE DATASET: lastpos @0x$hpos (curpos @0x$chpos)")
		npoints = max(0, npoints-1)
	end

	return resize!(d, npoints)
end

#Read in main sweep vector:
function readsweep(r::BlockReader, fmt::SpiceFormat, rowsize::Int)
	#Compute estimated signal length:
	curpos = position(r.io)
	sz = filesize(r.io)
	estimatedlen = div(sz - curpos, rowsize)

	data = Array{xtype(fmt)}(undef, estimatedlen)
	return readsignal!(r, data, 0, rowsize)
end

#Read in signal by number:
function readsignal(r::DataReader, signum::Int)
	nsigs = length(r.signalnames)
	if signum < 1 || signum > nsigs
		throw("Invalid signal number: $signum ∉ [1, $nsigs].")
	end
	blkreader = BlockReader(r.io, r.endianness, start=r.datastart)
	_xtype = xtype(r.format); _ytype = ytype(r.format)
	offset = sizeof(_xtype) + (signum-1) * sizeof(_ytype)
	data = Array{_ytype}(undef, length(r.sweep))
	return readsignal!(blkreader, data, offset, r.rowsize)
end

#Read in a SPICE file from path:
function _open(filepath::String)
	io = open(filepath, "r")
	endianness = _read(io, Endianness)
	blkreader = BlockReader(io, endianness, start=0)

	#Read in signal counts:
	count1 = _read(blkreader, String, 4)
	#What are other counts for? Are they counts?
	count2 = _read(blkreader, String, 4)
	count3 = _read(blkreader, String, 4)
	count4 = _read(blkreader, String, 4)

	try
		count1 = parse(Int, count1)
		count2 = parse(Int, count2)
	catch
		throw("Invalid signal count.")
	end
	datacolumns = Int(count1)+Int(count2)

	#Read in file format:
	format = _read(blkreader, SpiceFormat)

	#Read in "tags":
	header = SpiceTags(
		strip(_read(blkreader, String, 4*16)), #id
		strip(_read(blkreader, String, 16)), #date
		strip(_read(blkreader, String, 8)), #time
		strip(_read(blkreader, String, 4*16+8)) #comments
	)

	#Read in signal names:
	_skip(blkreader, 5*16) #Why? What is here?
	sweepname, signalnames = readnames(blkreader, datacolumns)

	#Compute row size:
	nsigs = length(signalnames)
	_xtype = xtype(format); _ytype = ytype(format)
	rowsize = sizeof(_xtype) + nsigs*sizeof(_ytype)

	#Move to start of first data block:
	nextblock(blkreader)
	datastart = position(blkreader.io) - WRITEBLOCK_HEADERSIZE

	sweep = readsweep(blkreader, format, rowsize)

	return DataReader(io, filepath, format,
		sweepname, signalnames, sweep,
		header, datastart, rowsize, endianness
	)
end

#==Higher-level interface
===============================================================================#
#_open(filepath::String)
Base.names(r::DataReader) = r.signalnames
Base.read(r::DataReader, signum::Int) = readsignal(r, signum)
function Base.read(r::DataReader, signame::String)
	signum = findfirst(isequal(signame), r.signalnames)
	if nothing == signum
		throw("Signal not found: $signame.")
	end
	return readsignal(r, signum)
end
Base.close(r::DataReader) = close(r.io)

#Last line


text:
using SciMLCovidForecasting
using Test

@testset "SciMLCovidForecasting.jl" begin
    # Write your own tests here.
end


text:
################################################################################

"sort hidden Markov model by amplitude."
function sortHMM!(self::HMM)

  ampVc = map(amplitude, self.data)
  sortedAmpVc = sort(ampVc)
  ampIxVc = map(sortedAmpVc) do ζ
   findall(χ -> χ == ζ, ampVc)
  end
  templ = Array{Array{T, 1} where T <: Number, 1}()

  for ι ∈ ampIxVc
    push!(templ, self.data[ι[1]])
  end
  self.data = templ
end

################################################################################

max_stars_repo_name: DanielRivasMD/HiddenMarkovModelReaders

max_issues_repo_name: DanielRivasMD/HiddenMarkovModelReaders


text:
# push!(LOAD_PATH,"../src/")

using LinearMapsAA
using Documenter

DocMeta.setdocmeta!(LinearMapsAA, :DocTestSetup, :(using LinearMapsAA); recursive=true)

makedocs(;
    modules = [LinearMapsAA],
    authors = "Jeff Fessler <fessler@umich.edu> and contributors",
    repo = "https://github.com/JeffFessler/LinearMapsAA.jl/blob/{commit}{path}#{line}",
    sitename = "LinearMapsAA.jl",
    format = Documenter.HTML(;
        prettyurls = get(ENV, "CI", "false") == "true",
#       canonical = "https://JeffFessler.github.io/LinearMapsAA.jl/stable",
#       assets = String[],
    ),
    pages = [
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo = "github.com/JeffFessler/LinearMapsAA.jl.git",
    devbranch = "main",
    devurl = "dev",
    versions = ["stable" => "v^", "dev" => "dev"]
#   push_preview = true,
)


text:
# Resolve

# 1. This is the (constant) scalar field the elevation measurements refer to
#
# struct ScalarField1D
#     # scalar field over one variable f:R->R (i.e. h= f(x))
#     # [a special case of ScalarField, h = f(·)]
#     # ...
# end

# 2. This is the actual scalar measurement sequence factor
# It takes two arguments:
#  - the elevation measurement sequence, and
#  - the odometry estimate at each elevation measurement
#
# addFactor!(fg, [:x1, :topography], ScalarFieldSequenceFactorNorthSouth(x,y))
struct ScalarFieldSequenceFactorNorthSouth{T} <: IIF.AbstractRelativeMinimize
  # FIXME: objective is to <: IIF.AbstractRelativeConstant (final name TBD)
  prevSequence::Vector{Float64} # sample location (odometry)
  nextSequence::Vector{Float64} # sample value (elevation meas)
  partial::T
end

# 3. Objective/fitness function evaluating match between a map and a template for different displacement values
# This is used within the factor to generate the likelihood over the terrain
#
# NOTE: this function assumes map and template are on an equivalent grid
# this function does not care for where/what the grid is
# NOTE: result will assume base pose is first in sequence
# for each i
#   energy[i] = sum( (map[i:(i+length(template))] - template).^2 )
# density = exp.(-energy)     # Woodward
# returns vector the size of length(map)
function _mySSDCorr(map, template)

    maplen = length(map)
    len = length(template)

    mnm = Statistics.mean(map)
    map_ = [map; mnm*ones(len-1)]
    density = zeros(maplen)
    for i in 1:maplen
      density[i] = -sum( (map_[i:(i+len-1)] .- template).^2 )
    end
    adaptive = abs(maximum(density))

    density ./= adaptive
    density .= exp.(density)
    density ./= sum(density)
    return density
end

function IIF.getSample(s::CalcFactor{<:ScalarFieldSequenceFactorNorthSouth}, N::Int=1)
  # locality if desired from current variable estimate,
  # X0 = getBelief(s.metadata.fullvariables[1]) # to consider only local map info

  # assuming user addConstant!(fg, :terrain, ScalarField1D(dem))

  # global terrain # should come from s.constants[:terrain] ?????

  # ensure equal spacing as terrain
  # assumptions:
  # - terrain.t is equally spaced
  # - s.factor.gridSequence increases monotonically
  # s.factor.gridSequence is the sequence of odometry estimates at which
  #  the elevation measurements (s.factor.scalarSequence) were taken
    # measurement = LinearInterpolation(s.factor.scalarSequence,s.factor.gridSequence)

    # upsample = 1 # avoid quantization effects from gridded sequence/terrain data
    # dx = diff(terrain.t)[1]/upsample
    # template = measurement.(s.factor.gridSequence[1]:dx:s.factor.gridSequence[end])
    # map_ = [terrain.u;]
    # map2_grid = range(terrain.t[1],terrain.t[end],length=upsample*length(map_))
    # map2 = terrain.(map2_grid)

  intensity_ = _mySSDCorr(prevSequence, nextSequence )

  # AliasingSS only works for 1D at this time
  # NOTE, if you stray out of region of support, things blow up!
  bss = AliasingScalarSampler([map2_grid;], intensity_) # [terrain.t;]

  # buffer with kde for continuous sampling (not grid spaced)
  stagedsmpls = rand(bss,N)
  smpls = rand(kde!(stagedsmpls, [3*dx;]),N)

  return (reshape(smpls,1,N),)
end


text:
module ParametricSurfaces

using StaticArrays
using LinearAlgebra
using ForwardDiff # for computing derivatives of parametric elements
using RecipesBase

using WavePropBase
using WavePropBase.Utils
using WavePropBase.Geometry
using WavePropBase.Interpolation
using WavePropBase.Mesh

WavePropBase.@import_interface

export
    # re-exported from WavePropBase
    clear_entities!,
    ElementaryEntity,
    Domain,
    skeleton,
    internal_boundary,
    external_boundary,
    HyperRectangle,
    ReferenceLine,
    ReferenceSquare,
    #types
    ParametricEntity,
    ParametricElement,
    #functions
    line,
    meshgen

include("parametricentity.jl")
include("parametricelement.jl")
include("meshgen.jl")
include("simpleshapes.jl")

WavePropBase.@export_interface

end # module

max_stars_repo_name: WaveProp/ParametricSurfaces

max_issues_repo_name: WaveProp/ParametricSurfaces

max_forks_repo_name: WaveProp/ParametricSurfaces


text:
# This file is part of Kpax3. License is MIT.

function test_prior_row_exceptions()
  @test_throws Kpax3.KDomainError Kpax3.EwensPitman(1.0, 0.0)
  @test_throws Kpax3.KDomainError Kpax3.EwensPitman(2.0, 0.0)

  @test_throws Kpax3.KDomainError Kpax3.EwensPitman(0.5, -1.0)
  @test_throws Kpax3.KDomainError Kpax3.EwensPitman(0.5, -0.5)

  @test_throws Kpax3.KDomainError Kpax3.EwensPitman(0.0, -1.0)
  @test_throws Kpax3.KDomainError Kpax3.EwensPitman(0.0, 0.0)

  @test_throws Kpax3.KDomainError Kpax3.EwensPitman(0.5, 0)
  @test_throws Kpax3.KDomainError Kpax3.EwensPitman(0.0, 0)
  @test_throws Kpax3.KDomainError Kpax3.EwensPitman(0.5, 1)
  @test_throws Kpax3.KDomainError Kpax3.EwensPitman(0.0, 1)

  @test_throws Kpax3.KDomainError Kpax3.EwensPitman(-1.0, 0)
  @test_throws Kpax3.KDomainError Kpax3.EwensPitman(-1.0, -1)
  @test_throws Kpax3.KDomainError Kpax3.EwensPitman(-1.0, 1.0)

  nothing
end

test_prior_row_exceptions()

function test_prior_row_constructor()
  priorR = Kpax3.EwensPitman(0.5, -0.1)
  @test isa(priorR, Kpax3.EwensPitmanPAUT)
  @test priorR.α == 0.5
  @test priorR.θ == -0.1

  priorR = Kpax3.EwensPitman(0.5, 0.0)
  @test isa(priorR, Kpax3.EwensPitmanPAZT)
  @test priorR.α == 0.5

  priorR = Kpax3.EwensPitman(0, 2.0)
  @test isa(priorR, Kpax3.EwensPitmanZAPT)
  @test priorR.θ == 2.0

  priorR = Kpax3.EwensPitman(-2, 5)
  @test isa(priorR, Kpax3.EwensPitmanNAPT)
  @test priorR.α == -2.0
  @test priorR.L == 5

  nothing
end

test_prior_row_constructor()


text:
using ADCME
using PyCall
using LinearAlgebra
using PyPlot
using Random
Random.seed!(233)

fint_comp = load_op_and_grad("./build/libFintComp", "fint_comp")
################## End Load Operator ##################

neqs_ = 20
θ = constant(rand(3,3))
A_ = rand(10,3)
fints = A_*θ
el = constant(rand(0:neqs_, 10, 3), dtype=Int32)
neqs = constant(neqs_, dtype=Int32)
# TODO: specify your input parameters
u = fint_comp(fints,el,neqs)
sess = tf.Session()
init(sess)
run(sess, u)
error()

# TODO: change your test parameter to `m`
# gradient check -- v
function scalar_function(θ)
    m = A_*θ
    return sum(tanh(fint_comp(m,el,neqs)))
end

# TODO: change `m_` and `v_` to appropriate values
m_ = constant(rand(3,3))
v_ = rand(3,3)
y_ = scalar_function(m_)
dy_ = gradients(y_, m_)
ms_ = Array{Any}(undef, 5)
ys_ = Array{Any}(undef, 5)
s_ = Array{Any}(undef, 5)
w_ = Array{Any}(undef, 5)
gs_ =  @. 1 / 10^(1:5)

for i = 1:5
    g_ = gs_[i]
    ms_[i] = m_ + g_*v_
    ys_[i] = scalar_function(ms_[i])
    s_[i] = ys_[i] - y_
    w_[i] = s_[i] - g_*sum(v_.*dy_)
end

sess = Session()
init(sess)
sval_ = run(sess, s_)
wval_ = run(sess, w_)
close("all")
loglog(gs_, abs.(sval_), "*-", label="finite difference")
loglog(gs_, abs.(wval_), "+-", label="automatic differentiation")
loglog(gs_, gs_.^2 * 0.5*abs(wval_[1])/gs_[1]^2, "--",label="\$\\mathcal{O}(\\gamma^2)\$")
loglog(gs_, gs_ * 0.5*abs(sval_[1])/gs_[1], "--",label="\$\\mathcal{O}(\\gamma)\$")

plt.gca().invert_xaxis()
legend()
xlabel("\$\\gamma\$")
ylabel("Error")


text:
using StatsBase, SnpArrays

function snparray_simulation(maf::Vector{Float64}, n::Int)
    genotype = genotype_sim.(maf, n)
    snparray = make_snparray(genotype)
    return snparray
end

# this is for generating genotypes given minor allele frequency under HWE

function genotype_sim(maf::Float64, n::Int)
    pAA = (1 - maf)^2 #0
    pAa = 2maf * (1 - maf) #1
    paa = maf^2     #2
    pvec = [pAA, pAa, paa]
    genes = [0, 1, 2]
    g = StatsBase.wsample(genes, pvec, n)
    return g
end

"""
Make a SnpArray
"""
function make_snparray(genotype)
    n = length(genotype[1])
    p = length(genotype)
    x = SnpArray(undef, n, p)
    for j in 1:p, i in 1:n
        j_locus = genotype[j]
        if j_locus[i] == 0
            x[i, j] = 0x00
        elseif j_locus[i] == 1
            x[i, j] = 0x02
        elseif j_locus[i] == 2
            x[i, j] = 0x03
        else
            throw(error("matrix shouldn't have missing values!"))
        end
    end
    return x
end

# Effect Size Simulation Conditional on MAF
function simulate_effect_size(maf)
    Simulated_ES = ones(length(maf))
    # Generating Effect Sizes where the lower the minor allele frequency the larger the effect size
    for i in eachindex(maf)
        Simulated_ES[i] = rand([-1, 1]) .* (0.01 / sqrt.(maf[i] .* (1 - maf[i])))
    end
    return Simulated_ES
end


text:
## Here is the code from unit 1

##
using DelimitedFiles
using SparseArrays
ei,ej = copy.(eachcol(Int.(readdlm("wiki-simple.edges"))))
xy = readdlm("wiki-simple.xy")
pages = readlines("wiki-simple.nodes")
n = size(xy,1)
A = sparse(ei,ej,1,n,n)

##  PageRank with coordinate relaxation
#
function relaxpagerank(A,α,v)
  @assert(0 ≤ α < 1, "needs probably α")
  @assert(all(vi -> vi ≥ 0, v), "needs non-negative v")
  v = v ./ sum(v) # we can normalize for them.
  d = vec(sum(A,dims=2)) # compute the degrees
  x = copy(v) # start of with v
  nsteps = 2*ceil(Int,log(eps(1.0))/log(α)) # upper bound on steps
  # we need row's of A', which are columns of A.
  # it's also easier to iterate over Di*x, instead of x
  d[d .== 0] .= 1 # switch 0's to 1's
  idx = x./d
  for i=1:nsteps
    for j=1:size(A,1)
      rhs = (1-α)*v[j]  + α*(A[:,j]'*idx)
      gamma = rhs/(1 - α*A[j,j])
      idx[j] = gamma/d[j]
    end
  end
  x =  (d.*idx)
  return x/sum(x) # renormalize to probability
end
x = relaxpagerank(A,0.85,ones(n)/n)

##
using MatrixNetworks
pr = pagerank(A,0.85)
norm(x-pr)


text:
module LTC

using Reexport

using LazyArrays
using Distributions
import NPZ: npzread
export npzread
using Juno
using LoopVectorization
using StructArrays
using MacroTools: @forward

#using Flux

using DiffEqBase
using OrdinaryDiffEq
using DiffEqSensitivity
using DiffEqFlux
using DiffEqFlux: initial_params, paramlength, FastChain, FastDense, sciml_train
import DiffEqFlux: initial_params, paramlength, FastChain, FastDense, sciml_train
export sciml_train
using GalacticOptim
#export initial_params
#using DiffRules
#using RecursiveArrayTools
#import RecursiveArrayTools: DiffEqArray
#using DiffEqCallbacks
#using DistributionsAD
# using ModelingToolkit
#using BenchmarkTools
#using Plots
#using Noise
#@reexport using Functors
#export Functors
#using ComponentArrays
#using Parameters: @unpack
using Tullio
using Zygote
using Zygote: @adjoint, Numeric, literal_getproperty, accum
export Zygote
using Flux: reset!, Zeros
import Flux: reset!
#export reset!
# using ReverseDiff
using ForwardDiff

#using ZygoteRules
#using ForwardDiff
#using Sundials
#using ODEInterfaceDiffEq

#@reexport using Statistics

rand_uniform(TYPE, lb,ub,dims...) = TYPE.(rand(Uniform(lb,ub),dims...))# |> f64

#Zygote.@nograd rand_uniform, reshape

#include("ltc-modelcleanup.jl")
include("wiring.jl")
#include("ncp_ensemble.jl")
#include("ncp.jl")
#include("_ncp.jl")
include("ltcnet.jl")
#include("_ncp_1.jl")
#include("new_ncp.jl")
#include("new_ncp_ensemble.jl")
include("custom_train.jl")
#include("tullio-sigmoid.jl")
# include("param_ranges.jl")

export get_bounds

export Wiring, NCPWiring
export LTCNet, Mapper, gNN, get_bounds, my_custom_train!
export reset!, reset_state!
export initial_params, paramlength
end


text:
@testset "Q3D plain text parsing" begin
    data(x) = joinpath(@__DIR__, "data", x)

    # file should satisfy some basic checks
    @test_throws AssertionError Circuit(data("empty.txt"))
    @test_throws AssertionError Circuit(data("one_block.txt"))

    # test expected output
    d = data("dummy_gc_1.txt")
    @test Circuit(d).vertices == ["net1", "net2", "net3", "net4"]

    design_variation, vertex_names, capacitance_matrix =
        AdmittanceModels.parse_q3d_txt(d, :capacitance)

    # test negative val + no units
    @test design_variation["dummy1"] == (-1.23, "")

    # test negative exponent + units
    @test design_variation["dummy2"] == (9.0e-7, "mm")

    # test negative val, positive exponent + units
    @test design_variation["dummy3"] == (-1.0e21, "pF")

    # test units that start with "e" (the parser should look for things like 'e09' or 'e-12')
    @test design_variation["dummy4"] == (1.2, "eunits")

    # unit not implemented
    @test_throws ErrorException try
        Circuit(data("dummy_gc_2.txt"); matrix_types = [:conductance])
    catch e
        buf = IOBuffer()
        showerror(buf, e)
        message = String(take!(buf))
        @test occursin("not implemented", message)
        rethrow(e)
    end

    @test_throws ErrorException try
        Circuit(data("dummy_gc_3.txt"); matrix_types = [:conductance])
    catch e
        buf = IOBuffer()
        showerror(buf, e)
        message = String(take!(buf))
        @test occursin("units not given", message)
        rethrow(e)
    end
end