text:
struct MetaGraph{
    T<:Integer,
    Label,
    Graph,
    VertexMeta,
    EdgeMeta,
    GraphMeta,
    WeightFunction,
    U<:Real,
} <: AbstractGraph{T}
    graph::Graph
    vprops::Dict{Label,Tuple{T,VertexMeta}}
    eprops::Dict{Tuple{Label,Label},EdgeMeta}
    gprops::GraphMeta
    weightfunction::WeightFunction
    defaultweight::U
    metaindex::Dict{T,Label}
end

#= MetaGraph(
    g::Graph,
    vprops::Dict{Label,Tuple{T,VertexMeta}},
    eprops::Dict{Tuple{Label,Label},EdgeMeta},
    gprops::GraphMeta,
    weightfunction::WeightFunction,
    defaultweight::U,
    metaindex::Dict{T,Label}
) where {
    Label,
    T,
    Graph,
    VertexMeta,
    EdgeMeta,
    GraphMeta,
    WeightFunction,
    U,
} = MetaGraph{
    T,
    Label,
    Graph,
    VertexMeta,
    EdgeMeta,
    GraphMeta,
    WeightFunction,
    U
}(
    g,
    vprops,
    metaindex,
    eprops,
    gprops,
    weightfunction,
    defaultweight
)
=#

"""
    MetaGraph(g;
        Label = Symbol,
        VertexMeta = nothing,
        EdgeMeta = nothing,
        gprops = nothing,
        weightfunction = eprops -> 1.0,
        defaultweight = 1.0
    )

Construct a new meta graph based on `g`, where `Label` is the type of the vertex labels, `VertexMeta` is the type
of the metadata at a vertex, and `EdgeMeta` is the type of the metadata at an
edge. You can also attach arbitrary graph level metadata as `gprops`. It is recommended not
to set `Label` to `Int` to avoid confusion between vertex labels and codes.

```jldoctest example
julia> using LightGraphs

julia> using MetaGraphsNext

julia> colors = MetaGraph(Graph(), VertexMeta = String, EdgeMeta = Symbol, gprops = "special")
Meta graph based on a {0, 0} undirected simple Int64 graph with vertices indexed by Symbol(s), String(s) vertex metadata, Symbol(s) edge metadata, "special" as graph metadata, and default weight 1.0
```

Use `setindex!` to add a new vertex with the given metadata. If a vertex with the given
index does not exist, it will be created automatically; otherwise, `setindex!` will modify
the metadata for the existing vertex.

```jldoctest example
julia> colors[:red] = "warm";

julia> colors[:yellow] = "warm";

julia> colors[:blue] = "cool";
```

You can access and change the metadata using indexing: zero arguments for graph metadata,
one label for vertex metadata, and two labels for edge metadata.

```jldoctest example
julia> colors[]
"special"

julia> colors[:blue] = "very cool";

julia> colors[:blue]
"very cool"

julia> colors[:red, :yellow] = :orange;

julia> colors[:red, :yellow]
:orange

julia> haskey(colors, :red, :yellow)
true

julia> haskey(colors, :yellow, :red) # undirected graph, so vertex order doesn't matter
true
```

You can delete vertices and edges with `delete!`.

```jldoctest example
julia> delete!(colors, :red, :yellow);

julia> delete!(colors, :blue);
```

You can use the `weightfunction` keyword to specify a function which will
transform vertex metadata into a weight. This weight must always be the same
type as the `defaultweight`.

```jldoctest example
julia> weighted = MetaGraph(Graph(), EdgeMeta = Float64, weightfunction = identity);

julia> weighted[:red] = nothing; weighted[:blue] = nothing; weighted[:yellow] = nothing;

julia> weighted[:red, :blue] = 1.0; weighted[:blue, :yellow] = 2.0;

julia> the_weights = LightGraphs.weights(weighted)
metaweights

julia> size(the_weights)
(3, 3)

julia> the_weights[1, 3]
1.0

julia> diameter(weighted)
3.0
```

MetaGraphs inherit many methods from LightGraphs. In general, inherited methods refer to
vertices by codes, not labels, for compatibility with AbstractGraph. Vertex codes get
reassigned after `rem_vertex!`, so I recommend using label indexing if possible.

```jldoctest example
julia> is_directed(colors)
false

julia> nv(zero(colors))
0

julia> ne(copy(colors))
0

julia> add_vertex!(colors, :white, "neutral")
true

julia> add_edge!(colors, 1, 3, :pink)
true

julia> rem_edge!(colors, 1, 3)
true

julia> rem_vertex!(colors, 3)
true

julia> rem_vertex!(colors, 3)
false

julia> eltype(colors) == Int
true

julia> edgetype(colors) == Edge{Int}
true

julia> vertices(colors)
Base.OneTo(2)

julia> has_edge(colors, 1, 2)
false

julia> has_vertex(colors, 1)
true

julia> LightGraphs.SimpleGraphs.fadj(colors, 1) == Int[]
true

julia> LightGraphs.SimpleGraphs.badj(colors, 1) == Int[]
true

julia> colors == colors
true

julia> issubset(colors, colors)
true

julia> SimpleGraph(colors)
{2, 0} undirected simple Int64 graph
```

You can seemlessly make MetaGraphs based on DiGraphs as well.

```jldoctest example
julia> rock_paper_scissors = MetaGraph(DiGraph(), Label = Symbol, EdgeMeta = Symbol);

julia> rock_paper_scissors[:rock] = nothing; rock_paper_scissors[:paper] = nothing; rock_paper_scissors[:scissors] = nothing;

julia> rock_paper_scissors[:rock, :scissors] = :rock_beats_scissors; rock_paper_scissors[:scissors, :paper] = :scissors_beats_paper; rock_paper_scissors[:paper, :rock] = :paper_beats_rock;

julia> is_directed(rock_paper_scissors)
true

julia> haskey(rock_paper_scissors, :scissors, :rock)
false

julia> haskey(reverse(rock_paper_scissors), :scissors, :rock)
true

julia> SimpleDiGraph(rock_paper_scissors)
{3, 3} directed simple Int64 graph

julia> sub_graph, _ = induced_subgraph(rock_paper_scissors, [1, 3]);

julia> haskey(sub_graph, :rock, :scissors)
true

julia> delete!(rock_paper_scissors, :paper);

julia> rock_paper_scissors[:rock, :scissors]
:rock_beats_scissors
```
"""
function MetaGraph(
    g::AbstractGraph{T};
    Label = Symbol,
    VertexMeta = Nothing,
    EdgeMeta = Nothing,
    gprops = nothing,
    weightfunction = eprops -> 1.0,
    defaultweight = 1.0,
) where {Vertex,T}
    MetaGraph(
        g,
        Dict{Label,Tuple{T,VertexMeta}}(),
        Dict{Tuple{Label,Label},EdgeMeta}(),
        gprops,
        weightfunction,
        defaultweight,
        Dict{T,Label}(),
    )
end

const MetaUndirectedGraph = MetaGraph{<:Any,<:Any,<:SimpleGraph}

SimpleGraph(g::MetaUndirectedGraph) = g.graph

is_directed(::Type{<:MetaUndirectedGraph}) = false

function arrange(g::MetaUndirectedGraph, label_1, label_2)
    vprops = g.vprops
    arrange(g::MetaUndirectedGraph, label_1, label_2, vprops[label_1], vprops[label_2])
end

function arrange(g::MetaUndirectedGraph, label_1, label_2, u, v)
    if u > v
        (label_2, label_1)
    else
        (label_1, label_2)
    end
end

zero(
    g::MetaGraph{T,Label,Graph,VertexMeta,EdgeMeta,GraphMeta},
) where {T,Label,Graph,VertexMeta,EdgeMeta,GraphMeta} = MetaGraph(
    Graph();
    Label = Label,
    VertexMeta = VertexMeta,
    EdgeMeta = EdgeMeta,
    gprops = g.gprops,
    weightfunction = g.weightfunction,
    defaultweight = g.defaultweight,
)


text:
using FileIO: FileIO

# Required for compatibility with pre 1.6 and post 1.6
# drop in the future if the minimum requirement is bumped
if isdefined(FileIO, :action)
    # FileIO >= 1.6
    _format_stream(format, io) = FileIO.Stream{format}(io)
else
    _format_stream(format, io) = FileIO.Stream(format, io)
end

ext: jl


text:
function LinSolver(k, l, U)
    e = [[1., 0.] [0., 1.]]
    function unitStrain(k,l)
        return e[k,:]*e[l,:]'
    end

    Chia, Chib = unitStrain(k,l)*L*e[1,:], unitStrain(k,l)*L*e[2,:]

    Chicorner = vcat(zeros(2, 1), Chia, Chia+Chib, Chib)
    Chiinternal = zeros(Float64, 2*Ns[2])
    Chimaster = zeros(Float64, 2*Ns[3])
    Chislave = Chimaster + D*Chicorner
    Chi = vcat(Chicorner, Chiinternal, Chimaster, Chislave)

    function Assembly(U::Matrix{Float64})
        gps = gaussPoints2D(2)
        K, Fint = zeros(Float64, Ns[6], Ns[6]), zeros(Float64, Ns[6], 1)
        for i = 1:Ns[5]
            Ce = vec(DOFs[i, :, :]')
            elK, elFint = zeros(Float64, 8, 8), zeros(Float64, 8, 1)
            for gp in eachrow(gps)
                g, h = gp[2:end]
                B, detJ, PK1_1D, A_2D = constitutiveEqs(U[Ce, :], elCoords[i, :, :], g, h, E, nu)
                elK += gp[1]*(B'*A_2D*B)*detJ
            end
            K[Ce, Ce] += elK
        end
        return sparse(K)
    end

    function solve(K::SparseMatrixCSC, Chicorner)
        Kic = K[2*Ns[1]+1:2*Ns[7], 1:2*Ns[1]]
        Kii = K[2*Ns[1]+1:2*Ns[7], 2*Ns[1]+1:2*Ns[7]]
        Kim = K[2*Ns[1]+1:2*Ns[7], 2*Ns[7]+1:2*Ns[9]]
        Kis = K[2*Ns[1]+1:2*Ns[7], 2*Ns[9]+1:2*Ns[11]]

        Kmc = K[2*Ns[7]+1:2*Ns[9], 1:2*Ns[1]]
        Kmi = K[2*Ns[7]+1:2*Ns[9], 2*Ns[1]+1:2*Ns[7]]
        Kmm = K[2*Ns[7]+1:2*Ns[9], 2*Ns[7]+1:2*Ns[9]]
        Kms = K[2*Ns[7]+1:2*Ns[9], 2*Ns[9]+1:2*Ns[11]]

        Ksc = K[2*Ns[9]+1:2*Ns[11], 1:2*Ns[1]]
        Ksi = K[2*Ns[9]+1:2*Ns[11], 2*Ns[1]+1:2*Ns[7]]
        Ksm = K[2*Ns[9]+1:2*Ns[11], 2*Ns[7]+1:2*Ns[9]]
        Kss = K[2*Ns[9]+1:2*Ns[11], 2*Ns[9]+1:2*Ns[11]]

        Qi = -(Kic+Kis*D)*Chicorner
        Qm = -(Kmc+Kms*D+Ksc+Kss*D)*Chicorner
        Kupper = hcat(Kii, Kim + Kis)
        Klower = hcat(Kmi + Ksi, Kmm + Kms + Ksm + Kss)

        A = vcat(Kupper, Klower)
        Q = vcat(Qi, Qm)

        sol = A \ Q

        return sol[1:2*Ns[2],:], sol[2*Ns[2]+1:end,:]
    end

    K = Assembly(U)
    Chiinternal, Chimaster = solve(K, Chicorner)
    Chislave = Chimaster + D*Chicorner
    Chi = vcat(Chicorner, Chiinternal, Chimaster, Chislave)

    return Chi
end

max_stars_repo_name: likemaestro/Nonlinear-Homogenization

max_issues_repo_name: likemaestro/Nonlinear-Homogenization

max_forks_repo_name: likemaestro/Nonlinear-Homogenization


text:
@testset "AlgAssAbsOrd" begin
  include("AlgAssAbsOrd/Order.jl")
  include("AlgAssAbsOrd/Ideal.jl")
  include("AlgAssAbsOrd/PicardGroup.jl")
  include("AlgAssAbsOrd/LocallyFreeClassGroup.jl")
  include("AlgAssAbsOrd/ICM.jl")
end


text:
module PoreMatMod_Test

using Test, Graphs, PoreMatMod, LinearAlgebra

NiPyC_fragment_trouble = Crystal("NiPyC_fragment_trouble.cif")
irmof1 = Crystal("IRMOF-1.cif")
PyC = moiety("PyC.xyz")
p_phenylene = moiety("p-phenylene.xyz")
tagged_p_phenylene = moiety("2-!-p-phenylene.xyz")
acetamido_p_phen = moiety("2-acetylamido-p-phenylene.xyz")

@testset "replacement split across PB" begin
    parent = deepcopy(NiPyC_fragment_trouble)
    query = deepcopy(PyC)
    replacement = moiety("PyC-CH3.xyz")

    @test_throws AssertionError replace(parent, query => replacement)
    infer_bonds!(parent, true)
    child = replace(parent, query => replacement)

    @test ne(child.bonds) == ne(parent.bonds) + 3 * 2 # added H -> CH3 on two PyC ligands

    # test conglomerate worked.
    remove_bonds!(child)
    infer_bonds!(child, true)
    @test ne(child.bonds) == ne(parent.bonds) + 3 * 2 # added H -> CH3 on two PyC ligands
end

@testset "split across PB, non-connected replacement" begin
    parent = Crystal("NiPyC_fragment_trouble.cif", convert_to_p1=false)
    infer_bonds!(parent, true)

    query = moiety("PyC_split.xyz")
    replacement = moiety("PyC_split_replacement.xyz")
    child = replace(parent, query => replacement)

    # ensure O atoms aligned
    parent_Os = parent[parent.atoms.species .== :O]
    child_Os  = child[child.atoms.species .== :O]
    nb_Os = parent_Os.atoms.n
    @test parent_Os.atoms.n == nb_Os

    min_distances = [Inf for _ = 1:nb_Os] # from parent
    for i = 1:nb_Os
        for j = 1:nb_Os
            r = norm(parent_Os.atoms.coords.xf[:, i] - child_Os.atoms.coords.xf[:, j])
            if r < min_distances[i]
                min_distances[i] = r
            end
        end
    end
    @test all(min_distances .< 0.01)
end

@testset "replacement spans twice the unit cell" begin
    parent = deepcopy(NiPyC_fragment_trouble)
    infer_bonds!(parent, true)

    query = deepcopy(PyC)
    replacement = moiety("PyC-long_chain.xyz")

    @test_throws Exception replace(parent, query => replacement)
end

@testset "substructure_search" begin
    xtal = deepcopy(irmof1)
    strip_numbers_from_atom_labels!(xtal)
    infer_bonds!(xtal, true)
    timil125 = Crystal("Ti-MIL-125.cif")
    strip_numbers_from_atom_labels!(timil125)
    infer_bonds!(timil125, true)
    query = deepcopy(p_phenylene)
    p_phenylene_w_R_grp = deepcopy(tagged_p_phenylene)
    search1 = query ∈ xtal

    @test nb_isomorphisms(search1) == 96

    @test nb_locations(search1) == 24

    @test nb_ori_at_loc(search1)[1] == 4

    @test search1.isomorphisms[1][1] == Dict([q => p for (q, p) in enumerate([233, 306, 318, 245, 185, 197, 414, 329, 402, 341])])

    search2 = query ∈ timil125

    @test nb_isomorphisms(search2) == 48

    @test nb_locations(search2) == 12

    @test nb_ori_at_loc(search2)[1] == 4

    @test search2.isomorphisms[1][1] == Dict([q => p for (q, p) in enumerate([8, 140, 144, 141, 7, 133, 186, 185, 190, 189])])

    search3 = p_phenylene_w_R_grp ∈ timil125

    @test nb_isomorphisms(search3) == 48

    @test nb_locations(search3) == 12

    @test nb_ori_at_loc(search3)[1] == 4

    @test search3.isomorphisms[1][1] == Dict([q => p for (q, p) in enumerate([8, 140, 144, 141, 7, 133, 185, 190, 189, 186])])

    query = moiety("!-S-bromochlorofluoromethane.xyz")
    parent = moiety("S-bromochlorofluoromethane.xyz")
    search = query ∈ parent

    @test search.isomorphisms[1][1] == Dict([q => p for (q, p) in enumerate([1, 2, 3, 4, 5])])

    @test [parent.atoms.species[search.isomorphisms[1][1][i]] for i in 1:4] == query.atoms.species[1:4] &&
        query.atoms.species[5] == :H! && parent.atoms.species[search.isomorphisms[1][1][5]] == :H

    query = deepcopy(p_phenylene)
    parent = Crystal("IRMOF-1_one_ring.cif")
    strip_numbers_from_atom_labels!(parent)
    infer_bonds!(parent, true)
    search = query ∈ parent

    @test search.isomorphisms[1][1] == Dict([q => p for (q, p) in enumerate([34, 38, 39, 36, 26, 27, 51, 46, 50, 48])])
end # test set: substructure_search

@testset "find_replace" begin
    parent = deepcopy(irmof1)
    strip_numbers_from_atom_labels!(parent)
    infer_bonds!(parent, true)
    query = deepcopy(tagged_p_phenylene)
    replacement = deepcopy(acetamido_p_phen)
    new_xtal = replace(parent, query => replacement)

    @test new_xtal.atoms.n == 592

    new_xtal = replace(parent, query => replacement, nb_loc=1)

    @test new_xtal.atoms.n == 431

    new_xtal = replace(parent, query => replacement, loc=[2,3])

    @test new_xtal.atoms.n == 438

    new_xtal = replace(parent, query => replacement, loc=[2,3,4], ori=[1,1,1])

    @test new_xtal.atoms.n == 445

    replacement = deepcopy(p_phenylene)
    new_xtal = replace(parent, query => replacement)

    @test ne(new_xtal.bonds) == ne(parent.bonds)

    replacement = deepcopy(acetamido_p_phen)
    new_xtal = replace(parent, query => replacement, nb_loc=1)

    @test ne(new_xtal.bonds) == (ne(parent.bonds) - ne(query.bonds) + ne(replacement.bonds))

    xtal = deepcopy(irmof1)
    strip_numbers_from_atom_labels!(xtal)
    infer_bonds!(xtal, true)
    query = deepcopy(tagged_p_phenylene)
    nb_bonds(xtal) = ne(xtal.bonds)
    # test that a "no-op" leaves the number of bonds unchanged
    replacement = deepcopy(p_phenylene)

    @test nb_bonds(replace(xtal, query => replacement)) == nb_bonds(xtal)

    # test that adding a new moiety increases the number of bonds correctly
    replacement = deepcopy(acetamido_p_phen)

    @test ne((replace(xtal, query => replacement, nb_loc=1)).bonds) ==
        (ne(xtal.bonds) - ne(query.bonds) + ne(replacement.bonds))
end

@testset "conglomerate test" begin
    xtal = Crystal("conglomerate_test.cif")
    infer_bonds!(xtal, false)
    @test ne(xtal.bonds) == 1
    remove_bonds!(xtal)
    @test ne(xtal.bonds) == 0
    infer_bonds!(xtal, true)
    @test ne(xtal.bonds) == 5
    PoreMatMod.conglomerate!(xtal)
    remove_bonds!(xtal)
    translate_by!(xtal.atoms.coords, Frac([0.5, 0.5, 0.5]))
    infer_bonds!(xtal, false)
    @test ne(xtal.bonds) == 5
end

@testset "remove duplicates" begin
    parent = moiety("ADC.xyz")
    new_box = replicate(unit_cube(), (10, 10, 10))
    new_atoms = Frac(Cart(parent.atoms, parent.box), new_box)
    new_charges = Frac(Cart(parent.charges, parent.box), new_box)
    parent = Crystal(parent.name, new_box, new_atoms, new_charges)
    infer_bonds!(parent, false)
    query = moiety("naphthyl_fragment.xyz")
    replacement = moiety("F_naphthyl_fragment.xyz")
    child = replace(parent, query => replacement, nb_loc=2, remove_duplicates=true, reinfer_bonds=true)

    @test child.atoms.n == parent.atoms.n

    @test nv(child.bonds) == nv(parent.bonds) && ne(child.bonds) == ne(parent.bonds)
end

end


text:
import .SpecialFunctions
import .SpecialFunctions: erf, erfc, gamma, lgamma,
    besselj0, besselj1, besselj, bessely0, bessely1, bessely

include("agm.jl")
include("bessel.jl")
include("elliptic.jl")
include("gamma_erf.jl")


text:
"""
    dp_evaluate_policy_textbook!(V, Q, fmdp, policy, γ; tol, maxiter)

Evaluates (deterministic or stochastic) decision policy `policy` for the given finite Markov
decision process `fmdp` and given discount factor `γ`.

The procedure is 'textbook', in the sense that unnecessary allocations and steps have
been made in order to make the procedure more explicit and readable.

By saying that a policy has been "evaluated", we imply that state values `V` and
state-action values `Q` are computed for this policy. This function evaluates policies in-place,
approximately, in an iterrative manner, until maximal difference between state values in consequtive
iterations is less than the specified `tol`, or until a maximal number of iterations `maxiter`
is performed.

The function returns a 2-tuple containing: the number of iterations performed, and a boolean
indicator of convergence (which is `true` if the vector of state values is evaluated with the
required precission, and `false` otherwise).
"""
function dp_evaluate_policy_textbook!(V::AbstractVector{<:Real}, Q::AbstractMatrix{<:Real}, fmdp::FiniteMDP, policy, γ::Real; tol = 1e-2, maxiter = 100)
    for i = 1:maxiter
        Q_from_V!(Q, V, fmdp, γ)
        Δ = V_from_Q!(V, Q, policy)
        if Δ < tol
            return i, true
        end
    end # for: iterations
    return maxiter, false
end

"""
    dp_update_V!(V, fmdp, [policy,] γ::Real)

Performs one-step in-place update of the vector of state values `V` for the given finite Markov Decision Process
`fmdp` (deterministic or stochastic) and given discount factor `γ`.

If decision policy `policy` is given, the update will be perform for that particular policy. Otherwise, the
update will be performed for the deterministic greedy policy with respect to the current state values.

Return the maximal absolute update of any element of `V`.
"""
function dp_update_V!(V::AbstractVector{<:Real}, fmdp::FiniteMDP, policy, γ::Real)
    Δ = -Inf
    for s = 1:states_no(fmdp)
        new_value = expected_return(fmdp, s, policy, V, γ)
        Δ = max(Δ, abs(V[s] - new_value))
        V[s] = new_value
    end # for: states
    return Δ
end

function dp_update_V!(V::AbstractVector{<:Real}, fmdp::FiniteMDP, γ::Real)
    Δ = -Inf
    for s = 1:states_no(fmdp)
        new_value = -Inf
        for a = 1:actions_no(fmdp)
            new_value = max(new_value, expected_return(fmdp, s, a, V, γ))
        end # for: actions
        Δ = max(Δ, abs(V[s] - new_value))
        V[s] = new_value
    end # for: states
    return Δ
end


text:
using FewShotAnomalyDetection
using Flux
using MLDataPattern
using FluxExtensions
using Adapt
using DataFrames
using CSV

import FewShotAnomalyDetection: loss, zparams
include("experimentalutils.jl")
include("vae.jl")

inputDim = 1
hiddenDim = 100
latentDim = 1
numLayers = 2
nonlinearity = "relu"
layerType = "Dense"

# Create data
T = Float64
train = T[10000]
train = randn(T, 1, 5000) .- 1
train = randn(T, 1, 5000) .+ 1

batchSize = 100
numBatches = 10000

encoder = Adapt.adapt(T, FluxExtensions.layerbuilder(inputDim, hiddenDim, latentDim * 2, numLayers, nonlinearity, "", layerType))
decoder = Adapt.adapt(T, FluxExtensions.layerbuilder(latentDim, hiddenDim, inputDim + 1, numLayers + 1, nonlinearity, "linear", layerType))
outerVAE = VAE(encoder, decoder, T(1), :scalarsigma)
opt = Flux.Optimise.ADAM(3e-4)
cb = Flux.throttle(() -> println("$datasetName outer VAE: $(printingloss(outerVAE, train[1]))"), 5)
Flux.train!((z, y) -> loss(outerVAE, z), Flux.params(outerVAE), RandomBatches((train, zero(train)), batchSize, numBatches), opt, cb = cb)


text:
using Revise, DiscreteEventsLite

f(a; k) = println(a, k)
g(scheduler, message) = register!(scheduler, f, 5.0, message; k=2)

scheduler = Scheduler(trace=true, store=true)
register!(scheduler, f, after, .99, "hi "; k=1, description="some event")
register!(scheduler, g, at, 2.0, scheduler, "I'm done "; id = "1")
register!(scheduler, stop!, at, 10.5, scheduler)
register!(scheduler, ()->(), every, 1.0; description="repeating")

run!(scheduler, 11)

replay_events(scheduler)


text:
"""
**SeisCopy**

*Copy a seis file*

**IN**
* in: input filename
* out: output filename

**OUT**

*Credits: AS, 2015*

"""
function SeisCopy(in::String,out::String)

	filename_data_in = ParseDataName(in)
	filename_headers_in = ParseHeaderName(in)
	DATAPATH = get(ENV,"DATAPATH",join([pwd(),"/"]))
	filename_data_out = join([DATAPATH out "@data@"])
	if filename_headers_in != "NULL"
		filename_headers_out = join([DATAPATH out "@headers@"])
	else
		filename_headers_out = "NULL"
	end
	extent = ReadTextHeader(in)
	WriteTextHeader(out,extent,"native_float",4,filename_data_out,filename_headers_out)
	run(`cp $filename_data_in $filename_data_out`)
	if filename_headers_in != "NULL"
		run(`cp $filename_headers_in $filename_headers_out`)
	end

end

function SeisCopy(in::Array{String,1},out::Array{String,1})

	for j = 1 : length(in)
		SeisCopy(in[j],out[j])
	end

end


text:
module SavitskyGolayFilters

import DSP: conv

vandermonde(x, k) = x .^ transpose(0:k)
fit_matrix(n, k) = vandermonde(-floor(n/2):floor(n/2), k)

"""
    savitskygolayfilter(n,k,d=0)

Compute the Savitsky-Golay filter of length `n`, degree `k` and derivative
order `d`.

See also: [`savitskygolay`](@ref)
"""
savitskygolayfilter(n,k,d=0) = begin
    M = fit_matrix(n, k)
    B = inv(transpose(M)*M)*transpose(M)
    B[d+1,end:-1:1]
end

"""
    savitskygolay(x, n, k, d=0)
    savitskygolay(x, filter)

Apply a Savitsky-Golay filter on `x`. The applied filter is either the already
calculated `filter` or a newly calculated filter of length `n`, degree `k`
and derivative order `d`.

# Examples
```jldoctest
julia> using Random

julia> x = -3π:0.1:3π;

julia> y = sin.(x) .+ randn(length(x));

julia> filtered = savitskygolay(y, 10, 2);

julia> size(filtered)
(158,)

julia> size(y)
(158,)

julia> filtered_derivative = savitskygolay(y, 10, 2, 1);

julia> size(filtered_derivative)
(158,)

```

See also: [`savitskygolayfilter`](@ref)
"""
function savitskygolay end

savitskygolay(x, n, k, d=0) = savitskygolay(x, savitskygolayfilter(n,k,d))

savitskygolay(x, filter::AbstractArray) = begin
    len_filter = length(filter)
    border_repeat = floor(Int, len_filter/2)
    conv(
         reverse(filter),
        [repeat([x[1]],border_repeat);x;repeat([x[end]],border_repeat)],
       )[len_filter:end-len_filter+1]
end

export savitskygolayfilter, savitskygolay

end


text:
using CalculusWithJulia

# hacky build script
CalculusWithJulia.WeaveSupport.weave_all(force=false, build_list=(:html,))
builddir = joinpath(@__DIR__, "build")
htmldir = joinpath(@__DIR__, "..", "html")
isdir(builddir) || mkdir(builddir)
for d in readdir(htmldir)
    dir = joinpath(htmldir, d)
    if isdir(dir)
        cp(dir, joinpath(builddir, basename(dir)), force=true)
    end
end
toc = joinpath(htmldir, "misc", "toc.html")
cp(toc, joinpath(builddir,"index.html"), force=true)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.

import Documenter
#Documenter.deploydocs(repo = "github.com/CalculusWithJulia/CalculusWithJulia.jl.git")
#=deploydocs(
    repo = "<repository url>"
)=#


text:
"""
    AbstractEpidemicModel

An abstract type representing an epidemic model.
"""
abstract type AbstractEpidemicModel end

"""
    init_state(em, t)

Setup the initial state for the Gillespie simulation.
"""
function init_state(em::AbstractEpidemicModel, state_0)
    return copy(state_0)
end

"""
    rate(em::AbstractEpidemicModel, k, state)

Return the reaction rate for reaction `k` for the epidemic model `em` at state `state`.
"""
function rate end

"""
    init_rates(em, state)

Compute the rates vector `a` for the state `state`.
"""
function init_rates end

"""
    init_output(em, state, n)

Initialize the output of the algorithm, allocating it and processing the initial `state`.
"""
function init_output(em::AbstractEpidemicModel, state::T, nmax) where T
    output = Vector{T}(undef, nmax)
    output[1] = copy(state)
    return output
end

"""
    update_state!(em, state, n)

Update the `state` of the simulation and the rates vector `a` with reaction `k` from the rates vector.
"""
function update_state_and_rates! end

"""
    update_output!(output, state, n, k, em)

Update the `output` of the simulation.

Arguments:
* `output`: the output array. `output[n]` gets mutated.
* `state`: the current state of the simulation
* `n`: the current time step.
* `k`: ...
"""
function update_output!(output, state, n, k, em::AbstractEpidemicModel)
    # NB. This is inefficient, so you'll want to overload it in most cases
    output[n] = copy(state)
end

"""
    finalize_output(em, output)

Finalize the output.

Overwrite this to massage your output for easier plotting.
"""
finalize_output(em::AbstractEpidemicModel, output) = output


text: compute(n::Int)::Int = (((4n + 3)n + 8)n - 9) ÷ 6

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
#module juleva

export @safe_cfunction

macro safe_cfunction(f, rt, args)
  return esc(:(@cfunction($f, $rt, $args), $rt, [$(args.args...)]))
end

# Add the C library to PATH
const mylib=joinpath(pwd(),"../C","libgen.so")
# Arguments
argv=["-c sc"]

fpi = convert(Float64,π)
# Parameters
x = Array{Float64,1}(undef,3)
y = Array{Float64,1}(undef,3)

xs = [0.; 0.; 0.]
xl = [0.; 0.; 0.]
xu = [fpi; fpi; fpi]

mutable struct GOpts
    algo :: Cint
    dim  :: Cint
    x_s  :: Ptr{Cdouble}
    x_l  :: Ptr{Cdouble}
    x_u  :: Ptr{Cdouble}
    iter :: Cint
    retRes :: Ptr{Cdouble}
    retErr :: Ptr{Cdouble}
    func  :: Ptr{Cvoid}
    GOpts() = new(-1,-1,C_NULL,C_NULL,C_NULL,-1,C_NULL,C_NULL,C_NULL)
end

function integrand(x::Ptr{Float64},dim::Ptr{Float64},params::Ptr{Float64})
    A = 1.0 / (pi * pi * pi)
    return 3*A / (1.0 - cos(unsafe_load(x,1))*cos(unsafe_load(x,2))*cos(unsafe_load(x,3)))::Cdouble
end

function get_ptr(f::Function)
    integrand_c = @cfunction( $f , Cdouble,(Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble}))
#    cf =  Base.cconvert(Ptr{Cvoid}, integrand_c)
#    GC.@preserve cf begin
#        fptr = Base.unsafe_convert(Ptr{Cvoid}, cf)
#    end
    #    return fptr
    return integrand_c
end

# Init options
g_opts = GOpts()
g_opts.algo = 2
g_opts.dim = 3
println(g_opts.dim)
g_opts.iter = 10
g_opts.retRes = pointer(x)
g_opts.retErr = pointer(y)
g_opts.x_s = pointer(xs)
g_opts.x_l = pointer(xl)
g_opts.x_u = pointer(xu)

function TestCallBackFunc(f::Ptr{Cvoid}, o::GOpts)
    g_opts.func = f
    dump(g_opts)
    return ccall((:g_optimize2,mylib),Int32,(Int32,Ptr{Ptr{UInt8}},Ref{GOpts}),length(argv), argv, Ref(g_opts))
end

function f end

function g(func::Function, o::GOpts)
    f = func
    return TestCallBackFunc(
        @cfunction(f,Cdouble,(Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble}))
        , o)
end

function f(x::Ptr{Float64},dim::Ptr{Float64},params::Ptr{Float64})
    A = 1.0 / (pi * pi * pi)
    return 3*A / (1.0 - cos(unsafe_load(x,1))*cos(unsafe_load(x,2))*cos(unsafe_load(x,3)))::Cdouble
end

g(f, g_opts)

#test(g_opts)

max_stars_repo_name: denisbertini/juleva

max_issues_repo_name: denisbertini/juleva

max_forks_repo_name: denisbertini/juleva


text:
function very_inefficient_solve(input)
	ret = 0
	conditions = ["shiny gold"]
	bags = [split(x, r" bags?", limit = 2) for x in input]
	i = 1
	while i <= length(bags)
		if any([contains(bags[i][2], condition) for condition in conditions])
			push!(conditions, bags[i][1])
			deleteat!(bags, i)
			ret += 1
			i = 1
		else
			i += 1
		end
	end
	ret
end

function solve(input)
	bags = split.(input, r" bags?")
	rules = Dict([ bag[1] => [ !contains(rule, "no other") ? match(r"\d [ \w]*", rule).match : nothing for rule in bag[2: end - 1] ] for bag in bags ])
	ret = 0
	queue = [("shiny gold", 1)]
	while (!isempty(queue))
		(bag, count) = pop!(queue)
		if !isnothing(rules[bag][1])
			children = split.(rules[bag], ' ', limit = 2)
			for child in children
				n = parse(Int64, child[1])
				push!(queue, (child[2], n * count))
				ret += n * count
			end
		end
	end
	ret
end

@assert very_inefficient_solve(readlines("input/test/test_day7.txt")) == 4 ["test case 1 failed"]
@assert solve(readlines("input/test/test_day7.txt")) == 32 ["test case 2 failed"]

println("First star: " * string(@time very_inefficient_solve(readlines("input/day7.txt"))))
println("Second star: " * string(@time solve(readlines("input/day7.txt"))))


text:
function angle_of_incidence(lat,lon,jd,time)
# angle_of_incidence calculates the cosine of the sun incident angle
#  y = ANGLE_OF_INCIDENCE(lat,lon,jd,time) calculates the cosine of the sun incident angle for
#   a given location (lat,lon), julian day (jd) and time in hours and fraction of hours.

# inclination angle
delta=(asin(-sin(23.45/180*pi)*cos(360/365.25*(jd+10)/180*pi)));
# lat to radians
lamda=lat/180*pi;
# time to radians
omega=time/24*2*pi+lon/360*2*pi+pi;

y=cos(delta)*cos(lamda)*cos(omega)+sin(delta)*sin(lamda);

return y

end


text:
"""

    calc_coeffs(X, Y, Z, XTX, ZTZ)

Calculates the the coefficient estimates

# Arguments

- X = 2d array of floats consisting of the row covariates, with all
  categorical variables coded in appropriate contrasts
- Y = 2d array of floats consisting of the multivariate response
- Z = 2d array of floats consisting of the column covariates, with all
  categorical variables coded in appropriate contrasts
- XTX = X*transpose(X) product as a 2d array of floats
- ZTZ = Z*transpose(Z) product as a 2d array of floats

# Value

2d array of floats

"""
function calc_coeffs(X::AbstractArray{Float64,2}, Y::AbstractArray{Float64,2},
                     Z::AbstractArray{Float64,2},
                     XTX::AbstractArray{Float64,2},
                     ZTZ::AbstractArray{Float64,2})

    return transpose(ZTZ\(transpose((XTX\transpose(X)*Y)*Z)))
end

"""
    calc_sigma(resid, targetType)

Estimates variance of errors and the shrinkage coefficient, without variance
shrinkage.

# Arguments

- resid = 2d array of floats consisting of the residuals
- targetType =`nothing`

# Value

Tuple
- sigma: 2d array of floats; estimated variance of errors
- lambda: 0.0

# Some notes

Since this version of `calc_sigma` does not implement variance shrinkage, the
shrinkage coefficient lambda is 0.

"""
function calc_sigma(resid::AbstractArray{Float64,2}, targetType::Nothing)

    # Residual sum of squares
    RSS = transpose(resid)*resid
    # Divide by the number of samples
    sigma = RSS ./ (size(resid,1) - 1)

    return sigma, 0.0
end

"""
    calc_sigma(resid, targetType)

Estimates variance of errors and the shrinkage coefficient, with variance
shrinkage.

# Arguments

- resid = 2d array of floats consisting of the residuals
- targetType = string indicating the target type toward which to shrink the
  variance. Acceptable inputs are "A", "B", "C", and "D".
    - "A": Target is identity matrix
    - "B": Target is diagonal matrix with constant diagonal
    - "C": Target is has same diagonal element, and same off-diagonal element
    - "D": Target is diagonal matrix with unequal entries

# Value

Tuple
- sigma: 2d array of floats; shrunk estimated variance of errors
- lambda: floating scalar; estimated shrinkage coefficient
  (0 = no shrinkage, 1 = complete shrinkage)

"""
function calc_sigma(resid::AbstractArray{Float64,2},
                    targetType::AbstractString)

    return shrink_sigma(resid, targetType)
end

"""
    calc_var(X, Z, XTX, ZTZ, sigma)

Calculate the variance (diagonal of the covariance matrix) of the coefficient
estimates.

# Arguments

- X = 2d array of floats consisting of the row covariates, with all
  categorical variables coded in appropriate contrasts
- Z = 2d array of floats consisting of the column covariates, with all
  categorical variables coded in appropriate contrasts
- XTX = X*transpose(X) product as a 2d array of floats
- ZTZ = Z*transpose(Z) product as a 2d array of floats
- sigma = 2d array of floats consisting of the estimated sigma

# Value

2d array of floats

"""
function calc_var(X::AbstractArray{Float64,2}, Z::AbstractArray{Float64,2},
                  XTX::AbstractArray{Float64,2},
                  ZTZ::AbstractArray{Float64,2},
                  sigma::AbstractArray{Float64,2})

    # LHS of covariance matrix
    varLeft = inv(XTX)
    # RHS of covariance matrix
    varRight = transpose(ZTZ\(transpose((ZTZ\(transpose(Z)*sigma))*Z)))

    # Diagonal of covariance matrix, aka the variance
    varDiag = transpose(kron_diag(varLeft, varRight))

    return varDiag
end


text:
# FibonacciAnyons
struct FibonacciAnyon <: Sector
    isone::Bool
    function FibonacciAnyon(s::Symbol)
        s == :I || s == :τ || throw(ArgumentError("Unknown FibonacciAnyon $s."))
        new(s === :I)
    end
end

Base.IteratorSize(::Type{SectorValues{FibonacciAnyon}}) = HasLength()
Base.length(::SectorValues{FibonacciAnyon}) = 2
Base.iterate(::SectorValues{FibonacciAnyon}, i = 0) =
    i == 0 ? (FibonacciAnyon(:I), 1) : (i == 1 ? (FibonacciAnyon(:τ), 2) : nothing)
function Base.getindex(S::SectorValues{FibonacciAnyon}, i)
    if i == 1
        return FibonacciAnyon(:I)
    elseif i == 2
        return FibonacciAnyon(:τ)
    else
        throw(BoundsError(S, i))
    end
end
findindex(::SectorValues{FibonacciAnyon}, s::FibonacciAnyon) = 2 - s.isone

Base.convert(::Type{FibonacciAnyon}, s::Symbol) = FibonacciAnyon(s)
Base.convert(::Type{Symbol}, a::FibonacciAnyon) = a.s ? (:I) : (:τ)
Base.one(::Type{FibonacciAnyon}) = FibonacciAnyon(:I)
Base.conj(s::FibonacciAnyon) = s

const _goldenratio = (1 + sqrt(5)) / 2
dim(a::FibonacciAnyon) = isone(a) ? one(_goldenratio) : _goldenratio

Base.@pure FusionStyle(::Type{FibonacciAnyon}) = SimpleNonAbelian()
Base.@pure BraidingStyle(::Type{FibonacciAnyon}) = Anyonic()
Base.isreal(::Type{FibonacciAnyon}) = false

⊗(a::FibonacciAnyon, b::FibonacciAnyon) = FibonacciIterator(a,b)

struct FibonacciIterator
    a::FibonacciAnyon
    b::FibonacciAnyon
end
Base.IteratorSize(::Type{FibonacciIterator}) = Base.HasLength()
Base.IteratorEltype(::Type{FibonacciIterator}) = Base.HasEltype()
Base.length(iter::FibonacciIterator) = (isone(iter.a) || isone(iter.b)) ? 1 : 2
Base.eltype(::Type{FibonacciIterator}) = FibonacciAnyon
function Base.iterate(iter::FibonacciIterator, state = 1)
    I = FibonacciAnyon(:I)
    τ = FibonacciAnyon(:τ)
    if state == 1 # first iteration
        iter.a == I && return (iter.b, 2)
        iter.b == I && return (iter.a, 2)
        return (I, 2)
    elseif state == 2
        (iter.a == iter.b == τ) && return (τ, 3)
        return nothing
    else
        return nothing
    end
end

Nsymbol(a::FibonacciAnyon, b::FibonacciAnyon, c::FibonacciAnyon) =
    isone(a) + isone(b) + isone(c) != 2 # zero if one tau and two ones

function Fsymbol(a::FibonacciAnyon, b::FibonacciAnyon, c::FibonacciAnyon,
                 d::FibonacciAnyon, e::FibonacciAnyon, f::FibonacciAnyon)
    Nsymbol(a, b, e) || return zero(_goldenratio)
    Nsymbol(e, c, d) || return zero(_goldenratio)
    Nsymbol(b, c, f) || return zero(_goldenratio)
    Nsymbol(a, f, d) || return zero(_goldenratio)

    I = FibonacciAnyon(:I)
    τ = FibonacciAnyon(:τ)
    if a == b == c == d == τ
        if e == f == I
            return +1/_goldenratio
        elseif e == f == τ
            return -1/_goldenratio
        else
            return +1/sqrt(_goldenratio)
        end
    else
        return one(_goldenratio)
    end
end

function Rsymbol(a::FibonacciAnyon, b::FibonacciAnyon, c::FibonacciAnyon)
    Nsymbol(a, b, c) || return 0*exp((0π/1)*im)
    if isone(a) || isone(b)
        return exp((0π/1)*im)
    else
        return isone(c) ? exp(+(4π/5)*im) : exp(-(3π/5)*im)
    end
end

Base.show(io::IO, ::Type{FibonacciAnyon}) = print(io, "FibonacciAnyon")
function Base.show(io::IO, a::FibonacciAnyon)
    s = isone(a) ? ":I" : ":τ"
    return get(io, :typeinfo, nothing) === FibonacciAnyon ?
        print(io, s) : print(io, "FibonacciAnyon(", s, ")")
end

Base.hash(a::FibonacciAnyon, h::UInt) = hash(a.isone, h)
Base.isless(a::FibonacciAnyon, b::FibonacciAnyon) = isless(!a.isone, !b.isone)


text:
using EasyFit
using Plots

x = sort(rand(10)) .- 0.5
y = sort(rand(10)) .- 0.5

plot(layout=(3,2))
plot!(framestyle=:box,grid=false)

sp=1
plot!(title="Linear fit",subplot=sp)
scatter!(x,y,label="",subplot=sp)
fit = fitlinear(x,y)
plot!(fit.x,fit.y,label="",linewidth=2,subplot=sp)

sp=2
y2 = y.^2
plot!(title="Quadratic fit",subplot=sp)
scatter!(x,y2,label="",subplot=sp)
fit = fitquad(x,y2)
plot!(fit.x,fit.y,label="",linewidth=2,subplot=sp)

sp=3
y2 = y.^3
plot!(title="Cubic fit",subplot=sp)
scatter!(x,y2,label="",subplot=sp)
fit = fitcubic(x,y2)
plot!(fit.x,fit.y,label="",linewidth=2,subplot=sp)

sp=4
y2 = 0.3*exp.(5*y)
plot!(title="Mono-exponential fit",subplot=sp)
scatter!(x,y2,label="",subplot=sp)
fit = fitexp(x,y2)
plot!(fit.x,fit.y,label="",linewidth=2,subplot=sp)

sp=5
y2 = 0.3*exp.(5*y) + 0.7*exp.(3*y2)
plot!(title="Bi-exponential fit",subplot=sp)
scatter!(x,y2,label="",subplot=sp)
fit = fitexp(x,y2,n=2)
plot!(fit.x,fit.y,label="",linewidth=2,subplot=sp)

sp=6
y2 = y.^5
plot!(title="Spline fit",subplot=sp)
scatter!(x,y2,label="",subplot=sp)
fit = fitexp(x,y2,n=2)
plot!(fit.x,fit.y,label="",linewidth=2,subplot=sp)

plot!(size=(600,600))


text:
using Requests
using DataFrames
################################################################################
# this function returns a DataFrame of articles (pmid, title, abstract)
function searchBreastCancerArticles(cancerType="breast neoplasms", researchType="diagnosis",
                                    minDate=1990,maxDate=2030, retmax=1000)

    search_terms = "\"$cancerType\"[mh] and \"$researchType\"[mh]"

    # Retrieve number of results - from ESearch results
    # define base search query for eutils
    base_search_query = "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi"

    # run esearch to get pmids
    search_result = readstring(post(base_search_query;
    data = Dict("db" => "pmc", "term" => "$search_terms", "retmax" => retmax,
                 "mindate"=>minDate, "maxdate"=>maxDate)))
    # println(search_result)

    # retrieve the pmids into an array
    pmids = []
    for result_line in split(search_result, "\n")
      # get list of pmids
      pmid = match(r"<Id>(\d+)<\/Id>", result_line)
      if pmid != nothing
        push!(pmids,pmid[1])
      end
    end

    # concatenate pmids into a single comma separated string
    id_string = join(collect(pmids), ",")

    # define base fetch query for eutils
    base_fetch_query = "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi"

    # retrieve metadata for pmid Set
    fetch_result_medline = readstring(post(base_fetch_query; data = Dict("db" => "pmc",
    "id" => id_string, "rettype" => "medline", "retmode" => "text")))

    fetch_result_xml = readstring(post(base_fetch_query; data = Dict("db" => "pmc",
      "id" => id_string, "rettype" => "xml", "retmode" => "text")))

    # save metadata in a file
    output_file = open("output/$researchType.txt", "w")
    write(output_file,fetch_result_medline)
    close(output_file)

    # field set
    pmcid_array =[] # pmcid set
    pmid_array = [] # pmid set
    title_array = [] # title set
    abstract_array=[] # article abstract
    date_published_array=[] # published online Date (DEP).
    # read through each line
    for fetch_article in split(fetch_result_medline, "PMC - ")
        # clear all field
        pmcid =""
        pmid=""
        date_published=""
        title=""

        for fetch_line in split(fetch_article, "\n")
            # extract PMCID
            pmcid_str = match(r"(PMC[0-9]+)", fetch_line)
            if pmcid_str != nothing
              pmcid = pmcid_str[1]
            end

            # extract PMID
            pmid_str = match(r"PMID- ([\w\W \r\n]+)", fetch_line)
            if pmid_str != nothing
              pmid = pmid_str[1]
            end

            # date published
            date_published_str = match(r"DEP - ([\w\W \r\n]+)", fetch_line)
            if date_published_str != nothing
                date_published = date_published_str[1]
            end

            # date published
            title_str = match(r"TI  - ([\w\W \r\n]+)", fetch_line)
            if title_str != nothing
                title = title_str[1]
            end

        end
        # push everything
        if pmid != ""
            push!(pmcid_array, pmcid)
            push!(pmid_array, pmid)
            push!(date_published_array, date_published)
            push!(title_array, title)

        end

    end

#=    for fetch_article in split(fetch_result_xml, "<front>")
        # clear all field
        full_title=""
        full_abstract=""

       # get Title
        title = match(r"<title-group>([\w\W \r\n]+)</title-group>", fetch_article)
        if title != nothing
            full_title=title[1]
        end
        # get abstract
        abstract_str = match(r"<abstract>([\w\W \r\n]+)<\/abstract>", fetch_article)
        if abstract_str != nothing
            full_abstract=abstract_str[1]
        end
            push!(title_array, full_title)
#            push!(abstract_array, full_abstract)
#        end
    end
=# # return dictionary of pmid and title
    return(DataFrame(pmcid = pmcid_array,
                     pmid = pmid_array,
                     date_published = date_published_array,
                     title = title_array
#                    main_abstract = abstract_array
                    )
            )
end

# usage
#df = fetchBreastCancerArticles("breast cancer", "diagnosis",2015,2018,10)

#output_file = open("output/diagnosis.csv", "w")
#close(output_file)
#writetable("output/diagnosis.csv",df)

# ref
#https://www.nlm.nih.gov/bsd/mms/medlineelements.html#ab

max_stars_repo_name: jmulloka/Breast-Cancer-Text-Mining

max_issues_repo_name: jmulloka/Breast-Cancer-Text-Mining

max_forks_repo_name: jmulloka/Breast-Cancer-Text-Mining


text:
# Minimum spanning tree in complete graphs.
# Only the upper part of the reward matrix is considered!

abstract type SpanningTreeSolver end
function build!(solver::SpanningTreeSolver, graph::SimpleGraph)
  nothing
end

struct SpanningTree <: CombinatorialInstance{Tuple{Int, Int}}
  # Mandatory properties.
  n_arms::Int
  optimal_average_reward::Float64

  # Probability distributions for the arm rewards.
  graph::SimpleGraph
  reward::Dict{Tuple{Int, Int}, Distribution} # Roughly a matrix, except that not all entries should be filled.
  # Even for a complete graph, the matrix has too many elements: the diagonal is not used.
  # For a very large number of reads into the data structure, the dictionary seems a good fit.
  # TODO: uniformise reward vs. rewards.

  # Internal solver.
  solver::SpanningTreeSolver

  function SpanningTree(graph::SimpleGraph, reward::Dict{Tuple{Int, Int}, Distribution}, solver::SpanningTreeSolver)
    if length(reward) != ne(graph)
      error("The edges and rewards do not perfectly match: there are $(length(reward)) rewards and $(ne(graph)) edges.")
    end
    for v in keys(reward)
      if ! has_edge(graph, v[1], v[2])
        error("The edges and rewards do not perfectly match: the edge $v has a reward, but does not exist in the graph.")
      end
    end

    # Prepare the solver to be used (if required).
    build!(solver, graph)

    avg_reward = Dict{Tuple{Int, Int}, Float64}(k => mean(v) for (k, v) in reward)
    opt_sol = solve_linear(solver, avg_reward)
    opt = sum(avg_reward[i] for i in opt_sol)

    # Done!
    return new(ne(graph), opt, graph, reward, solver)
  end
end

function initial_state(instance::SpanningTree)
  zero_counts = Dict(k => 0 for (k, _) in instance.reward)
  zero_rewards = Dict(k => 0.0 for (k, _) in instance.reward)
  return State{Tuple{Int, Int}}(0, 0.0, 0.0, zero_counts, zero_rewards, copy(zero_rewards))
end

solve_linear(instance::SpanningTree, rewards::Dict{Tuple{Int, Int}, Float64}) = solve_linear(instance.solver, rewards)
has_lp_formulation(instance::SpanningTree) = has_lp_formulation(instance.solver)
get_lp_formulation(instance::SpanningTree, rewards::Dict{Tuple{Int, Int}, Float64}) = has_lp_formulation(instance) ?
  get_lp_formulation(instance.solver, rewards) :
  error("The complete graph minimum spanning tree solver uses no LP formulation.")

function _mst_solution_normalise(reward::Dict{Tuple{Int, Int}, Float64}, solution::Vector{Tuple{Int, Int}})
  # Not all solvers return the right kind of edges: they are sometimes reversed with respect
  # to the input graph. Hence return those so that all arms in the output solution
  # correspond exactly to those in the input graph.
  normalised = Tuple{Int, Int}[]
  for e in solution
    if e in keys(reward)
      push!(normalised, e)
    else
      push!(normalised, (e[2], e[1]))
    end
  end
  return normalised
end

_mst_solution_normalise(instance::SpanningTree, solution::Vector{Tuple{Int, Int}}) =
  _mst_solution_normalise(instance.reward, solution)

function is_feasible(instance::SpanningTree, arms::Vector{Tuple{Int, Int}})
  n = nv(instance.graph)
  if length(arms) != n - 1
    return false
  end

  # An MST is feasible if all nodes are reachable from any given node (undirected graph).
  # However, not all edges are in order in the input vector.
  reachable_nodes = BitSet()
  sizehint!(reachable_nodes, n)
  push!(reachable_nodes, arms[1][1])

  arms_to_do = BitSet(1:(n - 1))

  while ! isempty(arms_to_do)
    has_changed = false

    for arm_id in arms_to_do
      arm = arms[arm_id]
      if arm[1] ∈ reachable_nodes
        push!(reachable_nodes, arm[2])
        delete!(arms_to_do, arm_id)
        has_changed = true
      end

      if arm[2] ∈ reachable_nodes
        push!(reachable_nodes, arm[1])
        delete!(arms_to_do, arm_id)
        has_changed = true
      end
    end

    if ! has_changed
      return false
    end
  end

  return true
end

function is_partially_acceptable(instance::SpanningTree, arms::Vector{Tuple{Int, Int}})
  if length(arms) == 0
    return true
  end

  if length(arms) > nv(instance.graph) - 1
    return false
  end

  # An MST is partially acceptable if it contains no cycle. However, not all edges are in order in the input vector.
  # Implement this with a graph search: if a previous node is found again (i.e. both ends of an edge have already been visited),
  # then there is a cycle.
  reached_nodes = BitSet()
  sizehint!(reached_nodes, length(arms) + 1)

  push!(reached_nodes, arms[1][1])
  push!(reached_nodes, arms[1][2])
  for i in 2:length(arms)
    arm1in = arms[i][1] ∈ reached_nodes
    arm2in = arms[i][2] ∈ reached_nodes

    if arm1in && arm2in
      return false
    end

    push!(reached_nodes, arms[i][1])
    push!(reached_nodes, arms[i][2])
  end

  return true
end


text:
function ∇ₛx end
function μⁱ∇ₛx end
function ∇x  end
function μⁱ∇x  end

function (sim::Simulation)(::typeof(∇x),gridtype::GridType)
    ∂x =   sim(∂, x̂,gridtype);
    ∂y =   sim(∂, ŷ,gridtype);
    ∂z =   sim(∂, ẑ,gridtype);
    O = spzeros(ncells(sim.grid),ncells(sim.grid));
     c = [O    -∂z     ∂y;
            ∂z     O     -∂x;
          -∂y    ∂x      O
        ]
end
function (sim::Simulation)(::typeof(∇ₛx),gridtype::GridType)
    ∂x =   sim(∂ₛ, x̂,gridtype);
    ∂y =   sim(∂ₛ, ŷ,gridtype);
    ∂z =   sim(∂ₛ, ẑ,gridtype);
    O = spzeros(ncells(sim.grid),ncells(sim.grid));
     c = [O    -∂z     ∂y;
            ∂z     O     -∂x;
          -∂y    ∂x      O
        ]
end

function (sim::Simulation)(::typeof(μⁱ∇x))
  μi =  sim(convert_to_diagonal_matrix,:μᵣ)
  μi  * sim(∇x, p̂)
end

function (sim::Simulation)(::typeof(μⁱ∇ₛx))
  μi =  sim(convert_to_diagonal_matrix,:μᵣ)
  μi  * sim(∇ₛx, p̂)
end


text:
using Test
using Bigsimr
using Distributions
using PDMats

@testset "Utilities" begin

    @testset "Special Functions" begin
        test_types = (Float64, Float32, Float16)
        for T in test_types
            # Must work for each type
            @test_nowarn Bigsimr._normpdf(T(1.0))
            @test_nowarn Bigsimr._normcdf(T(1.0))
            @test_nowarn Bigsimr._norminvcdf(T(0.5))

            # Must respect floating point types
            @test typeof(Bigsimr._normpdf(T(1.0))) === T
            @test typeof(Bigsimr._normcdf(T(1.0))) === T
            @test typeof(Bigsimr._norminvcdf(T(0.5))) === T
        end
    end

    @testset "PDCorMat" begin
        d = 10
        T = Float32
        p = cor_randPD(T, d)
        a = PDCorMat(p, Spearman)

        @test dim(a) == d
        @test size(a) == (d, d)
        @test size(a, 1) == size(a, 2) == d
        @test ndims(a) == 2
        @test eltype(a) == T
        @test Matrix(a) == p
        @test all(diag(a) .== one(T))

        @test_nowarn inv(a)
        @test_nowarn eigmax(a)
        @test_nowarn eigmin(a)
        @test_nowarn logdet(a)

        x = rand(T, d)
        @test_nowarn a * x
        @test_nowarn a \ x

        c = T(π)
        @test_nowarn a * c
        @test_nowarn c * a
        @test a * c == c * a
        @test eltype(a*c) == T

        b = PDCorMat(cor_randPD(Float32, 10), Spearman)
        @test_nowarn a + b
        @test_nowarn pdadd(a, b, c)

        m = rand(T, d, d)
        @test_nowarn pdadd(m, a)        # add `a` to a dense matrix `m` of the same size.
        @test_nowarn pdadd(m, a, c)     # add `a * c` to a dense matrix `m` of the same size.
        @test_nowarn pdadd!(m, a)       # add `a` to a dense matrix `m` of the same size inplace.
        @test_nowarn pdadd!(m, a, c)    # add `a * c` to a dense matrix `m` of the same size inplace.

        r = similar(m)
        @test_nowarn pdadd!(r, m, a)    # add `a` to a dense matrix `m` of the same size, and write the result to `r`.
        @test_nowarn pdadd!(r, m, a, c) # add `a * c` to a dense matrix `m` of the same size, and write the result to `r`.

        @test_nowarn quad(a, x)         # compute x' * a * x when `x` is a vector.
                                        # perform such computation in a column-wise fashion, when
                                        # `x` is a matrix, and return a vector of length `n`,
                                        # where `n` is the number of columns in `x`.
        @test_nowarn invquad(a, x)      # compute x' * inv(a) * x when `x` is a vector.
                                        # perform such computation in a column-wise fashion, when
                                        # `x` is a matrix, and return a vector of length `n`.

        r = similar(x)
        x = rand(T, d, d)
        @test_nowarn quad(a, x)
        @test_nowarn quad!(r, a, x)     # compute x' * a * x in a column-wise fashion, and write the results to `r`.
        @test_nowarn invquad(a, x)      # compute x' * inv(a) * x when `x` is a vector.
                                        # perform such computation in a column-wise fashion, when
                                        # `x` is a matrix, and return a vector of length `n`.
        @test_nowarn invquad!(r, a, x)  # compute x' * inv(a) * x in a column-wise fashion, and write the results to `r`.

        @test_nowarn X_A_Xt(a, x)       # compute `x * a * x'` for a matrix `x`.
        @test_nowarn Xt_A_X(a, x)       # compute `x' * a * x` for a matrix `x`.
        @test_nowarn X_invA_Xt(a, x)    # compute `x * inv(a) * x'` for a matrix `x`.
        @test_nowarn Xt_invA_X(a, x)    # compute `x' * inv(a) * x` for a matrix `x`.

        z = rmvn(1000, p)
        y = similar(z)
        r = similar(z)
        @test_nowarn whiten(a, z)       # whitening transform. `x` can be a vector or a matrix.
        @test_nowarn whiten!(a, y)      # whitening transform inplace, directly updating `y`.
        @test_nowarn whiten!(r, a, z)   # write the transformed result to `r`.

        x = randn(T, 1000, d)
        z = similar(x)
        r = similar(x)
        @test_nowarn unwhiten(a, x)     # inverse of whitening transform. `x` can be a vector or a matrix.
        @test_nowarn unwhiten!(a, z)    # un-whitening transform inplace, updating `z`.
        @test_nowarn unwhiten!(r, a, x) # write the transformed result to `r`.
    end

end


text:
@defcomp co2emissions begin
    region=Index()

    e_globalCO2emissions=Variable(index=[time],unit="Mtonne/year")
    e0_baselineCO2emissions=Parameter(index=[region],unit="Mtonne/year")
    e_regionalCO2emissions=Variable(index=[time,region],unit="Mtonne/year")
    er_CO2emissionsgrowth=Parameter(index=[time,region],unit="%")

    function run_timestep(p, v, d, t)

        #eq.4 in Hope (2006) - regional CO2 emissions as % change from baseline
        for r in d.region
            v.e_regionalCO2emissions[t,r]=p.er_CO2emissionsgrowth[t,r]*p.e0_baselineCO2emissions[r]/100
        end
        #eq. 5 in Hope (2006) - global CO2 emissions are sum of regional emissions
        v.e_globalCO2emissions[t]=sum(v.e_regionalCO2emissions[t,:])
    end
end

max_stars_repo_name: lrennels/PAGE-ICE

max_issues_repo_name: lrennels/PAGE-ICE

max_forks_repo_name: lrennels/PAGE-ICE


text:
using SnoopCompile

# Note: in order to generate a new precompile.jl file,
# you need to comment out the lines:
# ```
# include("../precompile/precompile.jl")
# _precompile_()
# ```
# in src/ITensors.jl (otherwise many functions will
# be precompiled and SnoopCompile won't think they
# need to get precompiled).
#
# Also, snooping on dmrg itself causes problems in
# precompilation. Explicitly add the command:
# ```
# precompile(Tuple{typeof(dmrg),MPO,MPS,Sweeps})
# ```
# to precompile.jl after it is generated.

inf_timing = @snoopi tmin=0.01 include("snoop/snoop.jl")
pc = SnoopCompile.parcel(inf_timing)
SnoopCompile.write("tmp", pc)
cp("tmp/precompile_ITensors.jl",
   "precompile.jl";
   force=true)


text:
export
    MultiStageCBS

module MultiStageCBS

using ..CRCBS
using Parameters, LightGraphs, DataStructures
using GraphUtils: get_graph, get_vtx

# NOTE MultiStageCBS Could also serve as the definition for "iterative MAPF" as
# defined in "Priority Inheritance with Backtracking for Iterative Multi-agent
# Path Finding", Okumura et al, IJCAI 2019

################################################################################
############################### ENVIRONMENT DEF ################################
################################################################################
# state
@with_kw struct State <: AbstractGraphState
    vtx::Int    = -1 # vertex of graph
    stage::Int  = -1 # which stage of the sequence
    t::Int      = -1
end
Base.convert(::Type{State},s::GraphState) = State(vtx=get_vtx(s),t=get_t(s))
# CRCBS.get_vtx(s::State) = s.vtx
CRCBS.get_t(s::State) = s.t
get_stage(s::State) = s.stage
Base.string(s::State) = "(v=$(get_vtx(s)),stage=$(get_stage(s)),t=$(get_t(s)))"
CRCBS.is_valid(state::State) = get_vtx(state) > 0
# action
@with_kw struct Action <: AbstractGraphAction
    e::Edge{Int}    = Edge(-1,-1)
    dt::Int         = 1
end
Base.convert(::Type{Action},a::GraphAction) = Action(e=get_e(a),dt=get_dt(a))
CRCBS.get_e(a::Action) = a.e
CRCBS.get_dt(a::Action) = a.dt
Base.string(a::Action) = "(e=$(get_e(a).src) → $(get_e(a).dst))"
@with_kw struct LowLevelEnv{C<:AbstractCostModel,H<:LowLevelSearchHeuristic,G<:AbstractGraph,T} <: GraphEnv{State,Action,C}
    graph::G                        = Graph()
    goal_sequence::Vector{State}    = Vector{State}()
    agent_idx::Int                  = -1
    constraints::T                  = discrete_constraint_table(nv(graph),nv(graph)^2,agent_idx)
    cost_model::C                   = SumOfTravelTime()
    heuristic::H                    = NullHeuristic() # MultiStagePerfectHeuristic(graph,Vector{Vector{Int}}())
end
# CRCBS.get_graph(env::LowLevelEnv)            = env.graph
# GraphUtils.get_graph(env::LowLevelEnv)            = env.graph
# CRCBS.get_cost_model(env::LowLevelEnv)       = env.cost_model
# CRCBS.get_agent_id(env::LowLevelEnv)         = env.agent_idx
# CRCBS.get_constraints(env::LowLevelEnv)      = env.constraints
# CRCBS.get_goal(env::LowLevelEnv)             = env.goal_sequence
# CRCBS.get_heuristic_model(env::LowLevelEnv)  = env.heuristic
CRCBS.base_env_type(env::LowLevelEnv)        = LowLevelEnv

function CRCBS.get_next_state(s::State,a::Action)
    @assert(is_valid(s))
    @assert(get_vtx(s) == get_e(a).src)
    State(get_e(a).dst, get_stage(s), get_t(s)+get_dt(a))
end
function CRCBS.get_next_state(env::E,s::State,a::Action) where {E<:LowLevelEnv}
    @assert(is_valid(s))
    @assert(get_stage(s) <= length(env.goal_sequence))
    stage = get_stage(s)
    if states_match(s, env.goal_sequence[get_stage(s)])
        stage = min(stage+1, length(env.goal_sequence))
    end
    return State(get_e(a).dst, stage, get_t(s)+get_dt(a))
end
CRCBS.get_possible_actions(env::LowLevelEnv,s)  = map(v->Action(e=Edge(get_vtx(s),v)),outneighbors(get_graph(env),get_vtx(s)))
CRCBS.wait(s::State) = Action(e=Edge(get_vtx(s),get_vtx(s)))
CRCBS.wait(env::E,s::State) where {E<:LowLevelEnv} = Action(e=Edge(get_vtx(s),get_vtx(s)))

function CRCBS.build_env(mapf::MAPF{E,S,G}, node::N, idx::Int)  where {S,G,E <: LowLevelEnv,N<:ConstraintTreeNode}
    goals = deepcopy(mapf.goals[idx])
    n = PathNode{State,Action}(sp=goals[end])
    s_constraints, _ = search_constraints(mapf.env,get_constraints(node,idx),n)
    t_goal = get_t(get_sp(n))
    for c in s_constraints
        t_goal = max(t_goal,get_time_of(c)+1)
    end
    goals[end] = State(get_sp(n),t=t_goal)
    E(
        graph = get_graph(mapf.env),
        constraints = get_constraints(node,idx),
        goal_sequence = goals,
        agent_idx = idx,
        cost_model = get_cost_model(mapf.env),
        heuristic = get_heuristic_model(mapf.env)
        )
end
function CRCBS.is_consistent(path::Path,start::State,goals::Vector{State})
    valid = true
    stage = 1
    if states_match(start,goals[stage])
        stage += 1
    end
    for k in 1:length(path)
        node = get_path_node(path,k)
        if states_match(get_sp(node),goals[stage])
            stage += 1
        end
        if stage > length(goals)
            return true
        end
    end
    return false
end
function CRCBS.is_goal(env::E,s::State) where {E<:LowLevelEnv}
    if get_stage(s) == length(env.goal_sequence)
        if states_match(s, env.goal_sequence[get_stage(s)])
            if get_t(s) >= get_t(env.goal_sequence[get_stage(s)])
                return true
            end
        end
    end
    return false
end
function CRCBS.get_heuristic_cost(h::MultiStagePerfectHeuristic,env::E,s::State) where {E<:LowLevelEnv}
# function CRCBS.get_heuristic_cost(env::E,h::MultiStagePerfectHeuristic,s::State) where {E<:LowLevelEnv}
    get_heuristic_cost(h, env.agent_idx, s.stage, s.vtx)
end

end # end module MultiStageCBS


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

"""
    ZScore()

The z-score (a.k.a. normal score) of `x` with mean `μ` and
standard deviation `σ` is the value `(x .- μ) ./ σ`.
"""
struct ZScore <: Colwise end

assertions(::Type{ZScore}) = [assert_continuous]

isrevertible(::Type{ZScore}) = true

function colcache(::ZScore, x)
  μ = mean(x)
  σ = std(x, mean=μ)
  (μ=μ, σ=σ)
end

colapply(::ZScore, x, c)  = @. (x - c.μ) / c.σ

colrevert(::ZScore, y, c) = @. c.σ * y + c.μ


text:
using Test
using PolynomialBasis
using ImplicitDomainQuadrature
# using Revise
using CutCellDG
include("../useful_routines.jl")

function displacement(alpha, x)
    u1 = alpha * x[2] * sin(pi * x[1])
    u2 = alpha * (x[1]^3 + cos(pi * x[2]))
    return [u1, u2]
end

function body_force(lambda, mu, alpha, x)
    b1 = alpha * (lambda + 2mu) * pi^2 * x[2] * sin(pi * x[1])
    b2 =
        -alpha * (6mu * x[1] + (lambda + mu) * pi * cos(pi * x[1])) +
        alpha * (lambda + 2mu) * pi^2 * cos(pi * x[2])
    return [b1, b2]
end

function onboundary(x, L, W)
    return x[2] ≈ 0.0 || x[1] ≈ L || x[2] ≈ W || x[1] ≈ 0.0
end

function error_for_plane_interface(
    x0,
    normal,
    nelmts,
    polyorder,
    numqp,
    penaltyfactor,
    eta,
)
    L = 1.0
    W = 1.0
    lambda, mu = 1.0, 2.0
    dx = 1.0 / nelmts
    penalty = penaltyfactor / dx * (lambda + mu)
    alpha = 0.1
    stiffness = CutCellDG.HookeStiffness(lambda, mu, lambda, mu)

    elasticitybasis = LagrangeTensorProductBasis(2, polyorder)
    levelsetbasis = HermiteTensorProductBasis(2)

    dim, nf = size(interpolation_points(levelsetbasis))
    refpoints = interpolation_points(elasticitybasis)

    cgmesh =
        CutCellDG.CGMesh([0.0, 0.0], [L, W], [nelmts, nelmts], nf)
    mesh = CutCellDG.DGMesh([0.0, 0.0], [L, W], [nelmts, nelmts], refpoints)

    levelset = CutCellDG.LevelSet(
        x -> plane_distance_function(x, normal, x0),
        cgmesh,
        levelsetbasis,
    )

    cutmesh = CutCellDG.CutMesh(mesh, levelset)
    cellquads = CutCellDG.CellQuadratures(cutmesh, levelset, numqp)
    interfacequads = CutCellDG.InterfaceQuadratures(cutmesh, levelset, numqp)
    facequads = CutCellDG.FaceQuadratures(cutmesh, levelset, numqp)

    mergedmesh =
        CutCellDG.MergedMesh!(cutmesh, cellquads, facequads, interfacequads)

    sysmatrix = CutCellDG.SystemMatrix()
    sysrhs = CutCellDG.SystemRHS()

    CutCellDG.assemble_displacement_bilinear_forms!(
        sysmatrix,
        elasticitybasis,
        cellquads,
        stiffness,
        mergedmesh,
    )
    CutCellDG.assemble_interelement_condition!(
        sysmatrix,
        elasticitybasis,
        facequads,
        stiffness,
        mergedmesh,
        penalty,
        eta,
    )
    CutCellDG.assemble_coherent_interface_condition!(
        sysmatrix,
        elasticitybasis,
        interfacequads,
        stiffness,
        mergedmesh,
        penalty,
        eta,
    )
    CutCellDG.assemble_penalty_displacement_bc!(
        sysmatrix,
        sysrhs,
        x -> displacement(alpha, x),
        elasticitybasis,
        facequads,
        stiffness,
        mergedmesh,
        x -> onboundary(x, L, W),
        penalty,
    )
    CutCellDG.assemble_body_force!(
        sysrhs,
        x -> body_force(lambda, mu, alpha, x),
        elasticitybasis,
        cellquads,
        mergedmesh,
    )

    matrix = CutCellDG.sparse_displacement_operator(sysmatrix, mergedmesh)
    rhs = CutCellDG.displacement_rhs_vector(sysrhs, mergedmesh)

    solution = matrix \ rhs
    nodaldisplacement = reshape(solution, 2, :)

    err = mesh_L2_error(
        nodaldisplacement,
        x -> displacement(alpha, x),
        elasticitybasis,
        cellquads,
        mergedmesh,
    )
end

x0 = [0.5, 0.0]
normal = [1.0, 0.0]
powers = [3, 4, 5]
nelmts = [2^p + 1 for p in powers]
polyorder = 1
numqp = required_quadrature_order(polyorder)
penaltyfactor = 1e2
eta = 1

err = [
    error_for_plane_interface(
        x0,
        normal,
        ne,
        polyorder,
        numqp,
        penaltyfactor,
        eta,
    ) for ne in nelmts
]

u1err = [er[1] for er in err]
u2err = [er[2] for er in err]

dx = 1.0 ./ nelmts
u1rate = convergence_rate(dx, u1err)
u2rate = convergence_rate(dx, u2err)

@test allapprox(u1rate, repeat([2.0], length(u1rate)), 0.1)
@test allapprox(u1rate, repeat([2.0], length(u2rate)), 0.1)

x0 = [0.5, 0.0]
normal = [1.0, 0.0]
powers = [3, 4, 5]
nelmts = [2^p + 1 for p in powers]
polyorder = 2
numqp = required_quadrature_order(polyorder)
penaltyfactor = 1e2
eta = 1

err = [
    error_for_plane_interface(
        x0,
        normal,
        ne,
        polyorder,
        numqp,
        penaltyfactor,
        eta,
    ) for ne in nelmts
]

u1err = [er[1] for er in err]
u2err = [er[2] for er in err]

dx = 1.0 ./ nelmts
u1rate = convergence_rate(dx, u1err)
u2rate = convergence_rate(dx, u2err)

@test allapprox(u1rate, repeat([3.0], length(u1rate)), 0.05)
@test allapprox(u1rate, repeat([3.0], length(u2rate)), 0.05)


text:
export Emoji

"""
An emoji.
More details [here](https://discordapp.com/developers/docs/resources/emoji#emoji-object).
"""
struct Emoji
    id::Union{Snowflake, Nothing}
    name::String
    roles::Union{Vector{Snowflake}, Missing}
    user::Union{User, Missing}
    require_colons::Union{Bool, Missing}
    managed::Union{Bool, Missing}
    animated::Union{Bool, Missing}
end
@boilerplate Emoji :dict :docs :lower :merge


text:
function inverseSide(side::HorizonSide)
    return HorizonSide(mod(Int(side)+2, 4))
end

function main(r::Robot)
    stepsToWest = moveSide(r, West)
    stepsToNord = moveSide(r, Sud)

    counter = moveSide(r, Ost)
    moveSide(r, West)

    while true
        putmarker!(r)
        moveSideByCount(r, Ost, counter)
        moveSide(r, West)
        if isborder(r, Nord) == true
            break
        end

        moveStep(r, Nord, false)
        counter -= 1
    end

    moveSide(r, West)
    moveSide(r, Sud)

    followByCount(r, Ost, stepsToWest)
    followByCount(r, Nord, stepsToNord)
end

function moveSide(r::Robot, side::HorizonSide, needMark=false)
    steps = 0
    while isborder(r, side) != true
        moveStep(r, side, needMark)
        steps += 1
    end
    return steps
end

function moveSideByCount(r::Robot, side::HorizonSide, count)
    for i in 1:count
        if isborder(r, side) != true
            moveStep(r, side, true)
        else
            break
        end
    end
end

function followByCount(r::Robot, side::HorizonSide, count::Int)
    for i = 0:count-1
        moveStep(r, side, false)
    end
end

function moveStep(r::Robot, side::HorizonSide, needMark::Bool)
    move!(r, side)
    if needMark
        putmarker!(r)
    end
end

max_stars_repo_name: TheDarkMoon/mirea-progs

max_issues_repo_name: TheDarkMoon/mirea-progs

max_forks_repo_name: TheDarkMoon/mirea-progs


text:
path = "/home/user/GraphEvolve.jl/src";
savepath = "/tmp";
# path = "/home/perot/julia/GraphEvolve.jl/src";
# savepath = "/home/perot/julia/scaling/data";
	push!(LOAD_PATH, path);
	using GraphEvolve;
	using DataFrames;
	using Dates;
	using ArgParse;
	using JLD;

# %%

function simulate_parallel_seeds(
	graph_type,
	evolution_method,
	system_size,
	n_steps,
	n_sims,
	starting_seed,
	t_data
	)
	row(g) = [
		replace(string(graph_type), "GraphEvolve." => ""),
		replace(replace(string(evolution_method), "GraphEvolve." => ""), "!" => ""),
		g.N,
		g.t,
		Int(g.rng.seed[1]),
		g.observables.delta[1],
		g.observables.delta[2],
		g.observables.delta[3],
		g.observables.delta[4],
	]

	rows = Array{Any, 1}(undef, n_sims)
	seeds = collect(starting_seed:starting_seed+n_sims-1)

	t₀ = now()
	Threads.@threads for i in 1:n_sims
		rows[i] = row(evolution_method(graph_type(system_size, seed=seeds[i]), n_steps, t_data, joinpath(savepath, "cluster_size_distribution")))
	end
	t₁ = now()

	runtime = Dates.value(t₁ - t₀) / 1000
	println("Time to run k = $(Int(log2(system_size))): $(runtime)s")

	return rows, runtime
end

function write_rows_to_csv(rows, file_name)
	open(file_name, "a") do f
    	for row in rows
			write(f, string(join(row, ","), "\n"))
		end
	end
end

function simulate_system_sizes(
	graph_type,
	evolution_method,
	system_sizes,
	n_sims,
	n_steps_ratio,
	starting_seed,
	t_data,
	file_name
	)
	# write the header row
	header = [
		"graph_type",
		"evolution_method",
		"N",
		"t",
		"seed",
		"r_0",
		"r_1",
		"m_0",
		"m_1"
	]

	open(file_name, "w") do f
		write(f, string(join(header, ","), "\n"))
	end

	# run the simulations for the given sizes
	runtimes = []
	for system_size in system_sizes
		rows, runtime = simulate_parallel_seeds(graph_type,
			evolution_method,
			system_size,
			Int(floor(n_steps_ratio * system_size))+1,
			n_sims,
			starting_seed,
			t_data,
		)
		push!(runtimes, (system_size, n_sims, Threads.nthreads(), runtime))
		write_rows_to_csv(rows, file_name)
	end

	# save the average run times
	open(joinpath(savepath, string("runtimes_$(Int(log2(system_sizes[1]))):$(Int(log2(system_sizes[end])))_$(starting_seed)_", Dates.now(), ".csv")), "w") do f
		write(f, "system_size,n_sims,n_threads,time\n")
		for runtime in runtimes
			write(f, string(join(runtime, ","), "\n"))
		end
	end
end

function parse_commandline()
	s = ArgParseSettings()

	@add_arg_table s begin
		"--graph_type"
			help = "Network, Lattice2D, or Lattice3D"
			arg_type = String
			default = "Network"
		"--evolution_method"
			help = "erdos_renyi, bohman_frieze, product_rule, or stochastic_edge_acceptance"
			arg_type = String
			default = "stochastic_edge_acceptance"
		"--n_sims"
			help = "Number of simulations to run for a given system size"
			arg_type = Int
			default = 100
		"--n_steps_ratio"
			help = "Relative (to n) number of edges to add to the graph"
			arg_type = Float64
			default = 0.8
		"--starting_seed"
			help = "Inital seed value with which to seed the rng"
			arg_type = Int
			default = 1
		"--k_min"
			help = "Minimum system size 2^k"
			arg_type = Int
			default = 15
		"--k_max"
			help = "Maximum system size 2^k"
			arg_type = Int
			default = 22
	end

	return parse_args(s)
end

# %%

graph_types = Dict(
	"Network" => Network,
	"Lattice2D" => Lattice2D,
	"Lattice3D" => Lattice3D
)

evolution_methods = Dict(
	"erdos_renyi" => erdos_renyi!,
	"bohman_frieze" => bohman_frieze!,
	"product_rule" => product_rule!,
	"stochastic_edge_acceptance" => stochastic_edge_acceptance!
)

parsed_args      = parse_commandline()
graph_type       = graph_types[parsed_args["graph_type"]]
evolution_method = evolution_methods[parsed_args["evolution_method"]]
n_sims           = parsed_args["n_sims"]
n_steps_ratio    = parsed_args["n_steps_ratio"]
starting_seed    = parsed_args["starting_seed"]
k_min            = parsed_args["k_min"]
k_max            = parsed_args["k_max"]

# %%

println(string("Available n_threads: ", Threads.nthreads()))
free_mem = Int(Sys.free_memory()) / 1024^3
println(string("Availabe memory [GB]: ", free_mem))

t_data = load("/home/user/thesis/code/scaling/t_data.jld")["t_data"]
# t_data = load("/home/perot/julia/scaling/t_data.jld")["t_data"]

# %%

simulate_system_sizes(
	graph_type,
	evolution_method,
	2 .^ collect(k_min:k_max),
	n_sims,
	n_steps_ratio,
	starting_seed,
	t_data,
	joinpath(
		savepath,
		"csv",
		string(
			parsed_args["graph_type"], "_",
			parsed_args["evolution_method"], "_",
			k_min, ":", k_max, "_",
			starting_seed, "_",
			Dates.today(),
			".csv"
		)
	)
)

# %%

# simulate_system_sizes(Network, stochastic_edge_acceptance!, 2 .^ [15, 16], 10, 0.8, 10, t_data, "/tmp/test.csv")

max_stars_repo_name: cameronperot/explosive-percolation

max_issues_repo_name: cameronperot/explosive-percolation

max_forks_repo_name: cameronperot/explosive-percolation


text:
@testset "num00a" begin
    rmatn = num00a(2400000.5, 53736.0)
    @test rmatn[1,1] ≈ 0.9999999999536227949 atol=1e-12
    @test rmatn[1,2] ≈ 0.8836238544090873336e-5 atol=1e-12
    @test rmatn[1,3] ≈ 0.3830835237722400669e-5 atol=1e-12
    @test rmatn[2,1] ≈ -0.8836082880798569274e-5 atol=1e-12
    @test rmatn[2,2] ≈ 0.9999999991354655028 atol=1e-12
    @test rmatn[2,3] ≈ -0.4063240865362499850e-4 atol=1e-12
    @test rmatn[3,1] ≈ -0.3831194272065995866e-5 atol=1e-12
    @test rmatn[3,2] ≈ 0.4063237480216291775e-4 atol=1e-12
    @test rmatn[3,3] ≈ 0.9999999991671660338 atol=1e-12
end

@testset "num00b" begin
    rmatn = num00b(2400000.5, 53736.0)
    @test rmatn[1,1] ≈ 0.9999999999536069682 atol=1e-12
    @test rmatn[1,2] ≈ 0.8837746144871248011e-5 atol=1e-12
    @test rmatn[1,3] ≈ 0.3831488838252202945e-5 atol=1e-12
    @test rmatn[2,1] ≈ -0.8837590456632304720e-5 atol=1e-12
    @test rmatn[2,2] ≈ 0.9999999991354692733 atol=1e-12
    @test rmatn[2,3] ≈ -0.4063198798559591654e-4 atol=1e-12
    @test rmatn[3,1] ≈ -0.3831847930134941271e-5 atol=1e-12
    @test rmatn[3,2] ≈ 0.4063195412258168380e-4 atol=1e-12
    @test rmatn[3,3] ≈ 0.9999999991671806225 atol=1e-12
end

@testset "num06a" begin
    rmatn = num06a(2400000.5, 53736.)
    @test rmatn[1,1] ≈ 0.9999999999536227668 atol=1e-12
    @test rmatn[1,2] ≈ 0.8836241998111535233e-5 atol=1e-12
    @test rmatn[1,3] ≈ 0.3830834608415287707e-5 atol=1e-12
    @test rmatn[2,1] ≈ -0.8836086334870740138e-5 atol=1e-12
    @test rmatn[2,2] ≈ 0.9999999991354657474 atol=1e-12
    @test rmatn[2,3] ≈ -0.4063240188248455065e-4 atol=1e-12
    @test rmatn[3,1] ≈ -0.3831193642839398128e-5 atol=1e-12
    @test rmatn[3,2] ≈ 0.4063236803101479770e-4 atol=1e-12
    @test rmatn[3,3] ≈ 0.9999999991671663114 atol=1e-12
end

@testset "numat" begin
    epsa =  0.4090789763356509900
    dpsi = -0.9630909107115582393e-5
    deps =  0.4063239174001678826e-4
    rmatn = numat(epsa, dpsi, deps)
    @test rmatn[1,1] ≈ 0.9999999999536227949 atol=1e-12
    @test rmatn[1,2] ≈ 0.8836239320236250577e-5 atol=1e-12
    @test rmatn[1,3] ≈ 0.3830833447458251908e-5 atol=1e-12
    @test rmatn[2,1] ≈ -0.8836083657016688588e-5 atol=1e-12
    @test rmatn[2,2] ≈ 0.9999999991354654959 atol=1e-12
    @test rmatn[2,3] ≈ -0.4063240865361857698e-4 atol=1e-12
    @test rmatn[3,1] ≈ -0.3831192481833385226e-5 atol=1e-12
    @test rmatn[3,2] ≈ 0.4063237480216934159e-4 atol=1e-12
    @test rmatn[3,3] ≈ 0.9999999991671660407 atol=1e-12
end

@testset "nut00a" begin
    dpsi, deps = nut00a(2400000.5, 53736.0)
    @test dpsi ≈ -0.9630909107115518431e-5 atol=1e-13
    @test deps ≈ 0.4063239174001678710e-4 atol=1e-13
end

@testset "nut00b" begin
    dpsi, deps = nut00b(2400000.5, 53736.0)
    @test dpsi ≈ -0.9632552291148362783e-5 atol=1e-13
    @test deps ≈ 0.4063197106621159367e-4 atol=1e-13
end

@testset "nut06a" begin
    dpsi, deps = nut06a(2400000.5, 53736.0)
    @test dpsi ≈ -0.9630912025820308797e-5 atol=1e-13
    @test deps ≈ 0.4063238496887249798e-4 atol=1e-13
end

@testset "nut80" begin
    dpsi, deps = nut80(2400000.5, 53736.0)
    @test dpsi ≈ -0.9643658353226563966e-5 atol=1e-13
    @test deps ≈ 0.4060051006879713322e-4 atol=1e-13
end

@testset "nutm80" begin
    rmatn = nutm80(2400000.5, 53736.)
    @test rmatn[1,1] ≈ 0.9999999999534999268 atol=1e-12
    @test rmatn[1,2] ≈ 0.8847935789636432161e-5 atol=1e-12
    @test rmatn[1,3] ≈ 0.3835906502164019142e-5 atol=1e-12
    @test rmatn[2,1] ≈ -0.8847780042583435924e-5 atol=1e-12
    @test rmatn[2,2] ≈ 0.9999999991366569963 atol=1e-12
    @test rmatn[2,3] ≈ -0.4060052702727130809e-4 atol=1e-12
    @test rmatn[3,1] ≈ -0.3836265729708478796e-5 atol=1e-12
    @test rmatn[3,2] ≈ 0.4060049308612638555e-4 atol=1e-12
    @test rmatn[3,3] ≈ 0.9999999991684415129 atol=1e-12
end


text:
using VulkanGen
const vg = VulkanGen

getfile(name) = dirname(dirname(pathof(VulkanCore))) *  "/gen/" * name
files = getfile.(["vk_common.jl", "vk_api.jl"])
api = API(files, sym -> @eval(vk, $(Meta.parse("$sym"))))

ignored_symbols = Dict(
    "structs" => [
        "VkBaseOutStructure",
        "VkBaseInStructure",
    ],
)

function filter_api(api)
    new_api = deepcopy(api)
    for (type, list) ∈ ignored_symbols, name ∈ list
        api_dict = getproperty(new_api, Symbol(type))
        haskey(api_dict, name) && delete!(api_dict, name)
    end
    new_api
end

f_api = filter_api(api)
w_api = wrap(f_api)
write(w_api, joinpath(dirname(dirname(@__DIR__)), "Vulkan", "generated", "_wrapped_api.jl"))
# write(w_api, joinpath(dirname(dirname(@__DIR__)), "Vulkan", "generated", "wrapped_api.jl"))


text:
using Test
using ProgressMeter
using LinearAlgebra
using Statistics
using SparseGaussianProcesses
# Random.seed!(0)

onfail(f, _::Test.Pass) = nothing
onfail(f, _::Tuple{Test.Fail,<:Any}) = f()

@testset "SparseGaussianProcesses" begin

  @testset "hyperprior" begin
    mean = [0.0,1.0]
    stddev = [1.0,1.0]
    hp = NormalHyperprior(mean, stddev)

    @test hp(mean) isa AbstractVector
    @test isapprox.(hp(mean), 0) |> all
    @test (hp(mean .+ 1) .> 0) |> all
  end

  @testset "kernel" begin
    for d in 1:2
      k = SquaredExponentialKernel(d)
      x = randn(d, 1)
      y = randn(d, 2)
      z = randn(d, 3)

      @test k.dims == (d,1)
      @test k(x,y) isa AbstractMatrix
      @test size(k(x,x)) == (1,1)
      @test (k(x,x) .== 1.0) |> all
      @test size(k(y,z)) == (2,3)
      @test (k(y,y)[1:1,1:1] .> k(y,y)[1:1,2:2]) |> all
      @test (k(y,y)[2:2,2:2] .> k(y,y)[2:2,1:1]) |> all
      @test isapprox.(k(z,z)', k(z,z)) |> all
      @test isapprox.(k(y .+ 1,y), [exp(-sqrt(d)^2) exp(-sum((y[:,2] .- (y[:,1] .+ 1)).^2));
                                     exp(-sum(((y[:,2] .+ 1) .- y[:,1]).^2)) exp(-sqrt(d)^2)]) |> all
    end
  end

  @testset "prior_basis" begin
    first_test = true
    for id in 1:2
      for od in 1:1
        scales = [0.1,1.,5.]
        for var in scales
          for ls in Iterators.product([scales for i in 1:id]...)
            l = 4096
            xd = 16
            s = 128
            t = 100

            k = SquaredExponentialKernel(id)
            k.log_variance .= log.(var)
            k.log_length_scales .= log.(ls)
            rf = EuclideanRandomFeatures(k, l)

            x = randn(id, xd); # sort!(x; dims=2)
            rf.weights = randn(l,s)
            K = zeros(xd, xd)

            @test size(rf(x,k)) == (1,xd,s)

            first_test ? first_test = false : begin @test_skip true; @info "Skipping slow RFF test"; continue end

            @showprogress for i in 1:t
              rf.weights = randn(l,s)
              out = rf(x,k)
              K .+= cov(reshape(out, (xd,s))')
            end
            K ./= t

            onfail(@test ((abs.(k(x,x) .- K) ./ var) .< 0.075) |> all) do
              @show od id var ls l xd s t;
              display(x); println("\n")
              display(k(x,x)); println("\n")
              display(K); println("\n")
              display(abs.(k(x,x) .- K) ./ var); println("\n")
              @show maximum(abs.(k(x,x) .- K) ./ var)
            end
          end
        end
      end
    end
  end

  @testset "inducing" begin
    for id in 1:2
      for od in 1:1
        m=10
        k = SquaredExponentialKernel(id)
        ip = MarginalInducingPoints(k, m)

        (z,mu,U,V) = ip()

        @test length(mu) == od*m
        @test size(U) == (m,m)
        @test size(V) == (m,m)

        @test isapprox.(mu, 0) |> all
        @test isapprox(U' * U, k(z,z) + V'*V)
        @test isapprox(V' * V, (0.00001*I(m)))

        x = randn(od, m)
        ip(x,k)

        (z,mu,U,V) = ip()

        @test size(z) == size(x)
        @test length(mu) == od*m
        @test size(U) == (m,m)
        @test size(V) == (m,m)

        @test isapprox(z, x)
        @test isapprox.(mu, 0) |> all
        @test isapprox(U' * U, k(z,z) + V'*V)
        @test isapprox(V' * V, (0.00001*I(m)))

        x = randn(od, m+1)
        ip(x,k)

        (z,mu,U,V) = ip()

        @test size(z) == size(x)
        @test length(mu) == od*(m+1)
        @test size(U) == (m+1,m+1)
        @test size(V) == (m+1,m+1)

        @test isapprox(z, x)
        @test isapprox.(mu, 0) |> all
        @test isapprox(U' * U, k(z,z) + V'*V)
        @test isapprox(V' * V, (0.00001*I(m+1)))
      end
    end
  end

  @testset "gp" begin
    for id in 1:1
      for od in 1:1
        scales = [0.5,1.,1.5]
        for var in scales
          for ls in Iterators.product([scales for i in 1:id]...)
            l = 4096
            zd = 8
            xd = 128
            ns = 16384

            k = SquaredExponentialKernel(id)
            k.log_variance .= log.(var)
            k.log_length_scales .= log.(ls)
            gp = SparseGaussianProcess(k)
            gp.prior_basis = EuclideanRandomFeatures(k, l)
            gp.inducing_points = MarginalInducingPoints(k, zd)

            x = 2 .* randn(id,xd); # sort!(x; dims=2)
            z = randn(id,zd)
            u = mapslices(x -> sum(sin.(x)), z; dims=1)
            gp.inducing_points(z,k)
            gp.inducing_points.log_jitter .= log.(0.01)./2
            (_,mu,U,V) = gp.inducing_points()
            W = V' * V

            @test size(gp(z)) == (od,zd,1)

            Q = cholesky(k(z,z) + W)
            mu .= Q \ vec(u)
            K_zz = k(z,z) - k(z,z) * (Q \ k(z,z))
            LinearAlgebra.copytri!(K_zz,'U')
            gp.inducing_points.covariance_triangle .= 0
            gp.inducing_points.covariance_triangle[diagind(gp.inducing_points.covariance_triangle)] .= log(0)

            mu_x = k(x,z) * (Q \ vec(u))
            K_xx = k(x,x) - k(x,z) * (Q \ k(z,x))

            rand!(gp; num_samples = ns)
            f = gp(x)[1,:,:]
            mu_f = mean(f; dims=2)
            K_ff = cov(f')

            onfail(@test (abs.(mu_f .- mu_x) .< 0.025) |> all) do
              @show od id var ls zd z u;
              display(mu_f'); println("\n")
              display(mu_x'); println("\n")
              @show maximum(abs.(mu_f .- mu_x))
            end

            onfail(@test ((abs.(K_xx .- K_ff) ./ var) .< 0.075) |> all) do
              @show od id var ls zd z u;
              display(x); println("\n")
              display(K_xx); println("\n")
              display(K_ff); println("\n")
              display(abs.(K_xx .- K_ff) ./ var); println("\n")
              @show maximum(abs.(K_xx .- K_ff) ./ var)
            end
          end
        end
      end
    end
  end

  @testset "loss" begin
    k = SquaredExponentialKernel(2)
    gp = SparseGaussianProcess(k)

    @test (SparseGaussianProcesses.prior_KL(gp) .>= -1e-10) |> all
    @test (SparseGaussianProcesses.prior_KL(gp) .<= 1) |> all

    gp.inducing_points.mean .= 1
    @test (SparseGaussianProcesses.prior_KL(gp) .>= 1) |> all

    gp.inducing_points.mean .= 0
    gp.inducing_points.covariance_triangle .= 0
    @test (SparseGaussianProcesses.prior_KL(gp) .>= 1) |> all

    x = randn(1,128)
    y = randn(1,128)

    @test loss(gp,x,y) isa AbstractArray

    import Flux.gradient
    import Flux.params

    g = gradient(() -> loss(gp,x,y)|>sum, params(gp))
    for p in params(gp)
      @test isfinite.(g[p]) |> all
    end
  end

  @testset "gpu" begin

  end
end


text:
using JuMP, Gurobi, OneClassSampling, SVDD, OneClassActiveLearning
include(joinpath(@__DIR__, "config.jl"))

experiment_name="warmup"

data_dirs = ["Parkinson"]

solver = with_optimizer(Gurobi.Optimizer; OutputFlag=0, Threads=1)
quality_metrics = Dict(:mcc => matthews_corr, :kappa => cohens_kappa, :f1 => f1_score)

threshold_strategies = [
    GroundTruthThresholdStrategy(),
    OutlierPercentageThresholdStrategy(),
    OutlierPercentageThresholdStrategy(0.05),
]

model_init_strategy = Dict{Symbol, Any}(
                            :svdd => Dict{Symbol, Any}(
                                :ROT_CTax => SVDD.SimpleCombinedStrategy(SVDD.RuleOfThumbScott(), SVDD.BoundedTaxErrorEstimate(0.05, 0.02, 0.98)),
                                :ROT_C1 => SVDD.SimpleCombinedStrategy(SVDD.RuleOfThumbScott(), SVDD.FixedCStrategy(1)),
                            ),
                            :sample_qe => Dict{Symbol, Any}(
                                :threshold_strategies => threshold_strategies
                        ))
NUM_RERUNS = 1
KDE_GAMMA = :scott

baselines = [
    RandomRatioSampler(0.1),
    HDS(0.1, KDE_GAMMA, OutlierPercentageThresholdStrategy())
]

competitors = [
    BPS(nothing, 0.05),
    DAEDS(30, 0.1, 0.3),
    DBSRSVDD(7, 0.3),
    HSR(5, 0.01),
    IESRSVDD(KDE_GAMMA, 0.5),
    FBPE(360),
    KFNCBD(100, 0.2),
    NDPSR(20, 10)
]

density_strategies = [
    RAPID(OutlierPercentageThresholdStrategy(), KDE_GAMMA),
]

sampling_strategies = Vector{Dict{Symbol, Any}}()

for b in baselines
    push!(sampling_strategies, Dict{Symbol, Any}(
        :method => b,
        :seed => 0
    ))
end

for c in competitors
    push!(sampling_strategies, Dict{Symbol, Any}(
        :method => c,
        :seed => 0
    ))
end

for t in threshold_strategies
    push!(sampling_strategies, Dict{Symbol, Any}(
        :method => PreFilteringWrapper(HSR(5, 0.1), KDE_GAMMA, t),
        :seed => 0
    ))
end

for d in density_strategies
    push!(sampling_strategies, Dict{Symbol, Any}(
        :method => d,
        :seed => 0
    ))
end

max_stars_repo_name: englhardt/ocs-evaluation

max_issues_repo_name: englhardt/ocs-evaluation

max_forks_repo_name: englhardt/ocs-evaluation


text:
include("euler/euler.jl")
using .Primes: prime_numbers
using .Numbers: multiplicative_order

function compute(n::Int)::Int
    max_cycle, result = 0, 0
    for prime ∈ prime_numbers(n)[3:end]
        cycle = multiplicative_order(10, prime)
        if cycle > max_cycle max_cycle, result = cycle, prime end
    end
    return result
end


text:
using OrdinaryDiffEq, DiffEqCallbacks, LinearAlgebra

# https://github.com/SciML/DiffEqBase.jl/issues/564 : Fixed
gravity = 9.8
stiffness = 500
equilibrium_length = 1
T = 5.0

f(u, p, t) = begin
    x1, x2, dx1, dx2 = u
    length = abs(x2 - x1)
    spring_force = stiffness * (equilibrium_length - length)
    ddx1 = -gravity - spring_force
    ddx2 = -gravity + spring_force
    if x1 <= 0
        ddx1 = max(0, ddx1)
    end
    if x2 <= 0
        ddx2 = max(0, ddx2)
    end
    [dx1, dx2, ddx1, ddx2]
end

sol = solve(
    ODEProblem(f, [5.0, 6.0, 0.0, 0.0], (0.0, T)),
    # Euler(),
    # dt=0.005,
    Rosenbrock23(),
    callback = ContinuousCallback((u, _, _) -> u[1], (integrator) -> (integrator.u[1] = 0; integrator.u[3] = 0)),
    # callback = ContinuousCallback((u, _, _) -> u[1], (integrator) -> (integrator.u[3] = 0)),
    reltol = 1e-3,
    abstol = 1e-3,
)

@show sol.destats

# https://github.com/SciML/DiffEqBase.jl/issues/553 : Floating point issue is resolved but some other error occurss
function model(du, u, p, t)
    du[1] = 0.
    for i in 2:(length(du)-1)
        du[i] = p[i] * (u[i-1] - u[i])
    end
    du[end] = p[end] * (p[1] * u[end-1] - u[end])
    return nothing
end

perror = [1.0, 0.02222434508140991, 0.017030281542289794, 0.015917011145559996, 0.1608874463597176, 0.13128016561792297, 0.11056834258380167, 0.5222141958458832, 1.0711942201995688, 0.2672878398678257, 8.900058706990183, 0.010760065201065117, 0.016319181296867765, 2.2693845639611925, 0.2152216345154439, 0.029186712540925457, 0.21419429135100806, 0.029177617589788596, 0.03064986043089549, 0.023280222517122397, 6.931251277770224]
y_max = 0.002604806609572015
u0 = [1, zeros(length(perror) - 1)...]
tspan = (0., 5000.)

condition(u, t, i) = (t == 1.)
affect!(i) = (i.u[1] = 0.)

condition2(u, t, i) = u[end] - y_max / 2.
t_half_1 = 0.
affect2!(i) = (t_half_1 = i.t)

prob = ODEProblem(model, u0, tspan, perror)
sol = solve(
    prob,
    Rosenbrock23();
    callback=CallbackSet(
        PositiveDomain(),
        DiscreteCallback(condition, affect!),
        ContinuousCallback(condition2, affect2!, terminate!),
    ),
    tstops = [1.],
    force_dtmin=true,
)

# https://github.com/SciML/DiffEqBase.jl/issues/515 : Fixed

using StaticArrays
using MultiScaleArrays

t_last = 0.0
function attactor(du, u, p, t)
    α, β = p
    n = length(u.nodes)
    return for k in 1:n
        du.nodes[k] = zero(du.nodes[k])
        for j in 1:n
            if (k == j)
                du.nodes[k] .+=
                    [u.nodes[k][3], u.nodes[k][4], -β * u.nodes[k][3], -β * u.nodes[k][4]]
            else
                du.nodes[k][3:4] .+= α * (u.nodes[j][1:2] - u.nodes[k][1:2])
            end
        end
    end
end

struct Thingy{B} <: AbstractMultiScaleArrayLeaf{B}
    values::Vector{B}
end

struct PhysicsLaw{T<:AbstractMultiScaleArray,B<:Number} <: AbstractMultiScaleArrayHead{B}
    nodes::Vector{T}
    values::Vector{B}
    end_idxs::Vector{Int}
end

Newton = construct(
    PhysicsLaw,
    [
        Thingy([-700.0, -350.0, 0.0, 0.0]),
        Thingy([-550.0, -150.0, 0.0, 0.0]),
        Thingy([-600.,15.,0.,10.]),
        Thingy([200.,-200.,5.,-5.])
    ],
)

parameters = [1e-2, 0.06]

function condition(out, u, t, integrator)
    i = 0
    n = length(u.nodes)
    for k in 1:n
        for l in (k+1):n
            i += 1
            out[i] = sum(abs2, u.nodes[k][1:2] .- u.nodes[l][1:2]) - 10000
        end
    end
end

function affect!(integrator, idx)
    i = 0
    u = integrator.u
    n = length(u.nodes)
    return for k in 1:n
        for l in (k+1):n
            i += 1
            if idx == i
                x₁ = u.nodes[k][1:2]
                v₁ = u.nodes[k][3:4]
                x₂ = u.nodes[l][1:2]
                v₂ = u.nodes[l][3:4]
                # https://stackoverflow.com/a/35212639
                v₁ = (
                    v₁ -
                    2 / (1 + 1) * (dot(v₁ - v₂, x₁ - x₂) / sum(abs2, x₁ - x₂) * (x₁ - x₂))
                )
                v₂ =
                    -(
                        v₂ -
                        2 / (1 + 1) *
                        (dot(v₂ - v₁, x₂ - x₁) / sum(abs2, x₂ - x₁) * (x₂ - x₁))
                    )

                println("Collision handeled.")

                m = (x₁ + x₂) / 2

                u.nodes[k][3:4] .= v₁
                u.nodes[l][3:4] .= v₂

                set_u!(integrator, u)
                println(sqrt(sum(abs2, x₁ .- x₂)) - 100, ":", v₁ ./ v₂)
                println(norm(v₁),":",norm(v₂),":",integrator.t,":",integrator.t - t_last)
                global t_last = integrator.t
                break
            end
        end
    end
end

cback = VectorContinuousCallback(
    condition,
    affect!,
    (x -> Int(((x - 1) * x) / 2))(length(Newton.nodes)), #Fix  was here
)

problemp = ODEProblem(attactor, Newton, (0.0, Inf), parameters)

world = init(problemp, AutoTsit5(Rosenbrock23()); save_everystep = false, callback = cback)

dt = 0.2

for i=1:1000
    step!(world, dt)
end

## https://github.com/SciML/OrdinaryDiffEq.jl/issues/1528

function f!(out, u, p, t)
    out[1] = 0
    out[2] = u[3]
    out[3] = -1.0*(u[2] - u[1])
end
u0 = [0, 0, 1.0]
function cond!(out, u, t, i)
    out[1] = u[3]
    nothing
end
function affect!(int, idx)
    terminate!(int)
end
cb = VectorContinuousCallback(cond!, affect!, nothing, 1)

u0 = [ 0.0, 0.0, 1.0 ]
prob = ODEProblem(f!, u0, (0.0, 10.0); callback=cb)
soln = solve(prob, Tsit5())
@test soln.t[end] ≈ 4.712347213360699


text:
struct Tag{Parent} end

@inline _find_dual(tag::T, l, i) where {T} = 0
@inline _find_dual(tag::T, l, i, x::Dual{T}, xs...) where {T} = i
@inline _find_dual(tag::T, l, i, x, xs...) where {T} = _find_dual(tag, l, i+1, xs...)

@inline oldertag(::Tag{Tag{T}}) where T = Tag{T}()
@inline oldertag(::Tag{T}) where T = nothing

@inline function find_dual(T::TT, xs...) where {TT<:Tag}
    _find_dual(T, length(xs), 1, xs...)
end

# Base case where T is not a Tag
@inline find_dual(T, xs...) = 0


text:
import Dierckx
import Calculus

using Formatting
using Crayons.Box

struct SplineFunction <: Function
    spline::Dierckx.Spline1D
    function_object::Function
    firstderivative_object::Function
    a::Float64
    b::Float64
    y_min::Float64
    y_max::Float64
    Δx_min::Float64
    Δx_max::Float64
    Δy_min::Float64
    Δy_max::Float64

    function SplineFunction(X::Vector{Float64}, Y::Vector{Float64})
        @assert(!isempty(X), "X must non be empty.")
        @assert(!isempty(Y), "Y must non be empty.")
        @assert(length(X) > 2, "The length of X must greater than 2.")
        @assert(length(Y) > 2, "The length of Y must greater than 2.")
        @assert(issorted(X), "The values in X must be sorted in the ascending order.")
        @assert(length(X) == length(Y), "X and Y must have equal sizes.")

        a = X[1]; b = X[end]
        y_min = minimum(Y); y_max = maximum(Y)
        Δx_min = minimum(X[2:end] .- X[1:end-1]); Δx_max = maximum(X[2:end] .- X[1:end-1])
        Δy_min = minimum(abs.(Y[2:end] .- Y[1:end-1])); Δy_max = maximum(abs.(Y[2:end] .- Y[1:end-1]))
        spline = Dierckx.Spline1D(X, Y; w=ones(length(X)), k=2, bc="extrapolate", s=0.0)

        this = new(spline,
            x::Number -> Dierckx.evaluate(spline, x),
            x::Number -> Dierckx.derivative(spline, x),
            a, b,
            y_min, y_max,
            Δx_min, Δx_max, Δy_min, Δy_max)

        return this
    end

    function (sfunction::SplineFunction)(x::Number)
        return sfunction.function_object(x)
    end
end
function Base.show(io::IO, sfunction::SplineFunction)
    compact = get(io, :compact, false)

    a = sfunction.a; b = sfunction.b
    y_min = sfunction.y_min; y_max = sfunction.y_max
    Δx_min = sfunction.Δx_min; Δx_max = sfunction.Δx_max
    Δy_min = sfunction.Δy_min; Δy_max = sfunction.Δy_max
    spline = sfunction.spline
    order = spline.k

    sb = IOBuffer()
    if compact
        printfmt(sb, "𝐒𝐩𝐥𝐢𝐧𝐞(𝐗 ∈ [{1:.6e}, {2:.6e}], 𝐘 ∈ [{3:.6e}, {4:.6e}];\n$(' '^15)𝚶 = {5}; Δxᵐⁱⁿ = {6:.6e}, Δxᵐᵃˣ = {7:.6e})",
            a, b, y_min, y_max, order, Δx_min, Δx_max)
    else
        printfmt(sb, "𝐒𝐩𝐥𝐢𝐧𝐞(𝐗 ∈ [{1:.6e}, {2:.6e}], 𝐘 ∈ [{3:.6e}, {4:.6e}];\n$(' '^15)𝚶 = {5}; Δxᵐⁱⁿ = {6:.6e}, Δxᵐᵃˣ = {7:.6e}; |Δyᵐⁱⁿ| = {8:.6e}, |Δyᵐᵃˣ| = {9:.6e})",
            a, b, y_min, y_max, order, Δx_min, Δx_max, Δy_min, Δy_max)
    end
    print(io, LIGHT_GRAY_FG(String(take!(sb))))
end
function Base.show(io::IO, ::MIME"text/plain", sfunction::SplineFunction)
    show(io, sfunction)
end

function ⋮(sfunction::SplineFunction, x::Number)
    return sfunction.firstderivative_object(x)
end
function d2(sfunction::SplineFunction, x::Number)
    return Calculus.derivative(sfunction.firstderivative_object, x)
end
function ∫(sfunction::SplineFunction, (a, b)::Tuple{Number, Number})
    return Dierckx.integrate(sfunction.spline, a, b)
end
function ⩪(sfunction::SplineFunction)
    return Dierckx.roots(sfunction.spline)
end