text:
"""
ObjectPool is a fixed-size one-dimensional array that does not store
any duplicate copies of the same object.  So the benefit is space-efficiency.
The tradeoff is the time used to maintain the index.
This is useful for denormalized data frames where string values
may be repeated many times.

An ObjectPool must be initialize with a default value and a fixed
array size.  If your requirement does not fit such assumptions,
you may want to look into using `PooledArrays` or
`CategoricalArrays` package instead.

The implementation is very primitive and is tailor for application
that knows exactly how much memory to allocate.
"""
mutable struct ObjectPool{T, S <: Unsigned} <: AbstractArray{T, 1}
    pool::Array{T}             # maintains the pool of unique things
    idx::Array{S}              # index references into `pool`
    indexcache::Dict{T, S}     # dict for fast lookups (K=object, V=index)
    uniqueitemscount::Int64    # how many items in `pool`, always start with 1
    itemscount::Int64          # how many items perceived in this array
    capacity::Int64            # max number of items in the pool
end

# Initially, there is only one item in the pool and the `idx` array has
# elements all pointing to that one default vaue.  The dictionary `indexcache`
# also has one item that points to that one value.  Hence `uniqueitemcount`
# would be 1 and `itemscount` would be `n`.
function ObjectPool{T, S}(val::T, n::Integer) where {T, S <: Unsigned}
    # Note: 64-bit case is constrainted by Int64 type (for convenience)
    maxsize = ifelse(S == UInt8, 2 << 7 - 1,
                ifelse(S == UInt16, 2 << 15 - 1,
                    ifelse(S == UInt32, 2 << 31 - 1,
                        2 << 62 - 1)))
    ObjectPool{T, S}([val], fill(1, n), Dict(val => 1), 1, n, maxsize)
end

# If the value already exist in the pool then just the index value is stored.
function Base.setindex!(op::ObjectPool, val::T, i::Integer) where T
    if haskey(op.indexcache, val)
        # The value `val` already exists in the cache.
        # Just set the array element to the index value from cache.
        op.idx[i] = op.indexcache[val]
    else
        if op.uniqueitemscount >= op.capacity
            throw(BoundsError("Exceeded pool capacity $(op.capacity). Consider using a larger pool size e.g. UInt32."))
        end
        # Encountered a new value `val`:
        # 1. add ot the object pool array
        # 2. increment the number of unique items
        # 3. store the new index in the cache
        # 4. set the array element with the new index value
        push!(op.pool, val)
        op.uniqueitemscount += 1
        op.indexcache[val] = op.uniqueitemscount
        op.idx[i] = op.uniqueitemscount
    end
    op
end

# AbstractArray trait
# Base.IndexStyle(::Type{<:ObjectPool}) = IndexLinear()

# single indexing
Base.getindex(op::ObjectPool, i::Integer) = op.pool[op.idx[convert(Int, i)]]

# general sizes
Base.size(op::ObjectPool)   = (op.itemscount, )
# Base.length(op::ObjectPool) = op.itemscount
# Base.endof(op::ObjectPool)  = op.itemscount

# typing
#Base.eltype(op::ObjectPool) = eltype(op.pool)

# make it iterable
# Base.start(op::ObjectPool)  = 1
# Base.next(op::ObjectPool, state) = (op.pool[op.idx[state]], state + 1)
# Base.done(op::ObjectPool, state) = state > op.itemscount

# custom printing
# function Base.show(io::IO, op::ObjectPool)
#     L = op.itemscount
#     print(io, "$L-element ObjectPool with $(op.uniqueitemscount) unique items:\n")
#     if L > 20
#         for i in 1:10  print(io, " ", op[i], "\n") end
#         print(io, " ⋮\n")
#         for i in L-9:L print(io, " ", op[i], "\n") end
#     else
#         for i in 1:L   print(io, " ", op[i], "\n") end
#     end
# end


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

module SimilaritySearch
abstract type AbstractSearchContext end

using Parameters

import Base: push!, append!
export AbstractSearchContext, SemiMetric, evaluate, search, searchbatch, getknnresult
include("distances/Distances.jl")

include("db.jl")
include("knnresult.jl")
@inline Base.length(searchctx::AbstractSearchContext) = length(searchctx.db)
@inline Base.getindex(searchctx::AbstractSearchContext, i::Integer) = searchctx.db[i]
@inline Base.eachindex(searchctx::AbstractSearchContext) = 1:length(searchctx)
@inline Base.eltype(searchctx::AbstractSearchContext) = eltype(searchctx.db)
include("perf.jl")
include("sequential-exhaustive.jl")
include("parallel-exhaustive.jl")
include("opt.jl")
include("searchgraph/SearchGraph.jl")
include("allknn.jl")
include("neardup.jl")
include("closestpair.jl")

const GlobalKnnResult = [KnnResult(32)]   # see __init__ function at the end of this file

"""
    getknnresult(k::Integer, pools=nothing) -> KnnResult

Generic function to obtain a shared result set for the same thread and avoid memory allocations.
This function should be specialized for indexes and pools that use shared results or threads in some special way.
"""
@inline function getknnresult(k::Integer, pools=nothing)
    res = @inbounds GlobalKnnResult[Threads.threadid()]
    reuse!(res, k)
end

"""
    searchbatch(index, Q, k::Integer; parallel=false, pools=GlobalKnnResult) -> indices, distances

Searches a batch of queries in the given index (searches for k neighbors).

- `parallel` specifies if the query should be solved in parallel at object level (each query is sequentially solved but many queries solved in different threads).
- `pool` relevant if `parallel=true`. If it is explicitly given it should be an array of `Threads.nthreads()` preallocated `KnnResult` objects used to reduce memory allocations.
    In most case uses the default is enought, but different pools should be used when indexes use internal indexes to solve queries (e.g., using index's proxies or database objects defined as indexes).

Note: The i-th column in indices and distances correspond to the i-th query in `Q`
Note: The final indices at each column can be `0` if the search process was unable to retrieve `k` neighbors.
"""
function searchbatch(index, Q, k::Integer; parallel=false, pools=getpools(index))
    m = length(Q)
    I = zeros(Int32, k, m)
    D = Matrix{Float32}(undef, k, m)
    searchbatch(index, Q, I, D; parallel, pools)
end

"""
    searchbatch(index, Q, I::AbstractMatrix{Int32}, D::AbstractMatrix{Float32}; parallel=false, pools=getpools(index)) -> indices, distances

Searches a batch of queries in the given index and `I` and `D` as output (searches for `k=size(I, 1)`)

"""
function searchbatch(index, Q, I::AbstractMatrix{Int32}, D::AbstractMatrix{Float32}; parallel=false, pools=getpools(index))
    k = size(I, 1)

    if parallel
        Threads.@threads for i in eachindex(Q)
            res, _ = search(index, Q[i], getknnresult(k, pools); pools)
            k_ = length(res)
            @inbounds I[1:k_, i] .= res.id
            @inbounds D[1:k_, i] .= res.dist
        end
    else
        @inbounds for i in eachindex(Q)
            res, _ = search(index, Q[i], getknnresult(k, pools); pools)
            k_ = length(res)
            I[1:k_, i] .= res.id
            D[1:k_, i] .= res.dist
        end
    end

    I, D
end

"""
    searchbatch(index, Q, KNN::AbstractVector{KnnResult}; parallel=false, pools=getpools(index)) -> indices, distances

Searches a batch of queries in the given index using an array of KnnResult's; each KnnResult object can specify different `k` values.

"""
function searchbatch(index, Q, KNN::AbstractVector{KnnResult}; parallel=false, pools=getpools(index))
    if parallel
        Threads.@threads for i in eachindex(Q)
            @inbounds search(index, Q[i], KNN[i]; pools)
        end
    else
        @inbounds for i in eachindex(Q)
            search(index, Q[i], KNN[i]; pools)
        end
    end

    KNN
end

function __init__()
    __init__visitedvertices()
    __init__beamsearch()
    __init__neighborhood()

    for _ in 2:Threads.nthreads()
        push!(GlobalKnnResult, KnnResult(32))
    end
end

# precompile as the final step of the module definition:

if ccall(:jl_generating_output, Cint, ()) == 1   # if we're precompiling the package
    @info "precompiling common combinations of indexes, distances, and databases"
    let
        # Note: something happens that the warming stage is not removed
        #=
        function run_functions(E, queries)
            I, D = searchbatch(E, queries, 3)
            @assert size(I) == (3, length(queries))
            I, D = allknn(E, 3)
            @assert size(I) == (3, length(queries))
            p1, p2, d = closestpair(E)
            @assert p1 != p2
        end

        # running common combinations for precompilation
        for X in [rand(Float32, 2, 16), rand(Float64, 2, 16)] # 32 to force calling triggers / callbacks
            for c in eachcol(X)
                normalize!(c)
            end

            # for dist in [L2Distance(), SqL2Distance(), L1Distance(), CosineDistance(), NormalizedCosineDistance(), AngleDistance(), NormalizedAngleDistance()]
            for dist in [L2Distance(), SqL2Distance(), CosineDistance(), NormalizedCosineDistance()]
                #for db in [MatrixDatabase(X), VectorDatabase(X)]
                for db in [MatrixDatabase(X)]
                    #for db in [db_, rand(db_, 15)]
                        E = ExhaustiveSearch(; db, dist)
                        run_functions(E, db)
                        G = SearchGraph(; db, dist, verbose=false)
                        index!(G)
                        run_functions(G, db)
                        optimize!(G, ParetoRecall())
                        optimize!(G, MinRecall(0.8))
                        optimize!(G, ParetoRadius())
                        try
                            # static databases will throw an error
                            push!(G, db[1])
                            append!(G, db)
                        finally
                            continue
                        end
                    #end
                end
            end
        end
        =#
        #=
        # Skip precompilation of the following combinations since they are not so used but increase significantly the compilation times
        for X in [
                rand(Int64(1):Int64(6), 4, 32), rand(Int32(1):Int32(6), 4, 32), rand(Int16(1):Int16(6), 4, 32), rand(Int8(1):Int8(6), 4, 32),
                rand(UInt64(1):UInt64(6), 4, 32), rand(UInt32(1):UInt32(6), 4, 32), rand(UInt16(1):UInt16(6), 4, 32), rand(UInt8(1):UInt8(6), 4, 32)
            ] # 32 to force calling triggers / callbacks

            for c in eachcol(X)
                sort!(c)
            end

            for dist in [StringHammingDistance(), LevenshteinDistance(), LcsDistance(), JaccardDistance(), DiceDistance(), CosineDistanceSet()]
                for db_ in [MatrixDatabase(X), VectorDatabase(X)]
                    for db in [db_, rand(db_, 30)]
                        E = ExhaustiveSearch(; db, dist)
                        run_functions(E, db)
                        G = SearchGraph(; db, dist, verbose=false)
                        index!(G)
                        run_functions(G, db)
                        optimize!(G, ParetoRecall())
                        optimize!(G, MinRecall(0.8))
                        optimize!(G, ParetoRadius())
                        try
                            # static databases will throw an error
                            push!(G, db[1])
                            append!(G, db)
                        finally
                            continue
                        end
                    end
                end
            end
        end
        =#
    end
end

end  # end SimilaritySearch module

ext: jl

max_stars_repo_name: sadit/SimilaritySearch

max_issues_repo_name: sadit/SimilaritySearch

max_forks_repo_name: sadit/SimilaritySearch


text:
#########################################################
# Elementary actions of Abstract permutations

Base.eachindex(p::AbstractPerm) = Base.OneTo(degree(p))

@inline fixes(p::GroupElement, pt, op=^) = op(pt, p) == pt
@inline fixes(p::AbstractPerm, v::AbstractVector, op=^) =
    all( i-> v[i] == v[op(i, p)], eachindex(v))

fixedpoints(p::AbstractPerm, range=eachindex(p), op=^) =
    [i for i in range if fixes(p, i, op)]
nfixedpoints(p::AbstractPerm, range=eachindex(p), op=^) = count(i->fixes(p, i, op), range)

for (fname, findname) in [(:firstmoved, :findfirst), (:lastmoved, :findlast)]
    @eval begin
        function $fname(p::AbstractPerm, op=^)
            k = $findname(i -> i != op(i, p), eachindex(p))
            k === nothing && return k
            return eltype(p)(k)
        end
    end
end

@inline Base.:^(v::Tuple, p::AbstractPerm) = ntuple(i->v[i^p], length(v))
@inline function Base.:^(v::AbstractVector, p::AbstractPerm)
    res = similar(v);
    @inbounds for i in eachindex(v)
        res[i] = v[i^p]
    end
    return res
end

# specific definitions

@inline @inbounds function Base.:^(n::Integer, p::Perm)
    if 1 <= n <= length(p.d)
        return oftype(n, p.d[n])
    end
    return n
end

Base.:^(n::Integer, p::Permutation) = n^perm(p)


text:
################################################################################
#
#  Base field
#
################################################################################

@doc Markdown.doc"""
    base_field(L::NumField) -> NumField

Given a number field $L/K$ this function returns the base field $K$.
For absolute extensions this returns $\mathbf{Q}$.
"""
base_field(::NumField)

_base_ring(K::NumField) = base_field(K)

_base_ring(::FlintRationalField) = FlintQQ

################################################################################
#
#  Predicates
#
################################################################################

export isabsolute

@doc Markdown.doc"""
    isabsolute(L::NumField) -> Bool

Returns whether $L$ is an absolute extension, that is, whether the base field
of $L$ is $\mathbf{Q}$.
"""
isabsolute(::NumField)

isabsolute(::NumField) = false

isabsolute(::NumField{fmpq}) = true

@doc Markdown.doc"""
    elem_type(L::NumField) -> Type

Returns the type of the elements of $L$.
"""
elem_type(::NumField)

################################################################################
#
#  Degree
#
################################################################################

@doc Markdown.doc"""
    degree(L::NumField) -> Int

Given a number field $L/K$, this function returns the degree of $L$ over $K$.
"""
degree(A::NumField)

dim(K::NumField) = degree(K)

################################################################################
#
#  Absolute degree
#
################################################################################

@doc Markdown.doc"""
    absolute_degree(L::NumField) -> Int

Given a number field $L/K$, this function returns the degree of $L$ over
$\mathbf Q$.
"""
function absolute_degree(A::NumField)
  return absolute_degree(base_field(A)) * degree(A)
end

################################################################################
#
#  Is simple extension
#
################################################################################

@doc Markdown.doc"""
    issimple(L::NumField) -> Bool

Given a number field $L/K$ this function returns whether $L$ is simple, that is,
whether $L/K$ is defined by a univariate polynomial$.
"""
issimple(a::NumField)

################################################################################
#
#  Number field creation
#
################################################################################

@doc Markdown.doc"""
    NumberField(f::Poly{NumFieldElem}, s::String;
                cached::Bool = false, check::Bool = false) -> NumField, NumFieldElem

Given an irreducible polynomial $f \in K[x]$ over some number field $K$, this
function creates the simple number field $L = K[x]/(x)$ and returns $(L, b)$,
where $b$ is the class of $x$ in $L$. The string `s` is used only for printing
the primitive element $b$.

Testing that $f$ is irreducible can be disabled by setting the keyword argument
`check` to `false`.
"""
NumberField(f::PolyElem{<:NumFieldElem}, s::String;
            cached::Bool = false, check::Bool = false)

################################################################################
#
#  Is commutative
#
################################################################################

iscommutative(K::NumField) = true
iscommutative(::FlintRationalField) = true

################################################################################
#
#  Absolute field
#
################################################################################

@doc Markdown.doc"""
    absolute_field(L::NumField) -> NumField, Map

Given a number field $L$, this function returns an absolute simple number field
$M/\mathbf{Q}$ together with a $\mathbf{Q}$-linear isomorphism $M \to K$.
"""
absolute_field(::NumField)

################################################################################
#
#  Normal basis
#
################################################################################

@doc Markdown.doc"""
    normal_basis(L::NumField) -> NumFieldElem

Given a normal number field $L/K$, this function returns an element $a$ of $L$,
such that the orbit of $a$ under the Galois group of $L/K$ is an $K$-basis
of $L$.
"""
function normal_basis(L::NumField)
  n = degree(L)
  K = base_field(L)
  Aut = automorphisms(L)

  length(Aut) != degree(L) && error("The field is not normal over the rationals!")

  A = zero_matrix(K, n, n)
  r = one(L)
  while true
    r = rand(basis(L), -n:n)
    for i = 1:n
      y = Aut[i](r)
      for j = 1:n
        A[i,j] = coeff(y, j - 1)
      end
    end
    if rank(A) == n
      break
    end
  end
  return r
end


text:
function gamesize(g::AirDefenseGame)
    ncity, ndp, nap, ndpdc, ndpac, ndc, nac = g.ncity, g.ndp, g.nap, g.ndpdc, g.ndpac, g.ndc, g.nac
    nIc1 = 1
    nIa2 = binomial(ndp, ndpac)
    nIa5 = nIa2 * binomial(ndp, nac) * binomial(ndp, ndc)
    nIc3 = 1
    nId4 = binomial(ndp, ndpdc) * (2 ^ ndp)
    nIa = @sprintf("%.2E", nIa2 + nIa5)
    nId = @sprintf("%.2E", nId4)
    na1c = binomial(ndp, ndpdc)
    na2c = binomial(ndp, ndpac)
    na3a = binomial(ndpac, nac)
    na4c = 2 ^ nac
    na5d = binomial(ndpdc, ndc)
    na6a = binomial(ncity, nap)
    nas = [na1c, na2c, na3a, na4c, na5d, na6a]
    leafnodes = @sprintf("%.2E", prod(nas))
    nnodes = @sprintf("%.2E", sum(cumprod(nas)))
    (nIa, nId, leafnodes, nnodes)
end

# N : player indices (nature not included)
getplayers(g::AirDefenseGame) = 1:2

# Player index of current node (note: this is distinct from which player is currently updating CFR)
function getplayer(depth::Int64)
    (depth == 6 || depth == 3) ? 2 : (depth == 5) ? 1 : 3
    # 1 = defender, 2 = attacker, 3 = nature
end

actions(depth::Int64, An::Vector{UnitRange{Int64}}) = An[depth]
actions(depth::Int64, A::Vector{Vector{Vector{Int64}}}) = 1:length(A[depth])

function numactions(I::Int64, ni_stage, na_stage)
    if I in 1:ni_stage[3]
        return na_stage[3] # na3a
    elseif I in (ni_stage[3] + 1):(ni_stage[3] + ni_stage[5])
        return na_stage[5] # na5d
    elseif I in (ni_stage[3] + ni_stage[5] + 1):(ni_stage[3] + ni_stage[5] + ni_stage[6])
        return na_stage[6] # na5a
    else
        @warn "Infoset $I not found."
        return 0
    end
end

terminal(depth::Int64) = depth == 7

function details(h::SVector, A)
    len = findfirst(x -> x == 0, h) - 1
    (len == nothing) && (len = length(h))
    df = DataFrame(Stage = Int[], Domain = String[], Player = String[], ActionNumber = Int[], Action = Array[])
    for i = 1:len
        stage = i
        domain = (i == 6) ? "Physical" : "Cyber"
        p = getplayer(h[1:i])
        pstring = p == 2 ? "Attacker" : p == 1 ? "Defender" : "Chance"
        push!(df, (stage, domain, pstring, h[i], A[i][h[i]]))
    end
    df
end

function infoset(h::SVector, depth::Int64, g::AirDefenseGame, gs::GameSet) # btime: ~150 ns
    # findfirst benchmarks faster than numerically calculating infoset (see earlier commit)
    A, ni_stage, na_stage, sensors, iadsdefenses = gs.A, gs.ni_stage, gs.na_stage, gs.sensors, gs.iadsdefenses
    cyber_attacks, cyber_defenses = A[3], A[5]
    ndp, nac = g.ndp, g.nac
    len = depth - 1
    if len == 5
        defended_iads = getindex(gs.A[1][h[1]], gs.A[5][h[5]])
        idc = findfirst(x -> x == defended_iads, iadsdefenses)
        niadsdefenses = length(gs.iadsdefenses)
        return sum(ni_stage[1:5]) + (h[2] - 1) * na_stage[3] * niadsdefenses + (h[3] - 1) * niadsdefenses + idc
    elseif len == 4
        nsensor = length(sensors) # number of combinations of sensed cyber attacks
        sensed = getindex(gs.A[2][h[2]], getindex(A[3][h[3]], A[4][h[4]] .== 2)) # ndp id's of sensed attacks
        isensor = findfirst(x -> x == sensed, sensors) # infoset id of that combination
        return sum(ni_stage[1:4]) + (h[1] - 1) * nsensor + isensor
    elseif len == 2
        return h[2]
    else
        player = getplayer(h)
        @warn "No infoset found. Chance node? h = $h, player = $player."
        return 0
    end
end

#TODO: Make these probabilities more rational?
function chance(depth::Int64, a::Int64, g::AirDefenseGame, gs::GameSet)
    len = depth - 1
    if len == 3 # nature determines defense cyber SA
        nsensed = sum(gs.A[4][a] .== 2) # number of sensed attacks
        return (g.psc ^ nsensed) * (1 - g.psc)^(g.nac - nsensed)
    elseif len == 1 # nature deals offensive capes
        return 1.0 / gs.na_stage[2]
    elseif len == 0 # nature deals defensive capes
        return 1.0 / gs.na_stage[1]
    else
        @warn "Node at depth $depth not a chance node."
        return 0.0
    end
end

function leafutility(h::SVector, player::Int64, πi::Float64, πo::Float64, g::AirDefenseGame, gs::GameSet)
    c1, c2, ac, c4, dc, ap = h
    attacked_cities = @view gs.A[6][h[6]][:]
    defended_iads = getindex(g.iads, getindex(gs.A[1][h[1]], gs.A[5][h[5]]))
    attacked_iads = getindex(g.iads, getindex(gs.A[2][h[2]], gs.A[3][h[3]]))
    leaf_value = utility(g, gs.C, defended_iads, attacked_iads, attacked_cities)
    u = player == 1 ? leaf_value * πo : -1.0 * leaf_value * πo
    return u
end

# convert behavioral strategy to realization strategy (i.e. sequence form?)
function real_strat_calc(σ, ni, ni_stage, ns1, ns2, seqI1, seqI2, nextseqI1, nextseqI2)
    length(σ) != ni && @warn "Unexpected behavioral strategy length: $(length(σ))."
    σ1 = σ[(ni_stage[3] + 1):sum(ni_stage[1:5])]
    σ2 = vcat(σ[1:ni_stage[3]], σ[(sum(ni_stage[1:5]) + 1):end])
    rp1 = zeros(Float64, ns1)
    rp2 = zeros(Float64, ns2)
    rp1[1], rp2[1] = 1.0, 1.0
    for (i, s1) in enumerate(σ1)
        rp1[nextseqI1[i]] .= rp1[seqI1[i]] .* s1
    end
    for (i, s2) in enumerate(σ2)
        rp2[nextseqI2[i]] .= rp2[seqI2[i]] .* s2
    end
    return (rp1, rp2)
end

function real_strat_struct_calc(σ, ni, ni_stage, ns1, ns2, seqI1, seqI2, nextseqI1, nextseqI2)
    length(σ) != ni && @warn "Unexpected behavioral strategy length: $(length(σ))."
    σ1 = σ[(ni_stage[3] + 1):sum(ni_stage[1:5])]
    σ2 = vcat(σ[1:ni_stage[3]], σ[(sum(ni_stage[1:5]) + 1):end])
    rp1_struct, rp2_struct = similar(σ1), similar(σ2)
    rp1, rp2 = real_strat_calc(σ, ni, ni_stage, ns1, ns2, seqI1, seqI2, nextseqI1, nextseqI2)
    for (i, s1) in enumerate(σ1)
        rp1_struct[i] = rp1[nextseqI1[i]]
    end
    for (i, s2) in enumerate(σ2)
        rp2_struct[i] = rp2[nextseqI2[i]]
    end
    return (rp1_struct, rp2_struct)
end

function real_strat_struct_combined_calc(σ, ni, ni_stage, ns1, ns2, seqI1, seqI2, nextseqI1, nextseqI2)
    σ1 = σ[(ni_stage[3] + 1):sum(ni_stage[1:5])]
    σ2 = vcat(σ[1:ni_stage[3]], σ[(sum(ni_stage[1:5]) + 1):end])
    rp1_struct, rp2_struct = similar(σ1), similar(σ2)
    rp1, rp2 = real_strat_calc(σ, ni, ni_stage, ns1, ns2, seqI1, seqI2, nextseqI1, nextseqI2)
    for (i, s1) in enumerate(σ1)
        rp1_struct[i] = rp1[nextseqI1[i]]
    end
    for (i, s2) in enumerate(σ2)
        rp2_struct[i] = rp2[nextseqI2[i]]
    end
    rp_struct = similar(σ)
    rp_struct[1:ni_stage[3]] .= rp2_struct[1:ni_stage[3]]
    rp_struct[(ni_stage[3] + 1):sum(ni_stage[1:5])] .= rp1_struct
    rp_struct[(sum(ni_stage[1:5]) + 1):end] .= rp2_struct[(ni_stage[3] + 1):end]
    return rp_struct
end

function real_strat(σ::Vector{Vector{Float64}}, gs::GameSet, gos::GameOptSet)
    ni, ni_stage, ns1, ns2 = gs.ni, gs.ni_stage, gos.ns1, gos.ns2
    seqI1, seqI2, nextseqI1, nextseqI2 = gos.seqI1, gos.seqI2, gos.nextseqI1, gos.nextseqI2
    return real_strat_calc(σ, ni, ni_stage, ns1, ns2, seqI1, seqI2, nextseqI1, nextseqI2)
end

function real_strat(σ::Vector{Vector{Float64}}, g::AirDefenseGame)
    A, _, na_stage = build_action(g)
    ni, ni1, ni2, ni_stage = build_info(g)
    ns1, ns2, ns1_stage, ns2_stage = build_nseq(g, na_stage, ni_stage)
    _, _, (seqI1, seqI2), (nextseqI1, nextseqI2) = build_Iset(g, A, na_stage, ns1_stage, ns2_stage, ni1, ni2, ns1, ns2)
    return real_strat_calc(σ, ni, ni_stage, ns1, ns2, seqI1, seqI2, nextseqI1, nextseqI2)
end

function real_strat_struct(σ::Vector{Vector{Float64}}, g::AirDefenseGame)
    A, _, na_stage = build_action(g)
    ni, ni1, ni2, ni_stage = build_info(g)
    ns1, ns2, ns1_stage, ns2_stage = build_nseq(g, na_stage, ni_stage)
    _, _, (seqI1, seqI2), (nextseqI1, nextseqI2) = build_Iset(g, A, na_stage, ns1_stage, ns2_stage, ni1, ni2, ns1, ns2)
    return real_strat_struct_calc(σ, ni, ni_stage, ns1, ns2, seqI1, seqI2, nextseqI1, nextseqI2)
end

function real_strat_struct_combined(σ::Vector{Vector{Float64}}, g::AirDefenseGame)
    A, _, na_stage = build_action(g)
    ni, ni1, ni2, ni_stage = build_info(g)
    ns1, ns2, ns1_stage, ns2_stage = build_nseq(g, na_stage, ni_stage)
    _, _, (seqI1, seqI2), (nextseqI1, nextseqI2) = build_Iset(g, A, na_stage, ns1_stage, ns2_stage, ni1, ni2, ns1, ns2)
    return real_strat_struct_combined_calc(σ, ni, ni_stage, ns1, ns2, seqI1, seqI2, nextseqI1, nextseqI2)
end

function real_strat_struct_combined(σ::Vector{Vector{Float64}}, gs::GameSet, seqvals)
    ni, ni_stage = gs.ni, gs.ni_stage
    ns1, ns2, seqI1, seqI2, nextseqI1, nextseqI2 = seqvals
    return real_strat_struct_combined_calc(σ, ni, ni_stage, ns1, ns2, seqI1, seqI2, nextseqI1, nextseqI2)
end

# build consraint set on sequence-form strategy for CCFR
function build_ccfr_constraints(g, σfix, conf::Float64, sd::Float64)
    noise = Normal(0, sd)
    hw = quantile(noise, 1 - (1 - conf) / 2)
    σlb = [max.(σi .- hw, 0.0) for σi in σfix]
    σub = [min.(σi .+ hw, 1.0) for σi in σfix]
    rp2fix = real_strat_struct(σfix, g)[2]
    rp2lb = real_strat_struct(σlb, g)[2]
    rp2ub = real_strat_struct(σub, g)[2]
    !all(rp2lb .<= rp2fix .<= rp2ub) && @warn "Sequence-form strategy out of bounds."
    constraint_coeff = [[1, -1 .* rp2lb[i][j]] for i in 1:length(rp2lb), j in 1:length(rep2lb[i])]
    constraint_coeff = vcat(hcat(ones(length(rp2fix)), -1 .* rp2lb), hcat(-1 .* ones(length(rp2fix)), rp2ub))
    @show length(constraint_coeff), length(constraint_coeff[:,1])
    grad_constraint_coeff = hcat(constraint_coeff[:,1], zeros(length(constraint_coeff[:,1])))
    return constraint_coeff, grad_constraint_coeff
end

max_stars_repo_name: ajkeith/Cyber-Air-Defense

max_issues_repo_name: ajkeith/Cyber-Air-Defense


text:
# Check if file is a data file
function is_solution_restart_file(filename::String)
  # Open file for reading
  h5open(filename, "r") do file
    # If attribute "mesh_file" exists, this must be a data file
    return haskey(attributes(file), "mesh_file")
  end
end

# Use data file to extract mesh filename from attributes
function extract_mesh_filename(filename::String)
  # Open file for reading
  h5open(filename, "r") do file
    # Extract filename relative to data file
    mesh_file = read(attributes(file)["mesh_file"])

    return joinpath(dirname(filename), mesh_file)
  end
end

# Read in mesh file and return relevant data
function read_meshfile(filename::String)
  # Open file for reading
  h5open(filename, "r") do file
    # Extract basic information
    if haskey(attributes(file), "ndims")
      ndims_ = read(attributes(file)["ndims"])
    else
      ndims_ = read(attributes(file)["ndim"]) # FIXME once Trixi's 3D branch is merged & released
    end
    n_cells = read(attributes(file)["n_cells"])
    n_leaf_cells = read(attributes(file)["n_leaf_cells"])
    center_level_0 = read(attributes(file)["center_level_0"])
    length_level_0 = read(attributes(file)["length_level_0"])

    # Extract coordinates, levels, child cells
    coordinates = Array{Float64}(undef, ndims_, n_cells)
    coordinates .= read(file["coordinates"])
    levels = Array{Int}(undef, n_cells)
    levels .= read(file["levels"])
    child_ids = Array{Int}(undef, 2^ndims_, n_cells)
    child_ids .= read(file["child_ids"])

    # Extract leaf cells (= cells to be plotted) and contract all other arrays accordingly
    leaf_cells = similar(levels)
    n_cells = 0
    for cell_id in 1:length(levels)
      if sum(child_ids[:, cell_id]) > 0
        continue
      end

      n_cells += 1
      leaf_cells[n_cells] = cell_id
    end
    leaf_cells = leaf_cells[1:n_cells]

    coordinates = coordinates[:, leaf_cells]
    levels = levels[leaf_cells]

    return center_level_0, length_level_0, leaf_cells, coordinates, levels
  end
end

# Read in data file and return all relevant information
function read_datafile(filename::String)
  # Open file for reading
  h5open(filename, "r") do file
    # Extract basic information
    if haskey(attributes(file), "ndims")
      ndims_ = read(attributes(file)["ndims"])
    else
      ndims_ = read(attributes(file)["ndim"])
    end
    if haskey(attributes(file), "polydeg")
      polydeg = read(attributes(file)["polydeg"])
    else
      polydeg = read(attributes(file)["N"])
    end
    n_elements = read(attributes(file)["n_elements"])
    n_variables = read(attributes(file)["n_vars"])
    time = read(attributes(file)["time"])

    # Extract labels for legend
    labels = Array{String}(undef, 1, n_variables)
    for v = 1:n_variables
      labels[1, v] = read(attributes(file["variables_$v"])["name"])
    end

    # Extract data arrays
    n_nodes = polydeg + 1

    if ndims_ == 2
      data = Array{Float64}(undef, n_nodes, n_nodes, n_elements, n_variables)
      for v = 1:n_variables
        vardata = read(file["variables_$v"])
        @views data[:, :, :, v][:] .= vardata
      end
    elseif ndims_ == 3
      data = Array{Float64}(undef, n_nodes, n_nodes, n_nodes, n_elements, n_variables)
      for v = 1:n_variables
        vardata = read(file["variables_$v"])
        @views data[:, :, :, :, v][:] .= vardata
      end
    else
      error("Unsupported number of spatial dimensions: ", ndims_)
    end

    # Extract element variable arrays
    element_variables = Dict{String, Union{Vector{Float64}, Vector{Int}}}()
    index = 1
    while haskey(file, "element_variables_$index")
      varname = read(attributes(file["element_variables_$index"])["name"])
      element_variables[varname] = read(file["element_variables_$index"])
      index +=1
    end

    return labels, data, n_elements, n_nodes, element_variables, time
  end
end


text:
# """
# The Gaussian components of a mixture of 2d normals.  The kth component
# is weight_vec[k] * N(mu_vec[k], sigma_vec[k]).  Each vec must have the
# same length, K.
# """
# type PSFGaussianComponents{T <: Number}
#   mu_vec::Vector{Vector{T}}
#   sigma_vec::Vector{Matrix{T}}
#   weight_vec::Vector{T}
#
#   PSFGaussianComponents(
#     mu_vec::Vector{Vector{Number}},
#     sigma_vec::Vector{Matrix{Number}},
#     weight_vec::Vector{Number}) = begin
#
#     @assert length(mu_vec) == length(sigma_vec) == length(weight_vec)
#     new(mu_vec, sigma_vec, weight_vec)
#   end
# end

using Base.Test
import SloanDigitalSkySurvey: SDSS, PSF

const field_dir =
  joinpath(Pkg.dir("SloanDigitalSkySurvey"), "dat", "sample_field")
const run_num = "003900"
const camcol_num = "6"
const field_num = "0269"

using ForwardDiff
using Optim
using PyPlot

using PSF.get_x_matrix_from_psf
using PSF.render_psf
using PSF.unwrap_parameters

raw_psf_comp =
  SDSS.load_psf_data(field_dir, run_num, camcol_num, field_num, 1);
pixel_keep = 10:42

psf = PSF.get_psf_at_point(10.0, 10.0, raw_psf_comp);
psf = psf[pixel_keep, pixel_keep];
psf_2 = PSF.get_psf_at_point(500.0, 500.0, raw_psf_comp);
psf_2 = psf_2[pixel_keep, pixel_keep];

x_mat = PSF.get_x_matrix_from_psf(psf);

verbose = true

psf_max = maximum(psf)
#matshow(psf, vmax=1.2 * psf_max); PyPlot.colorbar(); PyPlot.title("PSF")

K = 2

# Why is NM the best?
# K=3 is much slower and not much better than K=2 by the looks of it.
opt_result, mu_vec, sigma_vec, weight_vec = PSF.fit_psf_gaussians_least_squares(psf, K=K,
  optim_method=:nelder_mead, verbose=true);

initial_par = opt_result.minimum

function evaluate_fit{T <: Number}(par::Vector{T})
  gmm_psf = render_psf(par, x_mat)
  local fit = sum((psf .- gmm_psf) .^ 2)
  if verbose && T == Float64
    println("-------------------")
    println("Fit: $fit")
    mu_vec, sigma_vec, weight_vec = unwrap_parameters(par)
    println(mu_vec)
    println(sigma_vec)
    println(weight_vec)
  end
  fit
end

# ForwardDiff reverses the arguments relative to what Optim expects.
evaluate_fit_g_rev! = ForwardDiff.gradient(evaluate_fit, mutates=true)
function evaluate_fit_g!(x::Vector{Float64}, output::Vector{Float64})
  evaluate_fit_g_rev!(output, x)
end

evaluate_fit_h_rev! = ForwardDiff.hessian(evaluate_fit, mutates=true)
function evaluate_fit_h!(x::Vector{Float64}, output::Matrix{Float64})
  evaluate_fit_h_rev!(output, x)
end

@time optim_result_newton =
  Optim.optimize(evaluate_fit, evaluate_fit_g!, evaluate_fit_h!,
                initial_par, method=Optim.Newton())

# @time optim_result_bfgs =
#   Optim.optimize(evaluate_fit, evaluate_fit_g!, evaluate_fit_h!,
#                 initial_par, method=Optim.BFGS())

@time optim_result_nelder =
    Optim.optimize(evaluate_fit, evaluate_fit_g!, evaluate_fit_h!,
                    initial_par, method=:nelder_mead)
h = ForwardDiff.hessian(evaluate_fit, initial_par)
maximum(eig(h)[1]) / minimum(eig(h)[1])

opt_result_noop = PSF.fit_psf_gaussians_least_squares(psf, K=2, initial_par=opt_result.minimum,
  optim_method=:nelder_mead, verbose=true, iterations=1);

# Why isn't this better?  It's because NM has no derivative information and so
# doesn't stay near the optimum.
opt_result_2 = fit_psf_gaussians_least_squares(psf_2, K=2, initial_par=opt_result.minimum,
  optim_method=:nelder_mead, verbose=true);

# Coordinate ascent?

mu_vec = Array(Vector{Float64}, 2)
sigma_vec = Array(Matrix{Float64}, 2)
weight_vec = Array(Float64, 2)

opt_result1, mu1, sigma1, weight1 =
  PSF.fit_psf_gaussians_least_squares(psf, K=1,
    optim_method=:nelder_mead, verbose=true);
gmm_psf1 = render_psf(opt_result1.minimum, x_mat);

opt_result2, mu2, sigma2, weight2 =
  PSF.fit_psf_gaussians_least_squares(psf - gmm_psf1, K=1,
    optim_method=:nelder_mead, verbose=true);

mu_vec[1] = mu1[1]
mu_vec[2] = mu2[1]
sigma_vec[1] = sigma1[1]
sigma_vec[2] = sigma2[1]
weight_vec[1] = weight1[1]
weight_vec[2] = weight2[1]

initial_par = PSF.wrap_parameters(mu_vec, sigma_vec, weight_vec)

@time optim_result_newton =
  Optim.optimize(evaluate_fit, evaluate_fit_g!, evaluate_fit_h!,
                initial_par, method=Optim.Newton())

# opt_result = fit_psf_gaussians_least_squares(psf, initial_par=opt_result.minimum, K=2,
#   optim_method=Optim.AcceleratedGradientDescent(), verbose=true, iterations=20);
#

unwrap_parameters(opt_result.minimum)
gmm_psf = render_psf(opt_result.minimum, x_mat);
matshow(gmm_psf, vmax=1.2 * psf_max); PyPlot.colorbar(); PyPlot.title("fit1")
psf_residual = psf - gmm_psf;
resid_max = 1.5 * maximum(abs(psf_residual))
matshow(psf_residual, vmax=resid_max, vmin=-resid_max)
PyPlot.colorbar(); PyPlot.title("residual")
sum(psf_residual .^ 2)

# Compare to EM
gmm, scale = PSF.fit_psf_gaussians_em(psf);
em_psf = [ PSF.evaluate_gmm(gmm, x_mat[i, j]')[1] for
           i=1:size(x_mat, 1), j=1:size(x_mat, 2)];
em_psf_residual = psf - em_psf;
matshow(em_psf_residual, vmax=resid_max, vmin=-resid_max)
PyPlot.colorbar(); PyPlot.title("em residual")
sum(em_psf_residual .^ 2)

mu_vec = Array(Vector{T}, K)
sigma_vec = Array(Matrix{T}, K)
weight_vec = zeros(T, K)

for k = 1:K
  mu_vec[k] = gmm.μ[1, :][:]
  sigma_vec[k] = sigma_chol' * sigma_chol + sigma_min
  weight_vec[k] = exp(par[offset + 6]) + weight_min
end

PyPlot.close("all")

psf_max = maximum(psf)
matshow(psf, vmax=1.2 * psf_max); PyPlot.colorbar(); PyPlot.title("PSF")

opt_result_1 = fit_psf_gaussians_least_squares(psf);
unwrap_parameters(opt_result_1.minimum)
gmm_psf1 = render_psf(opt_result_1.minimum);
psf_residual1 = psf - gmm_psf1;
matshow(gmm_psf1, vmax=1.2 * psf_max); PyPlot.colorbar(); PyPlot.title("fit1")
resid_max = 1.5 * maximum(abs(psf_residual1))
matshow(psf_residual1, vmax=resid_max, vmin=-resid_max)
PyPlot.colorbar(); PyPlot.title("residual1")

opt_result_2 = fit_psf_gaussians_least_squares(psf_residual1);
unwrap_parameters(opt_result_2.minimum)
gmm_psf2 = gmm_psf1 + render_psf(opt_result_2.minimum);
psf_residual2 = psf - gmm_psf2;
#matshow(gmm_psf2); PyPlot.colorbar(); PyPlot.title("fit2")
matshow(psf_residual2, vmax=resid_max, vmin=-resid_max)
PyPlot.colorbar(); PyPlot.title("residual2")

opt_result_3 = fit_psf_gaussians_least_squares(psf_residual2);
unwrap_parameters(opt_result_3.minimum)
gmm_psf3 = gmm_psf2 + render_psf(opt_result_3.minimum);
psf_residual3 = psf - gmm_psf3;
matshow(psf_residual3, vmax=resid_max, vmin=-resid_max)
PyPlot.colorbar(); PyPlot.title("residual2")


text:
include( "initpop.jl" )
include( "transitionflux.jl" )
include( "nodepop.jl" )


text:
module Datasets
using Random
using DelimitedFiles: readdlm

include("load_datasets.jl")
export load_iris

include("generators.jl")
export make_spiral,
       make_poly,
       make_sin

end

max_stars_repo_name: JuliaML/MLBase

max_issues_repo_name: JuliaML/MLBase


text:
"""
    range(start => stop)

Create a [`PartialRange`] method with a `start` and `stop` but unspecified
`step` or `length`. The `PartialRange` is callable where `length` can be
specified as the sole argument or `step` can be specified as a keyword argument.
"""
function Base.range(r::Pair; kwargs...)
    r = PartialRange(r.first, r.second)
    if !isempty(kwargs)
        return r(;kwargs...)
    end
    r
end

"""
    range(start => stop, length)

Create a range from `start` to `stop` (inclusive) with the specified length.

Equivalent to range(start, stop; length)
"""
Base.range(r::Pair, length) = Base.range(r.first, r.second; length)

"""
    range(length, start => stop)

Create a range from `start` to `stop` (inclusive) with the specified length.

Equivalent to range(start, stop; length)
"""
Base.range(length, r::Pair) = Base.range(r.first, r.second; length)


text:
__precompile__()
module FlexDates

using Dates: Date, Day, days, value

using DiscreteRanges: DiscreteRange
import DiscreteRanges: isdiscrete, discrete_gap, discrete_next

import Base:
    convert, eltype, show,
    isless, typemin, typemax,
    ==, +, -, zero, oneunit,
    hash, promote, length

export FlexDate

"""
    FlexDate{epoch, T}(date::Date)
    FlexDate{epoch, T}(year, month, day)
    FlexDate{epoch}(days::T)

Construct an object representing a given date as the number of days (stored in
type `T`) since `epoch`, which is a `Date` object.

When `T` is to narrow to represent a particular date, an `InexactError` is
thrown.

For working with a particular dataset, it is convenient to define a type
constant and work with that instead:

```julia
const MyDate = FlexDate{Date(2001, 1, 1), Int16}

MyDate(2010, 4, 31)
```

Limited arithmetic and comparisons are supported, but are fastest when using the
same epoch (unless conversion to `Date` may occur). For timespans, use `..`
(DiscreteRange) from `DiscreteRanges`.
"""
struct FlexDate{E, T <: Integer}
    Δ::T
    function FlexDate{E}(Δ::T) where {E, T <: Integer}
        @assert E isa Date
        new{E, T}(Δ)
    end
end

FlexDate{E, T}(Δ::Integer) where {E, T <: Integer} = FlexDate{E}(T(Δ))

FlexDate{E, T}(year, month, day) where {E, T} =
    FlexDate{E, T}(Date(year, month, day))

_print_ET(io, E, T) = print(io, " [$(E) + $(T) days]")

function show(io::IO, flexdate::FlexDate{E, T}) where {E, T}
    print(io, convert(Date, flexdate))
    get(io, :compact, false) || _print_ET(io, E, T)
end

function promote(x::FlexDate{Ex, <: Integer},
                 y::FlexDate{Ey, <: Integer }) where {Ex, Ey}
    xΔ, yΔ = promote(x.Δ, y.Δ)
    FlexDate{Ex}(xΔ), FlexDate{Ex}(oftype(yΔ, yΔ + days(Ey-Ex)))
end

function promote(x::FlexDate{E, <: Integer}, y::FlexDate{E, <: Integer}) where E
    xΔ, yΔ = promote(x.Δ, y.Δ)
    FlexDate{E}(xΔ), FlexDate{E}(yΔ)
end

function convert(::Type{Date}, flexdate::FlexDate{E}) where E
    E + Day(flexdate.Δ)
end

function convert(::Type{FlexDate{E,T}}, date::Date) where {E, T}
    FlexDate{E}(T(days(date - E)))
end

FlexDate{E,T}(date::Date) where {E,T} = convert(FlexDate{E,T}, date)

isless(x::FlexDate{E}, y::FlexDate{E}) where E = isless(x.Δ, y.Δ)

isless(x::FlexDate{E1}, y::FlexDate{E2}) where {E1, E2} =
    isless(convert(Date, x), convert(Date, y))

isless(x::FlexDate{E}, y::Date) where E = isless(convert(Date, x), y)

isless(x::Date, y::FlexDate{E}) where E = isless(x, convert(Date, y))

typemin(::Type{FlexDate{E, T}}) where {E, T} = FlexDate{E, T}(typemin(T))

typemax(::Type{FlexDate{E, T}}) where {E, T} = FlexDate{E, T}(typemax(T))

(==)(x::FlexDate{E}, y::FlexDate{E}) where E = x.Δ == y.Δ

(==)(x::FlexDate{E1}, y::FlexDate{E2}) where {E1, E2} =
    convert(Date, x) == convert(Date, y)

(==)(x::FlexDate{E}, y::Date) where E = convert(Date, x) == y

(==)(x::Date, y::FlexDate{E}) where E = x == convert(Date, y)

eltype(x::FlexDate{E, T}) where {E, T} = T

hash(x::FlexDate, h::UInt) = hash(convert(Date, x), h)

(-)(x::FlexDate{E, T}, y::FlexDate{E, T}) where {E, T} = Day(x.Δ - y.Δ)

# NOTE: arithmetic remains in the given type T, over/underflow is an error
# convert to `Date` for calculations with large spans, FlexDate is a storage
# format.
(+)(x::FlexDate{E, T}, y::Day) where {E, T} = FlexDate{E}(T(x.Δ + value(y)))

(-)(x::FlexDate{E, T}, y::Day) where {E, T} = FlexDate{E}(T(x.Δ - value(y)))

# support for DiscreteRanges

function show(io::IO, D::DiscreteRange{FlexDate{E,T}}) where {E,T}
    print(io, convert(Date, D.left))
    print(io, "..")
    print(io, convert(Date, D.right))
    get(io, :compact, false) || _print_ET(io, E, T)
end

isdiscrete(::Type{<:FlexDate}) = true
discrete_gap(x::FlexDate{E,T}, y::FlexDate{E,T}) where {E,T} = x.Δ - y.Δ
discrete_next(x::FlexDate{E,T}, Δ) where {E,T} = FlexDate{E}(T(x.Δ + Δ))

end # module


text:
module GeometryModule
export Geometry, readGeometryXYZ, angstrom2bohr, angstrom2bohr!, computeEnergyInteratomicRepulsion

using ..BaseModule
using ..AtomModule

type Geometry
	atoms::Array{Atom,1}
end

function readGeometryXYZ(filename::AbstractString; unit::Symbol=:Bohr)
	geo = Geometry([])
	lines = open(readlines,filename)
	for line in lines[3:end]	# line 1 is n. of atoms, line 2 is comment
		columns = split(line)
		elem = Element(columns[1])
		pos = Position(map(parse,columns[2:4])...)	# float is not recommended for use on tuples
		append!(geo.atoms,[Atom(elem,pos)])
	end
	if (unit == :Bohr) angstrom2bohr!(geo) end	# Should be default as Bohr is an atomic unit
	return geo
end

angstrom2bohr(x::Real) = (return x/0.529177210)
angstrom2bohr(pos::Position) = (return Position(angstrom2bohr(pos.x),angstrom2bohr(pos.y),angstrom2bohr(pos.z)))
function angstrom2bohr!(geo::Geometry)
  for (atom in geo.atoms)
    atom.position = angstrom2bohr(atom.position)
  end
end

function computeEnergyInteratomicRepulsion(
  atom1::Atom,
  atom2::Atom)
  q1 = atom1.element.atomicNumber
  q2 = atom2.element.atomicNumber
  A = atom1.position
  B = atom2.position
  return q1*q2/(distance(A,B))
end

function computeEnergyInteratomicRepulsion(geo::Geometry)
  result = 0.::Float64
  for (A in 1:length(geo.atoms))
    for (B in 1:A-1)
      result += computeEnergyInteratomicRepulsion(geo.atoms[A],geo.atoms[B])
    end
  end
  return result
end

end #module


text:
include("common.jl")

facts("Optimize single objective problems in 5, 10, and 30 dimensions with DE") do
  simple_problems = ["Sphere", "Schwefel2.22", "Schwefel2.22"]
  for(problem in simple_problems)
    context(problem) do
      p = BlackBoxOptim.example_problems[problem]

      @fact fitness_for_opt(p, 5, 20, 5e3, :de_rand_1_bin) < 0.01 => true
      @fact fitness_for_opt(p, 5, 20,  5e3, :de_rand_1_bin_radiuslimited) < 0.01 => true

      @fact fitness_for_opt(p, 10, 20, 1e4, :de_rand_1_bin) < 0.01 => true
      @fact fitness_for_opt(p, 10, 20, 1e4, :de_rand_1_bin_radiuslimited) < 0.01 => true

      @fact fitness_for_opt(p, 30, 25, 3e4, :de_rand_1_bin) < 0.01 => true
      @fact fitness_for_opt(p, 30, 25, 3e4, :de_rand_1_bin_radiuslimited) < 0.01 => true
    end
  end

  context("Schwefel1.2") do
    problem = "Schwefel1.2"
    p = BlackBoxOptim.example_problems[problem]
    @fact fitness_for_opt(p, 5, 20,  5e3, :de_rand_1_bin_radiuslimited) < 0.01 => true
    @fact fitness_for_opt(p, 10, 50, 5e4, :de_rand_1_bin_radiuslimited) < 0.01 => true

    #DE/rand/1/bin seems to have troubles...
    @fact fitness_for_opt(p, 30, 50, 2e5, :de_rand_1_bin_radiuslimited) < 10.0 => true
    @fact fitness_for_opt(p, 30, 50, 2e5, :adaptive_de_rand_1_bin) < 10.0 => true
    @fact fitness_for_opt(p, 30, 50, 2e5, :adaptive_de_rand_1_bin_radiuslimited) < 10.0 => true
  end

  context("Rosenbrock") do
    problem = "Rosenbrock"
    p = BlackBoxOptim.example_problems[problem]
    @fact fitness_for_opt(p, 5, 20,   1e4, :de_rand_1_bin_radiuslimited) < 100.0 => true
    @fact fitness_for_opt(p, 10, 20,  5e4, :de_rand_1_bin_radiuslimited) < 100.0 => true
    @fact fitness_for_opt(p, 30, 40, 2e5, :de_rand_1_bin_radiuslimited) < 100.0 => true

    @fact fitness_for_opt(p, 30, 40, 2e5, :adaptive_de_rand_1_bin) < 100.0 => true
    @fact fitness_for_opt(p, 30, 40, 2e5, :adaptive_de_rand_1_bin_radiuslimited) < 100.0 => true

    @fact fitness_for_opt(p, 50, 40, 3e5, :adaptive_de_rand_1_bin_radiuslimited) < 100.0 => true
  end
end


text:
"""
The `HTTP.Servers` module provides HTTP server functionality.

The main entry point is `HTTP.listen(f, host, port; kw...)` which takes
a `f(::HTTP.Stream)::Nothing` function argument, a `host`, a `port` and
optional keyword arguments.  For full details, see `?HTTP.listen`.

For server functionality operating on full requests, see `?HTTP.Handlers`
module and `?HTTP.serve` function.
"""
module Servers

export listen

using ..IOExtras
using ..Streams
using ..Messages
using ..Parsers
using ..ConnectionPool
using Sockets
using MbedTLS: SSLContext, SSLConfig
import MbedTLS

using Dates

# rate limiting
mutable struct RateLimit
    allowance::Float64
    lastcheck::Dates.DateTime
end

function update!(rl::RateLimit, rate_limit)
    current = Dates.now()
    timepassed = float(Dates.value(current - rl.lastcheck)) / 1000.0
    rl.lastcheck = current
    rl.allowance += timepassed * rate_limit
    return nothing
end

const RATE_LIMITS = [Dict{IPAddr, RateLimit}()]
check_rate_limit(tcp::Base.PipeEndpoint, rate_limit::Rational{Int}) = true
check_rate_limit(tcp, ::Nothing) = true

"""
`check_rate_limit` takes a new connection (socket), and checks in
the global RATE_LIMITS store for the last time a connection was
seen for the same ip address. If the new connection has come too
soon, it is closed and discarded, otherwise, the timestamp for the
ip address is updated in the global cache.
"""
function check_rate_limit(tcp, rate_limit::Rational{Int})
    ip = Sockets.getsockname(tcp)[1]
    rate = Float64(rate_limit.num)
    rl = get!(RATE_LIMITS[Threads.threadid()], ip, RateLimit(rate, Dates.now()))
    update!(rl, rate_limit)
    if rl.allowance > rate
        @warn "throttling $ip"
        rl.allowance = rate
    end
    if rl.allowance < 1.0
        @warn "discarding connection from $ip due to rate limiting"
        return false
    else
        rl.allowance -= 1.0
    end
    return true
end

"Convenience object for passing around server details"
struct Server{S <: Union{SSLConfig, Nothing}, I <: Base.IOServer}
    ssl::S
    server::I
    hostname::String
    hostport::String
end

Base.isopen(s::Server) = isopen(s.server)
Base.close(s::Server) = close(s.server)

Sockets.accept(s::Server{Nothing}) = accept(s.server)::TCPSocket
Sockets.accept(s::Server{SSLConfig}) = getsslcontext(accept(s.server), s.ssl)

function getsslcontext(tcp, sslconfig)
    ssl = MbedTLS.SSLContext()
    MbedTLS.setup!(ssl, sslconfig)
    MbedTLS.associate!(ssl, tcp)
    MbedTLS.handshake!(ssl)
    return ssl
end

"""
    HTTP.listen([host=Sockets.localhost[, port=8081]]; kw...) do http::HTTP.Stream
        ...
    end

Listen for HTTP connections and execute the `do` function for each request.

The `do` function should be of the form `f(::HTTP.Stream)::Nothing`.

Optional keyword arguments:
 - `sslconfig=nothing`, Provide an `MbedTLS.SSLConfig` object to handle ssl
    connections. Pass `sslconfig=MbedTLS.SSLConfig(false)` to disable ssl
    verification (useful for testing).
 - `reuse_limit = nolimit`, number of times a connection is allowed to be
   reused after the first request.
 - `tcpisvalid = tcp->true`, function `f(::TCPSocket)::Bool` to, check accepted
    connection before processing requests. e.g. to do source IP filtering.
 - `readtimeout::Int=60`, close the connection if no data is received for this
    many seconds. Use readtimeout = 0 to disable.
 - `reuseaddr::Bool=false`, allow multiple servers to listen on the same port.
 - `server::Base.IOServer=nothing`, provide an `IOServer` object to listen on;
    allows closing the server.
 - `connection_count::Ref{Int}`, reference to track the number of currently
    open connections.
 - `rate_limit::Rational{Int}=nothing"`, number of `connections//second`
    allowed per client IP address; excess connections are immediately closed.
    e.g. 5//1.
 - `verbose::Bool=false`, log connection information to `stdout`.

e.g.
```julia
HTTP.listen("127.0.0.1", 8081) do http
    HTTP.setheader(http, "Content-Type" => "text/html")
    write(http, "target uri: \$(http.message.target)<BR>")
    write(http, "request body:<BR><PRE>")
    write(http, read(http))
    write(http, "</PRE>")
    return
end

HTTP.listen("127.0.0.1", 8081) do http
    @show http.message
    @show HTTP.header(http, "Content-Type")
    while !eof(http)
        println("body data: ", String(readavailable(http)))
    end
    HTTP.setstatus(http, 404)
    HTTP.setheader(http, "Foo-Header" => "bar")
    startwrite(http)
    write(http, "response body")
    write(http, "more response body")
end
```

The `server=` option can be used to pass an already listening socket to
`HTTP.listen`. This allows manual control of server shutdown.

e.g.
```julia
using Sockets
server = Sockets.listen(Sockets.InetAddr(parse(IPAddr, host), port))
@async HTTP.listen(f, host, port; server=server)

# Closing server will stop HTTP.listen.
close(server)
```

To run the following HTTP chat example, open two Julia REPL windows and paste
the example code into both of them. Then in one window run `chat_server()` and
in the other run `chat_client()`, then type `hello` and press return.
Whatever you type on the client will be displayed on the server and vis-versa.

```
using HTTP

function chat(io::HTTP.Stream)
    @async while !eof(io)
        write(stdout, readavailable(io), "\\n")
    end
    while isopen(io)
        write(io, readline(stdin))
    end
end

chat_server() = HTTP.listen("127.0.0.1", 8087) do io
    write(io, "HTTP.jl Chat Server. Welcome!")
    chat(io)
end

chat_client() = HTTP.open("POST", "http://127.0.0.1:8087") do io
    chat(io)
end
```
"""
function listen end

const nolimit = typemax(Int)

getinet(host::String, port::Integer) = Sockets.InetAddr(parse(IPAddr, host), port)
getinet(host::IPAddr, port::Integer) = Sockets.InetAddr(host, port)

function listen(f,
                host::Union{IPAddr, String}=Sockets.localhost,
                port::Integer=8081
                ;
                sslconfig::Union{MbedTLS.SSLConfig, Nothing}=nothing,
                tcpisvalid::Function=tcp->true,
                server::Union{Base.IOServer, Nothing}=nothing,
                reuseaddr::Bool=false,
                connection_count::Ref{Int}=Ref(0),
                rate_limit::Union{Rational{Int}, Nothing}=nothing,
                reuse_limit::Int=nolimit,
                readtimeout::Int=0,
                verbose::Bool=false)

    inet = getinet(host, port)
    if server !== nothing
        tcpserver = server
    elseif reuseaddr
        tcpserver = Sockets.TCPServer(; delay=false)
        if Sys.isunix()
            if Sys.isapple()
                verbose && @warn "note that `reuseaddr=true` allows multiple processes to bind to the same addr/port, but only one process will accept new connections (if that process exits, another process listening will start accepting)"
            end
            rc = ccall(:jl_tcp_reuseport, Int32, (Ptr{Cvoid},), tcpserver.handle)
            Sockets.bind(tcpserver, inet.host, inet.port; reuseaddr=true)
        else
            @warn "reuseaddr=true may not be supported on this platform: $(Sys.KERNEL)"
            Sockets.bind(tcpserver, inet.host, inet.port; reuseaddr=true)
        end
        Sockets.listen(tcpserver)
    else
        tcpserver = Sockets.listen(inet)
    end
    verbose && @info "Listening on: $host:$port"

    tcpisvalid = let f=tcpisvalid
        x -> f(x) && check_rate_limit(x, rate_limit)
    end

    s = Server(sslconfig, tcpserver, string(host), string(port))
    return listenloop(f, s, tcpisvalid, connection_count,
                         reuse_limit, readtimeout, verbose)
end

""""
Main server loop.
Accepts new tcp connections and spawns async tasks to handle them."
"""
function listenloop(f, server, tcpisvalid, connection_count,
                       reuse_limit, readtimeout, verbose)
    count = 1
    while isopen(server)
        try
            io = accept(server)
            if !tcpisvalid(io)
                verbose && @info "Accept-Reject:  $io"
                close(io)
                continue
            end
            connection_count[] += 1
            conn = Connection(io)
            conn.host, conn.port = server.hostname, server.hostport
            let io=io, count=count
                @async try
                    verbose && @info "Accept ($count):  $conn"
                    handle_connection(f, conn, server, reuse_limit, readtimeout)
                    verbose && @info "Closed ($count):  $conn"
                catch e
                    if e isa Base.IOError && e.code == -54
                        verbose && @warn "connection reset by peer (ECONNRESET)"
                    else
                        @error exception=(e, stacktrace(catch_backtrace()))
                    end
                finally
                    connection_count[] -= 1
                    # handle_connection is in charge of closing the underlying io
                end
            end
        catch e
            close(server)
            if e isa InterruptException
                @warn "Interrupted: listen($server)"
                break
            else
                rethrow(e)
            end
        end
        count += 1
    end
    return
end

"""
Start a `check_readtimeout` task to close the `Connection` if it is inactive.
Create a `Transaction` object for each HTTP Request received.
After `reuse_limit + 1` transactions, signal `final_transaction` to the
transaction handler.
"""
function handle_connection(f, c::Connection, server, reuse_limit, readtimeout)
    wait_for_timeout = Ref{Bool}(true)
    if readtimeout > 0
        @async check_readtimeout(c, readtimeout, wait_for_timeout)
    end
    try
        count = 0
        # if the connection socket or original server close, we stop taking requests
        while isopen(c) && isopen(server) && count <= reuse_limit
            handle_transaction(f, Transaction(c);
                               final_transaction=(count == reuse_limit))
            count += 1
        end
    finally
        wait_for_timeout[] = false
    end
    return
end

"""
If `c` is inactive for a more than `readtimeout` then close the `c`."
"""
function check_readtimeout(c, readtimeout, wait_for_timeout)
    while wait_for_timeout[]
        if inactiveseconds(c) > readtimeout
            @warn "Connection Timeout: $c"
            try
                writeheaders(c.io, Response(408, ["Connection" => "close"]))
            finally
                close(c)
            end
            break
        end
        sleep(8 + rand() * 4)
    end
    return
end

"""
Create a `HTTP.Stream` and parse the Request headers from a `HTTP.Transaction`
(by calling `startread(::Stream`).
If there is a parse error, send an error Response.
Otherwise, execute stream processing function `f`.
If `f` throws an exception, send an error Response and close the connection.
"""
function handle_transaction(f, t::Transaction; final_transaction::Bool=false)
    request = Request()
    http = Stream(request, t)

    try
        startread(http)
    catch e
        if e isa EOFError && isempty(request.method)
            return
        elseif e isa ParseError
            status = e.code == :HEADER_SIZE_EXCEEDS_LIMIT  ? 413 : 400
            write(t, Response(status, body = string(e.code)))
            close(t)
            return
        else
            rethrow(e)
        end
    end

    request.response.status = 200
    if final_transaction || hasheader(request, "Connection", "close")
        setheader(request.response, "Connection" => "close")
    end

    @async try
        f(http)
        closeread(http)
        closewrite(http)
    catch e
        @error "error handling request" exception=(e, stacktrace(catch_backtrace()))
        if isopen(http) && !iswritable(http)
            http.message.response.status = 500
            startwrite(http)
            write(http, sprint(showerror, e))
        end
        final_transaction = true
    finally
        final_transaction && close(t.c.io)
    end
    return
end

function __init__()
    Threads.resize_nthreads!(RATE_LIMITS)
    return
end

end # module


text:
@testset "Check" begin
    game = MatrixGame()

    #=
    Wrong input size:
    First dense layer has dimension 2, but should be 1
    =#
    strat = Chain(Dense(2,20, sigmoid), Dense(20,20, relu), Dense(20,3), softmax)
    v_net = Chain(Dense(2,20, sigmoid), Dense(20,20, relu), Dense(20,3))

    sol = DeepCFRSolver(
        game;
        buffer_size = 100*10^3,
        batch_size = 128,
        traversals = 10,
        strategy = strat,
        values = (v_net, deepcopy(v_net)),
    )

    @test_throws DimensionMismatch DeepCFR.checksolver(sol)

    #=
    Wrong output size:
    output layer has dimension 4, but should be 3
    =#
    strat = Chain(Dense(1,20, sigmoid), Dense(20,20, relu), Dense(20,4), softmax)
    v_net = Chain(Dense(1,20, sigmoid), Dense(20,20, relu), Dense(20,4))

    sol = DeepCFRSolver(
        game;
        buffer_size = 100*10^3,
        batch_size = 128,
        traversals = 10,
        strategy = strat,
        values = (v_net, deepcopy(v_net)),
    )

    @test_throws AssertionError DeepCFR.checksolver(sol)

    #=
    Unnormalized strategy network output:
    sum of strategy network outputs != 1.0
    =#
    strat = Chain(Dense(1,20, sigmoid), Dense(20,20, relu), Dense(20,3))
    sol = DeepCFRSolver(
        game;
        buffer_size = 100*10^3,
        batch_size = 128,
        traversals = 10,
        strategy = strat,
    )

    @test_throws AssertionError DeepCFR.checksolver(sol)

    strat = Chain(Dense(1,20, sigmoid), Dense(20,20, relu), Dense(20,3), softmax)
    v_net = Chain(Dense(1,20, sigmoid), Dense(20,20, relu), Dense(20,3))

    sol = DeepCFRSolver(
        game;
        buffer_size = 100*10^3,
        batch_size = 128,
        traversals = 10,
        strategy = strat,
        values = (v_net, deepcopy(v_net)),
    )

    @test DeepCFR.checksolver(sol) == "Tests Pass"
end


text:
#=
qpds:
- Julia version:
- Author: marcin
- Date: 2018-05-25
=#

using QI

star(A, B) = sqrtm(B)*A*sqrtm(B) # operator z (#2)

abstract type MCPD{D <: MatrixDistribution} <: CPD{D} end

struct DiscreteMCPD{D <: DiscreteMatrixDistribution} <: MCPD{D}
    # name of the node
    # TODO in the future, this should also be NodeNames,
    # as this class will ultimately represent a system of entangled variables
    # such system is now a single variable
    target::NodeName

    parents::NodeNames
    parental_ncategories::Vector{Int}

    # in the future, a vector of Ints for the same reason as in target
    ncategories::Int

    # conditional distribution on parents
    conditional_distribution::D

    #TODO inner constructor with argument checking
end

function (dmcpd::DiscreteMCPD)(a::Assignment=Assignment())
    if haskey(a, dmcpd.target)
        return a[dmcpd]
    else
        if isempty(dmcpd.parents)
            distribution = dmcpd.conditional_distribution.p
        else
            parents_assignments = [a[p] for p in dmcpd.parents]
            total_distribution = foldl(kron, 1, vcat(parents_assignments, [eye(dmcpd.ncategories)]))
            distribution = star(dmcpd.conditional_distribution.p, total_distribution)
        end

        all_ncategories = vcat(dmcpd.parental_ncategories, dmcpd.ncategories)
        trace_out_indices = [p for p in 1:length(dmcpd.parents)]
        return ptrace(distribution, all_ncategories, trace_out_indices)

        # when the target is a single variable, this is enough
        # in the future, this function will have to support cases when we want a distribution of some part of athe system
        # (of only some variables in the system)
        # i imagine it will return a list of HermitianMatrices then
        # especially an empty list, when none of the quried nodes are *in* the given system

        # also, in this case, should any parts of the system be in the assignment,
        # event() function must be applied on the distribution of the system before
        # the result has been traced out from it
    end
end

name(dmcpd::DiscreteMCPD) = dmcpd.target
parents(dmcpd::DiscreteMCPD) = dmcpd.parents

# TODO
# nparams(dmcpd::DiscreteMCPD)


text:
@testset "datacontainer" begin
    n_samples = 10
    n_dim = 4
    n_out = 3
    X = [rand(n_dim) for _ in 1:n_samples]
    y = rand(n_samples)
    ys = [rand(n_samples) for _ in 1:n_out]
    ybad = rand(n_samples + 1)
    ysbad = vcat(ys, ybad)
    @testset "DataContainer" begin
        data = AGP.wrap_data(X, y)
        @test AGP.nSamples(data) == n_samples
        @test AGP.nDim(data) == n_dim
        @test AGP.nOutput(data) == 1
        @test AGP.input(data) == X
        @test AGP.output(data) == y
        @test_throws ErrorException AGP.wrap_data(X, ybad)
        # Multiple ys case
        data = AGP.wrap_data(X, ys)
        @test AGP.nSamples(data) == n_samples
        @test AGP.nDim(data) == n_dim
        @test AGP.nOutput(data) == 1
        @test AGP.input(data) == X
        @test AGP.output(data) == ys
        @test_throws ErrorException AGP.wrap_data(X, ysbad)
    end

    @testset "MODataContainer" begin
        data = AGP.wrap_modata(X, ys)
        @test AGP.nSamples(data) == n_samples
        @test AGP.nDim(data) == n_dim
        @test AGP.nOutput(data) == n_out
        @test AGP.input(data) == X
        @test AGP.output(data) == ys
        @test_throws ErrorException AGP.wrap_data(X, ysbad)
    end

    @testset "OnlineDataContainer" begin
        data = AGP.OnlineDataContainer()
        AGP.wrap_data!(data, AGP.wrap_X(X)[1], y)
        @test AGP.nSamples(data) == n_samples
        @test AGP.nDim(data) == n_dim
        @test AGP.nOutput(data) == 1
        @test AGP.input(data) == X
        @test AGP.output(data) == y
    end

    X = rand(5, 2)
    @testset "wrap_X" begin
        @test AGP.wrap_X(X, 1)[1] isa AbstractVector
        @test AGP.wrap_X(X, 2)[1] isa AbstractVector
        @test AGP.wrap_X(vec(X))[1] isa AbstractVector
        @test AGP.wrap_X(vec(X))[2] == Float64
    end

end


text:
using Documenter
using TimeParallel

PAGES = ["Home" => "index.md"]

makedocs(;
    sitename = "TimeParallel.jl",
    format = Documenter.HTML(),
    modules = [TimeParallel],
    pages = PAGES,
    authors = "Giancarlo A. Antonucci <giancarlo.antonucci@icloud.com>"
)

deploydocs(;
    repo = "https://github.com/giancarloantonucci/TimeParallel.jl"
)


text:
function rand(rng::AbstractRNG, ::Random.SamplerTrivial{Random.CloseOpen01{ArbFloat{P}}},n) where P
    bfprec = precision(BigFloat)
    setprecision(BigFloat, P)
    bfrand = rand(BigFloat,n)
    result = ArbFloat{P}.(bfrand)
    setprecision(BigFloat, bfprec); println("ok")
    return result
end

function rand(::Type{ArbFloat{P}}, n::Integer) where P
    bfprec = precision(BigFloat)
    setprecision(BigFloat, P)
    bfrand = rand(BigFloat,n)
    result = ArbFloat{P}.(bfrand)
    setprecision(BigFloat, bfprec)
    return result
end

function rand(rng::AbstractRNG, ::Random.SamplerTrivial{Random.CloseOpen01{ArbReal{P}}},n) where P
    bfprec = precision(BigFloat)
    setprecision(BigFloat, P)
    bfrand = rand(BigFloat,n)
    result = ArbReal{P}.(bfrand)
    setprecision(BigFloat, bfprec); println("ok")
    return result
end

function rand(::Type{ArbReal{P}}, n::Integer) where P
    bfprec = precision(BigFloat)
    setprecision(BigFloat, P)
    bfrand = rand(BigFloat,n)
    result = ArbReal{P}.(bfrand)
    setprecision(BigFloat, bfprec)
    return result
end


text:
module VariableTemplates

export varsize, Vars, Grad, @vars, varsindex

using StaticArrays

function varsindex(::Type{S}, insym::Symbol) where {S <: NamedTuple}
    offset = 0
    for varsym in fieldnames(S)
        T = fieldtype(S, varsym)
        if T <: Real
            offset += 1
            varrange = offset:offset
        elseif T <: SHermitianCompact
            LT = StaticArrays.lowertriangletype(T)
            N = length(LT)
            varrange = offset .+ (1:N)
            offset += N
        elseif T <: StaticArray
            N = length(T)
            varrange = offset .+ (1:N)
            offset += N
        else
            varrange = offset .+ (1:varsize(T))
            offset += varsize(T)
        end
        if insym == varsym
            return varrange
        end
    end
    error("symbol '$insym' not found")
end

"""
    varsize(S)

The number of elements specified by the template type `S`.
"""
varsize(::Type{T}) where {T <: Real} = 1
varsize(::Type{Tuple{}}) = 0
varsize(::Type{NamedTuple{(), Tuple{}}}) = 0
varsize(::Type{SVector{N, T}}) where {N, T <: Real} = N

include("var_names.jl")

# TODO: should be possible to get rid of @generated
@generated function varsize(::Type{S}) where {S}
    types = fieldtypes(S)
    isempty(types) ? 0 : sum(varsize, types)
end

function process_vars!(syms, typs, expr)
    if expr isa LineNumberNode
        return
    elseif expr isa Expr && expr.head == :block
        for arg in expr.args
            process_vars!(syms, typs, arg)
        end
        return
    elseif expr.head == :(::)
        push!(syms, expr.args[1])
        push!(typs, expr.args[2])
        return
    else
        error("Invalid expression")
    end
end

macro vars(args...)
    syms = Any[]
    typs = Any[]
    for arg in args
        process_vars!(syms, typs, arg)
    end
    :(NamedTuple{$(tuple(syms...)), Tuple{$(esc.(typs)...)}})
end

struct GetVarError <: Exception
    sym::Symbol
end
struct SetVarError <: Exception
    sym::Symbol
end

"""
    Vars{S,A,offset}(array::A)

Defines property overloading for `array` using the type `S` as a template. `offset` is used to shift the starting element of the array.
"""
struct Vars{S, A, offset}
    array::A
end
Vars{S}(array) where {S} = Vars{S, typeof(array), 0}(array)

Base.parent(v::Vars) = getfield(v, :array)
Base.eltype(v::Vars) = eltype(parent(v))
Base.propertynames(::Vars{S}) where {S} = fieldnames(S)
Base.similar(v::Vars{S, A, offset}) where {S, A, offset} =
    Vars{S, A, offset}(similar(parent(v)))

@generated function Base.getproperty(
    v::Vars{S, A, offset},
    sym::Symbol,
) where {S, A, offset}
    expr = quote
        Base.@_inline_meta
        array = parent(v)
    end
    for k in fieldnames(S)
        T = fieldtype(S, k)
        if T <: Real
            retexpr = :($T(array[$(offset + 1)]))
            offset += 1
        elseif T <: SHermitianCompact
            LT = StaticArrays.lowertriangletype(T)
            N = length(LT)
            retexpr = :($T($LT($([:(array[$(offset + i)]) for i in 1:N]...))))
            offset += N
        elseif T <: StaticArray
            N = length(T)
            retexpr = :($T($([:(array[$(offset + i)]) for i in 1:N]...)))
            offset += N
        else
            retexpr = :(Vars{$T, A, $offset}(array))
            offset += varsize(T)
        end
        push!(expr.args, :(
            if sym == $(QuoteNode(k))
                return @inbounds $retexpr
            end
        ))
    end
    push!(expr.args, :(throw(GetVarError(sym))))
    expr
end

@generated function Base.setproperty!(
    v::Vars{S, A, offset},
    sym::Symbol,
    val,
) where {S, A, offset}
    expr = quote
        Base.@_inline_meta
        array = parent(v)
    end
    for k in fieldnames(S)
        T = fieldtype(S, k)
        if T <: Real
            retexpr = :(array[$(offset + 1)] = val)
            offset += 1
        elseif T <: SHermitianCompact
            LT = StaticArrays.lowertriangletype(T)
            N = length(LT)
            retexpr = :(
                array[($(offset + 1)):($(offset + N))] .= $T(val).lowertriangle
            )
            offset += N
        elseif T <: StaticArray
            N = length(T)
            retexpr = :(array[($(offset + 1)):($(offset + N))] .= val[:])
            offset += N
        else
            offset += varsize(T)
            continue
        end
        push!(expr.args, :(
            if sym == $(QuoteNode(k))
                return @inbounds $retexpr
            end
        ))
    end
    push!(expr.args, :(throw(SetVarError(sym))))
    expr
end

"""
    Grad{S,A,offset}(array::A)

Defines property overloading along slices of the second dimension of `array` using the type `S` as a template. `offset` is used to shift the starting element of the array.
"""
struct Grad{S, A, offset}
    array::A
end
Grad{S}(array) where {S} = Grad{S, typeof(array), 0}(array)

Base.parent(g::Grad) = getfield(g, :array)
Base.eltype(g::Grad) = eltype(parent(g))
Base.propertynames(::Grad{S}) where {S} = fieldnames(S)
Base.similar(g::Grad{S, A, offset}) where {S, A, offset} =
    Grad{S, A, offset}(similar(parent(g)))

@generated function Base.getproperty(
    v::Grad{S, A, offset},
    sym::Symbol,
) where {S, A, offset}
    M = size(A, 1)
    expr = quote
        Base.@_inline_meta
        array = parent(v)
    end
    for k in fieldnames(S)
        T = fieldtype(S, k)
        if T <: Real
            retexpr = :(SVector{$M, $T}(
                $([:(array[$i, $(offset + 1)]) for i in 1:M]...),
            ))
            offset += 1
        elseif T <: StaticArray
            N = length(T)
            retexpr = :(SMatrix{$M, $N, $(eltype(T))}(
                $([:(array[$i, $(offset + j)]) for i in 1:M, j in 1:N]...),
            ))
            offset += N
        else
            retexpr = :(Grad{$T, A, $offset}(array))
            offset += varsize(T)
        end
        push!(expr.args, :(
            if sym == $(QuoteNode(k))
                return @inbounds $retexpr
            end
        ))
    end
    push!(expr.args, :(throw(GetVarError(sym))))
    expr
end

@generated function Base.setproperty!(
    v::Grad{S, A, offset},
    sym::Symbol,
    val,
) where {S, A, offset}
    M = size(A, 1)
    expr = quote
        Base.@_inline_meta
        array = parent(v)
    end
    for k in fieldnames(S)
        T = fieldtype(S, k)
        if T <: Real
            retexpr = :(array[:, $(offset + 1)] .= val)
            offset += 1
        elseif T <: StaticArray
            N = length(T)
            retexpr = :(array[:, ($(offset + 1)):($(offset + N))] .= val)
            offset += N
        else
            offset += varsize(T)
            continue
        end
        push!(expr.args, :(
            if sym == $(QuoteNode(k))
                return @inbounds $retexpr
            end
        ))
    end
    push!(expr.args, :(throw(SetVarError(sym))))
    expr
end

end # module

max_stars_repo_name: yairchn/CLIMA

max_issues_repo_name: yairchn/CLIMA

max_forks_repo_name: yairchn/CLIMA


text:
# Revised on 4/29/21 to fix approximate comparisons with 0.0.   Works.
#include("../src/CGP.jl")
##include("../src/InfTheory.jl")
using Main.CGP
using Test

# Compare 2 floats for approximate equality where one of the floats might be zero
function approx_compare( x::Float64, y::Float64 )
  max( abs(x), abs(y) ) < 1.0e-12 || x ≈ y
end

nreps = 4
function test_degen_complexity( p::Main.CGP.Parameters, nreps::Integer )
  funcs = default_funcs(p.numinputs)
  for i = 1:nreps
    println("i: ",i)
    c = Main.CGP.random_chromosome(p,funcs)
    #println("c: ")
    #print_chromosome(c)
    try
      @test degeneracy(c) ≈ degeneracy1(c)
    catch
      println((degeneracy(c),degeneracy1(c)))
    end
    try
      @test approx_compare( complexity5(c), complexity6(c) )
    catch
      println((complexity5(c),complexity6(c)))
    end
    try
      @test approx_compare( complexity5(c), complexity4(c) )
    catch
      println((complexity5(c),complexity4(c)))
    end
    try
      @test approx_compare( complexity5(c), complexity7(c) )
    catch
      println((complexity5(c),complexity7(c)))
    end
  end
end

@testset "testing multiple versions of degeneracy and complexity" begin
#=
numinputs = 2
p = Main.CGP.Parameters( numinputs=numinputs, numoutputs=2, numinteriors=6, numlevelsback=6 )
test_degen_complexity( p, nreps )
=#
numinputs = 3
p =  Main.CGP.Parameters( numinputs=numinputs, numoutputs=3, numinteriors=7, numlevelsback=6 )
test_degen_complexity( p, nreps )
end  # @testset


text:
## Load necessary code

using Pkg
Pkg.activate(@__DIR__)
pkg"instantiate"
pkg"precompile"

using FileIO

include(joinpath(@__DIR__, "CommonUtils.jl"))
# include(joinpath(@__DIR__, "Plotting.jl"))

## load specific case

pargs = if !isdefined(Main, :parsed_args)
  parse_commandline()
else
  parsed_args
end

##

pargs["reportDir"] = if !haskey(pargs, "reportDir") && isdefined(Main, :fg)
  getLogPath(fg)
else
  # "/tmp/caesar/2020-02-23T01:43:32.222/fg_before_x291.tar.gz"
  "/tmp/caesar/2020-02-23T01:43:32.222/fg_after_x1391.tar.gz"
2020-02-24T04:28:41.955
  # pargs["reportDir"]
end

# @show pargs["reportDir"]
# @show splitpath(pargs["reportDir"])

# load the factor graph needed
fg = if true || !isdefined(Main, :fg)
  println("going to load fg")
  fg = LightDFG{SolverParams}(params=SolverParams())
  @show pathElem = splitpath(pargs["reportDir"])
  @show getSolverParams(fg).logpath = joinpath(pathElem[1:end-1]...)
  loadDFG(pargs["reportDir"], Main, fg)
  fg
else
  fg
end

# random temp data
drt_data, TTm, XXm, YYm, XXf, YYf = loadResultsDRT(fg)

# load here due to name conflict with Point2
using MakieLayout
using Makie

include(joinpath(@__DIR__, "MakiePlotsFG.jl"))

## get a portion of the graph

fg2 = deepcopy(fg)

# trim down the factor graph
klist = (ls(fg, r"x\d") |> sortDFG)[end-120:end]
rmlist = setdiff(ls(fg), klist)
rmfct = union((x->ls(fg, x)).(rmlist)...)
(x->deleteFactor!(fg, x)).(rmfct)
(x->deleteVariable!(fg, x)).(rmlist)

##

DRT = true
PPE = true
REF = false

include(joinpath(@__DIR__, "MakiePlotImgNoDock.jl"))

scene
##

@show scenepath = joinLogPath(fg, "DRTafter.png")
Base.rm(scenepath, force=true)
Makie.save(scenepath, scene)

##


text:
module OCV
    using Cxx
    using FixedSizeArrays
    using GeometryTypes
    import Distances: euclidean
    import Base: convert

    include("../deps/deps.jl")
    const cv_headerpath = joinpath(dirname(opencv_core),"../include")

    function __init__()
        addHeaderDir(cv_headerpath, kind = C_System)
        Libdl.dlopen(opencv_core,Libdl.RTLD_GLOBAL)
        Libdl.dlopen(opencv_highgui,Libdl.RTLD_GLOBAL)
        Libdl.dlopen(opencv_imgproc,Libdl.RTLD_GLOBAL)
        Libdl.dlopen(opencv_objdetect,Libdl.RTLD_GLOBAL)
        Libdl.dlopen(opencv_text,Libdl.RTLD_GLOBAL)
    end

    __init__()
    cxx"""
    #include <opencv2/highgui/highgui.hpp>
    #include <opencv2/imgproc/imgproc.hpp>
    #include <opencv2/core/core.hpp>
    """
    @Cxx.list cxxt"std::vector<cv::Vec4i>"

    function toline(line::Union{cxxt"cv::Vec4i",cxxt"cv::Vec4i&"})
        (Point{2,Int32}(icxx"$line[0];",icxx"$line[1];"),
         Point{2,Int32}(icxx"$line[2];",icxx"$line[3];"))
    end

    function tocircle(circle::Union{cxxt"cv::Vec3f",cxxt"cv::Vec3f&"})
        Circle{Float32}(Point{2,Float32}(icxx"$circle[0];",icxx"$circle[1];"),icxx"$circle[2];")
    end

    euclidean(p1::Point{2},p2::Point{2}) = sqrt((p1[1]-p2[1])^2+(p1[2]-p2[2])^2)

    Base.convert{T}(::Type{Point{2,T}},point::cxxt"cv::Point_<$T>") =
        Point{2,T}(icxx"$point.x;",icxx"$point.y;")
    function Base.convert{T}(::Type{HyperRectangle{2,T}},
            rect::Union{cxxt"cv::Rect_<$T>",cxxt"cv::Rect_<$T>&"})
        HyperRectangle{2,T}(icxx"$rect.tl();",
                          icxx"$rect.br();")
    end
    Base.convert{T}(::Type{HyperRectangle},rect::Union{cxxt"cv::Rect_<$T>",cxxt"cv::Rect_<$T>&"}) =
        convert(HyperRectangle{2,T},rect)

    using Images, ColorTypes

    function Base.convert(::Type{cxxt"cv::Mat"}, img::Image)
       src = icxx"""
         cv::Mat{(int)$(size(img)[2]),
                 (int)$(size(img)[1]),CV_8UC3,
                 $(convert(Ptr{Void},pointer(data(convert(Image{BGR},img)))))};
       """
   end

end # module


text:
function process(
    zipfilein::String, h5fileout::String;
    compressionlevel=1, strlen=128,
    timestampformat::DateFormat=DateFormat("d/m/y H:M:S"),
    sample::String="Mean")

    systemdata, resultvalues = open_plexoszip(zipfilein)

    targetsample = findsample(systemdata, sample)

    h5open(h5fileout, "w") do h5file::HDF5.File
        addconfigs!(h5file, systemdata)
        membership_idxs = addcollections!(h5file, systemdata, strlen, compressionlevel)
        addtimes!(h5file, systemdata, timestampformat, compressionlevel)
        length(systemdata.intervals) > 0 &&
            addblocks!(h5file, systemdata, timestampformat, compressionlevel)
        addvalues!(h5file, systemdata, membership_idxs, resultvalues,
                   compressionlevel, targetsample)
    end

end

function addconfigs!(f::HDF5.File, data::PLEXOSSolutionDataset)

    rootattrs = attributes(f)
    rootattrs["h5plexos"] = H5PLEXOS_VERSION

    for config in data.configs
        rootattrs[config.element] = config.value
    end

end

function addcollections!(
    f::HDF5.File, data::PLEXOSSolutionDataset,
    strlen::Int, compressionlevel::Int)

    counts = Dict{PLEXOSCollection,Int}()

    for membership in data.memberships
        if membership.collection in keys(counts)
            counts[membership.collection] += 1
        else
            counts[membership.collection] = 1
        end
    end

    collections = Dict(
         collection => (Vector{Tuple{String,String}}(undef, counts[collection]), 0)
         for collection in keys(counts))

    membership_idxs = Dict{PLEXOSMembership, Int}()

    for membership in data.memberships

        collection = membership.collection
        collection_memberships, collection_idx = collections[collection]
        collection_idx += 1

        membership_idxs[membership] = collection_idx

        collection_memberships[collection_idx] =
            isobjects(collection) ?
                (membership.childobject.name, # object collection
                 membership.childobject.category.name) :
                (membership.parentobject.name, # relation collection
                 membership.childobject.name)

        collections[collection] = (collection_memberships, collection_idx)

    end

    h5meta = create_group(f, "metadata")
    h5objects = create_group(h5meta, "objects")
    h5relations = create_group(h5meta, "relations")

    for collection in keys(collections)

        collection_memberships, _ = collections[collection]
        name = h5plexosname(collection)

        group, colnames = isobjects(collection) ?
            (h5objects, ("name", "category")) :
            (h5relations, ("parent", "child"))

        string_table!(group, name, strlen, colnames, collection_memberships)

    end

    return membership_idxs

end

function addvalues!(
    f::HDF5.File, data::PLEXOSSolutionDataset,
    membership_idxs::Dict{PLEXOSMembership,Int},
    resultvalues::Dict{Int,Vector{UInt8}},
    compressionlevel::Int, target_sample::PLEXOSSample)

    propertybands = Dict{PLEXOSProperty,Int}()

    # Advance pass through keys to determine number of bands per property
    for key in data.keys
        if key.property in keys(propertybands)
            propertybands[key.property] =
                max(propertybands[key.property], key.band)
        else
            propertybands[key.property] = key.band
        end
    end

    h5data = create_group(f, "data")

    for ki in data.keyindices

        ki.key.sample === target_sample || continue

        dset = dataset!(h5data, ki, propertybands, compressionlevel)
        member_idx = membership_idxs[ki.key.membership]

        start_idx = ki.position + 1
        end_idx = ki.position + 8*ki.length
        rawvalues = view(resultvalues[ki.periodtype], start_idx:end_idx)
        values = reinterpret(Float64, rawvalues)

        dset[ki.key.band, :, member_idx] = values

    end

end

function dataset!(h5data::HDF5.Group, ki::PLEXOSKeyIndex,
                  propertybands::Dict{PLEXOSProperty,Int},
                  compressionlevel::Int)

    collection = ki.key.membership.collection
    property = ki.key.property

    phase = phasenames[ki.key.phase]
    period = periodnames[ki.periodtype]
    coll = h5plexosname(collection)
    summarydata = property.issummary && (ki.periodtype != 0)
    prop = summarydata ? property.summaryname : property.name

    h5phase = haskey(h5data, phase) ? h5data[phase] : create_group(h5data, phase)
    h5period = haskey(h5phase, period) ? h5phase[period] : create_group(h5phase, period)
    h5coll = haskey(h5period, coll) ? h5period[coll] : create_group(h5period, coll)

    if haskey(h5coll, prop)

        dset = h5coll[prop]

    else

        nbands = propertybands[property]
        ntimes = ki.length

        collectiontype = isobjects(collection) ? "objects" : "relations"
        members = HDF5.root(h5data)["metadata/" * collectiontype * "/" * coll]
        nmembers = length(members)

        dset = create_dataset(h5coll, prop, HDF5.datatype(Float64),
                        HDF5.dataspace(nbands, ntimes, nmembers),
                        chunk=(nbands, ntimes, 1),
                        compress=compressionlevel)

        dset_attrs = attributes(dset)
        dset_attrs["period_offset"] = ki.periodoffset
        dset_attrs["units"] =
            summarydata ? property.summaryunit.value : property.unit.value

    end

    return dset

end

function addtimes!(f::HDF5.File, data::PLEXOSSolutionDataset,
                   localformat::DateFormat, compressionlevel::Int)

    # PLEXOS data format notes:
    # Period type 0 - phase-native interval / block data
    # Maps to ST periodtype 0 via relevant phase table
    # Store both block and interval results on disk (if not ST)?

    # Period type 1-7 - period-type-specific data
    # Direct mapping to period labels

    h5times = create_group(f["metadata"], "times")

    stringtype_id = HDF5.h5t_copy(HDF5.hdf5_type_id(String))
    HDF5.h5t_set_size(stringtype_id, 19)
    datestringtype = HDF5.Datatype(stringtype_id)

    for (dfield, pfield) in periodtypes

        uselocalformat = dfield == :intervals || dfield == :hours
        periodset = getfield(data, dfield)
        nperiods = length(periodset)
        nperiods == 0 && continue

        period_dts = DateTime.(getfield.(periodset, pfield),
                               uselocalformat ? localformat : stdformat)
        issorted(period_dts) || error("$(string(dfield)) not sorted")

        dset = create_dataset(h5times, periodname(dfield), datestringtype,
                              HDF5.dataspace((nperiods,)))

        timestamps = ascii.(format.(period_dts, stdformat))
        any(x -> length(x) != 19, timestamps) && error("Unexpected timestamp format")

        charlists = Vector{Char}.(timestamps)
        rawdata = UInt8.(vcat(charlists...))
        HDF5.h5d_write(
            dset, datestringtype, HDF5.H5S_ALL, HDF5.H5S_ALL,
            HDF5.H5P_DEFAULT, rawdata)

    end

end

function addblocks!(f::HDF5.File, data::PLEXOSSolutionDataset,
                    localformat::DateFormat, compressionlevel::Int)

    h5blocks = create_group(f["metadata"], "blocks")
    colnames = ("interval", "block")

    for phase in phasetypes

        blockmaps = getfield(data, phase)
        n_blockmaps = length(blockmaps)
        n_blockmaps == 0 && continue

        period_dts = DateTime.(
            getfield.(getfield.(blockmaps, :interval), :datetime), localformat)
        periods = format.(period_dts, stdformat)
        blocks = UInt32.(getproperty.(blockmaps, :period))

        string_uint32_table!(h5blocks, phasename(phase), 19,
                             colnames, collect(zip(periods, blocks)))

    end

end


text:
#!/usr/bin/env julia
#
# Jerzy Wroczyński (nr. indeksu 250075)
# Obliczenia Naukowe, Laboratorium
# Lista 1, Zadanie 2.
#

# Kahan’s ϵ finding method
function kahans_eps(fl)
    return abs(fl( 3 * fl(fl(4/3) - 1) - 1 ))
end

# print out the results for different float arithmetic types
for fl in [Float16, Float32, Float64]
    println(fl)
    println("  - Kahan’s method = ", kahans_eps(fl), "; actual ϵ = ", eps(fl), "; equal = ", (kahans_eps(fl) == eps(fl)))
end

max_stars_repo_name: jerry-sky/academic-notebook

max_issues_repo_name: jerry-sky/academic-notebook

max_forks_repo_name: jerry-sky/academic-notebook


text:
using Test, Random, LinearAlgebra, StatsBase, StatsPlots, Distributions
using Traceur, Parameters, Logging, Dates, GLM

using .Types, .Func, .Params

function runmodel(params=Params.default(); fundselector="probabilistic",
        doplot=false, boundstest=false)
    @unpack bigk, bigm, bign, bigt, impactrange, perfwindow, plotpath = params
    println("Run for ", bigt, " periods.")
    println("Sizes: ", bign, " investors, ", bigm, " funds, ", bigk, " stocks.")
    println("Stock impact range:", impactrange, ".")
    # Set up types
    market = Types.MarketIndex(
        zeros(bigt))
    stocks = Types.Equity(
        zeros(bigm, bigt),
        zeros(bigm),
        zeros(bigm),
        zeros(bigm),
        zeros(bigm, bigt))
    investors = Types.RetailInvestor(
        zeros(bign, bigk + 1),
        zeros(bign),
        zeros(bign))
    funds = Func.Types.EquityFund(
        zeros(bigk, bigm),
        zeros(bigk, bign),
        zeros(bigk, bigt))

    Func.initialise(market, stocks, investors, funds, params)

    if boundstest
        Func.boundstest(market, stocks, investors, funds)
    end

    liquidationlog = Func.modelrun(market, stocks, investors, funds, params,
        fundselector)
    println(liquidationlog)

    if boundstest
        Func.boundstest(market, stocks, investors, funds)
    end

    logit = Func.logisticregression(liquidationlog, params)
    println("\nLogistic regression results with whitened indep vars\n", logit)

    stylefacts = Func.calc_stylisedfacts(market.value, stocks.value,
        stocks.volume, params)

    if doplot
        Func.plot_stylisedfacts(market.value, stocks.value, stylefacts, params,
            liquidationlog)
        io = open(joinpath(plotpath, "parameters.txt"), "w+")
        logger = SimpleLogger(io)
        with_logger(logger) do
            @info "Parameters" Dates.today() Dates.Time(Dates.now()) fundselector params
            @info "zScoreKurtoses" stylefacts[10]
            @info "zScoreLossGain" stylefacts[11]
            @info "zScoreVolumeVola" stylefacts[12]
            @info "resultKStest" stylefacts[13]
            @info "logit1" logit[1]
            @info "logit2" logit[2]
            @info "logit3" logit[3]
            @info "liquidationslog" liquidationlog
        end
        close(io)
    end
end # runmodel()

# TODO: Think about how I should have tested the initialisation
# TODO: Think about how I should have tested the running cycle
# TODO: Tests for cases with empty reviewers, divestments


text:
mutable struct ADMM{U} <: AbstractOptimizer
    λ::U
    ρ::U
end

"""
    ADMM()
    ADMM(λ, ρ)

`ADMM` is an implementation of Lasso using the alternating direction methods of multipliers and loosely based on [this implementation](https://web.stanford.edu/~boyd/papers/admm/lasso/lasso.html).

`λ` is the sparsification parameter, `ρ` the augmented Lagrangian parameter.

# Example
```julia
opt = ADMM()
opt = ADMM(1e-1, 2.0)
```
"""
ADMM() = ADMM(0.1, 1.0)

function set_threshold!(opt::ADMM, threshold)
    opt.λ = threshold*opt.ρ
end

get_threshold(opt::ADMM) = opt.λ/opt.ρ

init(o::ADMM, A::AbstractArray, Y::AbstractArray) =  A \ Y
init!(X::AbstractArray, o::ADMM, A::AbstractArray, Y::AbstractArray) =  ldiv!(X, qr(A, Val(true)), Y)

#soft_thresholding(x::AbstractArray, t::T) where T <: Real = sign.(x) .* max.(abs.(x) .- t, zero(eltype(x)))

function fit!(X::AbstractArray, A::AbstractArray, Y::AbstractArray, opt::ADMM; maxiter::Int64 = 1, convergence_error::T = eps()) where T <: Real
    n, m = size(A)

    g = NormL1(get_threshold(opt))

    x̂ = deepcopy(X)
    ŷ = zero(X)

    P = I(m)/opt.ρ - (A' * pinv(opt.ρ*I(n) + A*A') *A)/opt.ρ
    c = P*(A'*Y)

    x_i = similar(X)
    x_i .= X

    iters = 0

    @inbounds for i in 1:maxiter
        iters += 1

        x̂ .= P*(opt.ρ*X - ŷ) + c
        prox!(X, g, x̂ + ŷ/opt.ρ)
        ŷ .= ŷ + opt.ρ*(x̂ - X)

        if norm(x_i - X, 2) < convergence_error
            break
        else
            x_i .= X
        end

    end

    X[abs.(X) .< get_threshold(opt)] .= zero(eltype(X))
    return iters
end


text:
"""
    Cnk(n, k, σ)

Calculation of the Cnk using the direct formula with arbitrary precision. This is unstable for large values of k, more so for small values of σ.

# Examples
```julia-repl
julia> GibbsTypePriors.Cnk(6, 5, 0.5)
[0.23437500000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000 +/- 6.66e-1506]
```
"""
function Cnk(n, k, σ::Float64)
    # = Cnk(n, k, σ, 0)
    σ_arb = RR(σ)
    return Cnk(n, k, σ_arb)
end
function Cnk(n, k, σ::arb)
    return 1//fac(RR(k)) * sum([ (-1)^i * binom(k, i, RR) * risingfac(-i * σ,n) for i in 0:k])
end

"""
    noncentral_generalised_factorial_coefficient(n, k, s, r)

Calculation of the noncentral generalised factorial coefficient (Charalambides, 2005) using a recursive relation. This function uses memoization to speed up the recursion.
It is called by Cnk_rec but follows the sign convention of (Charalambides, 2005), which differs from that of Cnk.

References:
- Charalambides, Charalambos A. Combinatorial methods in discrete distributions. Vol. 600. John Wiley & Sons, 2005.

# Examples
```julia-repl
julia> GibbsTypePriors.noncentral_generalised_factorial_coefficient(10, 5, 0.2, 0.)
-40.19355648000003
```
"""
@memoize function noncentral_generalised_factorial_coefficient(n, k, s::T, r)::T where T<:Number
    @assert n >= 0
    @assert k >= 0
    if k==0
        if n==0
            return 1
        else
            return risingfac(r, n)
        end
    else
        if k>n
            return 0
        else
            return (s * k + r - n + 1) * noncentral_generalised_factorial_coefficient(n - 1, k, s, r) + s * noncentral_generalised_factorial_coefficient(n - 1, k - 1, s, r)
        end
    end
end

# @memoize function noncentral_generalised_factorial_coefficient(n, k, s::arb, r::arb)::arb
#     #Getting ERROR: MethodError: Cannot `convert` an object of type Int64 to an object of type arb, which I am not able to fix. Indeed, there is no default convert method, the bits of precision must be chosen. This is a specialised function to fix this error
#     @assert n >= 0
#     @assert k >= 0
#     if k==0
#         if n==0
#             return arb_1
#         else
#             return risingfac(r, n)
#         end
#     else
#         if k>n
#             return arb_0
#         else
#             return (s * k + r - n + arb_1) * noncentral_generalised_factorial_coefficient(n - 1, k, s, r) + s * noncentral_generalised_factorial_coefficient(n - 1, k - 1, s, r)
#         end
#     end
# end

@memoize function noncentral_generalised_factorial_coefficient(n, k, s::arb, r::arb; arb_0::arb=arb_0, arb_1::arb=arb_1)::arb
    #Getting ERROR: MethodError: Cannot `convert` an object of type Int64 to an object of type arb, which I am not able to fix. This is a specialised function to fix this error
    @assert n >= 0
    @assert k >= 0
    if k==0
        if n==0
            return arb_1
        else
            return risingfac(r, n)
        end
    else
        if k>n
            return arb_0
        else
            return (s * k + r - n + 1) * noncentral_generalised_factorial_coefficient(n - 1, k, s, r; arb_0=arb_0, arb_1=arb_1) + s * noncentral_generalised_factorial_coefficient(n - 1, k - 1, s, r; arb_0=arb_0, arb_1=arb_1)
        end
    end
end

"""
    Cnk_rec(n, k, σ)

Calculation of the Cnk using the recursive formula with arbitrary precision. This seems stable for large values of k, whatever the range of σ, but slower than Cnk(n, k, σ).

# Examples
```julia-repl
julia> GibbsTypePriors.Cnk_rec(6, 5, 0.5)
0.234375
```
"""
function Cnk_rec(n, k, σ)
    # factor_k = factorial(k)
#   sum((-1)^(1:k) * gmp::chooseZ(n = k, k = 1:k) * Rmpfr::pochMpfr(-(1:k) * Gama, n) / factor_k)
  (-1)^(n - k) * noncentral_generalised_factorial_coefficient(n, k, σ, 0)
end
function Cnk_rec(n, k, σ::arb)
    #Same comment as  noncentral_generalised_factorial_coefficient(n, k, s::arb, r::arb)::arb

  (-1)^(n - k) * noncentral_generalised_factorial_coefficient(n, k, σ, arb_0)
end

"""
    Cnk_robust(n, k, σ; verbose = false)

Calculation of the Cnk using either the direct formula with arbitrary precision or the recursive formula. Accuracy is checked for the fast direct computation, and if it is too low (roughly lower than double precision), the function switches to recursive computation of the Cnk.

# Examples
```julia-repl
julia> GibbsTypePriors.Cnk_robust(6, 5, 0.5)
[0.23437500000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000 +/- 6.66e-1506]
```
"""
function Cnk_robust(n, k, σ; verbose = false)
    res = Cnk(n, k, σ)
    if has_reasonable_precision(res)
        return res
    else
        rec_res = Cnk_rec(n, k, RR(σ))
        if verbose
            println("switching to recursive formula for Cnk")
        end
        if has_reasonable_precision(rec_res)
            return rec_res
        else
            error("Unable to compute Cnk with reasonable precision (see function has_reasonable_precision() )")
        end
    end
end

max_stars_repo_name: konkam/GibbsTypePriors

max_issues_repo_name: konkam/GibbsTypePriors

max_forks_repo_name: konkam/GibbsTypePriors