text:
"""
    inputvertex(name, size)

Return an immutable input type vertex with the given `name` and `size`.

Typically used as "entry" point to a computation graph.

# Examples
```julia-repl
julia> using NaiveNASlib

julia> inputvertex("input", 5)
InputSizeVertex(InputVertex("input"), outputs=[], 5)

```
"""
inputvertex(name, size) = InputSizeVertex(name, size)

"""
    vertex(computation, outsize::Integer, trait::MutationTrait, inputs::AbstractVertex...; mutation=IoChange)

Return a mutable computation type vertex.

# Examples
```julia-repl
julia> using NaiveNASlib

julia> v = vertex(x -> 5x, 1, SizeInvariant(), inputvertex("input", 1));

julia> v(3)
15
```
"""
vertex(computation, outsize::Integer, trait::MutationTrait, inputs::AbstractVertex...; mutation=IoChange) = MutationVertex(CompVertex(computation, inputs...), mutation(collect(nout.(inputs)), outsize), trait)

"""
    immutablevertex(computation, outsize::Integer, inputs::AbstractVertex...; mutation=IoChange, traitdecoration=identity)

Return an immutable computation type vertex.

# Examples
```julia-repl
julia> using NaiveNASlib

julia> v = immutablevertex(x -> x * [1 2 3; 4 5 6], 3, inputvertex("input", 2));

julia> nin(v)
1-element Array{Int64,1}:
 2

julia> nout(v)
3

julia> v([1 2])
1×3 Array{Int64,2}:
 9  12  15

```
"""
immutablevertex(computation, outsize::Integer, inputs::AbstractVertex...; mutation=IoChange, traitdecoration=identity) = vertex(computation, outsize, traitdecoration(Immutable()), inputs..., mutation=mutation)

"""
    absorbvertex(computation, outsize::Integer, inputs::AbstractVertex...; mutation=IoChange, traitdecoration=identity)

Return a mutable computation type vertex which absorbs size changes. Typical example of this is a neural network layer.

# Examples
```julia-repl
julia> using NaiveNASlib

julia> v = absorbvertex(x -> x * [1 2 3; 4 5 6], 3, inputvertex("input", 2));

julia> nin(v)
1-element Array{Int64,1}:
 2

julia> nout(v)
3

julia> v([1 2])
1×3 Array{Int64,2}:
 9  12  15

```
"""
absorbvertex(computation, outsize::Integer, inputs::AbstractVertex...; mutation=IoChange, traitdecoration=identity) = vertex(computation, outsize, traitdecoration(SizeAbsorb()), inputs..., mutation=mutation)

"""
    invariantvertex(computation, input; mutation=IoChange, traitdecoration=identity)

Return a mutable computation type vertex which is size invariant, i.e nin == nout.

# Examples
```julia-repl
julia> using NaiveNASlib

julia> v = invariantvertex(x -> 2 .* x, inputvertex("input", 2));

julia> nin(v)
1-element Array{Int64,1}:
 2

julia> nout(v)
2

julia> v([1 2])
1×2 Array{Int64,2}:
 2  4
```
"""
invariantvertex(computation, input; mutation=IoChange, traitdecoration=identity) = vertex(computation, nout(input), traitdecoration(SizeInvariant()), input, mutation=mutation)

"""
    conc(v::AbstractVertex, vs::AbstractVertex...; dims, mutation=IoChange, traitdecoration=identity)

Return a mutable vertex which concatenates input along dimension `dim`.
# Examples
```julia-repl
julia> using NaiveNASlib

julia> v = conc(inputvertex.(["in1", "in2", "in3"], 1:3)..., dims=1);

julia> nin(v)
3-element Array{Int64,1}:
 1
 2
 3

julia> nout(v)
6

julia> v([1], [2, 3], [4, 5, 6])
6-element Array{Int64,1}:
 1
 2
 3
 4
 5
 6
```
"""
conc(v::AbstractVertex, vs::AbstractVertex...; dims, mutation=IoChange, traitdecoration=identity, outwrap=identity) = vertex(outwrap((x...) -> cat(x..., dims=dims)), sum(nout.((v, vs...))), traitdecoration(SizeStack()), v, vs..., mutation=mutation)

"""
    VertexConf
    VertexConf(;mutation = IoChange, traitdecoration = identity, outwrap = identity)

Config struct to be used with element wise op syntax (`+`, `-`, `*`, `/`).

`mutation` determines how size changes of the vertex are handled
`traitdecoration` allows for decorating the vertex trait with stuff like logging, validation etc.
`outwrap` is a function which returns a function which will be applied to the computed output. For example, the following `outwrap` scales output by a factor of 2: `outwrap = f ->  (x...) -> 2f((x...)`
"""
struct VertexConf
    mutation
    traitdecoration
    outwrap
end
mutationconf(m) = VertexConf(mutation=m)
traitconf(t) = VertexConf(traitdecoration=t)
outwrapconf(o) = VertexConf(outwrap=o)
VertexConf(;mutation = IoChange, traitdecoration = identity, outwrap = identity)= VertexConf(mutation, traitdecoration, outwrap)

function invariant_outsize(vs::AbstractVertex...)
    outsize = unique([nout.(vs)...])
    length(outsize) == 1 || throw(DimensionMismatch("Dimensions must match! Got $([nout.(vs)...])"))
    return outsize[]
end

# Common wiring for all elementwise operations
elemwise(op, conf::VertexConf, vs::AbstractVertex...) = vertex(conf.outwrap((x...) -> op.(x...)), invariant_outsize(vs...), conf.traitdecoration(SizeInvariant()), vs..., mutation=conf.mutation)

"""
    >>(conf::VertexConf, v::AbstractVertex)
    >>(name::String, v::AbstractVertex)
    >>(outwrap::Function, v::AbstractVertex)

Return inputs as a tuple. Only used to enable the `conf >> v1 op v2 ...` syntax.

"""
Base.:>>(conf::VertexConf, v::AbstractVertex) = (conf, v)
Base.:>>(name::String, v::AbstractVertex) = traitconf(t -> NamedTrait(t, name)) >> v
Base.:>>(outwrap::Function, v::AbstractVertex) = outwrapconf(outwrap) >> v

"""
    +(v::AbstractVertex, vs::AbstractVertex...)

Return a mutable vertex which performs (broadcasted) elementwise addition of its inputs.

A `VertexConf` with functions to create `MutationState` and `MutationTrait` can be supplied through the `>>` operator.

# Examples

```julia-repl
julia> using NaiveNASlib

julia> v = inputvertex("in1", 2) + inputvertex("in2", 2) + inputvertex("in3" ,2);

julia> nin(v)
3-element Array{Int64,1}:
 2
 2
 2

julia> nout(v)
2

julia> v([1, 2], [3, 4], [5, 6])
2-element Array{Int64,1}:
  9
 12

julia> name(v)
"MutationVertex::SizeInvariant"

julia> typeof(op(v))
IoChange

julia> conf = VertexConf(IoSize, t -> NamedTrait(t, "v"));

julia> v = conf >> inputvertex("in1", 3) + inputvertex("in2", 3);

julia> name(v)
"v"

julia> typeof(op(v))
IoSize
```
"""
Base.:+((conf, v)::Tuple{VertexConf, <:AbstractVertex}, vs::AbstractVertex...) = elemwise(+, conf, v, vs...)
Base.:+(v::AbstractVertex, vs::AbstractVertex...) = +(VertexConf() >> v, vs...)

"""
    *(v::AbstractVertex, vs::AbstractVertex...)
    *(conf::VertexConf, v::AbstractVertex, vs::AbstractVertex...)

Return a mutable vertex which performs (broadcasted) elementwise multiplication of its inputs.

A `VertexConf` with functions to create `MutationState` and `MutationTrait` can be supplied through the `>>` operator.

# Examples

```julia-repl
julia> using NaiveNASlib

julia> v = inputvertex("in1", 2) * inputvertex("in2", 2) * inputvertex("in3" ,2);

julia> nin(v)
3-element Array{Int64,1}:
 2
 2
 2

julia> nout(v)
2

julia> v([1, 2], [3, 4], [5, 6])
2-element Array{Int64,1}:
 15
 48

julia> name(v)
"MutationVertex::SizeInvariant"

julia> typeof(op(v))
IoChange

julia> conf = VertexConf(IoSize, t -> NamedTrait(t, "v"));

julia> v = conf >> inputvertex("in1", 3) * inputvertex("in2", 3);

julia> name(v)
"v"

julia> typeof(op(v))
IoSize

```
"""
Base.:*((conf, v)::Tuple{VertexConf, <:AbstractVertex}, vs::AbstractVertex...) = elemwise(*, conf, v, vs...)
Base.:*(v::AbstractVertex, vs::AbstractVertex...) = *(VertexConf() >> v, vs...)

"""
    -(v1::AbstractVertex, v2::AbstractVertex)

Return a mutable vertex which performs (broadcasted) elementwise subtraction of its inputs.

A `VertexConf` with functions to create `MutationState` and `MutationTrait` can be supplied through the `>>` operator.

# Examples

```julia-repl
julia> using NaiveNASlib

julia> v = inputvertex("in1", 2) - inputvertex("in2", 2)

julia> nin(v)
2-element Array{Int64,1}:
 2
 2

julia> nout(v)
2

julia> v([1, 2], [3, 4])
2-element Array{Int64,1}:
 -2
 -2

julia> name(v)
"MutationVertex::SizeInvariant"

julia> typeof(op(v))
IoChange

julia> conf = VertexConf(IoSize, t -> NamedTrait(t, "v"));

julia> v = conf >> inputvertex("in1", 3) - inputvertex("in2", 3);

julia> name(v)
"v"

julia> typeof(op(v))
IoSize
```
"""
Base.:-((conf, v1)::Tuple{VertexConf, <:AbstractVertex}, v2::AbstractVertex) = elemwise(-, conf, v1, v2)
Base.:-(v1::AbstractVertex, v2::AbstractVertex) = -(VertexConf() >> v1, v2)

"""
    -(v::AbstractVertex)

Return a mutable vertex which performs elementwise negation of its input.

A `VertexConf` with functions to create `MutationState` and `MutationTrait` can be supplied through the `>>` operator. Due to operator precedence, this has to be done in the following order: `-(conf >> v)`

#Examples

```julia-repl
julia> using NaiveNASlib

julia> v = -inputvertex("in", 2);

julia> v([1,2])
2-element Array{Int64,1}:
 -1
 -2

julia> conf = VertexConf(IoSize, t -> NamedTrait(t, "v"));

julia> v = -(conf >> inputvertex("in", 2));

julia> name(v)
"v"

julia> typeof(op(v))
IoSize
```
"""
Base.:-((conf, v)::Tuple{VertexConf, <:AbstractVertex}) = elemwise(-, conf, v)
Base.:-(v::AbstractVertex) = -(VertexConf() >> v)

"""
    /(v1::AbstractVertex, v2::AbstractVertex)
    /(conf::VertexConf, v1::AbstractVertex, v2::AbstractVertex)

Return a mutable vertex which performs (broadcasted) elementwise division of its inputs.

A `VertexConf` with functions to create `MutationState` and `MutationTrait` can be supplied through the `>>` operator.

# Examples

```julia-repl
julia> using NaiveNASlib

julia> v = inputvertex("in1", 2) / inputvertex("in2", 2)

julia> nin(v)
2/element Array{Int64,1}:
 2
 2

julia> nout(v)
2

julia> v([6, 8], [2, 4])
2/element Array{Int64,1}:
 3
 2

julia> name(v)
"MutationVertex::SizeInvariant"

julia> typeof(op(v))
IoChange

julia> conf = VertexConf(IoSize, t -> NamedTrait(t, "v"));

julia> v = conf >> inputvertex("in1", 3) / inputvertex("in2", 3);

julia> name(v)
"v"

julia> typeof(op(v))
IoSize
```
"""
Base.:/((conf, v1)::Tuple{VertexConf, <:AbstractVertex}, v2::AbstractVertex) = elemwise(/, conf, v1, v2)
Base.:/(v1::AbstractVertex, v2::AbstractVertex) = /(VertexConf() >> v1, v2)


text:
mutable struct kNNRecommender
    # hyper parameters
    topK::Int
    shrink::Float64
    weighting::Union{Nothing, Symbol}
    # parameters
    rating_matrix::SparseMatrixCSC
    similarity_matrix::SparseMatrixCSC

    kNNRecommender(topK, shrink) = new(topK, shrink)
    kNNRecommender(topK, shrink, weighting) = new(topK, shrink, weighting)
end

function tfidf(rating_matrix::SparseMatrixCSC)
    U, I, R = findnz(rating_matrix)
    n_users, n_items = size(rating_matrix)

    bincount = zeros(n_users) # to avoid zero devision
    for u in U
        bincount[u] += 1
    end

    idf = log.(n_items ./ (bincount .+ 1e-6)) .+ 1

    for j in 1:length(U)
        R[j] = R[j] * idf[U[j]]
    end

    return sparse(U, I, R)
end

function compute_similarity_matrix(rating_matrix::SparseMatrixCSC, topK::Int, shrink::Float64)
    # (user, item)^T * (user, item) -> (item, item)
    # Return S[i, j] where j is full items, and i is related items at topK

    simJ = Int[]
    simI = Int[]
    simS = Float64[]

    U, I, R = findnz(rating_matrix)
    n_users, n_items = size(rating_matrix)

    norms = sqrt.(sum(rating_matrix.^2, dims=1))

    for j in 1:n_items
        Uj, Rj = findnz(rating_matrix[:, j])
        simj = zeros(n_items)
        for (u, ruj) in zip(Uj, Rj)
            Iu, Ri = findnz(rating_matrix[u, :])
            for (i, rui) in zip(Iu, Ri)
                s = rui * ruj
                s /= norms[j] * norms[i] + shrink + 1e-6
                simj[i] += s
            end
        end

        arg_sort_i = sortperm(simj, rev=true)[2:topK+1]
        append!(simI, arg_sort_i)
        append!(simS, simj[arg_sort_i])
        append!(simJ, fill(j, length(arg_sort_i)))
    end

    return sparse(simI, simJ, simS)
end

function compute_similarity_matrix(knn::kNNRecommender)
    rating_matrix = knn.rating_matrix

    if knn.weighting == :tfidf
        rating_matrix = tfidf(rating_matrix)
    end

    knn.similarity_matrix = compute_similarity_matrix(rating_matrix, knn.topK, knn.shrink)

    return nothing
end

function predict(knn::kNNRecommender, uidx::Int, topN::Int)
    preds = sortperm((knn.rating_matrix[uidx, :]' * knn.similarity_matrix)', rev=true)

    I, R  = findnz(knn.rating_matrix[uidx, :])
    filter!(p->!(p in I), preds)

    return preds[1:min(topN, length(preds))]
end

function predict(knn::kNNRecommender, uidxs::Array{Int64,1}, topN::Int)
    return map(uidx->predict(knn, uidx, topN), uidxs)
end

function evaluate_recall(knn::kNNRecommender, topN::Int, rating_matrix_test::SparseMatrixCSC)
    recalls = Float64[]

    for uidx in 1:size(knn.rating_matrix)[1]
        gts, = findnz(rating_matrix_test[uidx, :])
        if length(gts) == 0 continue end

        preds = predict(knn, uidx, topN)
        tps = intersect(Set(preds), Set(gts))

        push!(recalls, length(tps) / length(gts))
    end

    return recalls
end

ext: jl


text:
using Test

@testset "getsymbols!" begin

    ex = :(a)

    s = Set{Symbol}()

    @test(getsymbols!(ex, s) == ex)
    @test(isempty(s))

    ex = :(tunninginput(a))

    s = Set{Symbol}()

    @test(getsymbols!(ex, s) == :(a))
    @test(ex == :(tunninginput(a)))
    @test(s == Set([:a]))

    ex = :(f(tunninginput(a), b))

    s = Set{Symbol}()

    @test(getsymbols!(ex, s) == :(f(a, b)))
    @test(ex == :(f(tunninginput(a), b)))
    @test(s == Set([:a]))

    ex = :(f(cos(sin(tunninginput(a))), [tunninginput(b), tunninginput(c), 5]))

    s = Set{Symbol}()

    @test(getsymbols!(ex, s) == :(f(cos(sin(a)), [b, c, 5])))
    @test(ex == :(f(cos(sin(tunninginput(a))), [tunninginput(b), tunninginput(c), 5])))
    @test(s == Set([:a, :b, :c]))

end


text:
"""
This struct contains all the data for the transmission lines between zones.
    name: Path to the directory where input data is stored.
    from_zone: Name of originating zone.
    to_zone: Name of destination zone.
    active_power_limit: Maximum flow capacity (MW).
"""
mutable struct ZonalLine
    name::String
    from_zone::String
    to_zone::String
    active_power_limit::Float64
end

get_name(line::ZonalLine) = line.name
get_from_zone(line::ZonalLine) = line.from_zone
get_to_zone(line::ZonalLine) = line.to_zone
get_active_power_limit(line::ZonalLine) = line.active_power_limit

function  set_active_power_limit!(line::ZonalLine, active_power_limit::Float64)
    line.active_power_limit = active_power_limit
    return
end


text:
__precompile__()

module SequentialLOB

using LinearAlgebra
using Statistics
using Distributions
using Random
using SharedArrays
using Distributed
using SpecialFunctions
using Logging
using IOLogging
using ProgressMeter
using TSSM

include("source_function.jl")
include("reaction_diffusion_path.jl")
include("parse_params.jl")
include("reaction_diffusion_spde.jl")
include("objective_surface.jl")

__version__ = "Sequential LOB"

export SLOB,
       SourceTerm,
       parse_commandline,
       ObjectiveSurface

end # module


text:
module Golgi

using Enzyme
using Zygote
using ChainRules
using CompilerPluginTools

@generated function rrule_call(f, args...)
    transform = find_heuristic_best_transform(f, args)
    return quote
        rrule_call($transform, f, args...)
    end
end

abstract type CompilerTransform end

struct NativeTransform <: CompilerTransform end
struct EnzymeTransform <: CompilerTransform end
struct YaoTransform <: CompilerTransform end
struct GenTransform <: CompilerTransform end
struct ChainRulesTransform <: CompilerTransform end

# defines the intrinsic/entry of the eDSL compiler transform
# when fallback to the glue compiler implementation
# it is always true, but of course this is gonna be
# generating slower code, but who cares
is_entry_of(::NativeTransform, f, tt) = true

function is_entry_of(::ChainRulesTransform, f, tt::Tuple)
    rrule_rettype = Core.Compiler.return_type(ChainRules.rrule, Tuple{f, tt...})
    return !(rrule_rettype === Nothing)
end

# TODO: check on Yao's quantum program entry
function is_entry_of(::YaoTransform, f, tt::Tuple)
    return false
end

# function is_entry_of(::GenTransform, f, tt)
# end

function is_entry_of(::EnzymeTransform, f, tt::Tuple)
    return f <: Function
end

function find_heuristic_best_transform(f, tt)
    is_entry_of(ChainRulesTransform(), f, tt) && return ChainRulesTransform()
    # is_entry_of(YaoTransform(), f, tt) && return YaoTransform()
    is_entry_of(EnzymeTransform(), f, tt) && return EnzymeTransform()
    return NativeTransform()
end

function rrule_call(::ChainRulesTransform, f, args...)
    return ChainRules.rrule(f, args...)
end

function rrule_call(::NativeTransform, f, args...)
    Zygote.pullback(f, args...)
end

# wrap enzyme to respect ChainRules convention
# NOTE: we don't handle Julia callable struct
# in Enzyme - until it works in Enzyme natively
function rrule_call(::EnzymeTransform, f::Function, args...)
    f(args...), function enzyme_pullback(Δ)
        grad = Enzyme.pullback(f, args...)(Δ)
        return (ZeroTangent(), grad...)
    end
end

end


text:
"""
    Site

Struct containing the coordinates, parametrized via primitive translations plus basis index and in real space form, of a lattice site.
* `int :: SVector{4, Int64}`   : site given in term of primitive translation and basis index
* `vec :: SVector{3, Float64}` : site given in cartesian coordinates
"""
struct Site
    int :: SVector{4, Int64}
    vec :: SVector{3, Float64}
end

# transform int to vec representation
function get_vec(
    int :: SVector{4, Int64},
    uc  :: Unitcell
    )   :: SVector{3, Float64}

    vec_x = int[1] * uc.vectors[1][1] + int[2] * uc.vectors[2][1] + int[3] * uc.vectors[3][1] + uc.basis[int[4]][1]
    vec_y = int[1] * uc.vectors[1][2] + int[2] * uc.vectors[2][2] + int[3] * uc.vectors[3][2] + uc.basis[int[4]][2]
    vec_z = int[1] * uc.vectors[1][3] + int[2] * uc.vectors[2][3] + int[3] * uc.vectors[3][3] + uc.basis[int[4]][3]
    vec   = SVector{3, Float64}(vec_x, vec_y, vec_z)

    return vec
end

# generate lattice sites from unitcell
function get_sites(
    size      :: Int64,
    uc        :: Unitcell,
    euclidean :: Bool
    )         :: Vector{Site}

    # init buffers
    ints    = SVector{4, Int64}[SVector{4, Int64}(0, 0, 0, 1)]
    current = copy(ints)
    touched = copy(ints)

    # iteratively add sites with bond distance 1 until required Euclidean norm is reached
    if euclidean
        # compute nearest neighbor distance
        nn_distance = norm(get_vec(ints[1] + uc.bonds[1][1], uc))

        while true
            # init list for new sites generated in this step
            # new_ints contains sites inside sphere with radius corresponding to required Euclidean norm
            # out_ints contains sites outside sphere with radius corresponding to required Euclidean norm
            new_ints = SVector{4, Int64}[]
            out_ints = SVector{4, Int64}[]

            # add sites with bond distance 1 to new sites from last step
            for int in current
                for i in eachindex(uc.bonds[int[4]])
                    new_int = int .+ uc.bonds[int[4]][i]

                    # check if site was generated already
                    if in(new_int, touched) == false
                        if in(new_int, ints) == false
                            # check if site is inside sphere with required Euclidean norm
                            if norm(get_vec(new_int, uc)) <= size * nn_distance
                                push!(new_ints, new_int)
                                push!(ints, new_int)
                            else
                                push!(out_ints, new_int)
                            end
                        end
                    end
                end
            end

            # if no site exceed Euclidean norm, update lists and proceed
            if length(out_ints) == 0
                touched = current
                current = new_ints
            # otherwise check if sites have reentrant neighbors
            else
                # init lists for possibles sites with reentrant neighbors
                re_ints = SVector{4, Int64}[]

                # add sites with bond distance 1 to sites outside sphere from last step
                for int in out_ints
                    for i in eachindex(uc.bonds[int[4]])
                        new_int = int .+ uc.bonds[int[4]][i]

                        # check if site was generated already
                        if in(new_int, current) == false
                            if in(new_int, ints) == false
                                # check if site is inside sphere with required Euclidean norm
                                if norm(get_vec(new_int, uc)) <= size * nn_distance
                                    push!(re_ints, int)
                                    break
                                end
                            end
                        end
                    end
                end

                # if sites with reentrant neighbors exist, update lists and proceed
                if length(re_ints) > 0
                    touched = current
                    current = vcat(new_ints, re_ints)
                # otherwise terminate, if there are no more sites to add
                else
                    if length(new_ints) > 0
                        touched = current
                        current = new_ints
                    else
                        break
                    end
                end
            end
        end
    # iteratively add sites with bond distance 1 until required bond distance is reached
    else
        # init metric
        metric = 0

        while metric < size
            # init list for new sites generated in this step
            new_ints = SVector{4, Int64}[]

            # add sites with bond distance 1 to new sites from last step
            for int in current
                for i in eachindex(uc.bonds[int[4]])
                    new_int = int .+ uc.bonds[int[4]][i]

                    # check if site was generated already
                    if in(new_int, touched) == false
                        if in(new_int, ints) == false
                            push!(new_ints, new_int)
                            push!(ints, new_int)
                        end
                    end
                end
            end

            # update lists and increment metric
            touched  = current
            current  = new_ints
            metric  += 1
        end
    end

    # build sites
    sites = Site[Site(int, get_vec(int, uc)) for int in ints]

    return sites
end

"""
    get_metric(
        s1 :: Site,
        s2 :: Site,
        uc :: Unitcell
        )  :: Int64

Compute the bond metric (i.e. minimal number of bonds required to connect two sites within the lattice graph) for a given unitcell.
"""
function get_metric(
    s1 :: Site,
    s2 :: Site,
    uc :: Unitcell
    )  :: Int64

    # init buffers
    current   = SVector{4, Int64}[s2.int]
    touched   = SVector{4, Int64}[s2.int]
    metric    = 0
    not_found = true

    # check if sites are identical
    if s1.int == s2.int
        not_found = false
    end

    # add sites with bond distance 1 around s2 until s1 is reached
    while not_found
        # increment metric
        metric += 1

        # init list for new sites generated in this step
        new_ints = SVector{4, Int64}[]

        # add sites with bond distance 1 to new sites from last step
        for int in current
            for i in eachindex(uc.bonds[int[4]])
                new_int = int .+ uc.bonds[int[4]][i]

                # check if site was touched already
                if in(new_int, touched) == false
                    push!(new_ints, new_int)
                end
            end
        end

        # if s1 is reached stop searching, otherwise update lists and continue
        if in(s1.int, new_ints)
            not_found = false
        else
            touched = current
            current = new_ints
        end
    end

    return metric
end

# check if Float64 item in list within numerical tolerance
function is_in(
    e    :: Float64,
    list :: Vector{Float64}
    )    :: Bool

    in = false

    for item in list
        if abs(item - e) <= 1e-8
            in = true
            break
        end
    end

    return in
end

"""
    get_nbs(
        n    :: Int64,
        s    :: Site,
        list :: Vector{Site}
        )    :: Vector{Int64}

Returns list of n-th nearest neigbors (Euclidean norm) for a given site s, assuming s in list.
"""
function get_nbs(
    n    :: Int64,
    s    :: Site,
    list :: Vector{Site}
    )    :: Vector{Int64}

    # init buffers
    dist        = Float64[]
    dist_unique = Float64[]

    # collect possible distances
    for item in list
        d = norm(item.vec - s.vec)

        if is_in(d, dist_unique) == false
            push!(dist_unique, d)
        end

        push!(dist, d)
    end

    # determine the n-th nearest neighbor distance
    @assert length(dist_unique) >= n + 1 "Could not find neighbors in list."
    dn = sort(dist_unique)[n + 1]

    # collect all sites with distance dn
    nbs = Int64[]

    for i in eachindex(dist)
        if abs(dn - dist[i]) <= 1e-8
            push!(nbs, i)
        end
    end

    return nbs
end

# check if site in list within numerical tolerance
function is_in(
    e    :: Site,
    list :: Vector{Site}
    )    :: Bool

    in = false

    for item in list
        if norm(item.vec - e.vec) <= 1e-8
            in = true
            break
        end
    end

    return in
end

# obtain minimal test set to verify symmetry transformations
function get_test_sites(
    uc        :: Unitcell,
    euclidean :: Bool
    )         :: Tuple{Vector{Site}, Union{Int64, Float64}}

    # init buffers
    test_sites = Site[]

    # add basis sites and connected neighbors to list
    for i in eachindex(uc.basis)
        b = Site(SVector{4, Int64}(0, 0, 0, i), uc.basis[i])

        if is_in(b, test_sites) == false
            push!(test_sites, b)
        end

        for j in eachindex(uc.bonds[i])
            int = b.int .+ uc.bonds[i][j]
            bp  = Site(int, get_vec(int, uc))

            if is_in(bp, test_sites) == false
                push!(test_sites, bp)
            end
        end
    end

    if euclidean
        # determine the maximum Euclidean distance
        metric = maximum(Float64[norm(s.vec) for s in test_sites])

        return test_sites, metric
    else
        # determine the maximum bond distance
        metric = maximum(Int64[get_metric(test_sites[1], s, uc) for s in test_sites])

        return test_sites, metric
    end
end


text:
struct Polygon{N, T <: Real, L <: LineString{N, T}, V <: AbstractVector{L}} <: AbstractRegion{T}
    exterior::L
    interiors::V
end

if VERSION < v"0.7"
    function Polygon(exterior::L, interiors::AbstractVector{L}) where {N, T, L <: LineString{N, T}}
        Polygon{N, T, L, typeof(interiors)}(exterior, interiors)
    end
end
function Polygon{N}(exterior::L, interiors::AbstractVector{L}) where {N, T, L <: LineString{N, T}}
    Polygon{N, T, L, typeof(interiors)}(exterior, interiors)
end
function Polygon{N,T}(exterior::L, interiors::AbstractVector{L}) where {N, T, L <: LineString{N, T}}
    Polygon{N, T, L, typeof(interiors)}(exterior, interiors)
end

Polygon(points::StaticVector{N,T}...) where {N, T<:Real} = Polygon(collect(points))
Polygon(points::AbstractVector{<:StaticVector{N,T}}) where {N,T<:Real} = Polygon{N}(LineString{N}(points))
Polygon{N}(points::AbstractVector{T}...) where {N, T<:Real} = Polygon{N}(collect(points))
Polygon{N}(points::AbstractVector{<:AbstractVector{T}}) where {N, T<:Real} = Polygon{N,T}(LineString{N,T}(points))
Polygon{N,T}(points::AbstractVector{<:Real}...) where {N, T<:Real} = Polygon{N,T}(collect(points))
Polygon{N,T}(points::AbstractVector{<:AbstractVector{<:Real}}) where {N, T<:Real} = Polygon{N,T}(LineString{N,T}(points))

Polygon(ls::LineString) = Polygon(ls, Vector{typeof(ls)}())
Polygon{N}(ls::LineString{N}) where {N} = Polygon{N}(ls, Vector{typeof(ls)}())
Polygon{N,T}(ls::LineString{N,T}) where {N,T<:Real} = Polygon{N,T}(ls, Vector{typeof(ls)}())

convert2d(p::Polygon{2}) = p
convert2d(p::Polygon{3}) = Polygon(convert2d(p.exterior), convert2d.(p.interiors))
convert3d(p::Polygon{2}) = Polygon(convert3d(p.exterior), convert3d.(p.interiors))
convert3d(p::Polygon{2}, z::Real) = Polygon(convert3d(p.exterior, z), convert3d.(p.interiors, z))
convert3d(p::Polygon{3}) = p
convert3d(p::Polygon{3}, z::Real) = Polygon(convert3d(p.exterior, z), convert3d.(p.interiors, z))

# Note: This doesn't compare the ordering of points
function Base.:(==)(p1::Polygon{N}, p2::Polygon{N}) where N
    p1.exterior == p2.exterior && p1.interiors == p2.interiors
end

function Base.isequal(p1::Polygon{N}, p2::Polygon{N}) where N
    isequal(p1.exterior, p2.exterior) && isequal(p1.interiors, p2.interiors)
end

function Base.hash(p::Polygon, h::UInt)
    hash(p.exterior, hash(p.interiors, hash(UInt === UInt64 ? 0xde95b490c51c55a5 : 0x0f7345a4, h)))
end

function show(io::IO, polygon::Polygon{N}) where N
    print(io, "Polygon{$N}([")
    for i in 1:length(polygon.exterior.points)
        print(io, polygon.exterior.points[i])
        if i < length(polygon.exterior.points)
            print(io, ", ")
        end
    end
    print(io, "]")
    for ls in polygon.interiors
    	print(io, ", [")
    	for i in 1:length(ls.points)
	        print(io, ls.points[i])
	        if i < length(ls.points)
	            print(io, ", ")
	        end
	    end
    	print(io, "]")
    end
    print(io, ")")
end

function lines(p::Polygon)
    out = lines(p.exterior)
    for ls in p.interiors
        append!(out, lines(ls))
    end
    return out
end

function area(p::Polygon{2, T}) where {T}
    return area(p.exterior) + mapreduce(area, +, p.interiors; init = zero(T))
end

function in(p::StaticVector{2, <:Real}, polygon::Polygon{2})
    # Winding number algorith, for example read
    # http://geomalgorithms.com/a03-_inclusion.html

    # Count number of times an edge of the polygon crosses a ray from
    # p to infinity (doesn't matter which direction but we choose +x direction)
    wn = winding_number(p, polygon.exterior)
    for linestring in polygon.interiors
        wn += winding_number(p, linestring)
    end

    return wn != 0 # Valid polygons can be oriented clockwise or anticlockwise
end

intersects(p::StaticVector{2, <:Real}, polygon::Polygon{2}) = p ∈ polygon
intersects(polygon::Polygon{2}, p::StaticVector{2, <:Real}) = p ∈ polygon

function intersects(l::Line{2}, polygon::Polygon{2})
    # The polygon is a surface. Either `l` intersects the edge of
    # the polygon or it is entirely *inside* the polygon (or else it
    # doesn't intersect at all)
    if intersects(l, polygon.exterior)
        return true
    end

    for linestring in polygon.interiors
        if intersects(l, linestring)
            return true
        end
    end

    # Either both endpoints are inside, or outside, so just test one
    return l.p1 ∈ polygon
end
intersects(polygon::Polygon{2}, l::Line{2}) = intersects(l, polygon)

function intersects(ls::LineString{2}, polygon::Polygon{2})
    # The polygon is a surface. Either `ls` intersects the edge of
    # the polygon, or it is entirely *inside* the polygon (or else it
    # doesn't intersect at all)
    if intersects(ls, polygon.exterior)
        return true
    end

    for linestring in polygon.interiors
        if intersects(ls, linestring)
            return true
        end
    end

    for p in ls.points
        if p ∉ polygon
            return false
        end
    end

    # `ls` must be entirely inside `polygon`
    return true
end
intersects(polygon::Polygon{2}, ls::LineString{2}) = intersects(ls, polygon)

function intersects(p1::Polygon{2}, p2::Polygon{2})
    # If any line (interior OR exterior) intersects with any other line,
    # then the two polygons intersect. If not, check if one polygon
    # is entirely inside the other

    # Check if edges intersect
    if intersects(p1.exterior, p2)
        return true
    end

    for linestring in p1.interiors
        if intersects(linestring, p2)
            return true
        end
    end

    # One polygon might be entirely inside another
    # Since there can't be multple "onion" rings in any
    # polygon (it can have holes but is otherwise contiguous)
    # we only need to check the exteriors
    if all(p -> p in p2, p1.exterior.points)
        return true
    end
    if all(p -> p in p1, p2.exterior.points)
        return true
    end

    # Or else they do not intersect
    return false
end


text:
"""
    deferredacceptance(students, schools, capacities; verbose, rev)

Given an array of student preferences, where `students[i, j]` indicates the
rank that student `j` gave to school `i`, and an array of the transposed
shape indicating the schools' rankings over the students, uses
student-proposing DA to compute a stable assignment. Returns a vector of
schools corresponding to each student (`m + 1` indicates unassigned) and a
list of the student rankings associated with each match. Both sets of
preferences must be strict; use `STB()` or similar to preprocess
if your data does not satisfy this.

Set `rev=true` to use school-proposing DA instead.
"""
function deferredacceptance(students        ::Union{AbstractArray{Int, 2}, AbstractArray{UInt, 2}},
            schools         ::Union{AbstractArray{Int, 2}, AbstractArray{UInt, 2}},
            capacities_in   ::Union{AbstractArray{Int, 1}, AbstractArray{UInt, 1}};
            verbose         ::Bool=false,
            rev             ::Bool=false)
    n, m = size(schools)
    @assert (m,) == size(capacities_in)  "Dim mismatch between schools and capacities"
    @assert (m, n) == size(students)     "Dim mismatch between schools and students"

    done = false
    nit = 0

    students_inv = mapslices(invperm, students, dims=1)
    schools_inv = mapslices(invperm, schools, dims=1)

    if rev == false
        capacities = vcat(capacities_in, n)  # For students who never get assigned
        curr_assn = students_inv[1, :]

        while !done
            nit += 1
            verbose ? println("Round $nit") : nothing
            done = true
            proposals = falses(n, m + 1)
            for (s, c) in enumerate(curr_assn)
                proposals[s, c] = true
            end
            for (c, S) in enumerate(eachcol(proposals[:, 1:m]))
                 rejections = filter(i->S[i], schools_inv[:, c])[(capacities[c] + 1):end]
                for s in rejections
                    done = false
                    verbose ? println("  School $c rejects student $s") : nothing
                    curr_assn[s] = get(students_inv, CartesianIndex(students[c, s] + 1, s), m + 1)
                end
            end
        end
        verbose ? println("DA terminated in $nit iterations") : nothing

        return curr_assn, [get(students, CartesianIndex(c, s), m + 1) for (s, c) in enumerate(curr_assn)]

    else
        not_yet_rejected = trues(n, m)
        curr_assn = [schools_inv[:, c][1:q] for (c, q) in enumerate(capacities_in)]

        while !done
            nit += 1
            verbose ? println("Round $nit") : nothing
            done = true
            proposals = falses(m, n) # May not need n+1 here
            for (c, S) in enumerate(curr_assn)
                proposals[c, S] .= true
            end
            for (s, C) in enumerate(eachcol(proposals))
                rejections = filter(i->C[i], students_inv[:, s])[2:end]
                for c in rejections
                    done = false
                    verbose ? println("  Student $s rejects school $c") : nothing
                    not_yet_rejected[s, c] = false
                    curr_assn[c] = filter(x -> not_yet_rejected[x, c],
                                          schools_inv[:, c])[1:min(end, capacities_in[c])]
                end
            end
        end
        verbose ? println("DA terminated in $nit iterations") : nothing

        # Compute the assignment from students' perspective
        students_assn = m .+ ones(Int, n)
        for (c, S) in enumerate(curr_assn)
            for s in S
                students_assn[s] = c
            end
        end

        return students_assn, [get(students, CartesianIndex(c, s), m + 1) for (s, c) in enumerate(students_assn)]
    end
end

"""
    cycle_DFS(edges)

Given the edges of a graph (in dictionary form), uses depth-first search
to find cycles. Assumes all cycles are node disjoint.
"""
function cycle_DFS(edges::Dict{Int, Int})::Set{Array{Int,1}}
    if isempty(edges)
        return Set{Vector{Int}}()
    end

    nodes = union(edges...) # try union(keys(edges), values(edges))
    visited = Dict{Int, Bool}((i, false) for i in nodes)
    out = Set{Vector{Int}}()
    for nd in nodes
        if !visited[nd]
            curr_node = nd
            stack = Vector{Int}()
            while true
                if curr_node == :deadend
                    break
                elseif curr_node in stack
                    visited[curr_node] = true
                    push!(out, stack[findfirst( x -> x == curr_node, stack):end])
                    break
                else
                    visited[curr_node] = true
                    push!(stack, curr_node)
                    curr_node = get(edges, curr_node, :deadend)
                end
            end
        end
    end

    return out
end

"""
    TTC(students_inv, assn; verbose)

Uses the top-trading cycles allocation to find the market core, given an initial
assignment. The implementation follows Nisan et al. (2007), §10.3.
"""
function TTC(students_inv   ::Union{AbstractArray{Int, 2}, AbstractArray{UInt, 2}},
             assn           ::Union{AbstractArray{Int, 1}, AbstractArray{UInt, 1}};
             verbose        ::Bool=false,
             )::Union{AbstractArray{Int, 1}, AbstractArray{UInt, 1}}
    (m, n) = size(students_inv)
    @assert (n, ) == size(assn) "Size mismatch between students_inv and assn"

    prev_assn = copy(assn)
    curr_assn = copy(assn)
    swaps_done = falses(n)

    for k in 1:(m - 1)
        verbose ? println("Searching for Pareto-improving cycles at rank $k") : nothing
        verbose ? println("Current assignment: $curr_assn") : nothing

        swap_requests = Dict{Int, Int}()
        for s in 1:n
            # In the symmetrical case, curr_assn is always a permutation, and you could
            # invert it at the outset for a marginal performance improvement.
            swap_targets = findall(curr_assn .== students_inv[k, s])
            if !isempty(swap_targets) && !swaps_done[s]
                # We can use any index here, but randomness seems fair.
                swap_requests[s] = rand(swap_targets)
            end
        end
        verbose ? println("  Swap requests: $swap_requests") : nothing

        cycles = cycle_DFS(swap_requests)

        for cyc in cycles
            verbose ? println("  Cycle involving students $cyc") : nothing
            for s in cyc
                curr_assn[s] = prev_assn[swap_requests[s]]
                swaps_done[s] = true
            end
        end
    end

    return curr_assn
end

"""
    RSD(students_inv, capacities)

Random serial dictatorship mechanism for one-sided matching (i.e. schools have neutral
preferences).
"""
function RSD(students_inv   ::Union{AbstractArray{Int, 2}, AbstractArray{UInt, 2}},
             capacities_in  ::Union{AbstractArray{Int, 1}, AbstractArray{UInt, 1}},
             )::Union{AbstractArray{Int, 1}, AbstractArray{UInt, 1}}
    (m, n) = size(students_inv)
    @assert (m, ) == size(capacities_in) "Size mismatch between students_inv and capacities"

    capacities = copy(capacities_in)
    assn = ones(Int, n)

    order = randperm(n)
    for s in order
        for k in 1:(m + 1)
            if k == m + 1
                assn[s] = m + 1
                break
            elseif capacities[students_inv[k, s]] > 0
                assn[s] = students_inv[k, s]
                capacities[assn[s]] -= 1
                break
            else
                k += 1
            end
        end
    end

    return assn
end

"""
    TTC_match(students, capacities; verbose)

Uses TTC to find a heuritically optimal one-sided school assignment. Seeds with RSD.
"""
function TTC_match(students     ::Union{AbstractArray{Int, 2}, AbstractArray{UInt, 2}},
                   capacities   ::Union{AbstractArray{Int, 1}, AbstractArray{UInt, 1}};
                   verbose      ::Bool=false,
                  )::Tuple{Union{AbstractArray{Int, 1}, AbstractArray{UInt, 1}},
                           Union{AbstractArray{Int, 1}, AbstractArray{UInt, 1}}}
    (m, n) = size(students)
    students_inv = mapslices(invperm, students, dims=1)
    assn_ = RSD(students_inv, capacities)
    assn = TTC(students_inv, assn_, verbose=verbose)

    return assn, [get(students, CartesianIndex(c, s), m + 1) for (s, c) in enumerate(assn)]
end

"""
    DA_rank_dist(students, schools, capacities; verbose, rev)

Convenience function that runs DA and outputs the cumulative rank
distribution data.
"""
function DA_rank_dist(students      ::Union{AbstractArray{Int, 2}, AbstractArray{UInt, 2}},
                      schools       ::Union{AbstractArray{Int, 2}, AbstractArray{UInt, 2}},
                      capacities    ::Union{AbstractArray{Int, 1}, AbstractArray{UInt, 1}};
                      verbose       ::Bool=false,
                      rev           ::Bool=false,
                     )::Union{AbstractArray{Int, 1}, AbstractArray{UInt, 1}}
    (n, m) = size(schools)
    dist_cmap = countmap(deferredacceptance(students, schools, capacities, verbose=verbose, rev=rev)[2])
    rank_hist = [get(dist_cmap, i, 0) for i in 1:m]

    return cumsum(rank_hist)
end

"""
    isstable(students, schools, capacities, assn; verbose)

Test if a discrete matching is stable. Allows for ties in school rankings.
"""
function isstable(students      ::Union{AbstractArray{Int, 2}, AbstractArray{UInt, 2}},
                  schools       ::Union{AbstractArray{Int, 2}, AbstractArray{UInt, 2}},
                  capacities    ::Union{AbstractArray{Int, 1}, AbstractArray{UInt, 1}},
                  assn          ::Union{AbstractArray{Int, 1}, AbstractArray{UInt, 1}};
                  verbose       ::Bool=false,
                 )::Bool
    crit = trues(3)
    (m, n) = size(students)
    @assert (n, m) == size(schools)   "Dim mismatch between students and schools"
    @assert (m,) == size(capacities)  "Dim mismatch between data and capacities"
    @assert (n,) == size(assn)        "Dim mismatch between data and assn"
    x = falses(m, n)

    for (s, c) in enumerate(assn)
        if c ≤ m
            x[c, s] = true
        end
    end

    crit[1] = all(sum(x[:, s]) ≤ 1 for s in 1:n)                # Student feas
    crit[2] = all(sum(x[c, :]) ≤ capacities[c] for c in 1:m)    # School feas
    crit[3] = all(capacities[c] * x[c, s] +                      # Stability
                  capacities[c] * x[:, s]' * (students[:, s] .≤ students[c, s]) +
                  x[c, :]' * (schools[:, c] .≤ schools[s, c]) ≥ capacities[c]
                  for c in 1:m, s in 1:n)

    res = all(crit)

    if verbose
        for (test, pass) in zip(["Student feas. ",
                                 "School feas.  ",
                                 "Stability     "], crit)
            println("$test:  $pass")
        end
    end

    return res
end

max_stars_repo_name: maxkapur/DeferredAcceptance

max_issues_repo_name: maxkapur/DeferredAcceptance

max_forks_repo_name: maxkapur/DeferredAcceptance


text:
module AsymmetricRisk

using Statistics, LinearAlgebra
using StatsBase
using KernelDensity
using HCubature
using Distributions
using EmpiricalCDFs
import QuadGK: quadgk
# using GLM
# using NamedTupleTools

export semivariance, lowerpartialmom, upperpartialmom, valueatrisk, expectedshortfall

export quantiledep, exceedancecor, downsidebeta, upsidebeta, coskewness, cokurtosis

include("univariate.jl")
include("bivariate.jl")

end # module


text:
# Performance comparison of various methods to sample from a truncated normal

# plot performance for each upper bound b over lower bounds a < b
const plot_results = true

using Distributions
using DataFrames

if plot_results
    using Gadfly
end

include("../src/randtn.jl")
include("../src/randtn_other.jl")

# helper functions for each algorithm
immutable RandTNChopin; end
immutable RandTNRobert; end
immutable RandTNInvCDF; end
const randfuntypes = [RandTNChopin(), RandTNRobert(), RandTNInvCDF()]

randfun(::RandTNChopin) = randtnstd
randfunargs(::RandTNChopin, d::Distribution) = (d.lower, d.upper, d.tp)
randfunid(::RandTNChopin) = :chopin

randfun(::RandTNRobert) = randtn_robert_base
randfunargs(::RandTNRobert, d::Distribution) = (d.lower, d.upper, d.tp)
randfunid(::RandTNRobert) = :robert

randfun(::RandTNInvCDF) = randtn_invcdf
randfunargs(::RandTNInvCDF, d::Distribution) = (d, )
randfunid(::RandTNInvCDF) = :invcdf

# timer
function randtime(f::Function, args, n::Int)
    f(args...)
    gc_disable()
    t = @elapsed for i = 1:n
        f(args...)
    end
    gc_enable()
    gc()
    t
end

# run tests for different upper bounds b for all lower bounds a < b
const anum = 100
const amin = -2.5
const amax = 4.0
const bs = [tnzig_xmin - 0.1, tnzig_xmin + 0.1,
            tnzig_xmax - 0.1, tnzig_xmax + 0.1, Inf]
const nsamples = 10^6

const testfns = [randtn, randtn_robert, randtn_invcdf]
const testids = [:chopin, :robert, :invcdf]

for i = 1:length(bs)
    b = bs[i]
    println("upper bound b = ", b)
    as = linspace(amin, min(b - 0.01, amax), anum)
    df = DataFrame([Float64, Float64, Symbol], [:a, :t, :algo], anum * length(testfns))
    for j = 1:length(testfns)
        randfuntype = randfuntypes[j]
        testfn = randfun(randfuntype)
        testid = randfunid(randfuntype)
        print("Testing ", string(testid))
        for k = 1:anum
            d = TruncatedNormal(0, 1, as[k], b)
            l = (j-1) * anum + k
            df[l, :a] = as[k]
            df[l, :t] = randtime(testfn, randfunargs(randfuntype, d), nsamples)
            df[l, :algo] = testid
        end
        println(", <t> = ", mean(df[((j-1)*anum+1):(j*anum), :t]))
    end
    if plot_results
        bfmt = @sprintf("b = %5.3f", b)
        p = plot(df, x="a", y="t", color="algo", Geom.line, Guide.title(bfmt))
        draw(PDF("randtn_perf_b$bfmt.pdf", 15cm, 10cm), p)
    end
end

max_stars_repo_name: jdrugo/fasttn

max_issues_repo_name: jdrugo/fasttn

max_forks_repo_name: jdrugo/fasttn


text:
function splitidentifier(identifier)
    parts = split(identifier, '.')
    modulename = join(parts[1:end-1], '.')
    bindingname = Symbol(parts[end])
    return modulename, bindingname
end

const MODNAMES = Ref(string.(Base.Docs.modules))

function getmodule(s::String)
    # Refresh `MODNAMES`
    if length(Base.Docs.modules) != length(MODNAMES[])
        MODNAMES[] = string.(Docs.modules)
    end
    idx = findfirst(MODNAMES[] .== s)
    return isnothing(idx) ? nothing : Docs.modules[idx]
end

function getdocs(m::Module, s::Symbol)
    bs = [(b, md) for (b, md) in Docs.meta(m)]
    bindingnames = [b.var for (b, _) in bs]
    idx = findfirst(bindingnames .== s)
    return isnothing(idx) ? nothing : bs[idx][2]
end

# Conversion to x-expressions

function Base.convert(::Type{XTree}, md::Base.Docs.MultiDoc)
    return XNode(:docs, XTree[convert(XTree, v) for v in values(md.docs)])
end

function Base.convert(::Type{XTree}, doc::Base.Docs.DocStr)
    s = collect(doc.text)[1]
    # Parse from docstring and filter out redundant line breaks
    xdoc = replace(Pollen.parse(s, Pollen.Markdown()), XLeaf("\n"), SelectTag(:br))
    return XNode(:doc, doc.data, children(xdoc))
end


text:
# exp_xform
SUITE["exp_xform"] = BenchmarkGroup()
let T=ComplexF32, N=100, M=10000, D=20, L=10
    rng = MersenneTwister(0)
    X = randn(rng, T, N, L)
    U = randn(rng, T, D, N)
    V = randn(rng, T, D, M)

    for mode in (:matrix, :row, :element, :column)
        # timing tests: :matrix is fastest, with :row a close 2nd
        fn(x, u, v) = exp_xform(x, u, v; mode=mode)
        SUITE["exp_xform"][string(mode)] = @benchmarkable $fn($X, $U, $V)
    end
end


text:
# Functions are typically designated by the function name and parentheses, with an argument in the parentheses
println("println is a function and this is the argument")

# If there is no input, a blank line is printed
println()

# Operators are exceptions to the parentheses rule
a = 12 + 14
println(a)

# Creating a basic function
function add_numbers(x, y)
    x + y
    end

# Using the function
println(add_numbers(1, 2))
println(add_numbers(2, 3))
println(add_numbers(3, 4))

max_stars_repo_name: KorfLab/julie

max_issues_repo_name: KorfLab/julie

max_forks_repo_name: KorfLab/julie


text:
module Markov

include("MarkovCategories.jl")
include("StochExpr.jl")

end


text:
@muladd function calc_twopoint_random!(_dW, integrator, dW)
    @.. _dW = ifelse(sign(dW) > 0.0, integrator.sqdt, -integrator.sqdt)
end

@muladd function calc_twopoint_random(integrator, dW)
  return sign(dW) > 0.0 ? integrator.sqdt : -integrator.sqdt
end


text:
@testset "Check file availability" for fn in rs.DataFiles
    @test isfile(fn)
    df = DataFrame(CSV.File(fn))
end

exp = rs.getExpLevels()
cond = rs.getConditions()

prNamesRNA = [:Rtg1, :Rtg2, :Rtg3, :Bmh, :Mks]
prNamesCond = ["Rtg1", "Rtg2", "Rtg3", "Mks"]

@testset "Protein name availability (RNA-Seq)" for pr in prNamesRNA
    @test pr in collect(keys(exp))
end
@testset "Protein name availability (Conditions)" for pr in prNamesCond
    @test pr in names(cond)
end

@test "s" in names(cond)
@test "gfp" in names(cond)
@test "Trans2Nuc" in names(cond)

# Solutions of rtgM4
@testset "Solutions to rtgM4" begin
    df = rs.readCSV(rs.DataFiles.solution_rtgM4)
    num_sols = size(df)[1]
    for i in 1:num_sols
        @test rs.isValid(rs.rtgM4(i))
    end
end


text:
# Create tensors with certain symmetry from medium parameters

"""
    iso(; vp=nothing, vs=nothing, lam=nothing, mu=nothing, K=nothing, G=nothing) -> C

Return an isotropic Voigt stiffness matrix `C` from a pair of the following:

`vp`, `vs` : Isotropic velocities in m/s

`lam`, `mu`: Lamé parameters **divided by density** in m^2/s^2

`K`, `G`   : Bulk and shear moduli **divided by density** in m^2/s^2
"""
function iso(; vp=nothing, vs=nothing, lam=nothing, mu=nothing, K=nothing, G=nothing)
    C = zero(EC)
    @inbounds begin
        if vp != nothing && vs != nothing
            C[1,1] = vp^2
            C[4,4] = vs^2
            C[1,2] = C[1,1] - 2C[4,4]
        elseif lam != nothing && mu != nothing
            C[1,1] = lam + 2mu
            C[4,4] = mu
            C[1,2] = lam
        elseif K != nothing && G != nothing
            C[1,1] = K + 4G/3
            C[4,4] = G
            C[1,2] = C[1,1] - 2C[4,4]
        else
            throw(ArgumentError("iso: Must request a pair of vp,vs; lam,mu; K,G"))
        end
        C[2,2] = C[3,3] = C[1,1]
        C[5,5] = C[6,6] = C[4,4]
        C[1,3] = C[2,3] = C[3,1] = C[3,2] = C[2,1] = C[1,2]
    end
    return C
end

"""
    vti(; vpv, vsv, vph, vsh, eta)

Return a density-normalised elastic tensor for a VTI medium defined using the velocities
Vpv, Vsv, Vph, Vsh and anisotropy parameter η.

    vti(; A, C, L, N, F)

Return an unnormalised elastic tensor for a VTI medium defined using Love's (1927) notation.

The constants are symmetric about the 3-axis.

### References

- Love AEH (1927) A Treatise on the Mathematical Theory of Elasticity. New York: Dover Publications.
"""
function vti(; vpv=nothing, vsv=nothing, vph=nothing, vsh=nothing, eta=nothing,
        A=nothing, C=nothing, L=nothing, N=nothing, F=nothing)
    if all(x->x!==nothing, (vpv, vsv, vph, vsh, eta))
        c = zero(EC)
        c[1,1] = c[2,2] = vph^2
        c[3,3] = vpv^2
        c[4,4] = c[5,5] = vsv^2
        c[6,6] = vsh^2
        c[1,3] = c[2,3] = c[3,1] = c[3,2] = eta*(c[1,1] - 2c[4,4])
        c[1,2] = c[2,1] = c[1,1] - 2c[6,6]
        c
    elseif all(x->x!==nothing, (A, C, L, N, F))
        EC{DEFAULT_FLOAT}((  A,  A-2N,  F,   0,   0,   0,
                           A-2N,   A,   F,   0,   0,   0,
                             F,    F,   C,   0,   0,   0,
                             0,    0,   0,   L,   0,   0,
                             0,    0,   0,   0,   L,   0,
                             0,    0,   0,   0,   0,   N))
    else
        throw(ArgumentError("Insufficient parameters defined."))
    end
end

"""
    thom(vp, vs, eps, gam, del) -> C

Return the 6x6 Voigt matrix defined `C` by the weak anisotropy parameters of
Thomsen (1986).

Output is density-normalised tensor.

### References

Thomsen, L. (1986).  Weak elastic anisotropy.  Geophysics, 51, 10, 1954-1966.
"""
function thom(vp, vs, eps, gam, del)
    if vp <= 0
        error("CIJ.thom: vp must be greater than 0")
    elseif vs < 0
        error("CIJ.thom: vs must be greater than or equal to 0")
    end

    C = zero(EC)
    @inbounds begin
        C[3,3] = vp^2
        C[1,1] = C[2,2] = C[3,3]*(2*eps + 1)
        C[4,4] = C[5,5] = vs^2
        C[6,6] = C[4,4]*(2*gam + 1)

        b = 2C[4,4]
        term = C[3,3] - C[4,4]
        c = C[4,4]^2 - (2*del*C[3,3]*term + term^2)
        d = b^2 - 4*c
        if (d < 0)
            error("CIJ.thom: S-velocity too high or delta too negative")
        end
        C[1,3] = C[2,3] = C[3,1] = C[3,2] = -b/2 + sqrt(d)/2
        C[1,2] = C[2,1] = C[1,1] - 2*C[6,6]
    end
    return C
end

"""
    thom_st(vp, vs, eps, gam, delst) -> C

Return the 6x6 Voigt matrix `C` defined by the general anisotropy parameters of
Thomsen (1986).

Output has same units as input.
"""
function thom_st(vp, vs, eps, gam, delst)
    if vp <= 0
        error("CIJ.thom: vp must be greater than 0")
    elseif vs < 0
        error("CIJ.thom: vs must be greater than or equal to 0")
    end

    C = zero(EC)
    @inbounds begin
        C[3,3] = vp^2
        C[1,1] = C[2,2] = C[3,3]*(2*eps + 1)
        C[4,4] = C[5,5] = vs^2
        C[6,6] = C[4,4]*(2*gam + 1)
        a = 2.0
        b = 4*C[4,4]
        c = C[4,4]^2 - 2*delst*C[3,3]^2 - (C[3,3] - C[4,4])*(C[1,1] + C[3,3] - 2*C[4,4])
        if b^2 - 4*a*c < 0
            error("CIJ.thom_st: S velocity too high or delta too negative")
        end
        C[1,3] = C[3,1] = C[2,3] = C[3,2] = (-b + sqrt(b^2 - 4*a*c))/(2*a)
        C[1,2] = C[2,1] = C[1,1] - 2C[6,6]
    end
    return C
end

"""
    global_VTI(vp, vs, xi, phi, eta) -> C

Return the 6x6 Voigt matrix `C` defined by the radial anisotropy parameters
as used typically in global seismology.
"""
function global_VTI(vp, vs, xi, phi, eta)
    if vp <= 0
        error("CIJ.global_VTI: vp must be greater than 0")
    elseif vs < 0
        error("CIJ.global_VTI: vs must be greater than or equal to 0")
    elseif xi <= 0 || phi <= 0 || eta <= 0
        error("CIJ.global_VTI: xi, phi and eta must be greater then 0")
    end

    L = 15*((3*phi + 8 + 4*eta)*vs^2 - (phi + 1 - 2*eta)*vp^2) /
        ((6 + 4*eta + 5*xi)*(3*phi + 8 +4*eta) - 8*(phi + 1 -2*eta)*(1 - eta))
    A = (15*vp^2 - 8*(1 - eta)*L) / (3*phi + 8 + 4*eta)
    F = eta*(A - 2L)
    C = phi*A
    N = xi*L
    C12 = A - 2N
    return EC{DEFAULT_FLOAT}(( A,  C12, F,  0,  0,  0,
                              C12,  A,  F,  0,  0,  0,
                               F,   F,  C,  0,  0,  0,
                               0,   0,  0,  L,  0,  0,
                               0,   0,  0,  0,  L,  0,
                               0,   0,  0,  0,  0,  N))
end

"""
    Panning_VTI(vp, vs, xi, phi, eta) -> C

Return the 6x6 Voigt matrix `C` defined by the radial anisotropy parameters
defined by Panning and Romanowicz as an approximation to general VTI.
"""
function Panning_VTI(vp, vs, xi, phi)
    L = 3*vs^2/(2 + xi)
    N = xi*L
    A = 5*vp^2/(4 + phi)
    C = phi*A
    F = A - 2*L
    C12 = A - 2*N
    return EC{DEFAULT_FLOAT}(( A,  C12, F,  0,  0,  0,
                              C12,  A,  F,  0,  0,  0,
                               F,   F,  C,  0,  0,  0,
                               0,   0,  0,  L,  0,  0,
                               0,   0,  0,  0,  L,  0,
                               0,   0,  0,  0,  0,  N))
end

"""
    lame(vp, vs, ρ) -> λ, μ, κ

Return the Lamé parameters `λ`, `μ` and `κ` from the velocities `vp` and `vs`
and density `ρ`.
"""
function lame(vp, vs, ρ)
    μ = ρ*vs^2
    λ = ρ*(vp^2 - 2vs^2)
    κ = ρ*vp^2 - 4/3*μ
    λ, μ, κ
end


text:
"""
    SoilHeatParameterizations

Functions for volumetric heat capacity, temperature as a function
of volumetric internal energy, saturated thermal conductivity, thermal
conductivity, relative saturation and the Kersten number are included.
Heat capacities denoted by `ρc_` are volumetric, while `cp_` denotes an isobaric
specific heat capacity.
"""
module SoilHeatParameterizations

using CLIMAParameters
using CLIMAParameters.Planet: ρ_cloud_liq, ρ_cloud_ice, cp_l, cp_i, T_0, LH_f0
using CLIMAParameters.Atmos.Microphysics: K_therm
using DocStringExtensions

export volumetric_heat_capacity,
    volumetric_internal_energy,
    saturated_thermal_conductivity,
    thermal_conductivity,
    relative_saturation,
    kersten_number,
    volumetric_internal_energy_liq,
    temperature_from_ρe_int,
    k_solid,
    k_dry,
    ksat_unfrozen,
    ksat_frozen

"""
    function temperature_from_ρe_int(
        ρe_int::FT,
        θ_i::FT,
        ρc_s::FT,
        param_set::AbstractParameterSet
    ) where {FT}

Computes the temperature of soil given `θ_i` and volumetric
internal energy `ρe_int`.
"""
function temperature_from_ρe_int(
    ρe_int::FT,
    θ_i::FT,
    ρc_s::FT,
    param_set::AbstractParameterSet,
) where {FT}
    _ρ_i = FT(ρ_cloud_ice(param_set))
    _T_ref = FT(T_0(param_set))
    _LH_f0 = FT(LH_f0(param_set))
    T = _T_ref + (ρe_int + θ_i * _ρ_i * _LH_f0) / ρc_s
    return T
end

"""
    volumetric_heat_capacity(
        θ_l::FT,
        θ_i::FT,
        ρc_ds::FT,
        param_set::AbstractParameterSet
    ) where {FT}

Compute the expression for volumetric heat capacity.
"""
function volumetric_heat_capacity(
    θ_l::FT,
    θ_i::FT,
    ρc_ds::FT,
    param_set::AbstractParameterSet,
) where {FT}
    _ρ_i = FT(ρ_cloud_ice(param_set))
    ρcp_i = FT(cp_i(param_set) * _ρ_i)

    _ρ_l = FT(ρ_cloud_liq(param_set))
    ρcp_l = FT(cp_l(param_set) * _ρ_l)

    ρc_s = ρc_ds + θ_l * ρcp_l + θ_i * ρcp_i
    return ρc_s
end

"""
    volumetric_internal_energy(
        θ_i::FT,
        ρc_s::FT,
        T::FT,
        param_set::AbstractParameterSet
    ) where {FT}

Compute the expression for volumetric internal energy.
"""
function volumetric_internal_energy(
    θ_i::FT,
    ρc_s::FT,
    T::FT,
    param_set::AbstractParameterSet,
) where {FT}
    _ρ_i = FT(ρ_cloud_ice(param_set))
    _LH_f0 = FT(LH_f0(param_set))
    _T_ref = FT(T_0(param_set))
    ρe_int = ρc_s * (T - _T_ref) - θ_i * _ρ_i * _LH_f0
    return ρe_int
end

"""
    saturated_thermal_conductivity(
        θ_l::FT,
        θ_i::FT,
        κ_sat_unfrozen::FT,
        κ_sat_frozen::FT
    ) where {FT}

Compute the expression for saturated thermal conductivity of soil matrix.
"""
function saturated_thermal_conductivity(
    θ_l::FT,
    θ_i::FT,
    κ_sat_unfrozen::FT,
    κ_sat_frozen::FT,
) where {FT}
    θ_w = θ_l + θ_i
    if θ_w < eps(FT)
        κ_sat = FT(0.0)
    else
        κ_sat = FT(κ_sat_unfrozen^(θ_l / θ_w) * κ_sat_frozen^(θ_i / θ_w))
    end

    return κ_sat
end

"""
    relative_saturation(
            θ_l::FT,
            θ_i::FT,
            porosity::FT
    ) where {FT}

Compute the expression for relative saturation.
"""
function relative_saturation(θ_l::FT, θ_i::FT, porosity::FT) where {FT}
    S_r = (θ_l + θ_i) / porosity
    return S_r
end
"""
    kersten_number(
        θ_i::FT,
        S_r::FT,
        soil_param_functions::PS
    ) where {FT, PS}

Compute the expression for the Kersten number.
"""
function kersten_number(
    θ_i::FT,
    S_r::FT,
    soil_param_functions::PS,
) where {FT, PS}
    a = soil_param_functions.a
    b = soil_param_functions.b
    ν_ss_om = soil_param_functions.ν_ss_om
    ν_ss_quartz = soil_param_functions.ν_ss_quartz
    ν_ss_gravel = soil_param_functions.ν_ss_gravel

    if θ_i < eps(FT)
        K_e =
            S_r^((FT(1) + ν_ss_om - a * ν_ss_quartz - ν_ss_gravel) / FT(2)) *
            (
                (FT(1) + exp(-b * S_r))^(-FT(3)) -
                ((FT(1) - S_r) / FT(2))^FT(3)
            )^(FT(1) - ν_ss_om)
    else
        K_e = S_r^(FT(1) + ν_ss_om)
    end
    return K_e
end

"""
    thermal_conductivity(
        κ_dry::FT,
        K_e::FT,
        κ_sat::FT
    ) where {FT}

Compute the expression for thermal conductivity of soil matrix.
"""
function thermal_conductivity(κ_dry::FT, K_e::FT, κ_sat::FT) where {FT}
    κ = K_e * κ_sat + (FT(1) - K_e) * κ_dry
    return κ
end

"""
    volumetric_internal_energy_liq(
        T::FT,
        T_ref::FT,
    ) where {FT}

Compute the expression for the volumetric internal energy of liquid water.
"""
function volumetric_internal_energy_liq(
    T::FT,
    param_set::AbstractParameterSet,
) where {FT}
    _T_ref = FT(T_0(param_set))
    _ρ_l = FT(ρ_cloud_liq(param_set))
    ρcp_l = FT(cp_l(param_set) * _ρ_l)
    ρe_int_l = ρcp_l * (T - _T_ref)
    return ρe_int_l
end

"""
    function k_solid(
        ν_ss_om::FT,
        ν_ss_quartz::FT,
        κ_quartz::FT,
        κ_minerals::FT,
        κ_om::FT,
    ) where {FT}

Computes the thermal conductivity of the solid material in soil.
The `_ss_` subscript denotes that the volumetric fractions of the soil
components are referred to the soil solid components, not including the pore
space.
"""
function k_solid(
    ν_ss_om::FT,
    ν_ss_quartz::FT,
    κ_quartz::FT,
    κ_minerals::FT,
    κ_om::FT,
) where {FT}
    return κ_om^ν_ss_om *
           κ_quartz^ν_ss_quartz *
           κ_minerals^(FT(1) - ν_ss_om - ν_ss_quartz)
end

"""
    function ksat_frozen(
        κ_solid::FT,
        porosity::FT,
        κ_ice::FT
    ) where {FT}

Computes the thermal conductivity for saturated frozen soil.
"""
function ksat_frozen(κ_solid::FT, porosity::FT, κ_ice::FT) where {FT}
    return κ_solid^(FT(1.0) - porosity) * κ_ice^(porosity)
end

"""
    function ksat_unfrozen(
        κ_solid::FT,
        porosity::FT,
        κ_l::FT
    ) where {FT}

Computes the thermal conductivity for saturated unfrozen soil.
"""
function ksat_unfrozen(κ_solid::FT, porosity::FT, κ_l::FT) where {FT}
    return κ_solid^(FT(1.0) - porosity) * κ_l^porosity
end

"""
    function ρb_ss(porosity::FT, ρp::FT) where {FT}
Computes the dry soil bulk density from the dry soil particle
density.
"""
function ρb_ss(porosity::FT, ρp::FT) where {FT}
    return (FT(1.0) - porosity) * ρp
end

"""
    function k_dry(
        param_set::AbstractParameterSet
        soil_param_functions::PS,
    ) where {PS}

Computes the thermal conductivity of dry soil.
"""
function k_dry(
    param_set::AbstractParameterSet,
    soil_param_functions::PS,
) where {PS}
    κ_dry_parameter = soil_param_functions.κ_dry_parameter
    FT = typeof(κ_dry_parameter)
    porosity = soil_param_functions.porosity
    ρp = soil_param_functions.ρp
    κ_solid = soil_param_functions.κ_solid
    κ_air = FT(K_therm(param_set))
    ρb_val = ρb_ss(porosity, ρp)
    numerator = (κ_dry_parameter * κ_solid - κ_air) * ρb_val + κ_air * ρp
    denom = ρp - (FT(1.0) - κ_dry_parameter) * ρb_val
    return numerator / denom
end

end # Module


text:
# Copyright (c) 2021 Idiap Research Institute, http://www.idiap.ch/
#  Niccolò Antonello <nantonel@idiap.ch>

export get_paths, collectpaths

struct PathIterator{W,A,D,F}
  bfs::BFS{W,A,D,F}
end

"""
  `get_paths(fst::FST, [i = get_initial(fst;single=true)])`

Returns an iterator that generates all of the possible paths of `fst` starting from `i`.

Notice that the `fst` must be acyclic for the algorithm to stop (see [`is_acyclic`](@ref)).
"""
get_paths(fst::WFST, i=get_initial(fst;single=true)) = PathIterator(BFS(fst,i))

function Base.iterate(it::PathIterator, q=init_bfs(it.bfs.fst, it.bfs.initial))
  if isempty(q)
    return nothing
  else
    while true
      (i,p), q = iterate(it.bfs,q)
      if isfinal(it.bfs.fst,i)
        w = get_final_weight(it.bfs.fst,i)
        pf = update_weight(p,w)
        return pf, q
      end
    end
  end
end

"""
  `collectpaths(fst::FST, [i = get_initial(fst;first=true)])`

Returns an array containing all of the possible paths of `fst` starting from `i`.

Notice that the `fst` must be acyclic for the algorithm to stop (see [`is_acyclic`](@ref)).
"""
function collectpaths(fst::WFST{W,A,D,I,O},
                      i::D=get_initial(fst;single=true)) where {W,D,I,O,
                                                                  A <: Arc{W,D}}
  paths = Vector{Path{W,D}}()
  if !isempty(fst)
    pathsIt = BFS(fst, i)
    for (i,p) in pathsIt
      if isfinal(fst,i)
        w = get_final_weight(fst,i)
        pf = update_weight(p,w)
        push!(paths,pf)
      end
    end
  end
  return paths
end


text:
using GameOfLife
using Test

@testset "world.jl" begin
    @test_nowarn WorldState()
    @test_nowarn WorldState(initial_state = rand(0:1, 10, 10))
    @test_nowarn WorldState(initial_state = rand(0:1, 3, 3))

    @test_throws ErrorException WorldState(initial_state = rand(0:1, 10, 10, 10))
    @test_throws ErrorException WorldState(initial_state = rand(0:1, 10, 10, 1))
    @test_throws ErrorException WorldState(initial_state = rand(0:1, 10))
    @test_throws ErrorException WorldState(initial_state = rand(0:1, 2, 2))
    @test_throws ErrorException WorldState(initial_state = rand(0:1, 3, 2))
    @test_throws ErrorException WorldState(initial_state = rand(0:2, 3, 3))

    @test GameOfLife.state_check([2]) == false
    @test GameOfLife.state_check([0]) == true
    @test GameOfLife.state_check([1]) == true
    @test GameOfLife.state_check([0 1 1; 1 0 1;]) == true
    @test GameOfLife.state_check([0 1 1; 1 0 -1;]) == false

    world = WorldState()
    @test length(world)  == 100
    @test length(WorldState(max_days = 999999)) == 999999

    @test iterate(world, 10)[2] == 11
    @test iterate(world, 100)[2] == 101
    @test isnothing(iterate(world, 101)) == true
end

max_stars_repo_name: VasanthManiVasi/Conways-Game-of-Life

max_issues_repo_name: VasanthManiVasi/Conways-Game-of-Life

max_forks_repo_name: VasanthManiVasi/Conways-Game-of-Life


text:
"""
Module to access nflscrapR-data files

See https://github.com/ryurko/nflscrapR-data/ for additional data documentation.
"""
module NFLScrapR

using  CSV
using  DataFrames
using  Downloads
using  Logging
using  Scratch: @get_scratch!
using  ..NFLTables: POST, PRE, REG, SeasonPart

export getgamedata, getplaydata

"""
seasons with nflscrapR data
"""
const SEASONS = tuple(2009:2019...)

"""
The root of the nflscrapR-data GitHub repo, after the raw redirect.
"""
const REPOROOT = "https://github.com/ryurko/nflscrapR-data/raw/master/"

function getlocalpath(datatype::Symbol, seasonpart::SeasonPart, season::Integer)::String
    datatype === :games || datatype === :play_by_play || error("Invalid datatype $datatype")
    season in SEASONS || error("Invalid season $season")
    shortdata = datatype === :play_by_play ? "pbp" : string(datatype)
    shortpart = seasonpart |> string |> lowercase
    return "$(shortpart)_$(shortdata)_$season.csv"
end

function getremotepath(datatype::Symbol, seasonpart::SeasonPart, season::Integer)::String
    seasonfolder = seasonpart === REG ? "regular_season" : "$(lowercase(string(seasonpart)))_season"
    return joinpath("$(datatype)_data", seasonfolder, getlocalpath(datatype, seasonpart, season))
end

function load_data_from_disk(datatype::Symbol, seasonpart::SeasonPart, season::Integer; redownload=false, reporoot=REPOROOT)
    dirpath = @get_scratch!("nflscrapR")
    localpath = joinpath(dirpath, getlocalpath(datatype, seasonpart, season))
    if redownload || !isfile(localpath)
        remotepath = joinpath(reporoot, getremotepath(datatype, seasonpart, season))
        @info "downloading..." remotepath localpath
        Downloads.download(remotepath, localpath)
    end
    return DataFrame(CSV.File(localpath, missingstring="NA"))
end

"""
    getplaydata(season, part)

Create a dataframe of play-by-play data for `part` of `season`.

# Examples
```jldoctest
julia> df = NFLScrapR.getplaydata(2019, "POST");

julia> show(first(df[:, [:game_id, :home_team, :away_team, :yardline_100, :half_seconds_remaining]], 5), eltypes=false)
5×5 DataFrame
 Row │ game_id     home_team  away_team  yardline_100  half_seconds_remaining
─────┼────────────────────────────────────────────────────────────────────────
   1 │ 2020012600  APR        NPR                  75                    1800
   2 │ 2020012600  APR        NPR                  75                    1776
   3 │ 2020012600  APR        NPR                  80                    1757
   4 │ 2020012600  APR        NPR                  71                    1715
   5 │ 2020012600  APR        NPR                  71                    1709

```
"""
getplaydata(season::Integer, part::SeasonPart) = load_data_from_disk(:play_by_play, part, season)
getplaydata(season::Integer, part::String) = getplaydata(season, parse(SeasonPart, part))

"""
    getgamedata(season, part)

Create a dataframe of game data for `part` of `season`.

# Examples
```jldoctest
julia> df = NFLScrapR.getgamedata(2019, "POST");

julia> show(first(df[:, [:game_id, :home_team, :away_team, :home_score, :away_score]], 5), eltypes=false)
5×5 DataFrame
 Row │ game_id     home_team  away_team  home_score  away_score
─────┼──────────────────────────────────────────────────────────
   1 │ 2020012600  APR        NPR           missing     missing
   2 │ 2020010400  HOU        BUF                22          19
   3 │ 2020010401  NE         TEN                13          20
   4 │ 2020010500  NO         MIN                20          26
   5 │ 2020010501  PHI        SEA                 9          17

```
"""
getgamedata(season::Integer, part::SeasonPart) = load_data_from_disk(:games, part, season)
getgamedata(season::Integer, part::String) = getgamedata(season, parse(SeasonPart, part))

end  # module NFLScrapR

# to support legacy calls
nflscraprRgame = NFLScrapR.getgamedata
nflscrapRplaybyplay = NFLScrapR.getplaydata


text:
# Autogenerated wrapper script for utf8proc_jll for aarch64-apple-darwin
export libutf8proc

JLLWrappers.@generate_wrapper_header("utf8proc")
JLLWrappers.@declare_library_product(libutf8proc, "@rpath/libutf8proc.2.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libutf8proc,
        "lib/libutf8proc.2.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module PartitioningWrapper

using JavaCall

const AbstractRegion = @jimport org.hipparchus.geometry.partitioning.AbstractRegion
const AbstractSubHyperplane = @jimport org.hipparchus.geometry.partitioning.AbstractSubHyperplane
const BSPTree = @jimport org.hipparchus.geometry.partitioning.BSPTree
const BSPTreeVisitor = @jimport org.hipparchus.geometry.partitioning.BSPTreeVisitor
const BSPTreeVisitor_Order = @jimport org.hipparchus.geometry.partitioning.BSPTreeVisitor$Order
const BSPTree_LeafMerger = @jimport org.hipparchus.geometry.partitioning.BSPTree$LeafMerger
const BSPTree_VanishingCutHandler = @jimport org.hipparchus.geometry.partitioning.BSPTree$VanishingCutHandler
const BoundaryAttribute = @jimport org.hipparchus.geometry.partitioning.BoundaryAttribute
const BoundaryBuilder = @jimport org.hipparchus.geometry.partitioning.BoundaryBuilder
const BoundaryProjection = @jimport org.hipparchus.geometry.partitioning.BoundaryProjection
const BoundaryProjector = @jimport org.hipparchus.geometry.partitioning.BoundaryProjector
const BoundarySizeVisitor = @jimport org.hipparchus.geometry.partitioning.BoundarySizeVisitor
const Characterization = @jimport org.hipparchus.geometry.partitioning.Characterization
const Class = @jimport java.lang.Class
const Consumer = @jimport java.util.function.Consumer
const Embedding = @jimport org.hipparchus.geometry.partitioning.Embedding
const Enum = @jimport java.lang.Enum
const Hyperplane = @jimport org.hipparchus.geometry.partitioning.Hyperplane
const Iterable = @jimport java.lang.Iterable
const JIterator = @jimport java.util.Iterator
const JString = @jimport java.lang.String
const JVector = @jimport org.hipparchus.geometry.Vector
const List = @jimport java.util.List
const NodesSet = @jimport org.hipparchus.geometry.partitioning.NodesSet
const Object = @jimport java.lang.Object
const Optional = @jimport java.util.Optional
const Point = @jimport org.hipparchus.geometry.Point
const Region = @jimport org.hipparchus.geometry.partitioning.Region
const RegionFactory = @jimport org.hipparchus.geometry.partitioning.RegionFactory
const Region_Location = @jimport org.hipparchus.geometry.partitioning.Region$Location
const Side = @jimport org.hipparchus.geometry.partitioning.Side
const Spliterator = @jimport java.util.Spliterator
const SubHyperplane = @jimport org.hipparchus.geometry.partitioning.SubHyperplane
const SubHyperplane_SplitSubHyperplane = @jimport org.hipparchus.geometry.partitioning.SubHyperplane$SplitSubHyperplane
const Transform = @jimport org.hipparchus.geometry.partitioning.Transform

export AbstractRegion
export AbstractSubHyperplane
export BSPTree
export BSPTreeVisitor
export BoundaryAttribute
export BoundaryBuilder
export BoundaryProjection
export BoundaryProjector
export BoundarySizeVisitor
export Characterization
export Embedding
export Hyperplane
export NodesSet
export Region
export RegionFactory
export Side
export SubHyperplane
export Transform
export add
export add_all
export apply
export apply_transform
export build_convex
export build_new
export check_point
export compare_to
export contains
export copy_self
export describe_constable
export difference
export empty_hyperplane
export equals
export for_each
export get_attribute
export get_barycenter
export get_boundary_size
export get_cell
export get_class
export get_close_cuts
export get_complement
export get_cut
export get_declaring_class
export get_hyperplane
export get_inside_splitters
export get_minus
export get_offset
export get_original
export get_outside_splitters
export get_parent
export get_plus
export get_plus_inside
export get_plus_outside
export get_projected
export get_projection
export get_remaining_region
export get_size
export get_splitters
export get_tolerance
export get_tree
export hash_code
export insert_cut
export insert_in_tree
export inside_touching
export intersection
export is_empty
export is_full
export iterator
export merge
export name
export notify
export notify_all
export ordinal
export outside_touching
export project
export project_to_boundary
export prune_around_convex_cell
export reunite
export same_orientation_as
export set_attribute
export split
export spliterator
export to_space
export to_string
export to_sub_space
export touch_inside
export touch_outside
export union
export value_of
export values
export visit
export visit_internal_node
export visit_leaf_node
export visit_order
export wait
export whole_hyperplane
export whole_space
export xor

include("abstract_region.jl")
include("abstract_sub_hyperplane.jl")
include("boundary_attribute.jl")
include("boundary_builder.jl")
include("boundary_projection.jl")
include("boundary_projector.jl")
include("boundary_size_visitor.jl")
include("bsp_tree.jl")
include("bsp_tree_visitor.jl")
include("characterization.jl")
include("embedding.jl")
include("hyperplane.jl")
include("nodes_set.jl")
include("region.jl")
include("region_factory.jl")
include("side.jl")
include("sub_hyperplane.jl")
include("transform.jl")

end


text:
@generated function precomputeB(win, x, n::NTuple{D,Int}, m, M, sigma, T, U1, U2) where D
    quote
        I = zeros(Int64, (2*m+1)^D, M)
        J = zeros(Int64, M+1)
        V = zeros(T, (2*m+1)^D, M)

        J[1] = 1
        @inbounds @simd for k in 1:M
            @nexprs $D d -> xscale_d = x[d,k] * n[d]
            @nexprs $D d -> c_d = floor(Int, xscale_d)

            @nloops $D l d -> (c_d-m):(c_d+m) d->begin
                # preexpr
                gidx_d = rem(l_d+n[d], n[d]) + 1
                Iidx_d = l_d - c_d + m + 1
                idx = abs( (xscale_d - l_d) )
                tmpWin_d =  win(idx / n[d], n[d], m, sigma)

            end begin
                # bodyexpr
                v = 1
                @nexprs $D d -> v *= tmpWin_d
                i1 = 1
                @nexprs $D d -> i1 += (Iidx_d-1) * U1[d]
                i2 = 1
                @nexprs $D d -> i2 += (gidx_d-1) * U2[d]

                I[i1,k] = i2
                V[i1,k] = v
            end
            J[k+1] = J[k] + (2*m+1)^D
        end

        S = SparseMatrixCSC(prod(n), M, J, vec(I), vec(V))
        return S
    end
end

"""
precompute LUT for the multidimensional interpolation window
"""
function precomp_windowHatInvLUT(T::Type, win_hat::Function, N::NTuple{D,Int64}, sigma::Real, m::Int64) where D
    # size of oversampled grid
    n = ntuple(d->round(Int,sigma*N[d]), D)
    # lookup tables for 1d interpolation kernels
    windowHatInvLUT1d = Vector{Vector{T}}(undef,D)
    for d=1:D
        windowHatInvLUT1d[d] = [1.0/win_hat(k-1-N[d]/2, n[d], m, sigma) for k=1:N[d]]
    end
    # lookup table for multi-dimensional kernels
    windowHatInvLUT = zeros(Complex{T},N)
    if D==1
        windowHatInvLUT .= windowHatInvLUT1d[1]
    elseif D==2
        windowHatInvLUT .= reshape( kron(windowHatInvLUT1d[1], windowHatInvLUT1d[2]), N )
    elseif D==3
        windowHatInvLUT2d = kron(windowHatInvLUT1d[1], windowHatInvLUT1d[2])
        windowHatInvLUT .= reshape( kron(windowHatInvLUT2d, windowHatInvLUT1d[3]), N )
    else
        error("CuNFFT does not yet support $(D) dimensions")
    end
    return windowHatInvLUT
end

"""
precompute indices of the apodized image in the oversampled grid
"""
@generated function precomp_apodIdx(N::NTuple{D,Int64}, n::NTuple{D,Int64}) where D
    quote
        # linear indices of the oversampled grid
        linIdx = LinearIndices(n)
        # offsets to central NxN-region of the oversampled grid
        @nexprs $D d -> offset_d = round(Int, n[d] - N[d]/2) - 1
        # apodization indices
        apodIdx = zeros(Int64,N)
        @nloops $D l apodIdx d->(gidx_d = rem(l_d+offset_d, n[d]) + 1) begin
            (@nref $D apodIdx l) = (@nref $D linIdx gidx)
        end
        return vec(apodIdx)
    end
end


text:
module GA
import Base.isless

using GeneticAlgorithms
using Persa

using Memento

mutable struct Param{T}
    name::AbstractString
    values::Array{T}
end

mutable struct Gene <: Entity
    values::Array{Int}
    params::Array{Param}
    fitness
end

function Gene(params::Array{Param})
    values = Array{Int}(length(params))

    for i=1:length(values)
        values[i] = rand([1:length(params[i].values)...])
    end

    gene = Gene(values, params, +Inf)

    debug(logger, "Creating $gene")

    return gene
end

function Gene(parents::Array{Gene})
    child = Gene(parents[1].params)

    num_parents = length(parents)

    for i=1:num_parents
        debug(logger, "Parent $i: $(parents[i])")
    end

    for i=1:length(child.values)
        parent = (rand(UInt) % num_parents) + 1
        child.values[i] = parents[parent].values[i]
    end

    debug(logger, "Child: $child")

    return child
end

function Base.show(io::IO, gene::Gene)
    text = "Gene: fitness ($(gene.fitness)), values ("

    for i=1:length(gene.params)
        text = string(text, "($(gene.params[i].name), $(gene.params[i].values[gene.values[i]]))")
    end

    text = string(text, ")")

    show(io, text)
end

function mutate(gene::Gene)
    # let's go crazy and mutate mutation_rate% of the time
    rand(Float64) < (1 - g_mutation_ratio) && return

    i = rand(UInt) % length(gene.values) + 1
    gene.values[i] = rand([1:length(gene.params[i].values)...])

    return nothing
end

function create_entity(num)
    return Gene(parameters)
end

function group_entities(pop)
    global generation += 1
    global best_genes
    global dontchange
    global g_population

    if pop[1].fitness == 0
        return
    end

    if pop[1].fitness < best_genes[1].fitness
        dontchange = 0
    else
        dontchange += 1
    end

    push!(pop, best_genes...)

    sort!(pop; lt = (x, y) -> x.fitness < y.fitness)

    best_genes = pop[1:g_elite]

    info(logger, "BEST $generation/$g_generations: $(pop[1])")

    if generation >= g_generations
        return
    end

    if dontchange > g_patience
        return
    end

    if (time() - initialtime) > g_timelimit
        info(logger, "Time Limit")
        return
    end

    # simple naive groupings that pair the best entitiy with every other
    for i in 1:g_population
        ranks = rankroulette(length(pop))
        parent1 = rand(ranks)
        deleteat!(ranks, findin(ranks, [parent1]))
        parent2 = rand(ranks)

        debug(logger, "Crossover: $parent1, $parent2")
        produce([parent1, parent2])
    end
end

function rankroulette(total::Int)
    elements = Array{Int}(convert(Int, (total .* (total + 1) ./ 2)))
    element = total
    j = 1
    for i=1:total
        for w=1:element
            elements[j] = element
            j += 1
        end
        element -= 1
    end
    return elements
end

function isless(lhs::Gene, rhs::Gene)
    abs(lhs.fitness) > abs(rhs.fitness)
end

function crossover(group)
    return Gene([group...])
end

function fitness(ent)
    debug(logger, "Training a model")

    values = params(ent)

    ent.fitness = evalfun(estfun(values...))
end

function params(gene::Gene)
    values = Array{Any}(length(gene.params))

    for i=1:length(values)
        values[i] = gene.params[i].values[gene.values[i]]
    end

    return values
end

function run(dataset::Persa.CFDatasetAbstract,
                createmodel::Function,
                generations::Int,
                population::Int,
                mutation_ratio::Float64,
                elite::Int,
                k::Float64,
                verbose::Bool,
                patience::Int,
                timelimit::Real,
                params::Tuple{AbstractString, Any}...)

    if verbose
        global logger = Memento.config("debug"; fmt="[{date} | {level}]: {msg}")
    else
        global logger = Memento.config("info"; fmt="[{date} | {level}]: {msg}")
    end

    info(logger, "Starting Genetic Algorithm")

    # Global Parameters
    global g_generations = generations
    global g_mutation_ratio = mutation_ratio
    global g_elite = elite
    global g_population = population
    global g_timelimit = timelimit
    global initialtime = time()

    # Global Variables
    global generation = 0

    global g_patience = patience
    global dontchange = 0

    global parameters = Array{Param}(length(params))

    for i=1:length(params)
        parameters[i] = Param(params[i][1], params[i][2])
    end

    global best_genes = Array{Gene}(elite)

    for i=1:elite
        best_genes[i] = Gene(parameters)
    end

    ds_train, ds_val = Persa.get(Persa.HoldOut(dataset, k))

    global estfun(args...) = createmodel(ds_train, args...)
    global evalfun(model) = Persa.aval(model, ds_val).mae

    gamodel = runga(GA; initial_pop_size = population)

    values = Array{Any}(length(best_genes[1].params))

    for i=1:length(values)
        values[i] = best_genes[1].params[i].values[best_genes[1].values[i]]
    end

    return (createmodel(dataset, values...), values, best_genes[1].fitness)
end
end

function modeltunega(dataset::Persa.CFDatasetAbstract,
                createmodel::Function,
                params::Tuple{AbstractString, Any}...;
                generations::Int = 10,
                population::Int = 3,
                mutation_ratio::Float64 = 0.2,
                elite::Int = 2,
                k::Float64 = 0.9,
                verbose::Bool = true,
                patience::Int = 3,
                timelimit::Real = Inf
                )
    GA.run(dataset, createmodel,
            generations, population, mutation_ratio,
            elite, k, verbose, patience, timelimit, params...)
end

max_stars_repo_name: JuliaRecsys/ModelTuneCF

max_issues_repo_name: JuliaRecsys/ModelTuneCF

max_forks_repo_name: JuliaRecsys/ModelTuneCF


text:
using ACEbase: @def

@def aceimports begin

   import ACE

   import ACE:   PIBasis,
                 get_spec,
                 degree,
                 gensparse,
                 add_into_A!, add_into_A_dA!,
                 scaling,
                 AbstractACEModel,
                 AbstractProperty,
                 params, set_params!, nparams
end

@def baseimports begin
   import ACEbase

   import ACEbase: ACEBasis, ScalarACEBasis, OneParticleBasis, Discrete1pBasis,
                   valtype, gradtype,
                   evaluate, evaluate_d, evaluate_dd, evaluate_ed,
                   evaluate!, evaluate_d!, evaluate_dd!, evaluate_ed!,
                   precon!, _allfieldsequal,
                   AbstractState,
                   AbstractConfiguration,
                   AbstractContinuousState,
                   AbstractDiscreteState

   import ACEbase: FIO
   import ACEbase.FIO: read_dict, write_dict, save_json, load_json
end

# modules external to our own eco-system, rigorously separate using and import

@def extimports begin

   using Random: shuffle

   import Base: ==, length, kron, filter

   # for some reason it seems we need to export this to be used outside?!?!

   using LinearAlgebra: norm, dot, mul!, I

   using StaticArrays
end

# old stuff to delete or re-introduce ...
# , PIPotential, PIBasisFcn, site_alloc_B, site_alloc_dB,
# site_evaluate!, site_evaluate_d!,
# set_Aindices!, VecOrTup


text:
using Luxor, Test, Random
Random.seed!(42)

function point_on_polygon_test(fname)

    Drawing(1200, 1200, fname)
    origin()
    background("grey20")
    polys = [ngon(Point(x, y), rand(20:100), rand(3:8), rand() * 2π, vertices=true) for x in -500:200:500, y in -500:200:500]
    for pol in polys
        randomhue()
        poly(pol, :fill)

        pt1 = rand(1:length(pol))
        pt2 = mod1(pt1 + 1, length(pol))

        # all these should lie on their polygons
        pt = between(pol[pt1], pol[pt2], rand())
        flag = ispointonpoly(pt, pol)
        if flag
            sethue("black")
            circle(pt, 10, :fill)
        else
            sethue("red")
            circle(pt, 20, :fill)
        end

        # move the point, then it shouldn't
        pt = pt + (0.01, 0.01)
        flag = ispointonpoly(pt, pol)
        if flag
            sethue("red")
            circle(pt, 6, :fill)
        else
            sethue("green")
            circle(pt, 20, :stroke)
        end

    end

    @test finish() == true
    println("...finished test: output in $(fname)")
end

point_on_polygon_test("point-on-polygon.png")


text:
# DofMap
abstract DofMap

"""
Degree of freedom map for 1d discretizations

This type assumes that all elements have the same internal
discretization and same basis function. It also assumes that the
basis function can be split into:

 * interior modes that are zero on the boundary
 * boundary modes that are nonzero at a sinlge boundary node.
"""
type DofMap1d <: DofMap
    "Number of boundary modes"
    nb::Int

    "Number of boundary modes that are not on Dirichilet BCs"
    nbslv::Int

    "Number of Dirichilet boundary modes"
    nbdir::Int

    "Number of elements"
    nel::Int

    "Local numbering system map for each element (the same for all elements)"
    lmap::LocalNumSys1d

    "Boundary map. Maps local boundary modes into global modes"
    bmap::Array{Int,2}

    "Maps local modes into global modes"
    map::Array{Int,2}

    "Identifies Dirichilet modes"
    idir::Dict{Int,Vector{Int}}

    "Inner constructor. Should not usually be called"
    function DofMap1d(nel, nb, nbdir, lmap::LocalNumSys1d,
                      bmap::Array{Int,2}, idir::Dict{Int,Vector{Int}})
        nbslv = nb - nbdir
        nloc = nmodes(lmap)
        mp = zeros(Int, nloc, nel)
        for e = 1:nel
            mp[1,e] = bmap[1,e]
            mp[2,e] = bmap[2,e]
            for i = 1:ninterior(lmap)
                mp[2+i,e] = nb + (e-1)*ninterior(lmap) + i
            end
        end
        new(nb, nbslv, nbdir, nel, lmap, bmap, mp, idir)
    end
end

"Number of boundary modes in the discretization"
nbmodes(dof::DofMap) = dof.nb
"Number of boundary modes in the discretization"
nbndry(dof::DofMap) = dof.nb

"Number of boundary modes that are not Dirichilet BCs"
nbslvmodes(dof::DofMap) = dof.nbslv
"Number of elements"
num_elems(dof::DofMap) = dof.nel

"Total number of interior modes"
ninodes(dof::DofMap1d) = dof.nel * ninterior(dof.lmap)
"Total number of interior modes"
ninterior(dof::DofMap1d) = dof.nel * ninterior(dof.lmap)

"Total number of modes not in Dirichilet BCs (interior and boundary modes)"
nslvmodes(dof::DofMap1d) = nbslvmodes(dof) + ninodes(dof)

"Returns the local numbering system"
locmap(dof::DofMap1d, e=1) = dof.lmap

"Is the element located on a Dirichilet boundary condition"
hasdirbc(dof::DofMap1d, e) = haskey(dof.idir, e)

"Returns the Dirichilet modes for a given element"
idirbc(dof::DofMap1d, e) = dof.idir[e]

export DofMap1d
export nbmodes, nbslvmodes, num_elems, ninodes, nslvmodes, locmap, hasdirbc, idirbc

"""
Local to global degree of freedom mapping

This function creates a `DofMap1d` object given a local mapping scheme, a the
number of nodes and the index of nodes on Dirichilet boundary conditions. It
also allows for the case where the boundary condition is periodic, such that the first node
equals to the last node.
"""
function DofMap1d(lmap, nnodes, idir, iper=false)
    nbe = nbndry(lmap)
    nie = ninterior(lmap)

    nel = nnodes - 1
    bmap = zeros(Int, nbe, nel)

    for e = 1:nel
        bmap[1,e] = e
        bmap[2,e] = e+1
    end
    nb = nnodes
    nd = length(idir)
    if iper
        bmap[2,nel] = 1
        nb = nnodes - 1
        nd = 0
    elseif nd > 0 # There is a Dirichilet BC.
        if nd == 1
            if idir[1] == 1
                ii = [nnodes;  1:(nnodes-1);]
            else
                ii = [1:nnodes;]
            end
        else
            ii = [nnodes; 1:(nnodes-1);]
        end

        for e = 1:nel
            bmap[1,e] = ii[e]
            bmap[2,e] = ii[e+1]
        end

    end
    idir = Dict{Int,Vector{Int}}()
    nbslv = nb - nd

    ib = bndidx(lmap)
    for e = 1:nel
        if bmap[1,e] > nbslv && bmap[2,e] > nbslv
            idir[e] = [1,2]
        elseif bmap[1,e] > nbslv
            idir[e] =  [1]
        elseif bmap[2,e] > nbslv
            idir[e] = [2]
        end
    end
    return DofMap1d(nel, nb, nd, lmap, bmap, idir)
end

num_be(dof::DofMap, e) = nbndry(dof.lmap)
num_ie(dof::DofMap, e) = ninterior(dof.lmap)

bmap(dof, e) = view(dof.bmap, :, e)

function global2local(dof::DofMap, xg::Array{Float64,1})

    xloc = zeros(dof.nloc, dof.nel)
    for e = 1:dof.nel
        xloc[:,e] = xg[dof.map[:,e]]
    end
    return xloc
end


text:
"""
    GenericCallbacks

This module defines interfaces and wrappers for callbacks to be used with an
`AbstractODESolver`.

A callback `cb` defines three methods:

- `GenericCallbacks.init!(cb, solver, Q, param, t)`, to be called at solver
  initialization.

- `GenericCallbacks.call!(cb, solver, Q, param, t)`, to be called after each
  time step: the return value dictates what action should be taken:

   * `0` or `nothing`: continue time stepping as usual
   * `1`: stop time stepping after all callbacks have been executed
   * `2`: stop time stepping immediately

- `GenericCallbacks.fini!(cb, solver, Q, param, t)`, to be called at solver
  finish.

Additionally, _wrapper_ callbacks can be used to execute the callbacks under
certain conditions:

 - [`AtInit`](@ref)
 - [`AtInitAndFini`](@ref)
 - [`EveryXWallTimeSeconds`](@ref)
 - [`EveryXSimulationTime`](@ref)
 - [`EveryXSimulationSteps`](@ref)

For convenience, the following objects can also be used as callbacks:

- A `Function` object `f`, `init!` and `fini!` are no-ops, and `call!` will
  call `f()`, and ignore the return value.
- A `Tuple` object will call `init!`, `call!` and `fini!` on each element
  of the tuple.
"""
module GenericCallbacks

export AtInit,
    AtInitAndFini,
    EveryXWallTimeSeconds,
    EveryXSimulationTime,
    EveryXSimulationSteps

using MPI

init!(f::Function, solver, Q, param, t) = nothing
function call!(f::Function, solver, Q, param, t)
    f()
    return nothing
end
fini!(f::Function, solver, Q, param, t) = nothing

function init!(callbacks::Tuple, solver, Q, param, t)
    for cb in callbacks
        GenericCallbacks.init!(cb, solver, Q, param, t)
    end
end
function call!(callbacks::Tuple, solver, Q, param, t)
    val = 0
    for cb in callbacks
        val_i = GenericCallbacks.call!(cb, solver, Q, param, t)
        val_i = (val_i === nothing) ? 0 : val_i
        val = max(val, val_i)
        if val == 2
            return val
        end
    end
    return val
end
function fini!(callbacks::Tuple, solver, Q, param, t)
    for cb in callbacks
        GenericCallbacks.fini!(cb, solver, Q, param, t)
    end
end

abstract type AbstractCallback end

"""
    AtInit(callback) <: AbstractCallback

A wrapper callback to execute `callback` at initialization as well as
after each interval.
"""
struct AtInit <: AbstractCallback
    callback
end
function init!(cb::AtInit, solver, Q, param, t)
    init!(cb.callback, solver, Q, param, t)
    call!(cb.callback, solver, Q, param, t)
end
function call!(cb::AtInit, solver, Q, param, t)
    call!(cb.callback, solver, Q, param, t)
end
function fini!(cb::AtInit, solver, Q, param, t)
    fini!(cb.callback, solver, Q, param, t)
end

"""
    AtInitAndFini(callback) <: AbstractCallback

A wrapper callback to execute `callback` at initialization and at
finish as well as after each interval.
"""
struct AtInitAndFini <: AbstractCallback
    callback
end
function init!(cb::AtInitAndFini, solver, Q, param, t)
    init!(cb.callback, solver, Q, param, t)
    call!(cb.callback, solver, Q, param, t)
end
function call!(cb::AtInitAndFini, solver, Q, param, t)
    call!(cb.callback, solver, Q, param, t)
end
function fini!(cb::AtInitAndFini, solver, Q, param, t)
    call!(cb.callback, solver, Q, param, t)
    fini!(cb.callback, solver, Q, param, t)
end

"""
    EveryXWallTimeSeconds(callback, Δtime, mpicomm)

A wrapper callback to execute `callback` every `Δtime` wallclock time seconds.
`mpicomm` is used to syncronize runtime across MPI ranks.
"""
mutable struct EveryXWallTimeSeconds <: AbstractCallback
    "callback to wrap"
    callback
    "wall time seconds between callbacks"
    Δtime::Real
    "MPI communicator"
    mpicomm::MPI.Comm
    "time of the last callback"
    lastcbtime_ns::UInt64
    function EveryXWallTimeSeconds(callback, Δtime, mpicomm)
        lastcbtime_ns = zero(UInt64)
        new(callback, Δtime, mpicomm, lastcbtime_ns)
    end
end

function init!(cb::EveryXWallTimeSeconds, solver, Q, param, t)
    cb.lastcbtime_ns = time_ns()
    init!(cb.callback, solver, Q, param, t)
end
function call!(cb::EveryXWallTimeSeconds, solver, Q, param, t)
    # Check whether we should do a callback
    currtime_ns = time_ns()
    runtime = (currtime_ns - cb.lastcbtime_ns) * 1e-9
    runtime = MPI.Allreduce(runtime, max, cb.mpicomm)
    if runtime < cb.Δtime
        return 0
    else
        # Compute the next time to do a callback
        cb.lastcbtime_ns = currtime_ns
        return call!(cb.callback, solver, Q, param, t)
    end
end
function fini!(cb::EveryXWallTimeSeconds, solver, Q, param, t)
    fini!(cb.callback, solver, Q, param, t)
end

"""
    EveryXSimulationTime(f, Δtime)

A wrapper callback to execute `callback` every `time` simulation time seconds.
"""
mutable struct EveryXSimulationTime <: AbstractCallback
    "callback to wrap"
    callback
    "simulation time seconds between callbacks"
    Δtime::Real
    "time of the last callback"
    lastcbtime::Real
    function EveryXSimulationTime(callback, Δtime)
        new(callback, Δtime, 0)
    end
end

function init!(cb::EveryXSimulationTime, solver, Q, param, t)
    cb.lastcbtime = t
    init!(cb.callback, solver, Q, param, t)
end
function call!(cb::EveryXSimulationTime, solver, Q, param, t)
    # Check whether we should do a callback
    if (t - cb.lastcbtime) < cb.Δtime
        return 0
    else
        # Compute the next time to do a callback
        cb.lastcbtime = t
        return call!(cb.callback, solver, Q, param, t)
    end
end
function fini!(cb::EveryXSimulationTime, solver, Q, param, t)
    fini!(cb.callback, solver, Q, param, t)
end

"""
    EveryXSimulationSteps(callback, Δsteps)

A wrapper callback to execute `callback` every `nsteps` of the time stepper.
"""
mutable struct EveryXSimulationSteps <: AbstractCallback
    "callback to wrap"
    callback
    "number of steps between callbacks"
    Δsteps::Int
    "number of steps since last callback"
    steps::Int
    function EveryXSimulationSteps(callback, Δsteps)
        new(callback, Δsteps, 0)
    end
end

function init!(cb::EveryXSimulationSteps, solver, Q, param, t)
    cb.steps = 0
    init!(cb.callback, solver, Q, param, t)
end
function call!(cb::EveryXSimulationSteps, solver, Q, param, t)
    cb.steps += 1
    if cb.steps < cb.Δsteps
        return 0
    else
        cb.steps = 0
        return call!(cb.callback, solver, Q, param, t)
    end
end
function fini!(cb::EveryXSimulationSteps, solver, Q, param, t)
    fini!(cb.callback, solver, Q, param, t)
end

end

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
module ProjectEulerSolutions

    include("Problems.jl")
    include("digit.jl")
    include("factorization.jl")
    include("prime.jl")
    include("pyramid.jl")
    include("sequence.jl")

end

max_stars_repo_name: gnujosh/julia-euler-project

max_issues_repo_name: gnujosh/julia-euler-project

max_forks_repo_name: gnujosh/julia-euler-project