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

using Dates

"""
    abstract type NeighborhoodReduction end

Overrides `Base.reduce(::NeighborhoodReduction, res::KnnResult, index::SearchGraph)` to postprocess `res` using some criteria.
Called from `find_neighborhood`, and returns a new KnnResult struct (perhaps a copy of res) since `push_neighborhood` captures
the reference of its output.
"""
abstract type NeighborhoodReduction end
abstract type LocalSearchAlgorithm end

"""
    abstract type Callback end

Abstract type to trigger callbacks after some number of insertions.
SearchGraph stores the callbacks in `callbacks` (a dictionary that associates symbols and callback objects);
A SearchGraph object controls when callbacks are fired using `callback_logbase` and `callback_starting`

"""
abstract type Callback end

### Basic operations on the index

"""
    struct SearchGraph <: AbstractSearchContext

SearchGraph index. It stores a set of points that can be compared through a distance function `dist`.
The performance is determined by the search algorithm `search_algo` and the neighborhood policy.
It supports callbacks to adjust parameters as insertions are made.

- `hints`: Initial points for exploration (empty hints imply using random points)

Note: Parallel insertions should be made through `append!` or `index!` function with `parallel_block > 1`

"""
@with_kw struct SearchGraph{DistType<:SemiMetric, DataType<:AbstractDatabase, SType<:LocalSearchAlgorithm}<:AbstractSearchContext
    dist::DistType = SqL2Distance()
    db::DataType = VectorDatabase()
    links::Vector{Vector{Int32}} = Vector{Int32}[]
    locks::Vector{Threads.SpinLock} = Threads.SpinLock[]
    hints::Vector{Int32} = Int32[]
    search_algo::SType = BeamSearch()
    verbose::Bool = true
end

@inline Base.length(g::SearchGraph) = length(g.locks)
include("visitedvertices.jl")

## search algorithms

"""
    SearchGraphPools(results=GlobalKnnResult, vstates=GlobalVisitedVertices, beams=GlobalBeamKnnResult)

A set of pools to alleviate memory allocations in `SearchGraph` construction and searching. Relevant on multithreading scenarious where distance functions, `evaluate`
can call other metric indexes that can use these shared resources (globally defined).

Each pool is a vector of `Threads.nthreads()` preallocated objects of the required type.
"""
struct SearchGraphPools{VisitedVerticesType}
    beams::Vector{KnnResult}
    satnears::Vector{KnnResult}
    vstates::VisitedVerticesType
end

@inline function getvstate(len, pools::SearchGraphPools)
    @inbounds v = pools.vstates[Threads.threadid()]
    _init_vv(v, len)
end

@inline function getbeam(bsize::Integer, pools::SearchGraphPools)
    @inbounds reuse!(pools.beams[Threads.threadid()], bsize)
end

@inline function getsatknnresult(pools::SearchGraphPools)
    reuse!(pools.satnears[Threads.threadid()], 1)
end

"""
    getpools(index::SearchGraph)

Creates or retrieve caches for the search graph.
"""
getpools(::SearchGraph; beams=GlobalBeamKnnResult, satnears=GlobalSatKnnResult, vstates=GlobalVisitedVertices) = SearchGraphPools(beams, satnears, vstates)

include("beamsearch.jl")
## parameter optimization and neighborhood definitions
include("optbs.jl")
include("neighborhood.jl")
include("hints.jl")

"""
    search(index::SearchGraph, q, res; hints=index.hints, pools=getpools(index))

Solves the specified query `res` for the query object `q`.
"""
function search(index::SearchGraph, q, res::KnnResult; hints=index.hints, pools=getpools(index))
    if length(index) > 0
        search(index.search_algo, index, q, res, hints, pools)
    else
        (res=res, cost=0)
    end
end

max_stars_repo_name: sadit/SimilaritySearch

max_issues_repo_name: sadit/SimilaritySearch

max_forks_repo_name: sadit/SimilaritySearch


text:
using Images, TestImages, Colors, ZernikePolynomials, FFTW
using NumberTheoreticTransforms

image_float = channelview(testimage("cameraman"))
image_int = map(x -> x.:i, image_float) .|> Int64

blur_float = evaluateZernike(LinRange(-41,41,512), [12, 4, 0], [1.0, -1.0, 2.0], index=:OSA)
blur_float ./= (sum(blur_float))
blur_float = blur_float
blur_int = blur_float .|> Normed{UInt8, 8} .|> x -> x.:i .|> Int64

blur_int = fftshift(blur_int)

@assert sum(blur_int) == 256 # to not change brightness

blurred_image_float = ifft(fft(image_float) .* fft(fftshift(blur_float))) |> real

t = 4
(g, q) = (314, 2^2^t+1)
X = fnt(image_int, g, q)
@assert ifnt(X, g, q) == image_int
H = fnt(blur_int, g, q)
@assert ifnt(H, g, q) == blur_int
Y = mod.((X .* H), q)
y = ifnt(Y, g, q)
blurred_image_int = y / 2^16

using Statistics
@show maximum(abs.(blurred_image_float .- blurred_image_int))
@show mean(abs.(blurred_image_float .- blurred_image_int))
@show var(abs.(blurred_image_float .- blurred_image_int))

save("doc/src/fnt/original.jpg", Images.clamp01nan.(image_float))
save("doc/src/fnt/blurred_fft.jpg", Images.clamp01nan.(Gray.(Normed{UInt8, 8}.(blurred_image_float))))
save("doc/src/fnt/blurred_fnt.jpg", Images.clamp01nan.(Gray.(Normed{UInt8, 8}.(blurred_image_int))))

ext: jl


text:
# ---
# title: 230. Kth Smallest Element in a BST
# id: problem230
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Binary Search, Tree
# link: <https://leetcode.com/problems/kth-smallest-element-in-a-bst/description/>
# hidden: true
# ---
#
# Given a binary search tree, write a function `kthSmallest` to find the **k**
# th smallest element in it.
#
#
#
# **Example 1:**
#
#
#
#     Input: root = [3,1,4,null,2], k = 1
#        3
#       / \
#      1   4
#       \
#        2
#     Output: 1
#
# **Example 2:**
#
#
#
#     Input: root = [5,3,6,2,4,null,null,1], k = 3
#            5
#           / \
#          3   6
#         / \
#        2   4
#       /
#      1
#     Output: 3
#
#
# **Follow up:**
# What if the BST is modified (insert/delete operations) often and you need to
# find the kth smallest frequently? How would you optimize the kthSmallest
# routine?
#
#
#
# **Constraints:**
#
#   * The number of elements of the BST is between `1` to `10^4`.
#   * You may assume `k` is always valid, `1 ≤ k ≤ BST's total elements`.
#
#
## @lc code=start
using LeetCode

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


text:
# boosting algorithm based on minimizing KL(q||p)
# ELBO
"approximated ELBO based on n samples"
function ELBO(approx_d::Distribution, ll::Function; n=100000)
    x = rand(approx_d, n)
    return ELBO(approx_d, ll, x)
end

"approximate ELBO using supplied samples x"
function ELBO{T <: Real}(approx_d::Distribution, ll::Function, x::Array{T})
    Z = ll(x) - logpdf(approx_d, x)
    return safemean(Z)
end

# mean, with Inf and NaN removed
function safemean{T <: Real}(x::AbstractVector{T})
    # @show mean(isfinite(x))
    x = x[isfinite.(x)]
    !isempty(x) || error("No finite value left for mean: x =", x)
    return mean(x)
end

# stochastic approximation of L as a function of α
"return approximate exclusive L as a function of α, with supplied samples"
function loss_exclusive{T <: Real}(q::Distribution, h::Distribution, ll::Function, samples_from_q::Array{T}, samples_from_h::Array{T})
    if length(q)==1
        @assert length(samples_from_q) == length(samples_from_h)
    else
        @assert size(samples_from_q) == size(samples_from_h)
    end
    f1 = _helper_dL_exclusive(q, h, ll, samples_from_h)
    f2 = _helper_dL_exclusive(q, h, ll, samples_from_q)
    return α -> α * f1(α) + (1-α) * f2(α)
end

"return approximate exclusive L as a function of α, with n samples"
function loss_exclusive(q::Distribution, h::Distribution, ll::Function; n=10000)
    samples_from_h = rand(h, n)
    samples_from_q = rand(q, n)
    return loss_exclusive(q, h, ll, samples_from_q, samples_from_h)
end

# stochastic approximations of dL/dα and d2L/dα2
"return ∫ p(x) [log ((1-α) q + α h) - log L] dx, as a function of α, with supplied samples x ~ p"
function _helper_dL_exclusive{T <: Real}(q::Distribution, h::Distribution, ll::Function, x::Array{T})
    z1 = pdf(q, x) + ϵ_lower_bound
    z2 = pdf(h, x) + ϵ_lower_bound
    z3 = logsumexp(ll(x) + log(ϵ_lower_bound))
    # @show mean(isfinite(z1)), mean(isfinite(z2)), mean(isfinite(z3))
    return α -> safemean( log.((1-α) * z1 + α * z2) - z3 )
end

"return approximate dL/dα as a function of α, with supplied samples"
function dL_exclusive{T <: Real}(q::Distribution, h::Distribution, ll::Function, samples_from_q::Array{T}, samples_from_h::Array{T})
    if length(q)==1
        @assert length(samples_from_q) == length(samples_from_h)
    else
        @assert size(samples_from_q) == size(samples_from_h)
    end
    f1 = _helper_dL_exclusive(q, h, ll, samples_from_h)
    f2 = _helper_dL_exclusive(q, h, ll, samples_from_q)
    return α -> f1(α) - f2(α)
end

"return approximate dL/dα as a function of α, with n samples"
function dL_exclusive(q::Distribution, h::Distribution, ll::Function; n=10000)
    samples_from_h = rand(h, n)
    samples_from_q = rand(q, n)
    return dL_exclusive(q, h, ll, samples_from_q, samples_from_h)
end

"return ∫ p(x) [(h - q) / ((1-α)q + αh)] dx, as a function of α, with supplied samples x ~ p"
function _helper_d2L_exclusive{T <: Real}(q::Distribution, h::Distribution, x::Array{T})
    hh = pdf(h, x) + ϵ_lower_bound
    qq = pdf(q, x) + ϵ_lower_bound
    # @show mean(isfinite(hh)), mean(isfinite(qq))
    return α -> safemean( (hh - qq) ./ ((1-α) * qq + α * hh) )
end

"return approximate d2L_exclusive/dα2 as a function of α, with supplied samples"
function d2L_exclusive{T <: Real}(q::Distribution, h::Distribution, samples_from_q::Array{T}, samples_from_h::Array{T})
    if length(q)==1
        @assert length(samples_from_q) == length(samples_from_h)
    else
        @assert size(samples_from_q) == size(samples_from_h)
    end
    f1 = _helper_d2L_exclusive(q, h, samples_from_h)
    f2 = _helper_d2L_exclusive(q, h, samples_from_q)
    return α -> f1(α) - f2(α)
end

"return approximate d2L/dα2 as a function of α, with n samples"
function d2L_exclusive(q::Distribution, h::Distribution; n=10000)
    samples_from_h = rand(h, n)
    samples_from_q = rand(q, n)
    return d2L_exclusive(q, h, samples_from_q, samples_from_h)
end

# minimization of L(α)
"search α with Newton's method"
function search_α_exclusive_newton(approx_d::Distribution, new_component::Distribution, ll::Function;
                  init_step_size=1., n=100, n_plot=1000, tol=1e-4, max_iter=1000, verbose=false, show_plot=false)
    # unicode plot
    if show_plot
        # get samples
        samples_from_q = rand(approx_d, n_plot)
        samples_from_h = rand(new_component, n_plot)
        approx_loss = loss_exclusive(approx_d, new_component, ll, samples_from_q, samples_from_h)
        α_vec = linspace(0.01, 0.99, 20)
        print(lineplot(α_vec, map(approx_loss, α_vec), title="loss"))
    end
    α = 0.
    iter = 0
    change = 1.
    while true
        iter += 1
        if iter > max_iter
            warn(@sprintf("change=%5g > %5g, max_iter=%2d reached", change, tol, max_iter))
            break
        end
        # get stochastic approximates
        samples_from_q = rand(approx_d, n)
        samples_from_h = rand(new_component, n)
        approx_dL = dL_exclusive(approx_d, new_component, ll, samples_from_q, samples_from_h)
        approx_d2L = d2L_exclusive(approx_d, new_component, samples_from_q, samples_from_h)
        # Newton's
        step_size =  init_step_size / iter
        α_new = α - step_size * approx_dL(α) / approx_d2L(α)
        α_new = min(α_new, 1.0)
        α_new = max(α_new, 0.0)
        change = abs(α_new - α)
        α = α_new
        if change < tol && iter > 200
            break
        end
    end
    if verbose
        println(@sprintf "α = %5g (change = %5g < tol = %5g) with %2d iterations." α change tol iter)
    end
    return α
end

# minimization of L(α)
"search α with SGD"
function search_α_exclusive_sgd(approx_d::Distribution, new_component::Distribution, ll::Function;
                  init_step_size=1.0, n=100, tol=1e-4, min_iter = 100, max_iter=10000, verbose=false)
    α = 0.
    iter = 0
    change = 1.
    while true
        iter += 1
        if iter > max_iter
            warn(@sprintf("change=%5g > %5g, max_iter=%2d reached", change, tol, max_iter))
            break
        end
        # get stochastic approximates
        samples_from_q = rand(approx_d, n)
        samples_from_h = rand(new_component, n)
        approx_dL = dL_exclusive(approx_d, new_component, ll, samples_from_q, samples_from_h)
        # Newton's
        step_size =  init_step_size / iter
        α_new = α - step_size * approx_dL(α)
        α_new = min(α_new, 1.0)
        α_new = max(α_new, 0.0)
        change = abs(α_new - α)
        α = α_new
        if change < tol && iter > min_iter
            break
        end
    end
    if verbose
        println(@sprintf "α = %5g (change = %5g < tol = %5g) with %2d iterations." α change tol iter)
    end
    return α
end

max_stars_repo_name: richardkwo/boosting-VB

max_issues_repo_name: richardkwo/boosting-VB

max_forks_repo_name: richardkwo/boosting-VB


text:
mutable struct FCGIRequest
    headers::Dict{String,String}
    in::IO
    out::IO
    err::IO
    keepconn::Bool
    isdone::Channel{Nothing}
    exception::Any
    function FCGIRequest(; headers::Dict{String,String}=Dict{String,String}(), in::IO=PipeBuffer(), out::IO=PipeBuffer(), err::IO=PipeBuffer(), keepconn::Bool=false)
        new(headers, in, out, err, keepconn, Channel{Nothing}(1), nothing)
    end
end
show(io::IO, request::FCGIRequest) = print(io, "FCGIRequest(", "complete=", isready(request.isdone), ", keepconn=", request.keepconn, ")")

mutable struct FCGIClient{T<:FCGISocket}
    addr::String
    csock::T
    nextreqid::UInt16
    processor::Union{Task,Nothing}
    requests::Dict{UInt16,FCGIRequest}
    lck::ReentrantLock
end

function FCGIClient(addr::String, csock::FCGISocket)
    client = FCGIClient(addr, csock, UInt16(0), nothing, Dict{UInt16,FCGIRequest}(), ReentrantLock())
    client.processor = @async process(client)
    client
end
FCGIClient(path::String) = FCGIClient(path, connect(path))
FCGIClient(ip::IPv4, port::Integer) = FCGIClient("$(ip):$(Int(port))", connect(ip, port))
show(io::IO, client::FCGIClient) = print(io, "FCGIClient(", client.addr, ", ", isopen(client.csock) ? "open" : "closed", ")")
close(client::FCGIClient) = close(client.csock)
isrunning(client::FCGIClient) = isopen(client.csock)
function fcgiwrite(client::FCGIClient{T}, rec::FCGIRecord) where {T <: FCGISocket}
    lock(client.lck) do
        fcgiwrite(client.csock, rec)
    end
end
function requestdone(client::FCGIClient, reqid::UInt16)
    lock(client.lck) do
        if reqid in keys(client.requests)
            req = client.requests[reqid]
            delete!(client.requests, reqid)
            put!(req.isdone, nothing)
            req.keepconn || close(client)
        end
    end
end
function process(client::FCGIClient)
    try
        while isopen(client.csock)
            record = FCGIRecord(client.csock)
            reqid = requestid(record.header)
            if record.header.type == FCGIHeaderType.END_REQUEST
                requestdone(client, reqid)
            elseif record.header.type == FCGIHeaderType.STDERR
                req = client.requests[reqid]
                writefully(req.err, record.content)
            elseif record.header.type == FCGIHeaderType.STDOUT
                req = client.requests[reqid]
                writefully(req.out, record.content)
            else
                @warn("unexpected message received", reqid, type=record.header.type)
            end
        end
    finally
        for reqid in collect(keys(client.requests))
            requestdone(client, reqid)
        end
        close(client.csock)
    end
end

function process(client::FCGIClient, request::FCGIRequest)
    reqid = lock(client.lck) do
        client.nextreqid += UInt16(1)
        client.requests[client.nextreqid] = request
        client.nextreqid
    end
    # begin request
    fcgiwrite(client, FCGIRecord(FCGIHeaderType.BEGIN_REQUEST, reqid, FCGIBeginRequest(FCGIRequestRole.RESPONDER, request.keepconn ? KEEP_CONN : UInt8(0))))
    # send params
    fcgiwrite(client, FCGIRecord(FCGIHeaderType.PARAMS, reqid, FCGIParams(request.headers)))
    # send params end
    fcgiwrite(client, FCGIRecord(FCGIHeaderType.PARAMS, reqid, FCGIParams()))

    try
        # send stdin
        buffin = BufferedOutput() do bytes
            isempty(bytes) || fcgiwrite(client, FCGIRecord(FCGIHeaderType.STDIN, reqid, bytes))
        end
        while !eof(request.in)
            write(buffin, readavailable(request.in))
        end
        flush(buffin)
        close(buffin)
        # send stdin end
        fcgiwrite(client, FCGIRecord(FCGIHeaderType.STDIN, reqid, UInt8[]))
    catch ex
        request.exception = ex
    end

    request
end


text:
include("DM1_1.jl")
include("DM2_1.jl")
include("DM1_2_corrige.jl")
#=========================================#
function ReactiveGRASP( #Utilisation du Reactive GRASP
    matrix, #Un matrice de taille n*m qui représente les conditions
    cost, #Une liste des couts de chaque variable
    n, #nombre de variables
    m, #nombre de contraintes
    ite, #on tierera ce nombre avant de recalculer des probabilités de chaque alpha
    coupe, #Le nombre de coupe du segment [0,1], nous permet de construire nos alpha
    alphaset,
    temps #La ressource en temps en seconde alloué
    )

    #initialisation
    (p,nb_iteration,z_cumul,zBest,zWorst) = intialiser(matrix,cost, n, m, ite, coupe, alphaset)

    evol_p=Float64[] #liste qui stokera l'évolution de nos probabilités au fur et à mesure des itérations.
    t=time()
    cpt = 1

    while (time()-t <= temps) #Condtion sur les ressources en temps
        append!(evol_p,p)
        while (cpt <= ite)

            for i in 1:length(p) #On s'assure que chaque probabilité se soit exprimé au moins une fois
                (SOL,z, crts) = GRASP(cost, matrix, n, m, p[i])
                z_cumul[i] += z
            end

            prob=rand(Float64)
            #prob = rand(1:length(p))

            alpha_choisit = choix_alpha(p,prob)
            (SOL, z, crts) = GRASP(cost, matrix, n, m, p[alpha_choisit])
            #(SOL, z, crts) = GRASP(cost, matrix, n, m, p[prob])

            #Amelioration
            (SOL, z) = exchange1_1(SOL,n,m,cost,crts,matrix)

            nb_iteration[alpha_choisit] += 1
            #nb_iteration[prob] += 1
            z_cumul[alpha_choisit] += z
            #z_cumul[prob] += z

            if z > zBest
                zBest = z
            end

            if z < zWorst
                zWorst = z
            end

            cpt = cpt+1 #Compteur du nombre d'itération
        end

        recalcul_p!(p,z_cumul,zBest,zWorst,nb_iteration,evol_p)
        nb_iteration = ones(Int64, length(p))
        z_cumul = zeros(Int64, length(p))

    end
    println("zBest: ", zBest, "    zWorst: ", zWorst)
    println(evol_p) #A commenter si on ne souhaite pas voir la liste d'évolution des probabilités
    return(evol_p)
end

#=======================================#
function intialiser(matrix,cost, n, m, ite, coupe, alphaset) #initialisationd des variables

    p = Float64[]
    if alphaset == 0
        for i in 1:coupe
            push!(p,i/coupe)
        end
    else
        p = alphaset
    end

    nb_iteration = ones(Int64, length(p))

    z_cumul = zeros(Int64, length(p))

    (SOL,z) =Glouton(cost, matrix, n, m)

    zBest = z

    zWorst = z

    return(p,nb_iteration,z_cumul,zBest,zWorst)
end

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

function recalcul_p!(p,z_cumul,zBest,zWorst,nb_iteration,evol_p)
    q=zeros(Float64,length(p))

    somme_q=0

    for i in 1:length(p)

        moyenne=(z_cumul[i]/nb_iteration[i])
        q[i]=(moyenne-zWorst)/(zBest-zWorst)
        somme_q+=q[i]

    end

    for i in 1:length(q)
        p[i]=q[i]/somme_q
    end

    for i in 2:length(p)-1
        p[i]+=p[i-1]
    end
    p[length(p)]=1

end

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

function choix_alpha(p,prob)
    for i in 1:length(p)
        if prob < p[i]
            return(i)
        elseif i==length(p)
            return(i)
        end
    end
end

max_stars_repo_path: Archives/Triste vie 10-08/DM2_2.jl

max_issues_repo_path: Archives/Triste vie 10-08/DM2_2.jl

max_forks_repo_path: Archives/Triste vie 10-08/DM2_2.jl


text:
struct Lexicon
    words::Dict{String, Set{String}}
    tags::Set{String}
end

"""
    Lexicon()
    Lexicon(io_or_filename)
    Lexicon(words)

Construct a lexicon. It can be empty (no parameters) or created from some IO-like object or a
sequence/iterable of words.

A lexicon is a list of (known) words, each of which can be tagged with various tags (e.g.
indicating how it is known, etc.).
"""
function Lexicon end

Lexicon(words) = Lexicon(Dict(word => Set{String}() for word in words), Set{String}())
Lexicon() = Lexicon(Dict{String, Set{String}}(), Set{String}())

function Lexicon(io::Union{IO, AbstractString})
    lex = Lexicon()
    for line in eachline(io)
        char, tag_string = split(line, ":")
        tags = split(tag_string, ",")
        push!(lex, char; tags=tags)
    end
    lex
end

Base.length(lex::Lexicon) = length(lex.words)
Base.in(item, lex::Lexicon) = haskey(lex.words, item)
Base.getindex(lex::Lexicon, index) = getindex(lex.words, index)
Base.print(io::IO, lex::Lexicon) =
    for (word, tags) in lex.words
        println(io, "$word:$(join(tags, ","))")
    end

function Base.push!(lex::Lexicon, words...; tags=())
    for word in words
        haskey(lex.words, word) || (lex.words[word] = Set{String}())
        for tag in tags
            push!(lex.words[word], tag)
        end
    end
    for tag in tags
        push!(lex.tags, tag)
    end
    nothing
end

tagged_with(lex::Lexicon, tag) =
    [word for (word, tags) in lex.words if tag in tags]

"""
    coverage(lexicon, charfreq)
    coverage(lexicon, text)

Compute a lexicon's coverage of a text (possibly via a precomputed character
frequency dictionary). Both token and type coverage are provided.
"""
function coverage(lex::Lexicon, cf::Accumulator)
    known_tokens, total_tokens, known_types, total_types = 0, 0, 0, length(lex)
    for (char, freq) in cf
        if char in lex
            known_tokens += freq
            known_types += 1
        end
        total_tokens += freq
        total_types += 1
    end
    (char_coverage=Float(known_tokens)/total_tokens, type_coverage=Float(known_types)/total_types)
end
coverage(lex::Lexicon, text) = coverage(charfreq(text), lex)


text:
#=
Copyright (c) 2015, Intel Corporation
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
- Redistributions of source code must retain the above copyright notice,
  this list of conditions and the following disclaimer.
- Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
THE POSSIBILITY OF SUCH DAMAGE.
=#

module APILibTest

using ParallelAccelerator
#importall ParallelAccelerator.API.Lib

#using CompilerTools
#CompilerTools.ReadWriteSet.set_debug_level(3)
#CompilerTools.AstWalker.set_debug_level(3)
#CompilerTools.OptFramework.set_debug_level(3)
#ParallelAccelerator.set_debug_level(3)
#ParallelAccelerator.DomainIR.set_debug_level(3)
#ParallelAccelerator.ParallelIR.set_debug_level(3)
#ParallelAccelerator.CGen.set_debug_level(3)

@acc minMax(A) = (indmin(A), indmax(A))

@acc vecLen(A) = sqrt(sumabs2(A))

@acc diagTest(A) = sum(diag(diagm(A)))

@acc traceTest(A) = trace(diagm(A))

@acc testEye(A) = eye(A) .* A

@acc testRepMat(A) = sum(repmat(A, 2, 3), 1)

@acc function normalizeW(W::Array{Float64,2}, k::Int)
# normalize each column
#=
    for j = 1:k
        w = view(W,:,j)
        scale!(w, 1.0 / sqrt(sumabs2(w)))
	#s = sqrt(sumabs2(W[:,j]))
	#W[:,j] /= s
    end
=#
    s::Array{Float64,1} = [sqrt(sumabs2(W[:,j])) for j in 1:k]
    scale(W,s)
end

@acc function testRepMat2(A, B)
   (l, m, n) = size(B)
   C = repmat(A, 1, 1, n) ./ B
   D = 1.0 ./ maximum(C, 2)
   return D
end

function test1()
  A = rand(5)
  minMax(A) == (@noacc minMax(A)) && minMax(A) == (Base.indmin(A), Base.indmax(A))
end

function test2()
  A = rand(5)
  abs(vecLen(A) - (@noacc vecLen(A))) < 1.0e-10 &&
  abs(vecLen(A) - sqrt(Base.sumabs2(A))) < 1.0e-10
end

function test3()
  A = rand(5)
  abs(diagTest(A) - sum(A)) < 1.0e-10
end

function test4()
  A = rand(5)
  abs(diagTest(A) - sum(A)) < 1.0e-10
end

function test5()
  A = rand(5, 5)
  abs(sum(testEye(A) .- Base.eye(A) .* A)) < 1.0e-10
end

function test6()
  A = rand(4, 3)
  abs(sum(testRepMat(A) .- sum(repmat(A, 2, 3), 1))) < 1.0e-10
end

function test7()
    m = 100
    k = 100
    W  = Array(Float64, m, k)
    fill!(W, 3)
    W0 = normalizeW(W, k)
    fill!(W, 3)
    W1 = @noacc normalizeW(W, k)
    abs(sum(W0) - sum(W1)) < 1.0e-10
end

function test8()
    A = testRepMat2(ones(3,3), ones(3,3,3))
    abs((sum(A) - 9.0)) < 1.0e-10
end

end

using Base.Test
println("Testing parallel library functions...")
@test APILibTest.test1()
@test APILibTest.test2()
@test APILibTest.test3()
@test APILibTest.test4()
@test APILibTest.test5()
@test APILibTest.test6()
@test APILibTest.test7()
@test APILibTest.test8()
println("Done testing parallel library functions.")


text:
julia> maximum([1,3,3,7])
7

julia> maximum([pi,e+2/5,cos(6)/5,sqrt(91/10)])
3.141592653589793

julia> maximum([1,6,Inf])
Inf

julia> maximum(Float64[])
maximum: argument is empty
at In[138]:1
 in maximum at abstractarray.jl:1591

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.8.0-DEV
#     language: julia
#     name: julia-1.8
# ---

# %%
using Plots
using LinearAlgebra
using Fmt

# %%
f(a) = -a + √(a^2 - 2)
a = range(√2, 100, 1000)
plot(a, f; ylim=(-1.5, 0), label="")

# %%
A(a) = Float64[
     0  1
    -2 -2a
]

U(a, t) = exp(t*A(a))

# %%
@gif for a in [0:0.02:1.4; fill(1.4, 40); 1.4:0.04:5]
    t = range(0, 4π; length=1000)
    u = (t -> U(a, t)[1, 1]).(t)
    v = (t -> U(a, t)[2, 1]).(t)
    P = plot(t, [u v]; label=["u1" "v1"], ylim=(-1, 1))
    u = (t -> U(a, t)[1, 2]).(t)
    v = (t -> U(a, t)[2, 2]).(t)
    plot!(t, [u v]; label=["u2" "v2"], ylim=(-1, 1), ls=:dash)
    plot!(; title=f"a = {$a:.2f}")
end

# %%
for a in 0.1:0.1:2
    t = reverse(range(0, 10π; length=5001))
    n = accumulate(max, norm.(U.(a, t)))
    k = findfirst(t -> norm(U(a, t)) > 0.1, t)
    @show a, t[k], norm(U(a, t[k]))
end

# %%
t = range(0, 4π; length=1000)
u = (t -> U(t)[1, 1]).(t)
v = (t -> U(t)[2, 1]).(t)
plot(t, [u v])

# %%
t = range(0, 10; length=1000)
u = (t -> U(t)[1, 2]).(t)
v = (t -> U(t)[2, 2]).(t)
plot(t, [u v])

# %%
a, V = eigen(A)

# %%
V/0.40824829046386296

# %%
p = [
    -1 + 1im
    2
]

A*p - (-1 - 1im)*p

# %%
q = [
    -1 - 1im
    2
]

A*q - (-1 + 1im)*q

# %%
[minimum(real.(eigvals(Float64[0 1; -2 -2a]))) for a in 0:0.1:2]

# %%
using SymPy

# %%
@vars a

# %%
A = [
    0 1
    -2 -2
]

# %%
A^2//2

# %%
A^3//6

# %%
A^4//24

# %%
B = [
    0 1
    -2 -2a
]

# %%
B^2/factorial(2)

# %%
expand.(B^3)/factorial(3)

# %%
expand.(B^4)/factorial(4)

# %%

max_stars_repo_path: 0014/absorbing condition.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0014/absorbing condition.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0014/absorbing condition.jl

max_forks_repo_name: genkuroki/public


text:
#=
CenterStrike
Number Of Strikes up and down
expiry

gender

market,

Vanna Volga Parameters

=#

function StrikeLadder(CallPutFlag, MarketPrice, CenterStrike, T, VVInput)

LowerStrike = CenterStrike - 5
UpperStrike = CenterStrike + 5
StrikePrices = collect(LowerStrike:0.5:UpperStrike)

SensitivityMatrix = zeros(Float64, length(StrikePrices), 7)

SensitivityMatrix[:,1] = StrikePrices
SensitivityMatrix[:,2] = 0.25

for i in 1:length(StrikePrices) SensitivityMatrix[i,3] = Black76(CallPutFlag, MarketPrice, StrikePrices[i], T, 0.05, SensitivityMatrix[i,2]) end

for i in 1:length(StrikePrices) SensitivityMatrix[i,4] = VannaVolga(VVInput.putStrike, VVInput.ATMStrike, VVInput.callStrike, StrikePrices[i], CallPutFlag, T, 0.05, 0, VVInput.putVol, VVInput.ATMVol, VVInput.callVol)[2] end

for i in 1:length(StrikePrices) SensitivityMatrix[i,5] = Black76(CallPutFlag, MarketPrice, StrikePrices[i], T, 0.05, SensitivityMatrix[i,4]) end

for i in 1:length(StrikePrices) SensitivityMatrix[i,6] = Delta(CallPutFlag, MarketPrice, StrikePrices[i], T, 0.05, 0, SensitivityMatrix[i,4]) end

for i in 1:length(StrikePrices) SensitivityMatrix[i,7] = Gamma(CallPutFlag, MarketPrice, StrikePrices[i], T, 0.05, 0, SensitivityMatrix[i,4]) end

SensitivityMatrix = DataFrame(SensitivityMatrix)
names!(SensitivityMatrix, [:Strike, :HVol, :ThPrice, :implVol, :implPrem, :implDelta, :implGamma])

SensitivityMatrix

end

#VVSep = VVParams(10, 11, 12, 0.25, 0.21, 0.3)

#StrikeLadder("c", 11, 11, 1/12, VVSep)


text:
using Slurm
import Slurm: api
import Parameters: reconstruct

type Job
  desc :: api.job_desc_msg_t
  refs :: Dict{Symbol, Any}

  function Job(;kwargs...)
    fields = Dict{Symbol, Any}()
    refs = Dict{Symbol, Any}()
    for (key, value) in kwargs
      if key ∉ fieldnames(api.job_desc_msg_t)
        error("$key is not a valid parameter for a slurm job.")
      end

      if isa(value, String)
        refs[key] = value
        fields[key] = Base.unsafe_convert(Ptr{Cchar}, value)
      elseif !isbits(value)
        error("$key => $value is not a bitstype")
      else
        fields[key] = value
      end
    end
    r_desc = Ref{api.job_desc_msg_t}()
    api.slurm_init_job_desc_msg(r_desc)
    new(reconstruct(r_desc[], fields), refs)
  end
end

function Base.unsafe_convert(P::Type{Ptr{api.job_desc_msg_t}}, j :: Job)
  return convert(P, Base.data_pointer_from_objref(j))
end

willrun(j :: Job) = api.slurm_job_will_run(j) == 1
# slurm_job_will_run2

# slurm_allocate_resources
# slurm_confirm_allocation
# slurm_free_resource_allocation_response_msg

# slurm_allocation_lookup

# slurm_allocate_resources_and_run
# slurm_free_resource_allocation_and_run_response_msg

# slurm_submit_batch_job
# slurm_free_submit_response_response_msg


text:
## *
function grad!(dy::TAny, ::Val{1},
               op::Call{typeof(*), Tuple{TArray{T1,N1}, TArray{T2,N2}}}) where {T1,N1,T2,N2}
    x, y = op.args
    yt = record!(dy.tape, Call, transpose, (y,))
    return record!(dy.tape, Call, *, (dy, yt))
end

function grad!(dy::TAny, ::Val{2},
               op::Call{typeof(*), Tuple{TArray{T1,N1}, TArray{T2,N2}}}) where {T1,N1,T2,N2}
    x, y = op.args
    xt = record!(dy.tape, Call, transpose, (x,))
    return record!(dy.tape, Call, *, (xt, dy))
end

## * with scalars
function grad!(dy::TAny, ::Val{1},
               op::Call{typeof(*), Tuple{TReal, TArray{T,N}}}) where {T,N}
    x, y = op.args
    return sum(y .* dy)
end
function grad!(dy::TAny, ::Val{2},
               op::Call{typeof(*), Tuple{TReal, TArray{T,N}}}) where {T,N}
    x, y = op.args
    return x .* dy
end

function grad!(dy::TAny, ::Val{1},
               op::Call{typeof(*), Tuple{TArray{T,N}, TReal}}) where {T,N}
    x, y = op.args
    return y .* dy
end
function grad!(dy::TAny, ::Val{2},
               op::Call{typeof(*), Tuple{TArray{T,N}, TReal}}) where {T,N}
    x, y = op.args
    return sum(x .* dy)
end

## +
grad!(dy::TAny, ::Val{1}, op::Call{typeof(+), Tuple{TArray{T,N}, TArray{T,N}}}) where {T,N} = dy
grad!(dy::TAny, ::Val{2}, op::Call{typeof(+), Tuple{TArray{T,N}, TArray{T,N}}}) where {T,N} = dy

## -
grad!(dy::TAny, ::Val{1}, op::Call{typeof(-), Tuple{TArray{T,N}, TArray{T,N}}}) where {T,N} = dy
grad!(dy::TAny, ::Val{2}, op::Call{typeof(-), Tuple{TArray{T,N}, TArray{T,N}}}) where {T,N} = -dy

## sum, mean
function grad!(dy::TAny, ::Val{1}, op::Call{typeof(sum), Tuple{TArray{T,N}}}) where {T,N}
    return record!(dy.tape, Call, sum_grad, (op.args[1], dy); kwargs=op.kwargs)
end
function grad!(dy::TAny, ::Val{1}, op::Call{typeof(mean), Tuple{TArray{T,N}}}) where {T,N}
    return record!(dy.tape, Call, mean_grad, (op.args[1], dy); kwargs=op.kwargs)
end

## transpose
grad!(dy::TAny, ::Val{1}, op::Call{typeof(transpose), Tuple{TArray{T,N}}}) where {T,N} =
    transpose(dy)
# TODO: .==
# TODO: minimum, maximum

## getindex
grad!(dy::TAny, ::Val{1}, op::Call{typeof(getindex), Tuple{TArray{T,N}, TReal}}) where {T,N} =
    ungetindex(dy, op.args[1], op.args[2])
grad!(dy::TAny, ::Val{2}, op::Call{typeof(getindex), Tuple{TArray{T,N}, TReal}}) where {T,N} =
    record!(dy.tape, Constant, 0)

## BROADCASTING

getfirst(x::TAny) = first(x.val)
getfirst(x) = first(x)

function merge_bcast_as_call!(tape::Tape, minitape::Tape, collect_from::Int, var_st::Dict)
    for i=collect_from:length(minitape)
        op = minitape[i]
        if op isa Call
            tape_args = map(x -> tape[var_st[x.id]].var, op.args)
            var = record!(tape, Bcast, op.fn, tape_args)
        elseif op isa Constant
            var = record!(tape, Constant, op.var.val)
        elseif op isa Assign
            # sometimes src size and dest size are different, e.g. in rand(3,4) .+ rand(3)
            # but we handle it in specific methods of grad!
            src_id_in_tape = var_st[op.src.id]
            var = record!(tape, Assign, tape[src_id_in_tape].var)
        else
            error("Unexpected operation in gradient minitape: $op")
        end
        var_st[op.var.id] = var.id  # add to mapping of minitape IDs to tape IDs
    end
end

function grad!(dy::Any, ::Val{Idx}, op::Bcast{FnType, TT}) where {Idx,FnType,TT}
    # record derivatives of first elements onto a minitape
    minitape = Tape()
    el_args = map(arg -> record!(minitape, Input, getfirst(arg)), op.args)
    # run grad! for first elements on the minitape
    el_y = op.fn(el_args...)
    el_dy = record!(minitape, Constant, getfirst(dy))
    el_op = minitape[el_y.id]
    collect_from = length(minitape) + 1
    el_dx = grad!(el_dy, Val(Idx), el_op)
    # merge minitape into tape, replacing Calls with Bcasts and changing op ids
    # initial variable substitution table (var_st) should contain mappings from
    # minitape input vars to corresponding original vars
    var_st = Dict(el_a.id => a.id for (el_a, a) in zip(el_args, op.args))
    var_st[el_dy.id] = dy.id
    tape = dy.tape
    merge_bcast_as_call!(tape, minitape, collect_from, var_st)
    dx = tape[var_st[el_dx.id]].var
    return dx
end

# function grad!(dy::Any, ::Val{Idx},
#                op::Bcast{FnType, Tuple{TArray{T1,N1}, TArray{T2,N2}}}) where {Idx,FnType,T1,N1,T2,N2}
# ...
# end

## SPECIAL BROADCASTING

# handle the case of e.g. rand(3, 4) .+ rand(3, 1)
function grad!(dy::Any, ::Val{1},
               op::Bcast{typeof(+), Tuple{TArray{T1,2}, TArray{T2,2}}}) where {T1,T2}
    x, y = op.args
    xsz = size(x.val)
    ysz = size(y.val)
    if xsz == ysz
         return record!(dy.tape, Assign, dy)
    else
        if xsz[1] == 1 && ysz[1] != 1
            return sum(dy; dims=1)
        elseif xsz[2] == 1 && ysz[2] != 1
            return sum(dy; dims=2)
        else
            return record!(dy.tape, Assign, dy)
        end
    end
end
function grad!(dy::Any, ::Val{2},
               op::Bcast{typeof(+), Tuple{TArray{T1,2}, TArray{T2,2}}}) where {T1,T2}
    x, y = op.args
    xsz = size(x.val)
    ysz = size(y.val)
    if xsz == ysz
         return record!(dy.tape, Assign, dy)
    else
        if xsz[1] != 1 && ysz[1] == 1
            return sum(dy; dims=1)
        elseif xsz[2] != 1 && ysz[2] == 1
            return sum(dy; dims=2)
        else
            return record!(dy.tape, Assign, dy)
        end
    end
end

function grad!(dy::Any, ::Val{1},
               op::Bcast{typeof(+), Tuple{TArray{T1,1}, TArray{T2,2}}}) where {T1,T2}
    return dropdims(sum(dy; dims=2); dims=2)
end
function grad!(dy::Any, ::Val{2},
               op::Bcast{typeof(+), Tuple{TArray{T1,2}, TArray{T2,1}}}) where {T1,T2}
    return dropdims(sum(dy; dims=2); dims=2)
end

function grad!(dy::Any, ::Val{1}, op::Bcast{typeof(+), Tuple{TReal, TArray{T,N}}}) where {T,N}
    return sum(dy)
end
function grad!(dy::Any, ::Val{2}, op::Bcast{typeof(+), Tuple{TArray{T,N}, TReal}}) where {T,N}
    return sum(dy)
end

function grad!(dy::Any, ::Val{1}, op::Bcast{typeof(-), Tuple{TReal, TArray{T,N}}}) where {T,N}
    return sum(dy)
end
function grad!(dy::Any, ::Val{2}, op::Bcast{typeof(-), Tuple{TArray{T,N}, TReal}}) where {T,N}
    return -sum(dy)
end

@require CuArrays="3a865a2d-5b23-5a0f-bc46-62713ec82fae" begin
    import CUDAnative

    function grad!(dy::Any, ::Val{1}, op::Bcast{typeof(CUDAnative.exp),
                                                Tuple{TArray{T,N}}}) where {T,N}
        w = record!(dy.tape, Bcast, CUDAnative.exp, (op.args[1],))
        return record!(dy.tape, Bcast, *, (w,))
    end

    function grad!(dy::Any, ::Val{1}, op::Bcast{typeof(CUDAnative.log),
                                                Tuple{TArray{T,N}}}) where {T,N}
        return record!(dy.tape, Bcast, /, (dy, op.args[1],))
    end

    function grad!(dy::Any, ::Val{1}, op::Bcast{typeof(CUDAnative.sqrt),
                                                Tuple{TArray{T,N}}}) where {T,N}
        # 0.5 .* x .^ (-0.5) .* dy
        w = record!(dy.tape, Bcast, CUDAnative.pow, (op.args[1], constant(dy.tape, -0.5)))
        return 0.5 .* w .* dy
    end

    function grad!(dy::Any, ::Val{1}, op::Bcast{typeof(CUDAnative.pow),
                                                Tuple{TArray{T,N}, R}}) where {T,N,R<:Real}
        # y .* x .^ (y-1) .* dy
        x, y = op.args
        w1 = record!(dy.tape, Bcast, CUDAnative.pow, (x, y-1.0f0))
        w2 = record!(dy.tape, Bcast, *, (y, w1))
        w3 = record!(dy.tape, Bcast, *, (w2, dy))
        return w3
    end

end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "normaliz"
version = v"3.8.5"

# Collection of sources required to build normaliz
sources = [
    ArchiveSource("https://github.com/Normaliz/Normaliz/releases/download/v$version/normaliz-$version.tar.gz",
                  "cf4fdaaa6ffcd8d268b1f16dd4b64cf86f1eab55177e611f8ef672e7365435a0")
]

# Bash recipe for building across all platforms
script = raw"""
cd normaliz-*
# avoid libtool problems
rm "${prefix}/lib/libgmpxx.la"
./configure --prefix=$prefix --host=$target --build=${MACHTYPE} --with-flint=$prefix --with-nauty=$prefix --with-gmp=$prefix CPPFLAGS=-I$prefix/include LDFLAGS=-L${libdir}
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
# windows build would require MPIR instead of GMP for 'long long'
platforms = filter(!Sys.iswindows, supported_platforms())
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("libnormaliz", :libnormaliz),
    ExecutableProduct("normaliz", :normaliz)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="GMP_jll", uuid="781609d7-10c4-51f6-84f2-b8444358ff6d")),
    Dependency(PackageSpec(name="FLINT_jll", uuid="e134572f-a0d5-539d-bddf-3cad8db41a82")),
    Dependency(PackageSpec(name="nauty_jll", uuid="55c6dc9b-343a-50ca-8ff2-b71adb3733d5")),
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version=v"6")

max_stars_repo_name: devmotion/Yggdrasil

max_issues_repo_name: devmotion/Yggdrasil

max_forks_repo_name: devmotion/Yggdrasil


text:
using ProjAssistant
@quickactivate

# ----------------------------------------------------------------------------
# ARGS
using ArgParse

set = ArgParseSettings()
@add_arg_table! set begin
    "--dry-run"
        help = "do all but change nothing on disk"
        action = :store_true
    "--redo"
        help = "ignore disk version (if exist)"
        action = :store_true
end

if isinteractive()
    # Dev values
    dry_run = false
    redo = true
else
    parsed_args = parse_args(set)
    dry_run = parsed_args["dry-run"]
    redo = parsed_args["redo"]
end

# ----------------------------------------------------------------------------
@time begin
    import Chemostat_EColi
    const ChE = Chemostat_EColi

    import Chemostat_Kayser2005
    const ChK = Chemostat_Kayser2005
    const KiJR = ChK.iJR904
    const Kd = ChK.KayserData

    import Chemostat_Nanchen2006
    const ChN = Chemostat_Nanchen2006
    const NiJR = ChN.iJR904
    const Nd = ChN.NanchenData

    import Chemostat_Folsom2014
    const ChF = Chemostat_Folsom2014
    const FiJR = ChF.iJR904
    const Fd = ChF.FolsomData

    import ChemostatPlots
    const ChP = ChemostatPlots

    import Chemostat
    import Chemostat.LP.MathProgBase
    const Ch = Chemostat
    const ChU = Ch.Utils
    const ChLP = Ch.LP
    const ChSS = Ch.SteadyState
    using ProgressMeter
    using Base.Threads
    using SparseArrays

    import SimTools
    const SimT = SimTools

    using Base.Threads
    using Serialization

    using MagicDicts

    using Statistics
end

# ----------------------------------------------------------------------------
function fix_Z_and_readj_model(iJR, Ed, exp;
        btol = 0.0,
        verbose = false, adjth = 0.001,
        proc_fun = (model) -> nothing
    )

    model = iJR.load_model("fva_models", exp) |> deepcopy
    biomider = iJR.BIOMASS_IDER
    biomidx = ChU.rxnindex(model, biomider)
    exglcider = iJR.EX_GLC_IDER

    # ----------------------------------------------------
    # fix z
    exp_growth = Ed.val("D", exp)

    verbose && @info("Before re-adjustement $("-"^50)")
    dflx = abs(exp_growth * btol)
    ChU.bounds!(model, biomider, exp_growth - dflx, exp_growth + dflx)

    # ----------------------------------------------------
    # find max Yzvg
    vgL, vgU = ChLP.fva(model, exglcider) .|> first
    min_vg = min(abs(vgL), abs(vgU)) # vg is negative
    verbose && @show min_vg

    exp_Yzvg = abs(Ed.val("D", exp) / Ed.val("uGLC", exp))
    model_Yzvg = abs(exp_growth/ min_vg)
    verbose && @show exp_Yzvg, model_Yzvg
    verbose && println()

    # ----------------------------------------------------
    # readjust biom reaction
    Sfac = 0.98
    last_Sfac = Sfac # in case of error
    biomS = collect(model.S[:, biomidx]) # backup
    x0 = Sfac
    x1 = 1.0
    gdth = adjth
    minΔx = 0.01
    maxΔx = 0.05
    target = exp_Yzvg
    gdmaxiter = 1000

    function readjust_S!(gdmodel)
        last_Sfac = Sfac
        Sfac = SimT.gd_value(gdmodel)

        verbose && println()
        verbose && @info("After re-adjustement $("-"^50)", Sfac)
        try
            # adjust biom
            model.S[:, biomidx] .= biomS .* Sfac

            # find max Yzvg
            vgL, vgU = ChLP.fva(model, exglcider) .|> first
            verbose && @show vgL, vgU

            min_vg = min(abs(vgL), abs(vgU)) # vg is negative
            verbose && @show min_vg

            model_Yzvg = abs(exp_growth/ min_vg)
            verbose &&  @show exp_Yzvg, model_Yzvg

            ret = proc_fun(model)
            ret === true &&

            # break
            last_Sfac == Sfac && return target

            return model_Yzvg
        catch err
            verbose && println()
            verbose && @warn("Error, using previous last_Sfac",
                last_Sfac, Sfac,
                (vgL, vgU),
                (exp_Yzvg, model_Yzvg),
                err
            )

            # use last one
            Sfac = last_Sfac
            model.S[:, biomidx] .= biomS .* Sfac

            return target
        end
    end

    gdmodel = SimT.grad_desc(readjust_S!;
        x0, x1, gdth, minΔx, maxΔx,
        target, maxiter = gdmaxiter,
        verbose = false
    )

    proc_fun(model)

    verbose && @show Sfac
    return Sfac
end;

## ----------------------------------------------------------------------------
const wlk = ReentrantLock()
thsv_set!(LPDAT, ks, v) = lock(() -> (LPDAT[ks...] = v), wlk)

## ----------------------------------------------------------------------------
function do_LP!(LPDAT, iJR, Ed, exp)

    biomider = iJR.BIOMASS_IDER
    costider = iJR.COST_IDER
    exglcider = iJR.EX_GLC_IDER
    atpmider = iJR.ATPM_IDER
    max_sense = ChLP.MAX_SENSE
    min_sense = ChLP.MIN_SENSE
    btol = 0.0

    # FBA_Z_FIX_MAX_ATP_MINMAX_COST
    let
        for (FBAmethod, atpm_sense, cost_sense) in [
            (:FBA_Z_FIX_MAX_ATP_MIN_COST, max_sense, min_sense),
            (:FBA_Z_FIX_MAX_ATP_MAX_COST, max_sense, max_sense),
        ]

            !redo && haskey(LPDAT, [FBAmethod, :model, exp]...) && continue

            model = iJR.load_model("base_model")
            exp_xi = Ed.val(:xi, exp)
            exp_intake_info = iJR.intake_info(exp)
            ChSS.apply_bound!(model, exp_xi, exp_intake_info; emptyfirst = true)

            # fix Z
            exp_growth = Ed.val("D", exp)
            dflx = abs(exp_growth * btol)
            ChU.bounds!(model, biomider, exp_growth - dflx, exp_growth + dflx)

            # max atpm
            # rxn[139]: ATPM (ATP maintenance requirement)
            # lb: 7.6, ub: 1000.0
            # (-1.0) atp_c + (-1.0) h2o_c ==> (1.0) adp_c + (1.0) h_c + (1.0) pi_c
            ChU.ub!(model, atpmider, iJR.ABS_MAX_BOUND)
            fbaout1 = ChLP.fba(model, atpmider; sense = atpm_sense)
            atpm = ChU.av(model, fbaout1, atpmider)
            dflx = abs(atpm * btol)
            ChU.bounds!(model, atpmider, atpm - dflx, atpm + dflx)

            # min cost
            fbaout = ChLP.fba(model, costider; sense = cost_sense)

            thsv_set!(LPDAT, [FBAmethod, :model, exp], ChU.compressed_model(model))
            thsv_set!(LPDAT, [FBAmethod, :fbaout, exp], fbaout)
        end
    end

    # FBA_Z_FIX_ADJ_MINMAX_VG_MINMAX_COST
    let
        for (FBAmethod, glc_sense, cost_sense) in [
            (:FBA_Z_FIX_ADJ_MIN_VG_MIN_COST, max_sense, min_sense),
            (:FBA_Z_FIX_ADJ_MIN_VG_MAX_COST, max_sense, max_sense),
            (:FBA_Z_FIX_ADJ_MAX_VG_MIN_COST, min_sense, min_sense),
            (:FBA_Z_FIX_ADJ_MAX_VG_MAX_COST, min_sense, max_sense),
        ]

            !redo && haskey(LPDAT, [FBAmethod, :model, exp]...) && continue

            # fix Z and readjust biom so exp and model yield (z/vg) match
            model = nothing
            fbaout = nothing
            function proc_fun(model_)
                model = model_

                # min vg (vg is negative)
                fbaout1 = ChLP.fba(model, exglcider; sense = glc_sense)
                exglc = ChU.av(model, fbaout1, exglcider)
                ChU.fixxing(model, exglcider, exglc; btol) do
                    # min cost
                    fbaout = ChLP.fba(model, costider; sense = cost_sense)
                end
            end
            # fix Z and readjust biom so exp and model yield (z/vg) match
            Sfac = fix_Z_and_readj_model(iJR, Ed, exp; proc_fun, verbose = false);

            # store
            thsv_set!(LPDAT, [FBAmethod, :Sfac, exp], Sfac)
            thsv_set!(LPDAT, [FBAmethod, :model, exp], ChU.compressed_model(model))
            thsv_set!(LPDAT, [FBAmethod, :fbaout, exp], fbaout)
        end
    end

    # FBA_Z_FIX_MINMAX_VG_MINMAX_COST
    let
        for (FBAmethod, glc_sense, cost_sense) in [
            (:FBA_Z_FIX_MIN_VG_MIN_COST, max_sense, min_sense),
            (:FBA_Z_FIX_MIN_VG_MAX_COST, max_sense, max_sense),
            (:FBA_Z_FIX_MAX_VG_MIN_COST, min_sense, min_sense),
            (:FBA_Z_FIX_MAX_VG_MAX_COST, min_sense, max_sense),
        ]

            !redo && haskey(LPDAT, [FBAmethod, :model, exp]...) && continue

            model = iJR.load_model("fva_models", exp)

            # fix Z
            exp_growth = Ed.val("D", exp)
            dflx = abs(exp_growth * btol)
            ChU.bounds!(model, biomider, exp_growth - dflx, exp_growth + dflx)

            # min vg (vg is negative)
            fbaout1 = ChLP.fba(model, exglcider; sense = glc_sense)
            exglc = ChU.av(model, fbaout1, exglcider)
            dflx = abs(exglc * btol)
            ChU.bounds!(model, exglcider, exglc - dflx, exglc + dflx)

            # min cost
            fbaout = ChLP.fba(model, costider; sense = cost_sense)

            thsv_set!(LPDAT, [FBAmethod, :model, exp], ChU.compressed_model(model))
            thsv_set!(LPDAT, [FBAmethod, :fbaout, exp], fbaout)
        end
    end

    # FBA_MAX_Z_MINMAX_COST
    let
        for (FBAmethod, biom_sense, cost_sense) in [
            (:FBA_MAX_Z_MIN_COST, max_sense, min_sense),
            (:FBA_MAX_Z_MAX_COST, max_sense, max_sense),
        ]
            !redo && haskey(LPDAT, [FBAmethod, :model, exp]...) && continue

            model = iJR.load_model("fva_models", exp)

            # max z
            fbaout1 = ChLP.fba(model, biomider; sense = biom_sense)
            objval = ChU.av(model, fbaout1, biomider)
            dflx = abs(objval * btol)
            ChU.bounds!(model, biomider, objval - dflx, objval + dflx)

            # min cost
            fbaout = ChLP.fba(model, costider; sense = cost_sense)

            thsv_set!(LPDAT, [FBAmethod, :model, exp], ChU.compressed_model(model))
            thsv_set!(LPDAT, [FBAmethod, :fbaout, exp], fbaout)
        end
    end

    # FBA_Z_FIX_MINMAX_COST
    let
        for (FBAmethod, cost_sense) in [
            (:FBA_Z_FIX_MIN_COST, min_sense),
            (:FBA_Z_FIX_MAX_COST, max_sense),
        ]
            !redo && haskey(LPDAT, [FBAmethod, :model, exp]...) && continue

            model = iJR.load_model("fva_models", exp)
            exp_growth = Ed.val("D", exp)

            # fix Z
            exp_growth = Ed.val("D", exp)
            dflx = abs(exp_growth * btol)
            ChU.bounds!(model, biomider, exp_growth - dflx, exp_growth + dflx)

            # min cost
            fbaout = ChLP.fba(model, costider; sense = cost_sense)

            thsv_set!(LPDAT, [FBAmethod, :model, exp], ChU.compressed_model(model))
            thsv_set!(LPDAT, [FBAmethod, :fbaout, exp], fbaout)
        end
    end

    # FBA_Z_VG_FIX_MINMAX_COST
    let
        for (FBAmethod, cost_sense) in [
            (:FBA_Z_VG_FIX_MIN_COST, min_sense),
            (:FBA_Z_VG_FIX_MAX_COST, max_sense),
        ]
            !redo && haskey(LPDAT, [FBAmethod, :model, exp]...) && continue

            model = iJR.load_model("fva_models", exp)

            # fix Z
            exp_growth = Ed.val("D", exp)
            dflx = abs(exp_growth * btol)
            ChU.bounds!(model, biomider, exp_growth - dflx, exp_growth + dflx)

            # fix vg
            exp_exglc = -abs(Ed.uval("GLC", exp))
            dflx = abs(exp_exglc * btol)
            ChU.bounds!(model, exglcider, exp_exglc - dflx, exp_exglc + dflx)

            # min cost
            fbaout = ChLP.fba(model, costider; sense = cost_sense)

            thsv_set!(LPDAT, [FBAmethod, :model, exp], ChU.compressed_model(model))
            thsv_set!(LPDAT, [FBAmethod, :fbaout, exp], fbaout)
        end
    end
    return LPDAT
end

## ----------------------------------------------------------------------------------
# create (iJR, Ed, exp) collection
function collect_experiments()
   collprod((iJR, Ed, EXPS)) = vec(collect(Iterators.product([iJR], [Ed], EXPS)))
   mod_dat = [
        (NiJR, Nd, 1:17),
        (KiJR, Kd, 5:13),
        (FiJR, Fd, 1:4),
    ]
    return vcat(map(collprod, mod_dat)...)
end

## ----------------------------------------------------------------------------------
let

    EXP_COLL = collect_experiments()
    PROGRESS = Dict() # to now when finished
    isfinished(iJR) = count((dat) -> isequal(first(dat), iJR), EXP_COLL) == PROGRESS[iJR]

    Ch = Channel(nthreads()) do Ch_
        for (iJR, Ed, exp) in EXP_COLL
            put!(Ch_, (iJR, Ed, exp))
        end
    end

    LPDATs = Dict()
    # @threads
    for _ in 1:nthreads()
        thid = threadid()

        for (iJR, Ed, exp) in Ch

            src = string(nameof(Ed))
            dfile = ChE.LP_DATfile(src)

            lock(wlk) do
                # init
                if !haskey(LPDATs, src)
                    LPDATs[src] = isfile(dfile) ? ldat(dfile) : MagicDict()
                end
                get!(PROGRESS, iJR, 0)

                # info
                @info(string("Doing ", "-"^30),
                    exp, src, thid
                ); println()
            end

            # --------------------------------------------
            do_LP!(LPDATs[src], iJR, Ed, exp)

            # --------------------------------------------
            lock(wlk) do
                # progress
                PROGRESS[iJR] += 1

                # finish
                if isfinished(iJR)
                    @info(string("Finished ", "-"^30),
                        src, thid
                    ); println()

                    # saving
                    dry_run || sdat(LPDATs[src], dfile)

                    # save memory
                    LPDATs[src] = nothing
                    GC.gc()
                end
            end
        end

    end # for thid
end;


text:
include("experiment_helpers.jl")

# Create Roadway
roadway_opts = MergingRoadwayOptions(lane_width=0.26)
roadway = build_roadway(roadway_opts)

# Create players
T = Float64
p = 2
model = DoubleIntegratorGame(p=p)
n = model.n
m = model.m

# Define the horizon of the problem
N = 20 # N time steps
dt = 0.1 # each step lasts 0.1 second
probsize = ProblemSize(N,model) # Structure holding the relevant sizes of the problem

# Define the objective of each player
# We use a LQR cost
Q = [Diagonal(3*ones(SVector{model.ni[i],T})) for i=1:p] # Quadratic state cost
R = [Diagonal(0.3*ones(SVector{model.mi[i],T})) for i=1:p] # Quadratic control cost
# Desrired state
xf = [SVector{model.ni[1],T}([1.0,+0.35,0.5,0]),
      SVector{model.ni[2],T}([1.0,-0.05,0.5,0]),
      ]
# Desired control
uf = [zeros(SVector{model.mi[i],T}) for i=1:p]
# Objectives of the game
game_obj = Algames.GameObjective(Q,R,xf,uf,N,model)
radius = 0.2*ones(p)
μ = 100.0*ones(p)
add_collision_cost!(game_obj, radius, μ)

# Define the constraints that each player must respect
game_con = Algames.GameConstraintValues(probsize)

# Define the initial state of the system
x0 = SVector{model.n,T}([
    0.2,   0.2,
   -0.45, -0.05,
    0.4,   0.5,
    0.1,   0.0,
    ])

# Define the Options of the solver
opts = Options()
opts.ls_iter = 15
opts.outer_iter = 20
opts.inner_iter = 20
opts.ρ_0 = 1e0
opts.reg_0 = 1e-8
opts.α_dual = 0.0
opts.λ_max = 1.0*1e7
opts.ϵ_dyn = 1e-6
opts.ϵ_sta = 1e-6
opts.ϵ_con = 1e-6
opts.ϵ_opt = 1e-6
opts.regularize = true
# Define the game problem
prob = Algames.GameProblem(N,dt,x0,model,opts,game_obj,game_con)

# Solve the problem
newton_solve!(prob)

plot!(prob.model, prob.pdtraj.pr)
plot!(prob.stats)

players = Vector{Player{T}}(undef, p)
players[1] = Player(model, roadway.lane[3])
players[2] = Player(model, roadway.lane[6])

# Create Scenario
sce = Scenario(model, roadway, players)

function get_figure(prob::Algames.GameProblem, sce::Scenario, vis_sol::Visualizer, vis_init::Visualizer)
    N = prob.probsize.N
    model = prob.model

    prob.opts.shift = 2^10
    prob.opts.α_dual = 0.0
    prob.opts.αx_dual = zeros(p)
    prob.opts.f_init = x -> (rand(x) .- 0.5)
    prob.opts.amplitude_init = 1.0

    # Set tolerances
    prob.opts.ϵ_opt = 1e-6
    prob.opts.ϵ_dyn = 1e-6
    prob.opts.ϵ_con = 1e-6
    prob.opts.ϵ_sta = 1e-6

    init = []
    sol = []
    plt_init = plot()
    plt_sol = plot()
    for k = 1:20
        # Randomize the initialization
        prob.opts.seed = k

        # Get the initial guess (random controls + dynamics rollout)
        prob.opts.outer_iter = 0
        newton_solve!(prob)
        init_traj = deepcopy(prob.pdtraj)

        # Get solution trajectory
        prob.opts.outer_iter = 20
        newton_solve!(prob)
        sol_traj = deepcopy(prob.pdtraj)

        # Only keep the solution that converged
        opt = mean(abs.(prob.stats.opt_vio[end].vio))
        sta = mean(abs.(prob.stats.sta_vio[end].vio))
        dyn = mean(abs.(prob.stats.dyn_vio[end].vio))
        con = mean(abs.(prob.stats.con_vio[end].vio))
        opts = prob.opts
        if opt < opts.ϵ_opt && sta < opts.ϵ_sta && dyn < opts.ϵ_dyn && con < opts.ϵ_con
            build_waypoint!(vis_init, sce.player, N, key=k)
            build_waypoint!(vis_sol,  sce.player, N, key=k)
            set_waypoint_traj!(vis_init, model, sce, init_traj.pr, key=k)
            set_waypoint_traj!(vis_sol,  model, sce, sol_traj.pr, key=k)

            plot!(model, init_traj.pr, plt=plt_init)
            plot!(model, sol_traj.pr, plt=plt_sol)
            push!(sol, sol_traj)
            push!(init, init_traj)
        end
    end
    set_state!(vis_init, model, sce, Algames.state(init[1].pr[1]))
    set_state!(vis_sol,  model, sce, Algames.state(sol[1].pr[1]))
    return init, sol
end

# Initialize visualizers
vis_init = Visualizer()
vis_sol = Visualizer()
open(vis_init)
open(vis_sol)

# Visualize trajectories
set_scenario!(vis_init, sce)
set_scenario!(vis_sol, sce)
set_env!(vis_init, VehicleState(1.0, 0.0, 0.0, 0.0))
set_env!(vis_sol, VehicleState(1.0, 0.0, 0.0, 0.0))
set_camera_birdseye!(vis_init, height=4.0)
set_camera_birdseye!(vis_sol, height=4.0)

init_traj, sol_traj = get_figure(prob, sce, vis_sol, vis_init)


text:
using BeamformerRecipes
using Test

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


text:
module CassetteIntroCar

using ..Ahorn, Maple

@mapdef Entity "brokemiahelper/cassetteIntroCar" CarCassette(x::Integer, y::Integer, index::Integer=0, tempo::Number=1)

colorNames = Dict{String, Int}(
    "Blue" => 0,
    "Rose" => 1,
    "Bright Sun" => 2,
    "Malachite" => 3
)

const placements = Ahorn.PlacementDict(
    "Cassette Intro Car ($index - $color)" => Ahorn.EntityPlacement(
        CarCassette,
        "point",
        Dict{String, Any}(
            "index" => index,
        )
    ) for (color, index) in colorNames
)

Ahorn.editingOptions(entity::CarCassette) = Dict{String, Any}(
    "index" => colorNames
)

colors = Dict{Integer, Ahorn.colorTupleType}(
    1 => (240, 73, 190, 255) ./ 255,
	2 => (252, 220, 58, 255) ./ 255,
	3 => (56, 224, 78, 255) ./ 255,
)

defaultColor = (73, 170, 240, 255) ./ 255
borderMultiplier = (0.9, 0.9, 0.9, 1)

barrierSprite = "scenery/car/barrier"
bodySprite = "scenery/car/body"
pavementSprite = "scenery/car/pavement"
wheelsSprite = "scenery/car/wheels"

function Ahorn.selection(entity::CarCassette)
    x, y = Ahorn.position(entity)

    hasRoadAndBarriers = get(entity.data, "hasRoadAndBarriers", false)

    rectangles = Ahorn.Rectangle[
        Ahorn.getSpriteRectangle(bodySprite, x, y, jx=0.5, jy=1.0),
        Ahorn.getSpriteRectangle(wheelsSprite, x, y, jx=0.5, jy=1.0),
    ]

    if hasRoadAndBarriers
        push!(rectangles, Ahorn.getSpriteRectangle(barrierSprite, x + 32, y, jx=0.0, jy=1.0))
        push!(rectangles, Ahorn.getSpriteRectangle(barrierSprite, x + 41, y, jx=0.0, jy=1.0))
    end

    return Ahorn.coverRectangles(rectangles)
end

function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::CarCassette, room::Maple.Room)
    x, y = Ahorn.position(entity)

    hasRoadAndBarriers = get(entity.data, "hasRoadAndBarriers", false)
    rng = Ahorn.getSimpleEntityRng(entity)

    pavementWidth = x - 48
    columns = floor(Int, pavementWidth / 8)

    Ahorn.drawSprite(ctx, wheelsSprite, x, y, jx=0.5, jy=1.0, tint=get(colors, Int(get(entity.data, "index", 0)), defaultColor))
    Ahorn.drawSprite(ctx, bodySprite, x, y, jx=0.5, jy=1.0, tint=get(colors, Int(get(entity.data, "index", 0)), defaultColor))

    if hasRoadAndBarriers
        Ahorn.drawSprite(ctx, barrierSprite, x + 32, y, jx=0.0, jy=1.0, tint=get(colors, Int(get(entity.data, "index", 0)), defaultColor))
        Ahorn.drawSprite(ctx, barrierSprite, x + 41, y, jx=0.0, jy=1.0, tint=get(colors, Int(get(entity.data, "index", 0)), defaultColor))

        for col in 0:columns - 1
            choice = col >= columns - 2 ? (col != columns - 2 ? 3 : 2) : rand(rng, 0:2)

            Ahorn.drawImage(ctx, pavementSprite, col * 8, y, choice * 8, 0, 8, 4, tint=get(colors, Int(get(entity.data, "index", 0)), defaultColor))
        end
    end
end

end

max_stars_repo_name: Brokemia/BrokemiaHelper

max_issues_repo_name: Brokemia/BrokemiaHelper

max_forks_repo_name: Brokemia/BrokemiaHelper


text:
module TestUtilities

# using Revise
using Test
using MLJBase

train, test = partition(1:100, 0.9)
@test collect(train) == collect(1:90)
@test collect(test) == collect(91:100)
train, test = partition(1:100, 0.9, shuffle=true)
@test length(train) == 90

y=[:a,:b,:a,:b,:c,:c]
part1,part2=collect(StratifiedKFold(y,2))
@test sort(y[part1])==sort(y[part2])
end # module
true


text:
"""
Basic Document-Term-Matrix (DTM) type.

# Fields
  * `dtm::SparseMatriCSC{T,Int}` the actual DTM; rows represent terms
and columns represent documents
  * `terms::Vector{String}` a list of terms that represent the lexicon of
the corpus associated with the DTM
  * `row_indices::OrderedDict{String, Int}` a map between the `terms` and the
rows of the `dtm`
"""
struct DocumentTermMatrix{T}
    dtm::SparseMatrixCSC{T, Int}
    terms::Vector{String}
    row_indices::OrderedDict{String, Int}
end

"""
    rowindices(terms)

Returns a dictionary that maps each term from the vector `terms`
to a integer idex.
"""
rowindices(terms::Vector{String}) =
    OrderedDict{String, Int}(term => i for (i,term) in enumerate(terms))

"""
    columnindices(terms)

Identical to `rowindices`. Returns a dictionary that maps
each term from the vector `terms` to a integer idex.
"""
columnindices = rowindices

"""
    DocumentTermMatrix{T}(crps::Corpus [,terms] [; ngram_complexity=DEFAULT_NGRAM_COMPLEXITY, tokenizer=DEFAULT_TOKENIZER])

Auxiliary constructor(s) of the `DocumentTermMatrix` type. The type `T` has to be
a subtype of `Real`. The constructor(s) requires a corpus `crps` and
a `terms` structure representing the lexicon of the corpus. The latter
can be a `Vector{String}`, an `AbstractDict` where the keys are the lexicon, or can
be missing, in which case the `lexicon` field of the corpus is used.
"""
function DocumentTermMatrix{T}(crps::Corpus,
                               terms::Vector{String};
                               ngram_complexity::Int=DEFAULT_NGRAM_COMPLEXITY,
                               tokenizer::Symbol=DEFAULT_TOKENIZER) where {T<:Real}
    row_indices = rowindices(terms)
    m = length(terms)
    n = length(crps)
    rows = Vector{Int}(undef, 0)
    columns = Vector{Int}(undef, 0)
    values = Vector{T}(undef, 0)
    for (i, doc) in enumerate(crps)
        ngs = ngrams(doc, ngram_complexity, tokenizer=tokenizer)
        for ngram in keys(ngs)
            j = get(row_indices, ngram, 0)
            v = ngs[ngram]
            if j != 0
                push!(columns, i)
                push!(rows, j)
                push!(values, v)
            end
        end
    end
    if length(rows) > 0
        dtm = sparse(rows, columns, values, m, n)
    else
        dtm = spzeros(T, m, n)
    end
    return DocumentTermMatrix(dtm, terms, row_indices)
end

DocumentTermMatrix(crps::Corpus,
                   terms::Vector{String};
                   ngram_complexity::Int=DEFAULT_NGRAM_COMPLEXITY,
                   tokenizer::Symbol=DEFAULT_TOKENIZER) =
    DocumentTermMatrix{DEFAULT_DTM_TYPE}(crps, terms, ngram_complexity=ngram_complexity, tokenizer=tokenizer)

DocumentTermMatrix{T}(crps::Corpus,
                      lex::AbstractDict;
                      ngram_complexity::Int=DEFAULT_NGRAM_COMPLEXITY,
                      tokenizer::Symbol=DEFAULT_TOKENIZER) where {T<:Real} =
    DocumentTermMatrix{T}(crps, collect(keys(lex)), ngram_complexity=ngram_complexity, tokenizer=tokenizer)

DocumentTermMatrix(crps::Corpus,
                   lex::AbstractDict;
                   ngram_complexity::Int=DEFAULT_NGRAM_COMPLEXITY,
                   tokenizer::Symbol=DEFAULT_TOKENIZER) =
    DocumentTermMatrix{DEFAULT_DTM_TYPE}(crps, lex, ngram_complexity=ngram_complexity, tokenizer=tokenizer)

DocumentTermMatrix{T}(crps::Corpus;
                      ngram_complexity::Int=DEFAULT_NGRAM_COMPLEXITY,
                      tokenizer::Symbol=DEFAULT_TOKENIZER) where {T<:Real} = begin
    DocumentTermMatrix{T}(crps, create_lexicon(crps, ngram_complexity),
                          ngram_complexity=ngram_complexity, tokenizer=tokenizer)
end

DocumentTermMatrix(crps::Corpus;
                   ngram_complexity::Int=DEFAULT_NGRAM_COMPLEXITY,
                   tokenizer::Symbol=DEFAULT_TOKENIZER) = begin
    DocumentTermMatrix{DEFAULT_DTM_TYPE}(crps, create_lexicon(crps, ngram_complexity),
                                         ngram_complexity=ngram_complexity, tokenizer=tokenizer)
end

DocumentTermMatrix(dtm::SparseMatrixCSC{T, Int}, terms::Vector{String}) where {T<:Real} =
    DocumentTermMatrix(dtm, terms, rowindices(terms))

"""
    dtm(d::DocumentTermMatrix)

Access the matrix of a `DocumentTermMatrix` `d`.
"""
dtm(d::DocumentTermMatrix) = d.dtm

"""
    dtm(crps::Corpus, eltype::Type{T}=DEFAULT_DTM_TYPE [; ngram_complexity=DEFAULT_NGRAM_COMPLEXITY, tokenizer=DEFAULT_TOKENIZER])

Access the matrix of the DTM associated with the corpus `crps`. The
`DocumentTermMatrix{T}` will first have to be created in order for
the actual matrix to be accessed.
"""
dtm(crps::Corpus,
    eltype::Type{T}=DEFAULT_DTM_TYPE;
    ngram_complexity::Int=DEFAULT_NGRAM_COMPLEXITY,
    tokenizer::Symbol=DEFAULT_TOKENIZER
   ) where {T<:Real} =
    dtm(DocumentTermMatrix{T}(crps, ngram_complexity=ngram_complexity, tokenizer=tokenizer))

# Produce the signature of a DTM entry for a document
function dtm_entries(d,
                     lex::OrderedDict{String, Int},
                     eltype::Type{T}=DEFAULT_DTM_TYPE;
                     ngram_complexity::Int=DEFAULT_NGRAM_COMPLEXITY,
                     tokenizer::Symbol=DEFAULT_TOKENIZER,
                     lex_is_row_indices::Bool=false) where {T<:Real}
    ngs = ngrams(d, ngram_complexity, tokenizer=tokenizer)
    indices = Vector{Int}(undef, 0)
    values = Vector{T}(undef, 0)
    local row_indices
    if lex_is_row_indices
        row_indices = lex
    else
        terms = collect(keys(lex))
        row_indices = rowindices(terms)
    end
    for ngram in keys(ngs)
        j = get(row_indices, ngram, 0)
        v = ngs[ngram]
        if j != 0
            push!(indices, j)
            push!(values, v)
        end
    end
    return (indices, values)
end

"""
    dtv(d, lex::OrderedDict{String,Int}, eltype::Type{T}=DEFAULT_DTM_TYPE [; ngram_complexity=DEFAULT_NGRAM_COMPLEXITY, tokenizer=DEFAULT_TOKENIZER])

Creates a document-term-vector with elements of type `T` for document `d`
using the lexicon `lex`. `d` can be an `AbstractString` or an `AbstractDocument`.
"""
function dtv(d,
             lex::OrderedDict{String, Int},
             eltype::Type{T}=DEFAULT_DTM_TYPE;
             ngram_complexity::Int=DEFAULT_NGRAM_COMPLEXITY,
             tokenizer::Symbol=DEFAULT_TOKENIZER,
             lex_is_row_indices::Bool=false) where {T<:Real}
    p = length(keys(lex))
    column = spzeros(T, p)
    indices, values = dtm_entries(d, lex, eltype,
                                  ngram_complexity=ngram_complexity,
                                  tokenizer=tokenizer,
                                  lex_is_row_indices=lex_is_row_indices)
    column[indices] = values
    return column
end

"""
    dtv(crps::Corpus, idx::Int, eltype::Type{T}=DEFAULT_DTM_TYPE [; ngram_complexity=DEFAULT_NGRAM_COMPLEXITY, tokenizer=DEFAULT_TOKENIZER])

Creates a document-term-vector with elements of type `T` for document `idx`
of the corpus `crps`.
"""
function dtv(crps::Corpus,
             idx::Int,
             eltype::Type{T}=DEFAULT_DTM_TYPE;
             ngram_complexity::Int=DEFAULT_NGRAM_COMPLEXITY,
             tokenizer::Symbol=DEFAULT_TOKENIZER,
             lex_is_row_indices::Bool=false) where {T<:Real}
    if idx >= length(crps.documents) || idx < 1
        error("DTV requires the document index in [1,$(length(crps.documents))]")
    end
    if isempty(crps.lexicon)
        lex = create_lexicon(crps, ngram_complexity)
    else
        lex = lexicon(crps)
    end
    return dtv(crps.documents[idx], lex, eltype,
               ngram_complexity=ngram_complexity,
               tokenizer=tokenizer,
               lex_is_row_indices=lex_is_row_indices)
end

function dtv(d;
             ngram_complexity::Int=DEFAULT_NGRAM_COMPLEXITY,
             tokenizer::Symbol=DEFAULT_TOKENIZER,
             lex_is_row_indices::Bool=false)
    throw(ErrorException("Cannot construct a DTV without a pre-existing lexicon"))
end

# Document is a list of regular expressions in text form
function dtm_regex_entries(d,
                           lex::OrderedDict{String, Int},
                           eltype::Type{T}=DEFAULT_DTM_TYPE;
                           ngram_complexity::Int=DEFAULT_NGRAM_COMPLEXITY,
                           tokenizer::Symbol=DEFAULT_TOKENIZER,
                           lex_is_row_indices::Bool=false) where {T<:Real}
    ngs = ngrams(d, ngram_complexity, tokenizer=tokenizer)
    patterns = Regex.(keys(ngs))
    indices = Vector{Int}(undef, 0)
    terms = collect(keys(lex))
    local row_indices
    if lex_is_row_indices
        row_indices = lex
    else
        row_indices = rowindices(terms)
    end
    for pattern in patterns
        for term in terms
            if occursin(pattern, term)
                j = row_indices[term]
                push!(indices, j)
            end
        end
    end
    values = ones(T, length(indices))
    return (indices, values)
end

"""
    dtv_regex(d, lex::OrderedDict{String,Int}, eltype::Type{T}=DEFAULT_DTM_TYPE [; ngram_complexity=DEFAULT_NGRAM_COMPLEXITY, tokenizer=DEFAULT_TOKENIZER])

Creates a document-term-vector with elements of type `T` for document `d`
using the lexicon `lex`. The tokens of document `d` are assumed to be regular
expressions in text format. `d` can be an `AbstractString` or an `AbstractDocument`.

# Examples
```
julia> dtv_regex(NGramDocument("a..b"), OrderedDict("aaa"=>1, "aaab"=>2, "accb"=>3, "bbb"=>4), Float32)
4-element Array{Float32,1}:
 0.0
 1.0
 1.0
 0.0
```
"""
function dtv_regex(d,
                   lex::OrderedDict{String, Int},
                   eltype::Type{T}=DEFAULT_DTM_TYPE;
                   ngram_complexity::Int=DEFAULT_NGRAM_COMPLEXITY,
                   tokenizer::Symbol=DEFAULT_TOKENIZER,
                   lex_is_row_indices::Bool=false) where {T<:Real}
    p = length(keys(lex))
    column = spzeros(T, p)
    indices, values = dtm_regex_entries(d, lex, eltype,
                                        ngram_complexity=ngram_complexity,
                                        tokenizer=tokenizer,
                                        lex_is_row_indices=lex_is_row_indices)
    column[indices] = values
    return column
end

"""
    hash_dtv(d, h::TextHashFunction, eltype::Type{T}=DEFAULT_DTM_TYPE [; ngram_complexity=DEFAULT_NGRAM_COMPLEXITY, tokenizer=DEFAULT_TOKENIZER])

Creates a hashed document-term-vector with elements of type `T` for document `d`
using the hashing function `h`. `d` can be an `AbstractString` or an `AbstractDocument`.
"""
function hash_dtv(d,
                  h::TextHashFunction,
                  eltype::Type{T}=DEFAULT_DTM_TYPE;
                  ngram_complexity::Int=DEFAULT_NGRAM_COMPLEXITY,
                  tokenizer::Symbol=DEFAULT_TOKENIZER) where {T<:Real}
    p = cardinality(h)
    res = spzeros(T, p)
    ngs = ngrams(d, ngram_complexity, tokenizer=tokenizer)
    for ng in keys(ngs)
        res[index_hash(ng, h)] += ngs[ng]
    end
    return res
end

function hash_dtv(d;
                  cardinality::Int=DEFAULT_CARDINALITY,
                  eltype::Type{T}=DEFAULT_DTM_TYPE,
                  ngram_complexity::Int=DEFAULT_NGRAM_COMPLEXITY,
                  tokenizer::Symbol=DEFAULT_TOKENIZER) where {T<:Real}
    hash_dtv(d, TextHashFunction(cardinality), eltype, ngram_complexity=ngram_complexity, tokenizer=tokenizer)
end

"""
    hash_dtm(crps::Corpus [,h::TextHashFunction], eltype::Type{T}=DEFAULT_DTM_TYPE [; ngram_complexity=DEFAULT_NGRAM_COMPLEXITY, tokenizer=DEFAULT_TOKENIZER])

Creates a hashed DTM with elements of type `T` for corpus `crps` using the
the hashing function `h`. If `h` is missing, the hash function of the `Corpus`
is used.
"""
function hash_dtm(crps::Corpus,
                  h::TextHashFunction,
                  eltype::Type{T}=DEFAULT_DTM_TYPE;
                  ngram_complexity::Int=DEFAULT_NGRAM_COMPLEXITY,
                  tokenizer::Symbol=DEFAULT_TOKENIZER) where {T<:Real}
    n, p = length(crps), cardinality(h)
    res = spzeros(T, p, n)
    for (i, doc) in enumerate(crps)
        res[:, i] = hash_dtv(doc, h, eltype, ngram_complexity=ngram_complexity, tokenizer=tokenizer)
    end
    return res
end

hash_dtm(crps::Corpus,
         eltype::Type{T}=DEFAULT_DTM_TYPE;
         ngram_complexity::Int=DEFAULT_NGRAM_COMPLEXITY,
         tokenizer::Symbol=DEFAULT_TOKENIZER) where {T<:Real} =
    hash_dtm(crps, hash_function(crps), eltype, ngram_complexity=ngram_complexity, tokenizer=tokenizer)

# Produce entries for on-line analysis when DTM would not fit in memory
struct EachDTV{U, S<:AbstractString, T<:AbstractDocument}
    corpus::Corpus{S,T}
    row_indices::OrderedDict{String, Int}
    ngram_complexity::Int
    tokenizer::Symbol
    function EachDTV{U,S,T}(corpus::Corpus{S,T},
                            row_indices::OrderedDict{String, Int},
                            ngram_complexity::Int=DEFAULT_NGRAM_COMPLEXITY,
                            tokenizer::Symbol=DEFAULT_TOKENIZER
                           ) where {U, S<:AbstractString, T<:AbstractDocument}
        @assert ngram_complexity >= 1 "Ngram complexity has to be >= 1"
        @assert tokenizer in [:default,
            :stringanalysis] "Tokenizer has to be either :default or :stringanalysis"
        new(corpus, row_indices, ngram_complexity, tokenizer)
    end
end

EachDTV{U}(crps::Corpus{S,T};
           ngram_complexity::Int=DEFAULT_NGRAM_COMPLEXITY,
           tokenizer::Symbol=DEFAULT_TOKENIZER) where {U,S,T} = begin
    row_indices = rowindices(collect(keys(create_lexicon(crps, ngram_complexity))))
    EachDTV{U,S,T}(crps, row_indices, ngram_complexity, tokenizer)
end

Base.iterate(edt::EachDTV, state=1) = begin
    if state > length(edt.corpus)
        return nothing
    else
        return next(edt, state)
    end
end

next(edt::EachDTV{U,S,T}, state::Int) where {U,S,T} =
    (dtv(edt.corpus.documents[state], edt.row_indices, U,
         ngram_complexity=edt.ngram_complexity,
         tokenizer=edt.tokenizer, lex_is_row_indices=true), state + 1)

"""
    each_dtv(crps::Corpus [; eltype::Type{U}=DEFAULT_DTM_TYPE, ngram_complexity=DEFAULT_NGRAM_COMPLEXITY, tokenizer=DEFAULT_TOKENIZER])

Iterates through the columns of the DTM of the corpus `crps` without
constructing it. Useful when the DTM would not fit in memory.
`eltype` specifies the element type of the generated vectors.
"""
each_dtv(crps::Corpus;
         eltype::Type{U}=DEFAULT_DTM_TYPE,
         ngram_complexity::Int=DEFAULT_NGRAM_COMPLEXITY,
         tokenizer::Symbol=DEFAULT_TOKENIZER) where U<:Real =
    EachDTV{U}(crps, ngram_complexity=ngram_complexity, tokenizer=tokenizer)

Base.eltype(::Type{EachDTV{U,S,T}}) where {U,S,T} = Vector{U}

Base.length(edt::EachDTV) = length(edt.corpus)

Base.size(edt::EachDTV) = (length(edt.corpus), edt.corpus.h.cardinality)

Base.show(io::IO, edt::EachDTV{U,S,T}) where {U,S,T} =
    print(io, "DTV iterator, tokenizer is $(edt.tokenizer), "*
          "$(length(edt)) elements of type $(eltype(edt)).")

struct EachHashDTV{U, S<:AbstractString, T<:AbstractDocument}
    corpus::Corpus{S,T}
    ngram_complexity::Int
    tokenizer::Symbol
    function EachHashDTV{U,S,T}(corpus::Corpus{S,T},
                                ngram_complexity::Int=DEFAULT_NGRAM_COMPLEXITY,
                                tokenizer::Symbol=DEFAULT_TOKENIZER) where
            {U, S<:AbstractString, T<:AbstractDocument}
        @assert ngram_complexity >= 1 "Ngram complexity has to be >= 1"
        @assert tokenizer in [:default,
            :stringanalysis] "Tokenizer has to be either :default or :stringanalysis"
        new(corpus, ngram_complexity, tokenizer)
    end
end

EachHashDTV{U}(crps::Corpus{S,T};
               ngram_complexity::Int=DEFAULT_NGRAM_COMPLEXITY,
               tokenizer::Symbol=DEFAULT_TOKENIZER) where {U,S,T} =
    EachHashDTV{U,S,T}(crps, ngram_complexity, tokenizer)

Base.iterate(edt::EachHashDTV, state=1) = begin
    if state > length(edt.corpus)
        return nothing
    else
        return next(edt, state)
    end
end

next(edt::EachHashDTV{U,S,T}, state::Int) where {U,S,T} =
    (hash_dtv(edt.corpus.documents[state], edt.corpus.h, U,
              ngram_complexity=edt.ngram_complexity,
              tokenizer=edt.tokenizer), state + 1)

"""
    each_hash_dtv(crps::Corpus [; eltype::Type{U}=DEFAULT_DTM_TYPE, ngram_complexity=DEFAULT_NGRAM_COMPLEXITY, tokenizer=DEFAULT_TOKENIZER])

Iterates through the columns of the hashed DTM of the corpus `crps` without
constructing it. Useful when the DTM would not fit in memory.
`eltype` specifies the element type of the generated vectors.
"""
each_hash_dtv(crps::Corpus;
              eltype::Type{U}=DEFAULT_DTM_TYPE,
              ngram_complexity::Int=DEFAULT_NGRAM_COMPLEXITY,
              tokenizer::Symbol=DEFAULT_TOKENIZER) where U<:Real =
    EachHashDTV{U}(crps, ngram_complexity=ngram_complexity, tokenizer=tokenizer)

Base.eltype(::Type{EachHashDTV{U,S,T}}) where {U,S,T} = Vector{U}

Base.length(edt::EachHashDTV) = length(edt.corpus)

Base.size(edt::EachHashDTV) = (length(edt.corpus), edt.corpus.h.cardinality)

Base.show(io::IO, edt::EachHashDTV{U,S,T}) where {U,S,T} =
    print(io, "Hash-DTV iterator, tokenizer is $(edt.tokenizer) "*
          "$(length(edt)) elements of type $(eltype(edt)).")

## getindex() methods
Base.getindex(dtm::DocumentTermMatrix, k::AbstractString) = dtm.dtm[dtm.row_indices[k], :]

Base.getindex(dtm::DocumentTermMatrix, i) = dtm.dtm[i]

Base.getindex(dtm::DocumentTermMatrix, i, j) = dtm.dtm[i, j]


text:
using Dates
using Mmap

using TeaFiles

function _read_all_items_mmap(
    ::Type{T},
    io::IOStream,
    metadata::TeaFiles.Header.TeaFileMetadata
)::Vector{T} where T
    # Memory map the file from the start of the item section
    seek(io, metadata.item_start)
    map = Mmap.mmap(io)

    num_items = TeaFiles.Body.get_num_items(sizeof(map), metadata)

    # Output buffer for the items that we're going to read.
    item_vector = Vector{T}(undef, num_items)

    # NOTE - one can gain some performance by using unsafe_copyto!
    item_vector .= reinterpret(T, map)
    return item_vector
end

function _read_all_items_syscall(
    ::Type{T},
    io::IOStream,
    metadata::TeaFiles.Header.TeaFileMetadata
)::Vector{T} where T
    seek(io, metadata.item_start)
    num_items = TeaFiles.Body.get_num_items(io, metadata)
    # In practice one may want to define "read" for type T
    item_vector = Vector{T}(undef, num_items)
    unsafe_read(io, Ref(item_vector, 1), sizeof(T) * num_items)
    return item_vector
end

function _test_write_read(data::Vector{T}; metadata_kwargs...) where T
    # Temporary file to use for experiments. We have to use a file to test memory mapping,
    # since this can't be used with a buffer.
    path = tempname()

    metadata = create_metadata(T; metadata_kwargs...)

    open(path; write=true) do io
        # Create the file and write the header and data to it.
        write(io, metadata)
        write(io, data)
    end

    result_mmap, result_syscall = open(path; read=true) do io
        # Read the file.
        read_metadata = read(io, TeaFiles.Header.TeaFileMetadata)
        @test read_metadata == metadata

        result_mmap = _read_all_items_mmap(Item, io, metadata)
        result_syscall = _read_all_items_syscall(Item, io, metadata)
        return result_mmap, result_syscall
    end

    @test result_mmap == data
    @test result_syscall == data

    # Necessary to remove any remnants of the memory mapping on windows, at least. Otherwise
    # we can find that the file remains open.
    GC.gc()
    GC.gc()
end

@testset "simple" begin
    struct Item
        time::Int64
        a::Int32
        b::Float64
    end

    some_data = [
        Item(1000 + i, i, 2 * i)
        for i in 1:100
    ]
    _test_write_read(some_data)
    _test_write_read(some_data; time_fields=[:time])
end


text:
using Graphs
using Test

function has_match(regex, str)
    match(regex, str) !== nothing
end

function xor(a, b)
    (a || b) && !(a && b)
end

# Vertex attributes get layed out correctly
let g=inclist(ExVertex, is_directed=false)
    add_vertex!(g, ExVertex(1, "label1"))
    add_vertex!(g, ExVertex(2, "label2"))
    add_edge!(g, vertices(g)[1], vertices(g)[2])
    attrs = attributes(vertices(g)[1], g)

    @test to_dot(attrs) == ""

    attrs["foo"] = "bar"
    @test to_dot(attrs) == "[\"foo\"=\"bar\"]"

    attrs["baz"] = "qux"
    # Below here we do some messing around so as not to assert on dict iteration order.

    @test to_dot(attrs) in ["[\"foo\"=\"bar\",\"baz\"=\"qux\"]",
                            "[\"baz\"=\"qux\",\"foo\"=\"bar\"]"]

    global sp = split(to_dot(g), "\n")
    @test ("1 [\"foo\"=\"bar\",\"baz\"=\"qux\"]" in sp) ||
          ("1 [\"baz\"=\"qux\",\"foo\"=\"bar\"]" in sp)
end

# Edge attributes get layed out correctly
let
    global g = inclist(ExVertex, ExEdge{ExVertex}, is_directed=false)
    add_vertex!(g, ExVertex(1, "label1"))
    add_vertex!(g, ExVertex(2, "label2"))

    add_edge!(g, vertices(g)[1], vertices(g)[2])

    global e = out_edges(vertices(g)[2], g)[1]

    attrs = attributes(e, g)
    attrs["foo"] = "bar"
    @test to_dot(attrs) == "[\"foo\"=\"bar\"]"

    attrs["baz"] = "qux"

    global sp = split(to_dot(g), "\n")
    @test ("1 -- 2 [\"foo\"=\"bar\",\"baz\"=\"qux\"]" in sp) ||
          ("1 -- 2 [\"baz\"=\"qux\",\"foo\"=\"bar\"]" in sp)
end

# Graph attributes get layed out correctly
let
    global g = inclist(ExVertex, ExEdge{ExVertex}, is_directed=false)
    add_vertex!(g, ExVertex(1, "label1"))
    add_vertex!(g, ExVertex(2, "label2"))
    add_edge!(g, vertices(g)[1], vertices(g)[2])
    global e = out_edges(vertices(g)[2], g)[1]
    graph_attrs = AttributeDict()

    graph_attrs["foo"]="bar"
    global sp = split(to_dot(g, graph_attrs), "\n")
    @test ("\"foo\"=\"bar\";" in sp)
    @test !("\"baz\"=\"qux\";" in sp)

    graph_attrs["baz"]="qux"
    global sp = split(to_dot(g, graph_attrs), "\n")
    @test ("\"foo\"=\"bar\";" in sp)
    @test ("\"baz\"=\"qux\";" in sp)
end

let g=simple_graph(0, is_directed=false)
    @test to_dot(g) == "graph graphname {\n}\n"
end

let g=simple_graph(0, is_directed=true)
    @test to_dot(g) == "digraph graphname {\n}\n"
end

let g=simple_adjlist(3)
    add_edge!(g, 1, 2)
    add_edge!(g, 1, 3)
    add_edge!(g, 2, 3)
    add_edge!(g, 3, 2)
    str = to_dot(g)
    @test has_match(r"1 -> 2", str)
    @test has_match(r"1 -> 3", str)
    @test has_match(r"2 -> 3", str)
    @test has_match(r"3 -> 2", str)
    @test !has_match(r"--", str)
end

let g=simple_adjlist(3, is_directed=false)
    add_edge!(g, 1, 2)
    add_edge!(g, 1, 3)
    add_edge!(g, 2, 3)
    str = to_dot(g)
    @test xor(has_match(r"1 -- 2", str), has_match(r"2 -- 1", str))
    @test xor(has_match(r"1 -- 3", str), has_match(r"3 -- 1", str))
    @test xor(has_match(r"2 -- 3", str), has_match(r"3 -- 2", str))
    @test !has_match(r"->", str)
end


text:
__precompile__(true)
module CSV

using DataStreams, WeakRefStrings, InternedStrings, Missings, CategoricalArrays, DataFrames

using Mmap, Dates

struct ParsingException <: Exception
    msg::String
end

const RETURN  = UInt8('\r')
const NEWLINE = UInt8('\n')
const COMMA   = UInt8(',')
const QUOTE   = UInt8('"')
const ESCAPE  = UInt8('\\')
const PERIOD  = UInt8('.')
const SPACE   = UInt8(' ')
const TAB     = UInt8('\t')
const MINUS   = UInt8('-')
const PLUS    = UInt8('+')
const NEG_ONE = UInt8('0')-UInt8(1)
const ZERO    = UInt8('0')
const TEN     = UInt8('9')+UInt8(1)
Base.isascii(c::UInt8) = c < 0x80

readbyte(from::IO) = Base.read(from, UInt8)
peekbyte(from::IO) = Base.peek(from)

@inline function readbyte(from::IOBuffer)
    @inbounds byte = from.data[from.ptr]
    from.ptr = from.ptr + 1
    return byte
end

@inline function peekbyte(from::IOBuffer)
    @inbounds byte = from.data[from.ptr]
    return byte
end

substitute(::Type{Union{T, Missing}}, ::Type{T1}) where {T, T1} = Union{T1, Missing}
substitute(::Type{T}, ::Type{T1}) where {T, T1} = T1
substitute(::Type{Missing}, ::Type{T1}) where {T1} = Missing

"""
Represents the various configuration settings for delimited text file parsing.

Keyword Arguments:

 * `delim::Union{Char,UInt8}`: how fields in the file are delimited; default `','`
 * `quotechar::Union{Char,UInt8}`: the character that indicates a quoted field that may contain the `delim` or newlines; default `'"'`
 * `escapechar::Union{Char,UInt8}`: the character that escapes a `quotechar` in a quoted field; default `'\\'`
 * `missingstring::String`: indicates how missing values are represented in the dataset; default `""`
 * `dateformat::Union{AbstractString,Dates.DateFormat}`: how dates/datetimes are represented in the dataset; default `Base.Dates.ISODateTimeFormat`
 * `decimal::Union{Char,UInt8}`: character to recognize as the decimal point in a float number, e.g. `3.14` or `3,14`; default `'.'`
 * `truestring`: string to represent `true::Bool` values in a csv file; default `"true"`. Note that `truestring` and `falsestring` cannot start with the same character.
 * `falsestring`: string to represent `false::Bool` values in a csv file; default `"false"`
 * `internstrings`: whether strings should be interned, rather than creating a new object for each string field; default `true`
"""
struct Options{D}
    delim::UInt8
    quotechar::UInt8
    escapechar::UInt8
    missingstring::Vector{UInt8}
    null::Union{Vector{UInt8},Nothing} # deprecated
    missingcheck::Bool
    dateformat::D
    decimal::UInt8
    truestring::Vector{UInt8}
    falsestring::Vector{UInt8}
    internstrings::Bool
    # non-public for now
    datarow::Int
    rows::Int
    header::Union{Integer,UnitRange{Int},Vector}
    types
end

Options(;delim=COMMA, quotechar=QUOTE, escapechar=ESCAPE, missingstring="", null=nothing, dateformat=nothing, decimal=PERIOD, truestring="true", falsestring="false", internstrings=true, datarow=-1, rows=0, header=1, types=Type[]) =
    Options(delim%UInt8, quotechar%UInt8, escapechar%UInt8,
            map(UInt8, collect(ascii(String(missingstring)))), null === nothing ? nothing : map(UInt8, collect(ascii(String(null)))), missingstring != "" || (null != "" && null != nothing),
            isa(dateformat, AbstractString) ? Dates.DateFormat(dateformat) : dateformat,
            decimal%UInt8, map(UInt8, collect(truestring)), map(UInt8, collect(falsestring)), internstrings, datarow, rows, header, types)
function Base.show(io::IO,op::Options)
    println(io, "    CSV.Options:")
    println(io, "        delim: '", Char(op.delim), "'")
    println(io, "        quotechar: '", Char(op.quotechar), "'")
    print(io, "        escapechar: '"); escape_string(io, string(Char(op.escapechar)), "\\"); println(io, "'")
    print(io, "        missingstring: \""); escape_string(io, isempty(op.missingstring) ? "" : String(collect(op.missingstring)), "\\"); println(io, "\"")
    println(io, "        dateformat: ", op.dateformat == nothing ? "" : op.dateformat)
    println(io, "        decimal: '", Char(op.decimal), "'")
    println(io, "        truestring: '$(String(op.truestring))'")
    println(io, "        falsestring: '$(String(op.falsestring))'")
    print(io, "        internstrings: $(op.internstrings)")
end

"""
A type that satisfies the `Data.Source` interface in the `DataStreams.jl` package.

A `CSV.Source` can be manually constructed in order to be re-used multiple times.

`CSV.Source(file_or_io; kwargs...) => CSV.Source`

Note that a filename string can be provided or any `IO` type. For the full list of supported
keyword arguments, see the docs for [`CSV.read`](@ref) or type `?CSV.read` at the repl

An example of re-using a `CSV.Source` is:
```julia
# manually construct a `CSV.Source` once, then stream its data to both a DataFrame
# and SQLite table `sqlite_table` in the SQLite database `db`
# note the use of `CSV.reset!` to ensure the `source` can be streamed from again
source = CSV.Source(file)
df1 = CSV.read(source, DataFrame)
CSV.reset!(source)
sq1 = CSV.read(source, SQLite.Sink, db, "sqlite_table")
```
"""
mutable struct Source{I, D} <: Data.Source
    schema::Data.Schema
    options::Options{D}
    io::I
    fullpath::String
    datapos::Int # the position in the IOBuffer where the rows of data begins
end

function Base.show(io::IO, f::Source)
    println(io, "CSV.Source: ", f.fullpath)
    println(io, f.options)
    show(io, f.schema)
end

"""
`CSV.TransposedSource(file_or_io; kwargs...) => CSV.TransposedSource`

Type that in all respects is identical to a `CSV.Source`, except when reading a csv file,
the data will be parsed "transposed", i.e. rows will become columns / columns will become rows.
This can be a huge convenience if the csv data happens to be transposed, as well as lead to huge
performance gains as the resulting data set can be more consistently typed.

Typical usage involves calling `CSV.read(file; transpose=true)`.
"""
mutable struct TransposedSource{I, D} <: Data.Source
    schema::Data.Schema
    options::Options{D}
    io::I
    fullpath::String
    datapos::Int # the position in the IOBuffer where the rows of data begins
    columnpositions::Vector{Int}
end

"""
A type that satisfies the `Data.Sink` interface in the `DataStreams.jl` package.

A `CSV.Sink` can be manually constructed in order to be re-used multiple times.

`CSV.Sink(file_or_io; kwargs...) => CSV.Sink`

Note that a filename string can be provided or any `IO` type. For the full list of supported
keyword arguments, see the docs for [`CSV.write`](@ref) or type `?CSV.write` at the repl

An example of re-using a `CSV.Sink` is:
```julia
# manually construct a `CSV.Source` once, then stream its data to both a DataFrame
# and SQLite table `sqlite_table` in the SQLite database `db`
# note the use of `CSV.reset!` to ensure the `source` can be streamed from again
source = CSV.Source(file)
df1 = CSV.read(source, DataFrame)
CSV.reset!(source)
sq1 = CSV.read(source, SQLite.Sink, db, "sqlite_table")
```
"""
mutable struct Sink{D, B} <: Data.Sink
    options::Options{D}
    io::IOBuffer
    fullpath::Union{String, IO}
    datapos::Int # the position in the IOBuffer where the rows of data begins
    header::Bool
    colnames::Vector{String}
    cols::Int
    append::Bool
    quotefields::B
end

include("parsefields.jl")
include("float.jl")
include("io.jl")
include("TransposedSource.jl")
include("Source.jl")
include("Sink.jl")
include("validate.jl")

end # module


text:
# using MultiScaleTreeGraph
# using Plots
# using ColorSchemes
# plotlyjs()

# file = joinpath(dirname(dirname(pathof(MultiScaleTreeGraph))), "test", "files", "simple_plant.mtg")
# mtg = read_mtg(file)
# plot(mtg)
@recipe function f(h::Node; mode = "2d")
    branching_order!(mtg)
    df_coordinates = MultiScaleTreeGraph.mtg_coordinates_df(mtg, force = true)
    df_coordinates[!,:branching_order] = descendants(mtg, :branching_order, self = true)

    x = df_coordinates.XX
    y = df_coordinates.YY
    z = df_coordinates.ZZ

    for i in 2:size(df_coordinates, 1)
        x2 = [df_coordinates.XX_from[i], df_coordinates.XX[i]]
        y2 = [df_coordinates.YY_from[i], df_coordinates.YY[i]]
        z2 = [df_coordinates.ZZ_from[i], df_coordinates.ZZ[i]]

        @series begin
            label := ""
            seriescolor := :black
            if mode == "2d"
                seriestype := :line
                x2, y2
            else
                seriestype := :line3d
                x2, z2, y2
            end
        end
    end

    @series begin
        label := ""
        seriescolor := :viridis
        marker_z := df_coordinates.branching_order
        colorbar_entry := false
        hover := string.(
            "name: `node_", df_coordinates.id,
            "`, link: `", df_coordinates.link,
            "`, symbol: `", df_coordinates.symbol,
            "`, index: `", df_coordinates.index, "`"
        )

        if mode == "2d"
            seriestype := :scatter
            x, y
        else
            seriestype := :scatter3d
            x, z, y
        end

    end
end

# p = plot()
    # for i in 2:size(df_coordinates, 1)
#     x2 = [df_coordinates.XX_from[i], df_coordinates.XX[i]]
#     y2 = [df_coordinates.YY_from[i], df_coordinates.YY[i]]
#     plot!(p, x2, y2, label = "", color = :black)
# end
# scatter!(p,x,y, marker_z = df_coordinates.branching_order)