text:
# Hyperparameter interface

"""
    HyperParameter{T}

A hyperparameter is any state that influences the
training and is not a parameter of the model.

Hyperparameters can be scheduled using the [`Scheduler`](#)
callback.
"""
abstract type HyperParameter{T} end

"""
    sethyperparameter!(learner, H, value)

Sets hyperparameter `H` to `value` on `learner`.
"""
function sethyperparameter! end

"""
    stateaccess(::Type{HyperParameter})

Defines what `Learner` state is accessed when calling
`sethyperparameter!` and `gethyperparameter`
"""
stateaccess(::Type{HyperParameter}) = ()

# Implementations

abstract type LearningRate <: HyperParameter{Float64} end

stateaccess(::Type{LearningRate}) = (optimizer = Write(),)
sethyperparameter!(learner, ::Type{LearningRate}, value) =
    setlearningrate!(learner.optimizer, value)

function setlearningrate!(optimizer, value)
    optimizer.eta = value
end


text:
using CarrierCapture
using DataFrames
using Test

include("test_tise.jl")
include("test_potential.jl")
include("test_rate.jl")
include("test_transfer_coord.jl")

ext: jl


text:
using GPUifyLoops
using Test
using InteractiveUtils

function kernel(A)
    @loop for i in (1:size(A,1);
                    threadIdx().x)
        A[i] = 2*A[i]
    end
    @synchronize
end

f2(x) = sin(x)
f3(x) = 1 + f2(x)
f4(x) = x^1.2
f5(x) = x^3

function kernel2!(A, B, h)
    @inbounds @loop for i in (1:size(A,1); threadIdx().x)
        A[i] = h(B[i])
    end
    nothing
end

@testset "Array" begin
    data = Array{Float32}(undef, 1024)
    kernel(data)
end

@static if Base.find_package("CuArrays") !== nothing
    using CuArrays
    using CUDAnative

    function kernel(A::CuArray)
        @launch CUDA() threads=length(A) kernel(A)
    end

    @testset "CuArray" begin
        data = CuArray{Float32}(undef, 1024)
        kernel(data)
    end

    @testset "contextualize" begin
        f(x) = 2*x
        g(x) = GPUifyLoops.contextualize(f)(x)
        @test g(3.0) == 6.0
        f(x) = 3*x

        if VERSION < v"1.3"
            @test_broken g(3.0) == 9.0
        else
            @test g(3.0) == 9.0
        end

        f1(x) = (sin(1.0 + x); return nothing)
        g1(x) = GPUifyLoops.contextualize(f1)(x)
        asm = sprint(io->CUDAnative.code_llvm(io, g1, Tuple{Float64}, kernel=true,
                                              optimize=false, dump_module=true))
        @test occursin(r"call .* double @__nv_sin", asm)
        @test occursin("fadd contract double", asm)

        @testset "don't overdub intrinsics" begin
            global simple_kernel, kernel
            simple_kernel(A, x) = (A[1] = 1 + x; return nothing)
            kernel(A, x) = GPUifyLoops.contextualize(simple_kernel)(A, x)
            CI, ret = CUDAnative.code_typed(kernel, Tuple{CUDAnative.CuDeviceArray{Int64,1, CUDAnative.AS.Global}, Int64}, debuginfo=:source)[1]

            intrinsics = findall(CI.code) do stmt
                if Base.Meta.isexpr(stmt, :call)
                    f = stmt.args[1]
                    if f isa GlobalRef
                        f = getfield(f.mod, f.name)
                        return f isa Core.IntrinsicFunction || f isa Core.Builtin
                    end
                end
                return false
            end

            for i in intrinsics
                lineinfo = CI.linetable[CI.codelocs[i]]
                @test !(lineinfo.method === :call ||
                        lineinfo.file === Symbol("context.jl"))
            end
        end

        begin
            global kernel2!
            data = rand(Float32, 1024)
            fdata = similar(data)

            kernel2!(fdata, data, f3)
            @test f3.(data) ≈ fdata

            kernel2!(fdata, data, f4)
            @test f4.(data) ≈ fdata

            kernel2!(fdata, data, f5)
            @test f5.(data) ≈ fdata

            function kernel2!(A::CuArray, B::CuArray, f)
                @launch CUDA() threads=length(A) kernel2!(A, B, f)
            end

            cudata = CuArray(data)
            cufdata = similar(cudata)

            kernel2!(cufdata, cudata, f3)
            @test f3.(data) ≈ cufdata

            kernel2!(cufdata, cudata, f4)
            @test f4.(data) ≈ cufdata

            kernel2!(cufdata, cudata, f5)
            @test f5.(data) ≈ cufdata
        end
    end

    @testset "read only load" begin
        global read_only_kernel, con_read_only_kernel
        read_only_kernel(A, B) = (B[1] = read_only_load(A, 1); return nothing)
        con_read_only_kernel(A, B) = GPUifyLoops.contextualize(read_only_kernel)(A, B)
        tt = Tuple{CUDAnative.CuDeviceArray{Float32, 1, CUDAnative.AS.Global},
                   CUDAnative.CuDeviceArray{Float32, 1, CUDAnative.AS.Global}}

        asm = sprint(io->CUDAnative.code_llvm(io, con_read_only_kernel, tt, kernel=true,
                                              optimize=false, dump_module=true))

        @test occursin("llvm.nvvm.ldg.global.f.f32", asm)

        A = rand(Float32, 1024)
        B = similar(A)
        @launch CPU() threads=length(A) read_only_kernel(A, B)
        @test A[1] == B[1]

        cuA = CuArray(A)
        cuB = similar(cuA)

        @launch CUDA() threads=length(cuA) read_only_kernel(cuA, cuB)
        @test A[1] == Array(cuB)[1]
    end
end

function kernel3!(A)
    s1 = @shmem eltype(A) (1024,)
    s2 = @shmem eltype(A) (1024,)

    @loop for i in (1:size(A,1); threadIdx().x)
        s1[i] = 2*A[i]
        s2[i] = 3*A[i]
    end
    @synchronize
    @loop for i in (1:size(A,1); threadIdx().x)
        A[i] = s1[i]
    end
    nothing
end

let
    function ker1!(::Val{Nq}) where Nq
      s_x = @shmem Float32 Nq
    end

    CI, rt = @code_typed ker1!(Val(10))
    @test Base.isconcretetype(rt)
end

@testset "shared memory" begin
    data = rand(Float32, 1024)
    cpudata = copy(data)

    @launch CPU() kernel3!(cpudata)
    @test cpudata ≈ 2 .* data

    @static if Base.find_package("CuArrays") !== nothing
        using CuArrays
        using CUDAnative

        cudata = CuArray(data)
        @launch CUDA() threads=length(cudata) kernel3!(cudata)
        @test Array(cudata) ≈ 2 .* data
    end
end

# Scratch arrays

function kernel_scratch(A, ::Val{N}) where N
    a = @scratch eltype(A) (N, N) 2
    b = @scratch eltype(A) (2, N, N) 2
    @loop for j in (1:size(A,2); threadIdx().y)
        @loop for i in (1:size(A,1); threadIdx().x)
            a[i, j] = A[i, j]
            b[1, i, j] = -A[i, j]
            b[2, i, j] = 2a[i, j]
        end
    end
end

function f1()
    A = @scratch Int64 (12, 3) 2
    @test A.data isa GPUifyLoops.MArray
    @test size(A.data) == (1,)
end

function f2()
    A = @scratch Int64 (12, 3) 1
    @test A.data isa GPUifyLoops.MArray
end

function f3()
    A = @scratch Int64 (12, 3) 1
    @test A isa GPUifyLoops.MArray
end

@testset "Scratch Arrays" begin
    contextualize(f1)()
    contextualize(f2)()
    f3()
    N = 10
    A = rand(N, N)
    @launch CPU() kernel_scratch(A, Val(N))

    @static if Base.find_package("CuArrays") !== nothing
        using CuArrays

        d_A = CuArray(A)
        @launch CUDA() kernel_scratch(d_A, Val(N))
    end
end

@testset "Loopinfo" begin
    # Right now test that we don't break things
    # Should probably test that codegen is correct.
    f(N) = @unroll 2 for i in 1:N
        @show i
    end
    f(10)

    f() = @unroll for i in 1:10
        @show i
    end
    f()
end

using StaticArrays
function kernel_MArray!(A)
  l_F = MArray{Tuple{3, 3}, eltype(A)}(undef)
  @inbounds for j = 1:3, i = 1:3
    l_F[i, j] = A[i, j]
  end
  nothing
end
function kernel_similar_MArray!(A)
  l_F = MArray{Tuple{3, 3}, eltype(A)}(undef)
  l_G = similar(l_F, Size(2,2))
  @inbounds for j = 1:2, i = 1:2
    l_G[i, j] = A[i, j]
  end
  nothing
end

@testset "StaticArrays" begin
  @static if Base.find_package("CuArrays") !== nothing
    using CuArrays
    using CUDAnative

    A = CuArray(rand(3,3))
    @launch CUDA() threads=(3,3) kernel_MArray!(A)
    @launch CUDA() threads=(3,3) kernel_similar_MArray!(A)
  end

  A = rand(3,3)
  @launch CPU() threads=(3,3) kernel_MArray!(A)
  @launch CPU() threads=(3,3) kernel_similar_MArray!(A)
end

# For the next three test we check that the _apply got inlined correctly
Base.@pure pure_f1(x, y) = x + y
let
    CI, rt = @code_typed GPUifyLoops.Cassette.overdub(GPUifyLoops.ctx, pure_f1, 1, 2)
    expr = CI.code[end-1]
    @test expr.head === :call || expr.head === :invoke
    @test expr.args[1] === GlobalRef(Base, :add_int)
end

Base.@pure pure_f2(x, ys...) = x + sum(ys)
let
    CI, rt = @code_typed GPUifyLoops.Cassette.overdub(GPUifyLoops.ctx, pure_f2, 1, 2, 3, 4)
    expr = CI.code[end-1]
    @test expr.head === :call || expr.head === :invoke
    @test expr.args[1] === GlobalRef(Base, :add_int)
end

Base.@pure pure_f3(ys...) = sum(ys)
let
    CI, rt = @code_typed GPUifyLoops.Cassette.overdub(GPUifyLoops.ctx, pure_f3, 1, 2, 3, 4)
    expr = CI.code[end-1]
    @test expr.head === :call || expr.head === :invoke
    @test expr.args[1] === GlobalRef(Base, :add_int)
end


text:
#
# script, or "quip", to run the CIRQ
#

push!(LOAD_PATH, homedir()*"/projects/QuIPS/src")

using QuIPS

M₁ = Measurement(1)

quip = [
    (:X, 1),
    (:H, 1)
]

QC = QCircuit(quip, 1)

run!(QC)

println(abs.(QC.wfn))

ψ, out = M₁(QC.wfn, QC.N)

println("Measurement Result: ", out)

max_stars_repo_name: aarontrowbridge/QuIPS

max_issues_repo_name: aarontrowbridge/QuIPS

max_forks_repo_name: aarontrowbridge/QuIPS


text:
using Pkg: Pkg
Pkg.activate(@__DIR__)
using Cascadia
using Gumbo
using CSV
using Markdown
using Downloads: download

include("piracy.jl")

function main()
    @info "starting"
    open(joinpath(dirname(@__DIR__), "README.md"), "w") do fh
        println(
            fh,
            """
            # 🐂 Hi Hi, I am Lyndon
            Feel encouraged to reach out to me. I like people.
            I am very contactable online. I am sure you can find me.

            GitHub doesn't make it easy to showcase all the projects I am involved in.
            Especially with so many being inside various github orgs.
            So I am using this profile README.md to experiment with providing something better.

            ***Note 1:** these lists are incomplete. When I find time I will remember the other 50 projects I am involved in.*\\
            ***Note 2:** this is just a list of projects I am involved in. A project listed here just means I think my involvement is in some sense significant. It doesn't mean I am running the project or even have commit rights.*
            """
        )
        write_section(fh, "main")
        println(fh)

        println(
            fh,
            """
            ## Dead-ends
            These are projects that I have spent a fair bit of time working on, but that I have now concluded are dead-ends.
            That they are not the way to continue to advance to solve this problem.
            They probably still function and work usefully.
            But they are deprecated and have other alternatives recommended.
            Not every project that I have abandonned is here, just ones that went a long way.
            You should ask me about why I think these are dead-ends, they are interesting things or they wouldn't have gotten this far.
            """
        )
        write_section(fh, "deadends")
        println(fh)

        println(
            fh,
            """

            ## Experiments and Early Ideas
            These are projects that are in an early stage, and exist more to prove a point.
            They may or may not be usable.
            I may or may not be actively working on them at this point in time.
            I might not have touched them in years.
            Regardless of this they are very cool ideas (in my very biased opinion).
            You should ask me about them.
            """
        )
        write_section(fh, "experiments")

        println(
            fh,
            """

            ## Conclusion
            I hope this has been useful to you to see what I have been up to.
            Do reach out to me, as I said, I love talking to other humans.

            You can find the script that generates this profile [here](https://github.com/oxinabox/oxinabox).
            It's pretty fun little webscrapy markdown generaty thing.
            """
        )
    end

    @info "done"
end

function write_section(fh, section_name)
    @info "gathering data: $section_name"
    repo_urls = CSV.File(joinpath(dirname(@__DIR__), "data", "$section_name.csv")).repo
    infos = asyncmap(get_info, repo_urls)
    @info "writing content: $section_name"
    for project_info in infos
        show(fh, MIME("text/markdown"), project_info)
    end
    @info "done: $section_name"
end

Base.@kwdef struct ProjectInfo
    url
    user
    name
    description
    icon
end

read_url(url) = parsehtml(String(take!(download(url, IOBuffer()))))

function get_info(url)
    doc = read_url(url)
    get_only(sel) = only(eachmatch(sel, doc.root))
    user = text(only(eachmatch(sel"[itemprop='author']", doc.root)))

    description = get_only(sel"meta[name='description']")."content"
    # strip github's cruft
    description = replace(description, r" - GitHub - .*"=>"")
    description = replace(description, r" Contribute to .*"=>"")
    description = strip(description)

    return ProjectInfo(;
        url,
        user,
        name = text(only(eachmatch(sel"[itemprop='name']", doc.root))),
        description,
        icon = get_avatar_url(user),
    )
end

const _AVATAR_URL_CACHE = Dict{String,String}()
function get_avatar_url(user)
    get!(_AVATAR_URL_CACHE, user) do
        doc = read_url(joinpath("https://github.com", user))
        eles = eachmatch(sel"img[itemprop='image']", doc.root)
        if isempty(eles)  # probably means is personal repo not org repo
            eles = eachmatch(sel"img.avatar-user", doc.root)
        end
        return first(eles)."src"
    end
end

function Base.show(io::IO, ::MIME"text/markdown", info::ProjectInfo)
    print(io, " - ")
    print(io, "  ")
    print(io, "[**$(info.user)/$(info.name)**]($(info.url)): ")
    print(io, "_$(info.description)_")
    println(io)
end

main()

max_stars_repo_name: oxinabox/oxinabox

max_issues_repo_name: oxinabox/oxinabox

max_forks_repo_name: oxinabox/oxinabox


text:
## Functions related to unstable nets (some vertices have the same equilibrium placement)

"""
    CollisionNode

Store the structure of a collision node through the subgraph `g` extracted with only the
edges bond to the vertices in the node.

The `num` field corresponds to the number of vertices in `g` that collide in the initial
graph, while the `neighs` field stores the indices of their neighbor out of the collision
site.

The colliding vertices are the `num` first vertices of `g`, the others are the neighbors.
In particular, `nv(g) == num + length(neighs)` and `g[(num+1):(nv(g))]` has no edge.
"""
struct CollisionNode
    g::PeriodicGraph{0}
    num::Int
    neighs::Vector{Int}
end
(Base.:(==))(c1::CollisionNode, c2::CollisionNode) = c1.num == c2.num && c1.g == c2.g
Base.hash(c::CollisionNode, h::UInt) = hash(c.g, hash(c.num, hash(c.neighs, h)))
Base.length(c::CollisionNode) = c.num

"""
    shrink_collisions(net::CrystalNet, collisions)

Remove all colliding vertices and replace them by one new vertex per `CollisionNode`, whose
neighbours are that of the vertices within.
"""
function shrink_collisions(net::CrystalNet{D,T}, collisions::Vector{UnitRange{Int}}) where {D,T}
    first_colliding = first(first(collisions))
    m = length(collisions)
    n = first_colliding + m - 1

    collision_vmap = collect(1:length(net.types))
    for (i, node) in enumerate(collisions)
        @toggleassert all(==(net.pge.pos[first(node)]), @view net.pge.pos[node])
        collision_vmap[node] .= first_colliding + i - 1
    end
    edgs = PeriodicEdge{D}[]
    for e in edges(net.pge.g)
        src = collision_vmap[e.src]
        dst = collision_vmap[e.dst.v]
        if src < first_colliding || dst < first_colliding
            push!(edgs, PeriodicEdge{D}(src, dst, e.dst.ofs))
        end
    end
    newgraph = PeriodicGraph{D}(edgs)

    @toggleassert nv(newgraph) == n

    newpos = net.pge.pos[1:first_colliding]
    append!(newpos, net.pge.pos[first(collisions[i])] for i in 2:m)
    newtypes = net.types[1:first_colliding-1]
    append!(newtypes, :* for i in 1:m) # type of collision nodes will be Symbol('*')

    return CrystalNet{D,T}(net.pge.cell, newtypes, newpos, newgraph, net.options)
end

function _order_collision(subgraph, subnodes)
    m = nv(subgraph)
    if m == 2 # 2
        return [1, 2]
    end

    if m == 3
        length(subnodes) == 1 && return sortperm(degree(subgraph))
        @toggleassert length(subnodes) == 2 # 1-2 or 2-1
        singled = length(subnodes[1]) == 1 ? subnodes[1][1] : subnodes[2][1]
        v = get(neighbors(subgraph, singled), 1, singled == 1 ? 2 : 1)
        remaining = 6 - singled - v
        return length(subnodes[1]) == 1 ? [singled, v, remaining] : [v, remaining, singled]
    end

    # m == 4
    degs = degree(subgraph)
    @toggleassert all(<(4), degs)
    J = sortperm(degs)
    num = length(subnodes)

    if num == 1
        if degs[J[1]] == degs[J[4]]
            if degs[1] == 1 # 1-1-1-1
                neigh1 = only(neighbors(subgraph, 1))
                other1 = ifelse(neigh1 == 2, 3, 2)
                other2 = 10 - neigh1 - other1 - 1
                @toggleassert only(neighbors(subgraph, other1)) == other2
                return [1, neigh1, other1, other2]
            end
            if degs[1] == 2 # 2-2-2-2
                neighs1 = neighbors(subgraph, 1)
                neighl = neighs1[1]
                neighr = neighs1[2]
                neigho = 10 - 1 - neighl - neighr
                @toggleassert issetequal(neighbors(subgraph, neigho), [neighl, neighr])
                return [neighl, 1, neighr, neigho]
            end
        end
        if degs[J] == [1, 1, 2, 2]
            neigha = only(neighbors(subgraph, J[1]))
            neighb = only(neighbors(subgraph, J[2]))
            @toggleassert issetequal(neighbors(subgraph, neigha), [J[1], neighb])
            return [J[1], J[2], neigha, neighb]
        end
        return J
    end

    if num == 2
        if length(subnodes[1]) == 2 # 2-2
            a, b = subnodes[1]
            rem_edge!(subgraph, a, b)
            c, d = subnodes[2]
            rem_edge!(subgraph, c, d)
            if degree(subgraph, b) < degree(subgraph, a)
                a, b = b, a
            end
            if (degree(subgraph, b) == 1 && only(neighbors(subgraph, b)) == c) ||
                degree(subgraph, d) < degree(subgraph, c)
                c, d = d, c
            end
            return [a, b, c, d]
        end
        # 1-3 or 3-1
        i1 = length(subnodes[1]) == 1 ? 1 : 2
        (t, others) = only(subnodes[i1]), subnodes[3-i1]
        u, v, w = others[sortperm(degs[others])]
        @toggleassert issetequal([u,v,w,t], 1:4)
        nts = neighbors(subgraph, t)
        dt = length(nts)
        if dt == 1
            nt = only(nts)
            if degs[v] == degs[w]
                if degs[w] == 1
                    @toggleassert degs[u] == 1
                    if nt == w
                        w, v = v, w
                    end
                    if nt == v
                        u, v = v, u
                    end
                elseif nt == v
                    w, v = v, w
                end
            end
        elseif dt == 2
            if degs[u] == degs[v]
                if degs[u] == 1
                    @toggleassert degs[w] == degs[v]+1 == 2
                    @toggleassert w ∈ nts
                    if v ∈ nts
                        u, v = v, u
                    end
                else
                    @toggleassert degs[u] == 2
                    if degs[v] == 2
                        @toggleassert degs[w] == 2
                        ww = 10 - t - nts[1] - nts[2]
                        if ww == u
                            u, w = w, u
                        elseif ww == v
                            v, w = w, v
                        end
                    end
                end
            end
        end
        return i1 == 1 ? [t, u, v, w] : [u, v, w, t]
    end

    @toggleassert num == 3 # 1-1-2 or 1-2-1 or 2-1-1
    n2 = length(subnodes[1]) == 2 ? 1 : length(subnodes[2]) == 2 ? 2 : 3
    xi, xj = subnodes[n2]
    rem_edge!(subgraph, xi, xj)
    xk = only(subnodes[n2 == 1 ? 2 : 1])
    xl = only(subnodes[n2 == 1 ? 3 : 5-n2])
    rem_edge!(subgraph, xk, xl)
    neighk = neighbors(subgraph, xk)
    if length(neighk) == 1
        if only(neighk) == xi
            xi, xj = xj, xi
        end
    else
        neighls = neighbors(subgraph, xl)
        if length(neighls) == 1 && neighls[1] == xi
            xi, xj = xj, xi
        end
    end
    return n2 == 1 ? [xi, xj, only(subnodes[2]), only(subnodes[3])] :
           n2 == 2 ? [only(subnodes[1]), xi, xj, only(subnodes[3])] :
                     [only(subnodes[1]), only(subnodes[2]), xi, xj]
end

"""
    order_collision(graph::PeriodicGraph, colliding)

Given collision node (in the form of the corresponding list of colliding vertices), find
an ordering of them which is independent of the current ordering of these vertices and of
vertices which are neither in the collision node nor any of its neighbours.
Return an empty list if it fails.

This function assumes that no vertex in the node has a neighbour in another collision node
and that there are no two representatives of the same vertex that are neighbour to some
vertices in the range.
"""
function order_collision(graph::PeriodicGraph{D}, colliding) where D
    m = length(colliding)
    neigh_per_v = [Int[] for _ in 1:m]
    subgraph = SimpleGraph(m)
    for (i, (nlist, u)) in enumerate(zip(neigh_per_v, colliding))
        for x in neighbors(graph, u)
            j = findfirst(==(x.v), colliding)
            if j isa Int
                add_edge!(subgraph, i, j)
            else
                push!(nlist, x.v)
            end
        end
        sort!(nlist)
    end
    I = sortperm(neigh_per_v)
    subnodes = Vector{Int}[[I[1]]]
    last_list = neigh_per_v[I[1]]
    for _i in 2:length(I)
        i = I[_i]
        nlist = neigh_per_v[i]
        if nlist == last_list
            push!(last(subnodes), i)
        else
            push!(subnodes, [i])
            last_list = nlist
        end
    end
    length(subnodes) == m && return I

    m > 4 && return Int[]
    return _order_collision(subgraph, subnodes)
end

function collision_utils(collisions::Vector{CollisionNode}, num_withcolliding::Int, rev_vmap=nothing)
    m = length(collisions)
    num_nocolliding = num_withcolliding - m
    n = num_nocolliding + sum(length, collisions)
    new_collisions = collisions

    if rev_vmap isa Vector{Int}
        colliding_nodes = rev_vmap[(num_nocolliding+1):num_withcolliding]
        new_collisions = [CollisionNode(collisions[i], rev_vmap) for i in sortperm(colliding_nodes)]
        # The order of collisions and is made independent of the initial representation
    end

    return n, m, num_withcolliding, num_nocolliding, new_collisions
end

function collision_utils(collisions::Vector{CollisionNode}, vmap::Vector{Int})
    num_withcolliding = length(vmap)
    rev_vmap = Vector{Int}(undef, num_withcolliding)
    for (i, j) in enumerate(vmap)
        rev_vmap[j] = i
    end
    return collision_utils(collisions, num_withcolliding, rev_vmap)
end

"""
    expand_collisions(collisions::Vector{CollisionNode}, graph::PeriodicGraph, vmap)

Expand each collision node into the appropriate number of vertices so that the resulting
graph is isomorphic to the initial one, in a manner that only depends on the current graph.
Return the resulting graph.

`vmap` is the map of vertices between `initial_graph` (with collapsed collision nodes) and
`graph`
"""
function expand_collisions(collisions::Vector{CollisionNode}, graph::PeriodicGraph{D}, vmap) where D
    n, m, num_withcolliding, num_nocolliding, collisions = collision_utils(collisions, vmap)
    # collisions are now sorted according to vmap, both in the list and for each subgraph

    @toggleassert nv(graph) == num_withcolliding
    @toggleassert length(vmap) == nv(graph)

    # We now push all the collision nodes to the end
    perm = Vector{Int}(undef, num_withcolliding)
    current_idx_colliding = 1
    for (i, j) in enumerate(vmap)
        if j > num_nocolliding # vertex i in the current graph is a collision node
            perm[num_nocolliding + current_idx_colliding] = i
            current_idx_colliding += 1
        else
            idx = i + 1 - current_idx_colliding
            perm[idx] = i
        end
    end
    graph = graph[perm]
    newgraph = graph[1:num_nocolliding] # subgraph of the non-colliding nodes keeping the order

    # At this point, the graph has all the collision nodes stacked at its end, but kept in
    # the same order as that given by the topological genome of the graph.

    add_vertices!(newgraph, n - num_nocolliding)
    offsetcounter = num_nocolliding
    for (i, node) in enumerate(collisions)
        k = length(node)
        # bonds internal to the collision site
        for j in 1:k
            for x in neighbors(node.g, j)
                x.v ≤ k || continue
                add_edge!(newgraph, offsetcounter + j, PeriodicVertex{D}(offsetcounter + x.v))
            end
        end
        # other bonds
        for (j, x) in enumerate(neighbors(graph, num_nocolliding + i))
            # x is a (non-colliding) neighbour of the collision node
            for y in neighbors(node.g, length(node) + j)
                add_edge!(newgraph, offsetcounter + y.v, x)
            end
        end
        offsetcounter += length(node)
    end

    # newgraph contains all vertices from initial_graph
    # The non-colliding ones are ordered like in graph, the others are stacked in the order
    # of their collision node in graph, each colliding subgraph sorted with order_collision.

    return newgraph
end

"""
    CollisionNode(graph::PeriodicGraph, node::UnitRange{Int}, vmap=nothing)

Return a canonical `CollisionNode` identifying the collision node.
This depends on the order of its neighbours, which is given by their order in the graph or
in `vmap` if provided.
"""
function CollisionNode(graph::PeriodicGraph, node::UnitRange{Int}, rev_vmap=nothing)
    return CollisionNode(unsorted_node(graph, node), rev_vmap)
end

function unsorted_node(graph::PeriodicGraph, node::UnitRange{Int})
    neighbors_of_node_set = Set{Int}()
    for u in node
        for x in neighbors(graph, u)
            push!(neighbors_of_node_set, x.v)
        end
    end
    setdiff!(neighbors_of_node_set, node)
    neighbors_of_node = collect(neighbors_of_node_set)

    newmap = Vector{Int}(undef, nv(graph))
    for (i, u) in enumerate(node)
        newmap[u] = i
    end
    l = length(node)
    for (i, u) in enumerate(neighbors_of_node)
        newmap[u] = l + i
    end

    edgs = PeriodicEdge{0}[]
    for u in node
        v = newmap[u]
        for x in neighbors(graph, u)
            push!(edgs, PeriodicEdge{0}(v, newmap[x.v], ()))
        end
    end

    newgraph = PeriodicGraph{0}(edgs)

    return CollisionNode(newgraph, length(node), neighbors_of_node)
end

function CollisionNode(c::CollisionNode, rev_vmap)
    sorted_neighbors = rev_vmap isa Vector{Int} ? rev_vmap[c.neighs] : c.neighs
    I = sortperm(sorted_neighbors)

    newmap = collect(1:nv(c.g))
    l = length(c)
    for (j, i) in enumerate(I)
        newmap[l+j] = l+i
    end
    sorted_neighbors = sorted_neighbors[I]

    newgraph = c.g[newmap]
    n = nv(newgraph)
    @toggleassert n == length(sorted_neighbors) + l
    neworder = order_collision(newgraph, 1:l)
    @toggleassert !isempty(neworder)
    perm = collect(1:nv(newgraph))
    perm[1:l] = neworder

    return CollisionNode(newgraph[perm], l, sorted_neighbors)
end

function collect_collisions(net::CrystalNet)
    sortedpos = issorted(net.pge.pos) ? net.pge.pos : sort(net.pge.pos)
    @toggleassert all(pos -> all(x -> 0 ≤ x < 1, pos), net.pge.pos)
    collision_sites = Vector{Int}[]
    collision_vertices = Int[]
    flag_collision = false
    for i in 2:length(sortedpos)
        if sortedpos[i-1] == sortedpos[i]
            if flag_collision
                push!(collision_sites[end], i)
                push!(collision_vertices, i)
            else
                push!(collision_sites, [i-1, i])
                push!(collision_vertices, i-1, i)
                flag_collision = true
            end
        else
            flag_collision = false
        end
    end
    return collision_sites, collision_vertices
end

"""
    collision_nodes(c::CrystalNet)

Check that the net is stable, i.e. that no two vertices have the same equilibrium placement.

A net is still considered stable if the collisions in equilibrium placement cannot lead to
different topological genomes. In practice, this happens when:
A) there is no edge between two collision sites and
B) there is no edge between a collision site and two representatives of the same vertex and
C) for each collision site either:
   - the site is made of at most 4 vertices, or
   - no 2 vertices on the site share the same exact set of neighbours out of the site.

In this case, return the list of corresponding `CollisionNode`, the list being empty if the
net is truly stable. Otherwise, return `nothing`.

Also return an updated net where the vertices in a `CollisionNode` are collapsed into a new
vertex, appearing after the non-colliding vertices.
"""
function collision_nodes(net::CrystalNet{D,T}) where {D,T}
    collision_sites, collision_vertices = collect_collisions(net)

    isempty(collision_sites) && return net, CollisionNode[] # no collision at all

    # Check that conditions A, B and C are respected
    collision_vertices_set = BitSet(collision_vertices)
    for site in collision_sites
        known_neighbors = Dict{Int,SVector{D,Int}}()
        known_nlist = Set{Vector{PeriodicVertex{D}}}()
        uniquenlist = true
        for u in site
            this_nlist = PeriodicVertex{D}[]
            for x in neighbors(net.pge.g, u)
                get!(known_neighbors, x.v, x.ofs) == x.ofs || return net, nothing # condition B
                x.v ∈ site && continue
                x.v ∈ collision_vertices_set && return net, nothing # condition A
                if uniquenlist
                    push!(this_nlist, x)
                end
            end
            if uniquenlist
                if this_nlist ∈ known_nlist
                    uniquenlist = false
                else
                    push!(known_nlist, this_nlist)
                end
            end
        end
        !uniquenlist && length(site) > 4 && return net, nothing # condition C
    end

    # Reorder the vertices of the net so that those belonging to a collision node are
    # contiguous, all collision nodes are sorted according to length and they appear after
    # non-colliding vertices.
    sort!(collision_sites; by=length)
    n = length(net.types)
    vmap = Vector{Int}(undef, n)
    idx = 1
    for i in 1:n
        if i ∉ collision_vertices_set
            vmap[idx] = i
            idx += 1
        end
    end
    @toggleassert idx == n - length(collision_vertices) + 1
    collisions = Vector{UnitRange{Int}}(undef, length(collision_sites))
    for (k, node) in enumerate(collision_sites)
        collisions[k] = idx:(idx+length(node)-1)
        for x in node
            vmap[idx] = x
            idx += 1
        end
    end
    @toggleassert idx == n+1

    newgraph = net.pge.g[vmap]
    newpos = net.pge.pos[vmap]
    ofs = newpos[1]
    offsets = Vector{SVector{D,Int}}(undef, n)
    for (i, x) in enumerate(newpos)
        offsets[i] = floor.(Int, x .- ofs)
        newpos[i] = x .- ofs .- offsets[i]
    end
    offset_representatives!(newgraph, .-offsets)
    # net.pge.pos[1] should always be [0,0,0]

    newnet = CrystalNet{D,T}(net.pge.cell, net.types[vmap], newpos, newgraph, net.options)
    newnodes = [CollisionNode(newnet.pge.g, node) for node in collisions]
    return shrink_collisions(newnet, collisions), newnodes
end


text:
abstract type AbstractFault end
abstract type AbstractActuatorFault <: AbstractFault end

# API
affect!(fault::AbstractActuatorFault, Λ, t) = error("Not implemented fault")

# FaultSet
FaultSet(args...) = AbstractFault[args...]

function Affect!(faults::Vector{AbstractActuatorFault})
    actuator_faults_groups = SplitApplyCombine.group(fault -> fault.index, faults)  # actuator_faults groups classified by fault index
    return function (Λ, t)
        for actuator_faults_group in actuator_faults_groups
            last_actuator_fault = select_last_before_t(actuator_faults_group, t)
            affect!(last_actuator_fault, Λ, t)
        end
    end
end

# select_last
"""
Select the last fault (concerning the applied time) among given faults.
"""
function select_last_before_t(faults::Vector{T} where T <: AbstractFault, t)
    fault_times_less_than_t = faults |> Map(fault -> fault.time) |> Filter(<=(t)) |> collect
    if length(fault_times_less_than_t) == 0
        # case 1: if there is no faults occured before t
        # return (t, u) -> u
        return NoFault()
    else
        # case 2: if there are some faults occured before t
        fault_last_indices = findall(fault -> fault.time == maximum(fault_times_less_than_t), faults)
        fault_last_idx = length(fault_last_indices) == 1 ? fault_last_indices[1] : error("Among given faults, more than one faults occured at the same time")
        return faults[fault_last_idx]
    end
end

struct NoFault <: AbstractFault
end
function affect!(fault::NoFault, Λ, t)
end

"""
Loss of effectiveness (LoE).
"""
struct LoE <: AbstractActuatorFault
    time::Real
    index
    level::Real
    function LoE(time, index, level)
        @assert level >= 0.0 && level <= 1.0
        new(time, index, level)
    end
end

"""
# Variables
Λ ∈ R^(n×n): effectiveness matrix
_Λ ∈ R^n: effectiveness vector
"""
function affect!(fault::LoE, Λ, t)
    @unpack time, index, level = fault
    _Λ = ones(size(diag(Λ)))
    if t >= time
        _Λ[index] = level
    end
    Λ .= Λ * Diagonal(_Λ) |> Matrix
end


text:
@testset "calibrate with distances: 3-cycle, non-identifiable example" begin
net2 = readTopology("((((D:0.1,C:0.2):1.5,(B:0.1)#H1:0.9::0.7):0.1,(#H1:0.01::0.3,A:0.3):0.8):0.1);")
# plot(net2, :RCall, useEdgeLength=true, showNodeNumber=true, showGamma=true);
dAll = pairwiseTaxonDistanceMatrix(net2, keepInternal=true)
@test dAll ≈ [0.0 .8 1.1 .1 .943 1.043 1.6 1.8 1.7
 0.8   0.0   0.3   0.9   1.263 1.363 2.4   2.6   2.5
 1.1   0.3   0.0   1.2   1.563 1.663 2.7   2.9   2.8
 0.1   0.9   1.2   0.0   0.903 1.003 1.5   1.7   1.6
 0.943 1.263 1.563 0.903 0.0   0.1   2.403 2.603 2.503
 1.043 1.363 1.663 1.003 0.1   0.0   2.503 2.703 2.603
 1.6   2.4   2.7   1.5   2.403 2.503 0.0   0.2   0.1
 1.8   2.6   2.9   1.7   2.603 2.703 0.2   0.0   0.3
 1.7   2.5   2.8   1.6   2.503 2.603 0.1   0.3   0.0]
taxa = [l.name for l in net2.leaf]
net2distances = pairwiseTaxonDistanceMatrix(net2)
@test net2distances ≈ [.0 .3 2.603 2.8
 0.3   0.0   2.703 2.9
 2.603 2.703 0.0   1.663
 2.8   2.9   1.663 0.0]
g = PhyloNetworks.pairwiseTaxonDistanceGrad(net2);
@test size(g) == (9,9,9)
# plot(net2, :RCall, useEdgeLength=true, showEdgeNumber=true);
na = getNodeAges(net2);
@test na ≈ [1.7,0.11,0.0,1.6,0.1,0.0,0.1,0.0,0.0]
g = PhyloNetworks.pairwiseTaxonDistanceGrad(net2, nodeAges=na);
@test size(g) == (9,9,9)
@test g[:,:,1] ≈ [0 1 1 1 1 1 1 1 1;1 0 0 2 1.4 1.4 2 2 2; 1 0 0 2 1.4 1.4 2 2 2
 1 2 2 0 .6 .6 0 0 0; 1 1.4 1.4 .6 0 0 .6 .6 .6; 1 1.4 1.4 .6 0 0 .6 .6 .6
 1 2 2 0 .6 .6 0 0 0; 1 2 2 0 .6 .6 0 0 0; 1 2 2 0 .6 .6 0 0 0]
# plot(net2, :RCall, useEdgeLength=true, showEdgeNumber=true);

net  = readTopology("((#H1:0.06::0.3,A:0.6):1.3,(B:0.1)#H1:0.7::0.7,(C,D):1.4);");
# same topology, unrooted, different BL, leaves ordered differently
# here: branch lengths not identifiable, even if minor fixed to 0
# plot(net, :RCall, showEdgeLength=true);
calibrateFromPairwiseDistances!(net, net2distances, taxa, #verbose=true,
  forceMinorLength0=true, ultrametric=false)
# plot(net, :RCall, showEdgeLength=true, useEdgeLength=true);
est = pairwiseTaxonDistanceMatrix(net, checkPreorder=false)
@test est ≈ [0 1.663 2.9 2.8; 1.663 0 2.703 2.603; 2.9 2.703 0 .3;
             2.8 2.603 .3 0] atol=1e-5
# deep    LN_BOBYQA: 146 evaluations
# got 0.0 at [0.00536, 1.32091, 0.14155, 0.84494, 0.2, 0.1, 1.37373] after 152 iterations (returned FTOL_REACHED)
# shallow LN_BOBYQA: 79 evaluations
# deep    LD_MMA: 47 evaluations
# got 0.0 at [0.20465, 1.03622, 0.19367, 0.77048, 0.2, 0.1, 1.45914] after 47 iterations (returned FTOL_REACHED)
# second call, starting from solution:
# got 0.0 at [0.00536, 1.32092, 0.14155, 0.84493, 0.2, 0.1, 1.37373] after 29 iterations (returned FTOL_REACHED)

net = readTopology("((#H1:0.06::0.3,A:0.6):1.3,(B:0.1)#H1:0.7::0.7,(C,D):1.4);");
calibrateFromPairwiseDistances!(net, net2distances, taxa,
  verbose=false, forceMinorLength0=false, ultrametric=false)
# plot(net, :RCall, showEdgeLength=true, useEdgeLength=true);
est = pairwiseTaxonDistanceMatrix(net)
@test est ≈ [0 1.663 2.9 2.8; 1.663 0 2.703 2.603; 2.9 2.703 0 .3;
             2.8 2.603 .3 0] atol=.00002

# Test symbol input
net = readTopology("((#H1:0.06::0.3,A:0.6):1.3,(B:0.1)#H1:0.7::0.7,(C,D):1.4);");
taxa = [Symbol(t) for t in taxa]
calibrateFromPairwiseDistances!(net, net2distances, taxa,
  verbose=false, forceMinorLength0=false, ultrametric=false)
estSymb = pairwiseTaxonDistanceMatrix(net)
@test estSymb ≈ est atol=.00002
end

@testset "calibrate with distances: 5-cycle example (when unrooted)" begin
# most lengths identifiable when we force
# minor hybrid to length 0 and ultrametric network
# not identifiable otherwise
net = readTopology("((Ag:3.0,(#H1:0.0::0.2,Ak:2.5):0.5):0.5,(((((Az:0.2,Ag2:0.2):1.3,As:1.5):1.0)#H1:0.3::0.8,Ap:2.8):0.2,Ar:3.0):0.5);");
#plot(net, :RCall, useEdgeLength=true, showNodeNumber=true, showGamma=true);
taxa = [l.name for l in net.leaf];
netdist = pairwiseTaxonDistanceMatrix(net);
@test getNodeAges(net) ≈ [3.5,3,0,2.8,0,3,2.5,0,2.5,1.5,0,.2,0,0,0]
net = readTopology("((((((Ag2,Az),As))#H1:::0.8,Ap),Ar),(Ag,(#H1:::0.2,Ak)));");
# same topology, no BL, leaves ordered differently
#plot(net, :RCall, showGamma=true);
o = [4,3,5,6,7,1,2]; # to get leaves in same order as in taxa. inverse: [6,7,2,1,3,4,5]
calibrateFromPairwiseDistances!(net, netdist, taxa,
  forceMinorLength0=true, ultrametric=false)
# got 0.0 at [0.19998, 0.19998, 1.30003, 1.5, 0.6802, 0.69964, 2.79995, 0.20002, 2.99995, 0.50006, 2.9999, 2.49952, 0.50049, 0.50006] after 1000 iterations (returned MAXEVAL_REACHED)
# plot(net, :RCall, showEdgeLength=true, useEdgeLength=true);
est = pairwiseTaxonDistanceMatrix(net, checkPreorder=false)
@test est ≈ netdist[o,o] atol=.0005
# most edge lengths are correct, except for 5 edges:
# - the 2 edges at the root: expected (the network should be unrooted)
#   when the optimization uses LN_BOBYQA
# - the 3 edges connecting to the hybrid node: lack of identifiability?

net = readTopology("((((((Ag2,Az),As))#H1:::0.8,Ap),Ar),(Ag,(#H1:::0.2,Ak)));");
# same topology, no BL, leaf ordered differently
for e in net.edge e.length=1.0; end
preorder!(net)
na = getNodeAges(net);
@test na ≈ [6.,1,4,0,0,5,0,4,0,3,2,0,1,0,0]
#plot(net, :RCall, useEdgeLength=true, showEdgeNumber=true);
pairwiseTaxonDistanceMatrix(net, nodeAges=na); # update edge lengths in net
# plot(net, :RCall, showEdgeLength=true, useEdgeLength=true);
@test [e.length for e in net.edge] ≈ [1.,1,1,2,1,1,4,1,5,1,1,1,4,-3,5]
g = PhyloNetworks.pairwiseTaxonDistanceGrad(net, nodeAges=na);
#plot(net, :RCall, useEdgeLength=true, showEdgeNumber=true);
net.edge[11].length = 4.; # to respect constraints
net.edge[14].length = 0.;
net.edge[15].length = 2.;
net1 = deepcopy(net);
fmin, na = calibrateFromPairwiseDistances!(net1, netdist, taxa, #verbose=true,
  forceMinorLength0=false, ultrametric=true)
# got 0.0 at [3.5, 3.0, 2.50001, 3.0, 2.8, 1.95797, 1.5, 0.2] after 58 iterations (returned FTOL_REACHED)
est = pairwiseTaxonDistanceMatrix(net1)
@test est ≈ netdist[o,o] atol=.0005
# plot(net1, :RCall, useEdgeLength=true, showEdgeLength=true);
# still not identifiable: all good but the 3 edges adjacent to hybrid node
# only the age of the hybrid node is wrong: zipped down nearer tips

calibrateFromPairwiseDistances!(net, netdist, taxa, #verbose=true,
  forceMinorLength0=true, ultrametric=true, NLoptMethod=:LN_COBYLA)
est = pairwiseTaxonDistanceMatrix(net)
@test est ≈ netdist[o,o] atol=.06
# hard problem? reached max # iterations
calibrateFromPairwiseDistances!(net, netdist, taxa, #verbose=true,
  forceMinorLength0=true, ultrametric=true)
est = pairwiseTaxonDistanceMatrix(net)
@test est ≈ netdist[o,o] atol=.06
# plot(net, :RCall, showEdgeLength=true, useEdgeLength=true);
est = [e.length for e in net.edge]
el = [.2,.2,1.3,1.5,1,.3,2.8,.2,3,.5,3,0,2.5,.5,.5]
@test est ≈ el atol=0.06
end


text:
using ForwardDiff
using DiffResults

struct GradObj{F,G}
    f::F
    cfg::G
end

function GradObj(f::Function, x::AbstractArray; chunksize::Int=5)
    cfg = ForwardDiff.GradientConfig(f, x, ForwardDiff.Chunk{chunksize}())
    return GradObj(f, cfg)
end

struct LogProbDeriv{T,C,F,G}
    transform::T
    cm::C
    ℓπ::F
    ∇ℓπ::G
end
function LogProbDeriv(cm::Soss.ConditionalModel, chunksize::Int)
    transform = xform(cm)
    function ℓπ(x)
        p, logjac = TV.transform_and_logjac(transform, x)
        return logdensity(cm, p) + logjac
    end
    cfg = ForwardDiff.GradientConfig(ℓπ, rand(transform.dimension), ForwardDiff.Chunk{chunksize}())

    function ∇ℓπ(x)
        ForwardDiff.gradient(ℓπ, x, cfg)
    end
    ∇ℓπ(rand(transform.dimension))
    @time ∇ℓπ(rand(transform.dimension))
    return LogProbDeriv(transform, cm, ℓπ, ∇ℓπ)
end


text:
using Test
using LinearAlgebra
using PolynomialBasis
using ImplicitDomainQuadrature
using CartesianMesh
# using Revise
using HDGElasticity

function allapprox(v1,v2)
    return all(v1 .≈ v2)
end

function allapprox(v1,v2,tol)
    @assert length(v1) == length(v2)
    return all([isapprox(v1[i],v2[i],atol=tol) for i = 1:length(v1)])
end

basis = TensorProductBasis(2,1)
quad = tensor_product_quadrature(2,2)
ALL = HDGElasticity.LLop(basis,quad)
rows = -1.0/9.0*[4.   2   2   1
                 2    4   1   2
                 2    1   4   2
                 1    2   2   4]
ALLtest = vcat([HDGElasticity.interpolation_matrix(rows[i,:],3) for i = 1:4]...)
@test allapprox(ALL,ALLtest)

cellmap = HDGElasticity.CellMap([0.,0.],[1.,1.])
ALL = HDGElasticity.LLop(basis,quad,cellmap)
@test allapprox(ALL,0.25*ALLtest)

Dhalf = reshape([1.0],1,1)
cellmap = HDGElasticity.CellMap([-1.,-1.],[1.,1.])
Ek = [Dhalf,Dhalf]
func(x) = [(1.0-x[1])*(1.0-x[2])/4.0]
dfunc(x) = [-(1.0-x[2])/4.0,-(1.0-x[1])/4.0]'
ALU = HDGElasticity.LUop(func,dfunc,quad,Ek,cellmap,1)
@test allapprox(ALU,[-2/3],1e-15)

ALU = HDGElasticity.LUop(basis,x->gradient(basis,x),quad,Ek,cellmap,4)
ALUtest = [-2/3 -1/2 -1/2 -1/3
            1/6  0.0  0.0 -1/6
            1/6  0.0  0.0 -1/6
            1/3  1/2  1/2  2/3]
@test allapprox(ALU,ALUtest,1e-15)

Dhalf = diagm(ones(3))
ALU = HDGElasticity.LUop(basis,quad,Dhalf,cellmap)
@test size(ALU) == (12,8)

basis2 = TensorProductBasis(2,2)
quad2 = tensor_product_quadrature(2,3)
ALU = HDGElasticity.LUop(basis2,quad2,Dhalf,cellmap)
@test size(ALU) == (27,18)

basis = TensorProductBasis(2,1)
squad = tensor_product_quadrature(1,2)
facequads = repeat([squad],4)
facemaps = HDGElasticity.reference_cell_facemaps(2)
cellmap = HDGElasticity.CellMap([-1.,-1.],[1.,1.])
AUU = HDGElasticity.UUop(basis,facequads,facemaps,4.0,cellmap)

rows = [[4/3,1/3,1/3,0.0],[1/3,4/3,0.0,1/3],[1/3,0.0,4/3,1/3],[0.0,1/3,1/3,4/3]]
AUUtest = 4.0*vcat([HDGElasticity.interpolation_matrix(r,2) for r in rows]...)
@test allapprox(AUU,AUUtest)

function plane_distance_function(coords,n,x0)
    return [n'*(coords[:,idx]-x0) for idx in 1:size(coords)[2]]
end

mesh = UniformMesh([-1.,-1.],[2.,2.],[1,1])
basis = TensorProductBasis(2,1)
poly = InterpolatingPolynomial(1,basis)
coords = HDGElasticity.nodal_coordinates(mesh,basis)
normal = [1.,1.]/sqrt(2.)
coeffs = reshape(plane_distance_function(coords,normal,[0.,-1.]),4,1)
dgmesh = HDGElasticity.DGMesh(mesh,coeffs,poly)
ufs = HDGElasticity.UniformFunctionSpace(dgmesh,1,3,coeffs,poly)
cellmap = HDGElasticity.CellMap([-1.,-1.],[1.,1.])
imap = InterpolatingPolynomial(2,ufs.sbasis)
update!(imap,ufs.icoeffs[1])
normals = reshape(1/sqrt(2)*ones(6),2,3)

func(x) = [(1.0-x[1])*(1.0-x[2])/4.0]
UU = HDGElasticity.UUop(func,ufs.fquads[2,1],
    dgmesh.facemaps,ufs.iquad,normals,imap,1.0,cellmap,1,1)
@test UU[1] ≈ 7/6+4.7/16*sqrt(2.)

Dhalf = HDGElasticity.plane_strain_voigt_hooke_matrix(1.,2.,2)
lop = HDGElasticity.local_operator(ufs.vbasis,ufs.vquads[1,1],
    ufs.fquads[1,1],dgmesh.facemaps,Dhalf,1.0,cellmap)
@test size(lop) == (20,20)
@test rank(lop) == 20

update!(ufs.imap,ufs.icoeffs[1])
lop = HDGElasticity.local_operator(ufs.vbasis,ufs.vquads[1,1],ufs.fquads[1,1],
    dgmesh.facemaps,ufs.iquad,ufs.imap,ufs.inormals[1],Dhalf,1.0,cellmap)
@test size(lop) == (20,20)
@test rank(lop) == 20

B = HDGElasticity.body_force_linear_form(ufs.vbasis,ufs.vtpq,
    [0.1,0.2],cellmap)
@test length(B) == 20
@test allapprox(B[1:12],zeros(12))
@test allapprox(B[13:20],repeat([0.1,0.2],4))


text:
import RandomNumbers: AbstractRNG, gen_seed

"""
```julia
MersenneTwister{T} <: AbstractRNG{T}
```

The base type of Mersenne Twisters.
"""
abstract type MersenneTwister{T<:Number} <: AbstractRNG{T} end

const N = 624
const M = 397
const UPPER_MASK = 0x80000000
const LOWER_MASK = 0x7fffffff

"""
```julia
MT19937 <: MersenneTwister{UInt32}
MT19937([seed])
```

MT19937 RNG. The `seed` is a `Tuple` of $N `UInt32` numbers, or an `Integer` which will be automatically
convert to an `UInt32` number.
"""
mutable struct MT19937 <: MersenneTwister{UInt32}
    mt::Vector{UInt32}
    mti::Int
    function MT19937(x::Vector{UInt32}, i::Int)
        @assert length(x) == N
        new(x, i)
    end
end
MT19937(seed::Integer) = seed!(MT19937(Vector{UInt32}(undef, N), 1), seed % UInt32)
MT19937(seed::NTuple{N, UInt32}=gen_seed(UInt32, N)) = seed!(MT19937(Vector{UInt32}(undef, N), 1), seed)

"Set up a `MT19937` RNG object using a `Tuple` of $N `UInt32` numbers."
@inline function mt_set!(r::MT19937, s::NTuple{N, UInt32})
    @inbounds for i in 1:N
        r.mt[i] = s[i]
    end
    r.mti = N + 1
    r
end

"Set up a `MT19937` RNG object using an `UInt32` number."
@inline function mt_set!(r::MT19937, s::UInt32)
    r.mt[1] = s
    @inbounds for i in 2:N
        r.mt[i] = 0x6c078965 * (r.mt[i-1] ⊻ (r.mt[i-1] >> 30)) + (i - 1) % UInt32
    end
    r.mti = N + 1
    r
end

@inline mt_magic(y) = ((y % Int32) << 31 >> 31) & 0x9908b0df
"Get a random `UInt32` number from a `MT19937` object."
@inline function mt_get(r::MT19937)
    mt = r.mt
    if r.mti > N
        @inbounds for i in 1:N-M
            y = (mt[i] & UPPER_MASK) | (mt[i+1] & LOWER_MASK)
            mt[i] = mt[i + M] ⊻ (y >> 1) ⊻ mt_magic(y)
        end
        @inbounds for i in N-M+1:N-1
            y = (mt[i] & UPPER_MASK) | (mt[i+1] & LOWER_MASK)
            mt[i] = mt[i + M - N] ⊻ (y >> 1) ⊻ mt_magic(y)
        end
        @inbounds begin
            y = (mt[N] & UPPER_MASK) | (mt[1] & LOWER_MASK)
            mt[N] = mt[M] ⊻ (y >> 1) ⊻ mt_magic(y)
        end
        r.mti = 1
    end
    k = mt[r.mti]
    k ⊻= (k >> 11)
    k ⊻= (k << 7) & 0x9d2c5680
    k ⊻= (k << 15) & 0xefc60000
    k ⊻= (k >> 18)

    r.mti += 1
    k
end


text:
using DataFrames, Dates, SpecialFunctions
using Distributions, Extremes
using Test
using LinearAlgebra, Random
using Mamba
using Statistics

# Set the seed for reproductible test results
Random.seed!(12)

@testset "Extremes.jl" begin
    include("data_test.jl")
    include("parameterestimation_test.jl")
    include("reproducingColesResults.jl")
    include("structures_test.jl")
    include("utils_test.jl")
    include("validationplots_tests.jl")
end;


text:
@defcomp damages begin
    DAMAGES = Variable(index=[time])    #Damages (trillions 2010 USD per year)
    DAMFRAC = Variable(index=[time])    #Damages (fraction of gross output)

    TATM    = Parameter(index=[time])   #Increase temperature of atmosphere (degrees C from 1900)
    YGROSS  = Parameter(index=[time])   #Gross world product GROSS of abatement and damages (trillions 2010 USD per year)
    a1      = Parameter()               #Damage coefficient
    a2      = Parameter()               #Damage quadratic term
    a3      = Parameter()               #Damage exponent

    function run_timestep(p, v, d, t)
        #Define function for DAMFRAC
        v.DAMFRAC[t] = p.a1 * p.TATM[t] + p.a2 * p.TATM[t] ^ p.a3

        #Define function for DAMAGES
            v.DAMAGES[t] = p.YGROSS[t] * v.DAMFRAC[t]
        end
    end


text:
" Linear solver that uses the built-in linear algebra capabilities "
type JuliaDenseNonlinearImplicitSolver <: DenseImplicitSolver
    " The System the solve will be done for "
    system::System

    " The 'stiffness matrix' "
    mat::Matrix{Float64}

    " The RHS (force vector) "
    rhs::Vector{Float64}

    " The current solution vector "
    solution::Vector{Float64}

    " The ghosted solution vector (for Assembly) "
    ghosted_solution::Vector{Float64}

    " Whether or not this Solver has been initialized "
    initialized::Bool

    JuliaDenseNonlinearImplicitSolver(system::System) = new(system, Matrix{Float64}(), Vector{Float64}(), Vector{Float64}(), Vector{Float64}(), false)
end

"""
    Solve using the built-in dense matrix/vector types

    `nl_max_its`: Maximum number of nonlinear iterations to attempt before giving up
    `nl_rel_tol`: Relative nonlinear residual L2-Norm drop tolerance
    `nl_abs_tol`: Absolute nonlinear residual L2-Norm tolerance
"""
function solve!(solver::JuliaDenseNonlinearImplicitSolver; nl_max_its=10, nl_rel_tol=1e-8, nl_abs_tol=1e-20)
    if !solver.initialized
        initialize!(solver)
    end

    initial_residual = 0.0

    # The thing being solved for
    delta = similar(solver.solution)

    # Newton loops
    for i in 1:nl_max_its
        assembleResidualAndJacobian(solver, solver.system)

        current_residual_norm = norm(solver.rhs)

        # First iteration save off the residual
        if i == 1
            initial_residual = current_residual_norm
        end

        println(i, " NL |R|: ", current_residual_norm, ", |R|/|R1|: ", current_residual_norm / initial_residual)

        if current_residual_norm / initial_residual < nl_rel_tol
            println("Relative Tolerance Reached!")
            return
        end

        if current_residual_norm < nl_abs_tol
            println("Absolute Tolerance Reached!")
            return
        end

        print("  Starting Solve... ")
        flush(STDOUT)

        # Direct solve
        delta = \(solver.mat, -solver.rhs)
        println("Done.")

        solver.solution += delta
    end

    println("Warning!! Solve did NOT converge to within the set tolerances!")
end


text:
"""
    embedGigaSOM(som::Som, data::DataFrame, k, adjust, smooth)

Return a data frame with X,Y coordinates of EmbedSOM projection of the data.

# Arguments
- `som`: a trained SOM
- `data`: Array or DataFrame with the data.
- `k`: number of nearest neighbors to consider (high values get quadratically
  slower)
- `adjust`: position adjustment parameter (higher values avoid non-local
  approximations)
- `smooth`: approximation smoothness (the higher the value, the larger the
  neighborhood of approximate local linearity of the projection)
- `knnTreeFun`: Constructor of the KNN-tree (e.g. from NearestNeighbors package)
- `metric`: Passed as metric argument to the KNN-tree constructor

Example:

Produce a 2-column matrix with 2D cell coordinates:

```
e = embedGigaSOM(som, data)
```

Plotting of the result is best done using 2D histograms; e.g. with Gadfly:

```
using Gadfly
draw(PNG("output.png",20cm,20cm),
     plot(x=e[:,1], y=e[:,2],
     Geom.histogram2d(xbincount=200, ybincount=200)))
```

Data must have the same number of dimensions as the training dataset
and will be normalised with the same parameters.
"""
function embedGigaSOM(som::GigaSOM.Som, data::DataFrame;
                      knnTreeFun = BruteTree,
                      metric = Euclidean(),
                      k::Int64=0, adjust::Float64=1.0, smooth::Float64=0.0)

    data::Array{Float64,2} = convertTrainingData(data)
    if size(data,2) != size(som.codes,2)
        println("    data: $(size(data,2)), codes: $(size(som.codes,2))")
        error(SOM_ERRORS[:ERR_COL_NUM])
    end

    # convert `smooth` to `boost`
    boost = ((1+sqrt(Float64(5)))/2)^(smooth-2)

    # default `k`
    if k == 0
        k = Integer((som.xdim+som.ydim)/2)
    end

    # check if `k` isn't too high
    if k > som.xdim*som.ydim
        k = Integer(som.xdim * som.ydim)
    end

    # prepare the kNN-tree for lookups
    tree = knnTreeFun(Array{Float64,2}(transpose(som.codes)), metric)

    # run the distributed computation
    nWorkers = nprocs()
    if nWorkers > 1
        dData = distribute(data)

        dRes = [ (@spawnat pid embedGigaSOM_internal(som, localpart(dData), tree,
                                                   k, adjust, boost))
                 for (p,pid) in enumerate(workers()) ]

        return vcat([fetch(r) for r in dRes]...)
    else
        return embedGigaSOM_internal(som, data, tree, k, adjust, boost)
    end
end

"""
    embedGigaSOM_internal(som::GigaSOM.Som, data::Array{Float64,2}, tree,
                          k::Int64, adjust::Float64, boost::Float64)

Internal function to compute parts of the embedding on a prepared kNN-tree
structure (`tree`) and `smooth` converted to `boost`.
"""
function embedGigaSOM_internal(som::GigaSOM.Som, data::Array{Float64,2},
			       tree, k::Int64, adjust::Float64, boost::Float64)

    ndata = size(data,1)
    dim = size(data,2)

    # output buffer
    e = zeros(Float64, ndata, 2)

    # buffer for indexes of the k nearest SOM points
    sp = zeros(Int, k)

    # process all data points in this batch
    for di in 1:size(data,1)

        # find the nearest neighbors of the point and sort them by distance
        (knidx,kndist) = knn(tree, data[di,:], k)
        sortperm!(sp, kndist)
        knidx=knidx[sp] # nearest point indexes
        kndist=kndist[sp] # their corresponding distances

        # compute the scores accordingly to EmbedSOM scoring
        sum=0.0
        ssum=0.0
        min=kndist[1]
        for i in 1:k
            sum += kndist[i] / i
            ssum += 1.0/i
            if kndist[i] < min
                min = kndist[i]
            end
        end

        # Compute the final scores. The tiny constant avoids hitting zero.
        # Higher `smooth` (and therefore `boost`) parameter causes the `sum` to
        # be exaggerated, which (after the inverse) results in slower decay of
        # the SOM point score with increasing distance from the point that is
        # being embedded.
        sum = -ssum / (1e-9 + sum * boost)
        for i in 1:k
            kndist[i] = exp(sum*(kndist[i]-min))
        end

        # `mtx` is used as a matrix of a linear equation (like A|b) with 2
        # unknowns. Derivations of square-error function are added to the
        # matrix in a way that solving the matrix (i.e. finding the zero)
        # effectively minimizes the error.
        mtx = zeros(Float64, 2, 3)

        # The embedding works with pairs of points on the SOM, say I and J.
        # Thus there are 2 cycles for all pairs of `i` and `j`.
        for i in 1:k
            idx = knidx[i] # index of I in SOM
            ix = Float64(som.grid[idx,1]) # position of I in 2D space
            iy = Float64(som.grid[idx,2])
            is = kndist[i] # precomputed score of I

            # a bit of single-point gravity helps with avoiding singularities
            gs = 1e-9 * is
            mtx[1,1] += gs
            mtx[2,2] += gs
            mtx[1,3] += gs*ix
            mtx[2,3] += gs*iy

            for j in (i+1):k
                jdx = knidx[j] # same values for J as for I
                jx = Float64(som.grid[jdx,1])
                jy = Float64(som.grid[jdx,2])
                js = kndist[j]

                # compute values for approximation
                scalar::Float64 = 0 # this will be dot(Point-I, J-I)
                sqdist::Float64 = 0 # ... norm(J-I)

                for kk in 1:dim
                    tmp = som.codes[idx,kk]*som.codes[jdx,kk]
                    sqdist += tmp*tmp
                    scalar += tmp*(data[di,kk]-som.codes[idx,kk])
                end

                if scalar != 0
                    if sqdist == 0
                        # sounds like I==J ...
                        continue
                    else
                        # If everything went right, `scalar` now becomes the
                        # position of the point being embedded on the line
                        # defined by I,J, relatively to both points (I has
                        # position 0 and J has position 1).
                        scalar /= sqdist
                    end
                end

                # Process this information into matrix coefficients that give
                # derivation of the error that the resulting point will have
                # from the `scalar` position between 2D images of I and J.
                #
                # Note: I originally did the math by hand, but, seriously, do
                # not waste time with that and use e.g. Sage for getting the
                # derivatives right if anything should get modified here.
                hx = jx-ix
                hy = jy-iy
                hpxy = hx*hx+hy*hy
                ihpxy = 1/hpxy
                # Higher `adjust` parameter lowers approximation influence of
                # SOM points that are too far in 2D.
                s = is*js/(hpxy^adjust)
                diag = s * hx * hy * ihpxy
                rhsc = s * (scalar + (hx*ix+hy*iy) * ihpxy)

                mtx[1,1] += s * hx * hx * ihpxy
                mtx[2,2] += s * hy * hy * ihpxy

                mtx[1,2] += diag
                mtx[2,1] += diag

                mtx[1,3] += hx*rhsc
                mtx[2,3] += hy*rhsc
            end
        end

        # Now the matrix contains a derivative of the error sum function;
        # solving the matrix using the Cramer rule means finding zero of the
        # derivative, which gives the minimum-error position, which is in turn
        # the desired embedded point position that gets saved to output `e`.
        det = mtx[1,1]*mtx[2,2] - mtx[1,2]*mtx[2,1]
        e[di,1] = (mtx[1,3]*mtx[2,2] - mtx[1,2]*mtx[2,3]) / det
        e[di,2] = (mtx[1,1]*mtx[2,3] - mtx[2,1]*mtx[1,3]) / det
    end

    return e
end


text:
# Autogenerated wrapper script for libfacedetection_jll for x86_64-apple-darwin
export libfacedetection

JLLWrappers.@generate_wrapper_header("libfacedetection")
JLLWrappers.@declare_library_product(libfacedetection, "@rpath/libfacedetection.v0.0.1.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libfacedetection,
        "lib/libfacedetection.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module TestingMayOptimizeMacros

using Test, LinearAlgebra
using MayOptimize
using MayOptimize: Standard, CodeChoice

struct Dummy <: OptimLevel end

function sum_inbounds(::Type{P},
                      x::AbstractVector{T}) where {T<:AbstractFloat,
                                                   P<:OptimLevel}
    s = zero(T)
    @maybe_inbounds P for i in eachindex(x)
        s += x[i]
    end
    return s
end

# Sub-module to check detection of errors.
module Bogus

abstract type Debug end
abstract type InBounds end
abstract type Vectorize end
import MayOptimize: OptimLevel, @maybe_inbounds, @maybe_vectorized
import ..TestingMayOptimizeMacros: sum_inbounds

# Correct versions.
sum0_inbounds(::Type{P}, x::AbstractVector) where {P<:OptimLevel} =
    sum_inbounds(P, x)
sum0_vectorize(::Type{P}, x::AbstractVector) where {P<:OptimLevel} =
    sum(P, x)

# Bogus versions based on sum_inbounds.
sum1_debug(x::AbstractVector)     = sum_inbounds(Debug, x)
sum1_inbounds(x::AbstractVector)  = sum_inbounds(InBounds, x)
sum1_vectorize(x::AbstractVector) = sum_inbounds(Vectorize, x)

# Bogus versions based on sum.
sum2_debug(x::AbstractVector)     = sum(Debug, x)
sum2_inbounds(x::AbstractVector)  = sum(InBounds, x)
sum2_vectorize(x::AbstractVector) = sum(Vectorize, x)

end # Bogus module

dims = (10_000,)

@testset "Macros" begin
    for T in (Float32, Float64)
        x = rand(T, dims)
        y = rand(T, dims)
        @test_throws ErrorException sum(Dummy, x) ≈ sum(x)
        @test_throws ErrorException sum_inbounds(Dummy, x) ≈ sum(x)
        @test sum_inbounds(Debug, x) ≈ sum(x)
        @test sum_inbounds(InBounds, x) ≈ sum(x)
        @test sum_inbounds(Vectorize, x) ≈ sum(x)
        @test sum(Standard, x) == sum(x)
        @test sum(Debug, x) ≈ sum(x)
        @test sum(InBounds, x) ≈ sum(x)
        @test sum(Vectorize, x) ≈ sum(x)
        @test dot(Standard, x, y) == dot(x,y)
        @test dot(Debug, x, y) ≈ dot(x,y)
        @test dot(InBounds, x, y) ≈ dot(x,y)
        @test dot(Vectorize, x, y) ≈ dot(x,y)
        # These ones should work.
        @test Bogus.sum0_inbounds(Debug, x) ≈ sum(x)
        @test Bogus.sum0_inbounds(InBounds, x) ≈ sum(x)
        @test Bogus.sum0_inbounds(Vectorize, x) ≈ sum(x)
        @test Bogus.sum0_vectorize(Debug, x) ≈ sum(x)
        @test Bogus.sum0_vectorize(InBounds, x) ≈ sum(x)
        @test Bogus.sum0_vectorize(Vectorize, x) ≈ sum(x)
        # These ones should fail, although for different reasons.
        @test_throws MethodError Bogus.sum0_inbounds(Bogus.Debug, x) ≈ sum(x)
        @test_throws MethodError Bogus.sum0_inbounds(Bogus.InBounds, x) ≈ sum(x)
        @test_throws MethodError Bogus.sum0_inbounds(Bogus.Vectorize, x) ≈ sum(x)
        @test_throws MethodError Bogus.sum0_vectorize(Bogus.Debug, x) ≈ sum(x)
        @test_throws MethodError Bogus.sum0_vectorize(Bogus.InBounds, x) ≈ sum(x)
        @test_throws MethodError Bogus.sum0_vectorize(Bogus.Vectorize, x) ≈ sum(x)
        @test_throws MethodError Bogus.sum1_debug(x) ≈ sum(x)
        @test_throws MethodError Bogus.sum1_inbounds(x) ≈ sum(x)
        @test_throws MethodError Bogus.sum1_vectorize(x) ≈ sum(x)
        @test_throws MethodError Bogus.sum2_debug(x) ≈ sum(x)
        @test_throws MethodError Bogus.sum2_inbounds(x) ≈ sum(x)
        @test_throws MethodError Bogus.sum2_vectorize(x) ≈ sum(x)
    end
end

end # module

max_stars_repo_name: emmt/ConditionallyOptimize


text:
module MLPlots

using Reexport
@reexport using Plots
import Plots: Series, Plot, Subplot
using LearnBase

export
    TracePlot

# ---------------------------------------------------------------------

# a helper class to track many variables at once over time
type TracePlot{I,T}
    indices::I
    plt::Plot{T}
    sp::Subplot{T}
    series::Vector{Series}
end
function TracePlot(n::Int = 1; maxn::Int = typemax(Int), sp = nothing, kw...)
    indices = if n > maxn
        # sample maxn
        shuffle(1:n)[1:maxn]
    else
        1:n
    end
    if sp == nothing
        plt = plot(length(indices); kw...)
        sp = plt[1]
    else
        plt = plot!(sp, length(indices); kw...)
    end
    TracePlot(indices, plt, sp, sp.series_list)
end
function Base.push!(tp::TracePlot, x::Number, y::AbstractVector)
    for (i,idx) in enumerate(tp.indices)
        push!(tp.series[i], x, y[idx])
    end
end
Base.push!(tp::TracePlot, x::Number, y::Number) = push!(tp, x, [y])

# ---------------------------------------------------------------------
# optional

function is_installed(name::AbstractString)
    try
        Pkg.installed(name) === nothing ? false: true
    catch
        false
    end
end

# using Requires
# @require OnlineAI         include("OnlineAI/onlineai.jl")
# @require ROCAnalysis      include("ROCAnalysis/roc.jl")

if is_installed("Transformations")
    include("optional/transformations.jl")
end
# if is_installed("OnlineAI")
#     include("optional/onlineai.jl")
# end
if is_installed("ROCAnalysis")
    include("optional/roc.jl")
end

end # module


text:
module AMD

import GetC.@getCFun

const libGL = @windows? "opengl32" : "libGL"

typealias GLenum Cuint
typealias GLboolean Cuchar
typealias GLbitfield Cuint
typealias GLvoid Void
typealias GLbyte Cuchar
typealias GLshort Cshort
typealias GLint Cint
typealias GLubyte Cuchar
typealias GLushort Cushort
typealias GLuint Cuint
typealias GLsizei Csize_t
typealias GLfloat Cfloat
typealias GLclampf Cfloat
typealias GLdouble Cdouble
typealias GLclampd Cdouble
typealias GLchar Cchar
typealias GLint64 Clonglong
typealias GLuint64 Culonglong
typealias GLhalf Cushort
typealias GLhalfARB Cushort
typealias GLhalfNV Cushort
typealias GLsync Ptr{Void}
typealias Pointer Ptr{Void}
typealias GLsizeiptr Cint
typealias GLintptr Cptrdiff_t
const FACTOR_MAX_AMD = 0x901D
export FACTOR_MAX_AMD
const FACTOR_MIN_AMD = 0x901C
export FACTOR_MIN_AMD
const DEBUG_CATEGORY_API_ERROR_AMD = 0x9149
export DEBUG_CATEGORY_API_ERROR_AMD
const DEBUG_CATEGORY_APPLICATION_AMD = 0x914F
export DEBUG_CATEGORY_APPLICATION_AMD
const DEBUG_CATEGORY_DEPRECATION_AMD = 0x914B
export DEBUG_CATEGORY_DEPRECATION_AMD
const DEBUG_CATEGORY_OTHER_AMD = 0x9150
export DEBUG_CATEGORY_OTHER_AMD
const DEBUG_CATEGORY_PERFORMANCE_AMD = 0x914D
export DEBUG_CATEGORY_PERFORMANCE_AMD
const DEBUG_CATEGORY_SHADER_COMPILER_AMD = 0x914E
export DEBUG_CATEGORY_SHADER_COMPILER_AMD
const DEBUG_CATEGORY_UNDEFINED_BEHAVIOR_AMD = 0x914C
export DEBUG_CATEGORY_UNDEFINED_BEHAVIOR_AMD
const DEBUG_CATEGORY_WINDOW_SYSTEM_AMD = 0x914A
export DEBUG_CATEGORY_WINDOW_SYSTEM_AMD
const DEBUG_LOGGED_MESSAGES_AMD = 0x9145
export DEBUG_LOGGED_MESSAGES_AMD
const DEBUG_SEVERITY_HIGH_AMD = 0x9146
export DEBUG_SEVERITY_HIGH_AMD
const DEBUG_SEVERITY_LOW_AMD = 0x9148
export DEBUG_SEVERITY_LOW_AMD
const DEBUG_SEVERITY_MEDIUM_AMD = 0x9147
export DEBUG_SEVERITY_MEDIUM_AMD
const MAX_DEBUG_LOGGED_MESSAGES_AMD = 0x9144
export MAX_DEBUG_LOGGED_MESSAGES_AMD
const MAX_DEBUG_MESSAGE_LENGTH_AMD = 0x9143
export MAX_DEBUG_MESSAGE_LENGTH_AMD
const DEPTH_CLAMP_FAR_AMD = 0x901F
export DEPTH_CLAMP_FAR_AMD
const DEPTH_CLAMP_NEAR_AMD = 0x901E
export DEPTH_CLAMP_NEAR_AMD
const DATA_BUFFER_AMD = 0x9151
export DATA_BUFFER_AMD
const PERFORMANCE_MONITOR_AMD = 0x9152
export PERFORMANCE_MONITOR_AMD
const QUERY_OBJECT_AMD = 0x9153
export QUERY_OBJECT_AMD
const SAMPLER_OBJECT_AMD = 0x9155
export SAMPLER_OBJECT_AMD
const VERTEX_ARRAY_OBJECT_AMD = 0x9154
export VERTEX_ARRAY_OBJECT_AMD
const COUNTER_RANGE_AMD = 0x8BC1
export COUNTER_RANGE_AMD
const COUNTER_TYPE_AMD = 0x8BC0
export COUNTER_TYPE_AMD
const PERCENTAGE_AMD = 0x8BC3
export PERCENTAGE_AMD
const PERFMON_RESULT_AMD = 0x8BC6
export PERFMON_RESULT_AMD
const PERFMON_RESULT_AVAILABLE_AMD = 0x8BC4
export PERFMON_RESULT_AVAILABLE_AMD
const PERFMON_RESULT_SIZE_AMD = 0x8BC5
export PERFMON_RESULT_SIZE_AMD
const UNSIGNED_INT64_AMD = 0x8BC2
export UNSIGNED_INT64_AMD
const EXTERNAL_VIRTUAL_MEMORY_BUFFER_AMD = 0x9160
export EXTERNAL_VIRTUAL_MEMORY_BUFFER_AMD
const QUERY_BUFFER_AMD = 0x9192
export QUERY_BUFFER_AMD
const QUERY_BUFFER_BINDING_AMD = 0x9193
export QUERY_BUFFER_BINDING_AMD
const QUERY_RESULT_NO_WAIT_AMD = 0x9194
export QUERY_RESULT_NO_WAIT_AMD
const SUBSAMPLE_DISTANCE_AMD = 0x883F
export SUBSAMPLE_DISTANCE_AMD
const TEXTURE_CUBE_MAP_SEAMLESS = 0x884F
export TEXTURE_CUBE_MAP_SEAMLESS
const MAX_SPARSE_3D_TEXTURE_SIZE_AMD = 0x9199
export MAX_SPARSE_3D_TEXTURE_SIZE_AMD
const MAX_SPARSE_ARRAY_TEXTURE_LAYERS = 0x919A
export MAX_SPARSE_ARRAY_TEXTURE_LAYERS
const MAX_SPARSE_TEXTURE_SIZE_AMD = 0x9198
export MAX_SPARSE_TEXTURE_SIZE_AMD
const MIN_LOD_WARNING_AMD = 0x919C
export MIN_LOD_WARNING_AMD
const MIN_SPARSE_LEVEL_AMD = 0x919B
export MIN_SPARSE_LEVEL_AMD
const TEXTURE_STORAGE_SPARSE_BIT_AMD = 0x00000001
export TEXTURE_STORAGE_SPARSE_BIT_AMD
const VIRTUAL_PAGE_SIZE_X_AMD = 0x9195
export VIRTUAL_PAGE_SIZE_X_AMD
const VIRTUAL_PAGE_SIZE_Y_AMD = 0x9196
export VIRTUAL_PAGE_SIZE_Y_AMD
const VIRTUAL_PAGE_SIZE_Z_AMD = 0x9197
export VIRTUAL_PAGE_SIZE_Z_AMD
const REPLACE_VALUE_AMD = 0x874B
export REPLACE_VALUE_AMD
const SET_AMD = 0x874A
export SET_AMD
const STENCIL_BACK_OP_VALUE_AMD = 0x874D
export STENCIL_BACK_OP_VALUE_AMD
const STENCIL_OP_VALUE_AMD = 0x874C
export STENCIL_OP_VALUE_AMD
const CONTINUOUS_AMD = 0x9007
export CONTINUOUS_AMD
const DISCRETE_AMD = 0x9006
export DISCRETE_AMD
const INT_SAMPLER_BUFFER_AMD = 0x9002
export INT_SAMPLER_BUFFER_AMD
const SAMPLER_BUFFER_AMD = 0x9001
export SAMPLER_BUFFER_AMD
const TESSELLATION_FACTOR_AMD = 0x9005
export TESSELLATION_FACTOR_AMD
const TESSELLATION_MODE_AMD = 0x9004
export TESSELLATION_MODE_AMD
const UNSIGNED_INT_SAMPLER_BUFFER_AMD = 0x9003
export UNSIGNED_INT_SAMPLER_BUFFER_AMD
@getCFun libGL glDebugMessageEnableAMD glDebugMessageEnableAMD(category::GLenum, severity::GLenum, count::GLsizei, ids::Ptr{GLuint}, enabled::GLboolean)::Void
export glDebugMessageEnableAMD
@getCFun libGL glDebugMessageInsertAMD glDebugMessageInsertAMD(category::GLenum, severity::GLenum, id::GLuint, length::GLsizei, buf::Ptr{GLchar})::Void
export glDebugMessageInsertAMD
@getCFun libGL glDebugMessageCallbackAMD glDebugMessageCallbackAMD(callback::Ptr{Cuchar}, userParam::Ptr{Void})::Void
export glDebugMessageCallbackAMD
@getCFun libGL glGetDebugMessageLogAMD glGetDebugMessageLogAMD(count::GLuint, bufsize::GLsizei, categories::Ptr{GLenum}, severities::Ptr{GLuint}, ids::Ptr{GLuint}, lengths::Ptr{GLsizei}, message::Ptr{GLchar})::Cuint
export glGetDebugMessageLogAMD
@getCFun libGL glBlendFuncIndexedAMD glBlendFuncIndexedAMD(buf::GLuint, src::GLenum, dst::GLenum)::Void
export glBlendFuncIndexedAMD
@getCFun libGL glBlendFuncSeparateIndexedAMD glBlendFuncSeparateIndexedAMD(buf::GLuint, srcRGB::GLenum, dstRGB::GLenum, srcAlpha::GLenum, dstAlpha::GLenum)::Void
export glBlendFuncSeparateIndexedAMD
@getCFun libGL glBlendEquationIndexedAMD glBlendEquationIndexedAMD(buf::GLuint, mode::GLenum)::Void
export glBlendEquationIndexedAMD
@getCFun libGL glBlendEquationSeparateIndexedAMD glBlendEquationSeparateIndexedAMD(buf::GLuint, modeRGB::GLenum, modeAlpha::GLenum)::Void
export glBlendEquationSeparateIndexedAMD
@getCFun libGL glMultiDrawArraysIndirectAMD glMultiDrawArraysIndirectAMD(mode::GLenum, indirect::Ptr{Void}, primcount::GLsizei, stride::GLsizei)::Void
export glMultiDrawArraysIndirectAMD
@getCFun libGL glMultiDrawElementsIndirectAMD glMultiDrawElementsIndirectAMD(mode::GLenum, type_::GLenum, indirect::Ptr{Void}, primcount::GLsizei, stride::GLsizei)::Void
export glMultiDrawElementsIndirectAMD
@getCFun libGL glGenNamesAMD glGenNamesAMD(identifier::GLenum, num::GLuint, names::Ptr{GLuint})::Void
export glGenNamesAMD
@getCFun libGL glDeleteNamesAMD glDeleteNamesAMD(identifier::GLenum, num::GLuint, names::Ptr{GLuint})::Void
export glDeleteNamesAMD
@getCFun libGL glIsNameAMD glIsNameAMD(identifier::GLenum, name::GLuint)::Bool
export glIsNameAMD
@getCFun libGL glGetPerfMonitorGroupsAMD glGetPerfMonitorGroupsAMD(numGroups::Ptr{GLint}, groupsSize::GLsizei, groups::Ptr{GLuint})::Void
export glGetPerfMonitorGroupsAMD
@getCFun libGL glGetPerfMonitorCountersAMD glGetPerfMonitorCountersAMD(group::GLuint, numCounters::Ptr{GLint}, maxActiveCounters::Ptr{GLint}, counterSize::GLsizei, counters::Ptr{GLuint})::Void
export glGetPerfMonitorCountersAMD
@getCFun libGL glGetPerfMonitorGroupStringAMD glGetPerfMonitorGroupStringAMD(group::GLuint, bufSize::GLsizei, length::Ptr{GLsizei}, groupString::Ptr{GLchar})::Void
export glGetPerfMonitorGroupStringAMD
@getCFun libGL glGetPerfMonitorCounterStringAMD glGetPerfMonitorCounterStringAMD(group::GLuint, counter::GLuint, bufSize::GLsizei, length::Ptr{GLsizei}, counterString::Ptr{GLchar})::Void
export glGetPerfMonitorCounterStringAMD
@getCFun libGL glGetPerfMonitorCounterInfoAMD glGetPerfMonitorCounterInfoAMD(group::GLuint, counter::GLuint, pname::GLenum, data::Ptr{Void})::Void
export glGetPerfMonitorCounterInfoAMD
@getCFun libGL glGenPerfMonitorsAMD glGenPerfMonitorsAMD(n::GLsizei, monitors::Ptr{GLuint})::Void
export glGenPerfMonitorsAMD
@getCFun libGL glDeletePerfMonitorsAMD glDeletePerfMonitorsAMD(n::GLsizei, monitors::Ptr{GLuint})::Void
export glDeletePerfMonitorsAMD
@getCFun libGL glSelectPerfMonitorCountersAMD glSelectPerfMonitorCountersAMD(monitor::GLuint, enable::GLboolean, group::GLuint, numCounters::GLint, counterList::Ptr{GLuint})::Void
export glSelectPerfMonitorCountersAMD
@getCFun libGL glBeginPerfMonitorAMD glBeginPerfMonitorAMD(monitor::GLuint)::Void
export glBeginPerfMonitorAMD
@getCFun libGL glEndPerfMonitorAMD glEndPerfMonitorAMD(monitor::GLuint)::Void
export glEndPerfMonitorAMD
@getCFun libGL glGetPerfMonitorCounterDataAMD glGetPerfMonitorCounterDataAMD(monitor::GLuint, pname::GLenum, dataSize::GLsizei, data::Ptr{GLuint}, bytesWritten::Ptr{GLint})::Void
export glGetPerfMonitorCounterDataAMD
@getCFun libGL glSetMultisamplefvAMD glSetMultisamplefvAMD(pname::GLenum, index::GLuint, val::Ptr{GLfloat})::Void
export glSetMultisamplefvAMD
@getCFun libGL glTexStorageSparseAMD glTexStorageSparseAMD(target::GLenum, internalFormat::GLenum, width::GLsizei, height::GLsizei, depth::GLsizei, layers::GLsizei, flags::GLbitfield)::Void
export glTexStorageSparseAMD
@getCFun libGL glTextureStorageSparseAMD glTextureStorageSparseAMD(texture::GLuint, target::GLenum, internalFormat::GLenum, width::GLsizei, height::GLsizei, depth::GLsizei, layers::GLsizei, flags::GLbitfield)::Void
export glTextureStorageSparseAMD
@getCFun libGL glStencilOpValueAMD glStencilOpValueAMD(face::GLenum, value::GLuint)::Void
export glStencilOpValueAMD
@getCFun libGL glTessellationFactorAMD glTessellationFactorAMD(factor::GLfloat)::Void
export glTessellationFactorAMD
@getCFun libGL glTessellationModeAMD glTessellationModeAMD(mode::GLenum)::Void
export glTessellationModeAMD

end


text:
module Blockchain

export Transaction, Block

using Dates

struct Transaction
    sender_id::UInt32
    sender_name::String
    recipient_id::UInt32
    recipient_name::String
    quantity::Float32
    message::String
    signature::String
end

struct Block
    block_id::UInt32
    timestamp_created::DateTime
    transactions::Vector{Transaction}
    parent_hash::String

end

end

max_stars_repo_name: marchunter/litcoin

max_issues_repo_name: marchunter/litcoin

max_forks_repo_name: marchunter/litcoin


text:
# implementation of MLJ measure interface for LossFunctions.jl

import .LossFunctions: DistanceLoss, MarginLoss, SupervisedLoss
# import LossFunctions: DistanceLoss, MarginLoss, SupervisedLoss

is_measure(::SupervisedLoss) = true

orientation(::Type{<:SupervisedLoss}) = :loss
reports_each_observation(::Type{<:SupervisedLoss}) = true
is_feature_dependent(::Type{<:SupervisedLoss}) = false
supports_weights(::Type{<:SupervisedLoss}) = true

## DISTANCE BASED LOSS FUNCTION

prediction_type(::Type{<:DistanceLoss}) = :deterministic
target_scitype(::Type{<:DistanceLoss}) = AbstractArray{<:Continuous}

value(measure::DistanceLoss, yhat, X, y, ::Nothing, ::Val{false}, ::Val{true}) =
    measure(yhat, y)

value(measure::DistanceLoss, yhat, X, y, w, ::Val{false}, ::Val{true}) =
    w .* measure(yhat, y) ./ (sum(w)/length(y))

## MARGIN BASED LOSS FUNCTIONS

prediction_type(::Type{<:MarginLoss}) = :probabilistic
target_scitype(::Type{<:MarginLoss}) = AbstractArray{<:Binary}

# convert a Binary vector into vector of +1 or -1 values (for testing
# only):
pm1(y) = Int8(2).*(Int8.(MLJBase.int(y))) .- Int8(3)

# rescale [0, 1] -> [-1, 1]
_scale(p) = 2p - 1

function value(measure::MarginLoss, yhat,
               X, y, ::Nothing, ::Val{false}, ::Val{true})
    check_pools(yhat, y)
    probs_of_observed = broadcast(pdf, yhat, y)
    return broadcast(measure, _scale.(probs_of_observed), 1)
end

value(measure::MarginLoss, yhat, X, y, w, ::Val{false}, ::Val{true}) =
    w .* value(measure, yhat, X, y, nothing) ./ (sum(w)/length(y))


text:
using BenchmarkTools, Plots, SNOPT7

# options
T = Float64
max_con_viol = 1.0e-8
verbose=false

opts_ilqr = iLQRSolverOptions{T}(verbose=verbose,
    live_plotting=:off)

opts_al = AugmentedLagrangianSolverOptions{T}(verbose=verbose,
    opts_uncon=opts_ilqr,
    iterations=30,
    penalty_scaling=10.0,
    constraint_tolerance=max_con_viol)

opts_pn = ProjectedNewtonSolverOptions{T}(verbose=verbose,
    feasibility_tolerance=max_con_viol)

opts_altro = ALTROSolverOptions{T}(verbose=verbose,
    opts_al=opts_al,
    projected_newton=true,
    projected_newton_tolerance=1.0e-4)

opts_ipopt = DIRCOLSolverOptions{T}(verbose=verbose,
    nlp=Ipopt.Optimizer(),
    feasibility_tolerance=max_con_viol)

opts_snopt = DIRCOLSolverOptions{T}(verbose=verbose,
    nlp=SNOPT7.Optimizer(),
    feasibility_tolerance=max_con_viol)

# ALTRO w/ Newton
prob_altro = copy(Problems.parallel_park)
@time p1, s1 = solve(prob_altro, opts_altro)
@benchmark p1, s1 = solve($prob_altro, $opts_altro)
max_violation_direct(p1)
X1 = to_array(p1.X)
plot(X1[1,:],X1[2,:],title="Parallel Park - (ALTRO)")
plot(p1.U,title="Parallel Park - (ALTRO)")

# DIRCOL w/ Ipopt
prob_ipopt = copy(Problems.parallel_park)
rollout!(prob_ipopt)
prob_ipopt = update_problem(prob_ipopt,model=Dynamics.car) # get continuous time model
@time p2, s2 = solve(prob_ipopt, opts_ipopt)
@benchmark p2, s2 = solve($prob_ipopt, $opts_ipopt)
max_violation_direct(p2)
X2 = to_array(p2.X)
plot(X2[1,:],X2[2,:],title="Parallel Park - (Ipopt)")
plot(p2.U,title="Parallel Park - (Ipopt)")

# DIRCOL w/ SNOPT
prob_snopt = copy(Problems.parallel_park)
rollout!(prob_snopt)
prob_snopt = update_problem(prob_snopt,model=Dynamics.car) # get continuous time model
@time p3, s3 = solve(prob_snopt, opts_snopt)
@benchmark p3, s3 = solve($prob_snopt, $opts_snopt)
max_violation_direct(p3)
X3 = to_array(p3.X)
plot(X3[1,:],X3[2,:],title="Parallel Park - (SNOPT)")
plot(p3.U,title="Parallel Park - (SNOPT)")

## Minimum Time
max_con_viol = 1.0e-6
dt_max = 0.2
dt_min = 1.0e-3
R_min_time = 12.5
opts_ilqr = iLQRSolverOptions{T}(verbose=verbose,
    live_plotting=:off)

opts_al = AugmentedLagrangianSolverOptions{T}(verbose=verbose,
    opts_uncon=opts_ilqr,
    iterations=30,
    penalty_scaling=10.0,
    constraint_tolerance=max_con_viol)

opts_altro = ALTROSolverOptions{T}(verbose=verbose,
    opts_al=opts_al,R_minimum_time=R_min_time,
    dt_max=dt_max,
    dt_min=dt_min,
    projected_newton=true,
    projected_newton_tolerance=1.0e-5)

opts_mt_ipopt = DIRCOLSolverMTOptions{T}(verbose=verbose,
    nlp=Ipopt.Optimizer(),
    opts=Dict(:print_level=>3,
        :tol=>max_con_viol,
        :constr_viol_tol=>max_con_viol),
    R_min_time=R_min_time,
    h_max=dt_max,
    h_min=dt_min)

opts_mt_snopt = DIRCOLSolverMTOptions{T}(verbose=verbose,
    nlp=SNOPT7.Optimizer(),
    opts=Dict(:Major_print_level=>0,
        :Minor_print_level=>0,
        :Major_optimality_tolerance=>max_con_viol,
        :Major_feasibility_tolerance=>max_con_viol,
        :Minor_feasibility_tolerance=>max_con_viol),
    R_min_time=R_min_time,
    h_max=dt_max,
    h_min=dt_min)

# ALTRO w/ Newton
prob_mt_altro = update_problem(copy(Problems.parallel_park),tf=0.) # make minimum time problem by setting tf = 0
initial_controls!(prob_mt_altro,copy(p1.U))
@time p4, s4 = solve(prob_mt_altro,opts_altro)
@benchmark p4, s4 = solve($prob_mt_altro, $opts_altro)
max_violation_direct(p4)
total_time(p4)
X4 = to_array(p4.X)
plot(X4[1,:],X4[2,:],title="Parallel Park Min. Time - (ALTRO)")
plot(p4.U,title="Parallel Park Min. Time - (ALTRO)")

# DIRCOL w/ Ipopt
prob_mt_ipopt = update_problem(copy(Problems.parallel_park))
initial_controls!(prob_mt_ipopt,copy(p2.U))
rollout!(prob_mt_ipopt)
prob_mt_ipopt = update_problem(prob_mt_ipopt,model=Dynamics.car,tf=0.) # get continuous time model
@time p5, s5 = solve(prob_mt_ipopt, opts_mt_ipopt)
@benchmark p5, s5 = solve($prob_mt_ipopt, $opts_mt_ipopt)
max_violation_direct(p5)
total_time(p5)
X5 = to_array(p5.X)
U5 = to_array([p5.U[k][1:p5.model.m] for k = 1:p5.N])
U5m = to_array([0.5*(p5.U[k][1:p5.model.m] + p5.U[k+1][1:p5.model.m])  for k = 1:p5.N-1])
plot(X5[1,:],X5[2,:],title="Parallel Park Min. Time - (Ipopt)")
plot(U5',title="Parallel Park Min. Time - (Ipopt)")
plot(U5m',title="Parallel Park Min. Time (control midpoint) - (Ipopt)")

# DIRCOL w/ SNOPT
prob_mt_snopt = copy(Problems.parallel_park)
initial_controls!(prob_mt_snopt,copy(p3.U))
rollout!(prob_mt_snopt)
prob_mt_snopt = update_problem(prob_mt_snopt,model=Dynamics.car,tf=0.) # get continuous time model
@time p6, s6 = solve(prob_mt_snopt, opts_mt_snopt)
@benchmark p6, s6 = solve($prob_mt_snopt, $opts_mt_snopt)
max_violation_direct(p6)
total_time(p6)
X6 = to_array(p6.X)
U6 = to_array([p6.U[k][1:p6.model.m] for k = 1:p6.N])
U6m = to_array([0.5*(p6.U[k][1:p6.model.m] + p6.U[k+1][1:p6.model.m]) for k = 1:p6.N-1])
plot(X6[1,:],X6[2,:],title="Parallel Park Min. Time - (SNOPT)")
plot(U6',title="Parallel Park Min. Time - (SNOPT)")
plot(U6m',title="Parallel Park Min. Time - (control midpoint) (SNOPT)")

###############################################
#            Create PGF Plot                  #
###############################################
include("vars.jl")

t1 = trajectory_plot_flip(p1, mark="none", legendentry="original", style="very thick, color=$col_altro, mark options={fill=$col_altro}");
t2 = trajectory_plot_flip(p2, mark="none", legendentry="Ipopt", style="very thick, color=$col_ipopt, mark options={fill=$col_ipopt}");
t3 = trajectory_plot_flip(p3, mark="none", legendentry="SNOPT", style="very thick, color=$col_snopt, mark options={fill=$col_snopt}");
t4 = trajectory_plot_flip(p4, mark="none", legendentry="mintime", style="very thick, color=black, mark options={fill=$col_altro}, dashed");
t5 = trajectory_plot_flip(p5, mark="none", legendentry="Ipopt", style="very thick, color=black, mark options={fill=$col_ipopt}, dashed");
t6 = trajectory_plot_flip(p6, mark="none", legendentry="SNOPT", style="very thick, color=$col_snopt, mark options={fill=$col_snopt}, dashed");

x0 = prob_altro.x0
xf = prob_altro.xf
goal = ([x0[2], xf[2]],
        [x0[1], xf[1]])
z = ["start","end"]
g = PGF.Plots.Scatter(goal[1], goal[2], z,
    scatterClasses="{start={yellow, mark=*, yellow, scale=2},
        end={mark=square*, red, scale=2}}",
    );

a = Axis([t3; t2; t1; t6; t5; t4; g],
    xmin=-0.1, ymin=-0.3, xmax=1.1, ymax=0.3,
    axisEqualImage=true,
    legendPos="outer north east",
    hideAxis=true)

a = Axis([t1; t4; g],
    xmin=-0.1, ymin=-0.3, xmax=1.1, ymax=0.3,
    axisEqualImage=true,
    legendPos="outer north east",
    hideAxis=true)

# Save to tikz format
# NOTE: To fix the problem with the legend for the start and goal points, replace \addplot+ with \addplot in the tikz file
save(joinpath(paper,"ppark_traj.tikz"), a, include_preamble=false)

max_issues_repo_name: MehrdadZareian/GMARA

max_forks_repo_name: MehrdadZareian/GMARA


text:
using ITensorBenchmarks

# Run all of the benchmarks with a single thread and write the results
# By default, a verbose output is written to a log file
runbenchmarks(write_results = true)


text:
effective_particles(pf) = effective_particles(expweights(pf))
effective_particles(we::AbstractVector) = 1/sum(abs2, we)

function shouldresample(pf::AbstractParticleFilter)
    resample_threshold(pf) == 1 && (return true)
    th      = num_particles(pf)*resample_threshold(pf)
    ne = effective_particles(pf)
    return ne < th
end

resample(pf::AbstractParticleFilter, M::Int=num_particles(pf)) = resample(resampling_strategy(pf), expweights(pf), state(pf).j, state(pf).bins, M)
resample(T::Type{<:ResamplingStrategy}, s::PFstate, M::Int=num_particles(s)) = resample(T, s.we, s.j, s.bins, M)
resample(T::Type{<:ResamplingStrategy}, we::AbstractVector, M::Int=length(we)) = resample(T, we, zeros(Int,M), zeros(length(we)), M)
resample(we::AbstractArray) = resample(ResampleSystematic,we)

function resample(::Type{ResampleSystematic}, we, j, bins, M = length(we))
    N = length(we)
    bins[1] = we[1]
    for i = 2:N
        bins[i] = bins[i-1] + we[i]
    end
    r = rand()*bins[end]/N
    s = r:(1/M):(bins[N]+r) # Added r in the end to ensure correct length (r < 1/N)
    bo = 1
    for i = 1:M
        @inbounds for b = bo:N
            if s[i] < bins[b]
                j[i] = b
                bo = b
                break
            end
        end
    end
    return j
end

# """
# There is probably lots of room for improvement here. All bins need not be formed in the beginning.
# One only has to keep 1 values, the current upper limit, no array needed.
# """
"""
    draw_one_categorical(pf,w)
Obs! This function expects log-weights
"""
function draw_one_categorical(pf,w)
    bins = state(pf).bins
    logsumexp!(w,bins)
    for i = 2:length(w)
        bins[i] += bins[i-1]
    end
    @assert bins[end] ≈ 1 "All expweigths 0"
    s = rand()*bins[end]
    # ind = findfirst(x->bins[x]>=s, 1:length(bins))
    midpoint = length(bins)÷2
    if s < bins[midpoint]
        for b = 1:midpoint
            if s <= bins[b]
                return b
            end
        end
    else
        for b = midpoint:length(bins)
            if s <= bins[b]
                return b
            end
        end
    end
    length(bins)
end


text:
"""
    struct CoreModel <: MetabolicModel

A "bare bones" core linear optimization problem of the form, with reaction and
metabolite names.
```
min c^T x
s.t. S x = b
      xₗ ≤ x ≤ xᵤ
```
"""
mutable struct CoreModel <: MetabolicModel
    S::SparseMat
    b::SparseVec
    c::SparseVec
    xl::SparseVec
    xu::SparseVec
    rxns::Vector{String}
    mets::Vector{String}

    function CoreModel(
        S::M,
        b::V,
        c::V,
        xl::V,
        xu::V,
        rxns::K,
        mets::K,
    ) where {V<:VecType,M<:MatType,K<:StringVecType}

        all([length(b), length(mets)] .== size(S, 1)) ||
            throw(DimensionMismatch("inconsistent number of metabolites"))

        all([length(c), length(xl), length(xu), length(rxns)] .== size(S, 2)) ||
            throw(DimensionMismatch("inconsistent number of reactions"))

        new(sparse(S), sparse(b), sparse(c), sparse(xl), sparse(xu), rxns, mets)
    end
end

"""
    reactions(a::CoreModel)::Vector{String}

Get the reactions in a `CoreModel`.
"""
function reactions(a::CoreModel)::Vector{String}
    a.rxns
end

"""
    metabolites(a::CoreModel)::Vector{String}

Metabolites in a `CoreModel`.
"""
function metabolites(a::CoreModel)::Vector{String}
    a.mets
end

"""
    stoichiometry(a::CoreModel)::SparseMat

`CoreModel` stoichiometry matrix.
"""
function stoichiometry(a::CoreModel)::SparseMat
    a.S
end

"""
    bounds(a::CoreModel)::Tuple{SparseVec,SparseVec}

`CoreModel` flux bounds.
"""
function bounds(a::CoreModel)::Tuple{SparseVec,SparseVec}
    (a.xl, a.xu)
end

"""
    balance(a::CoreModel)::SparseVec

`CoreModel` target flux balance.
"""
function balance(a::CoreModel)::SparseVec
    a.b
end

"""
    objective(a::CoreModel)::SparseVec

`CoreModel` objective vector.
"""
function objective(a::CoreModel)::SparseVec
    a.c
end

"""
    reaction_stoichiometry(model::CoreModel, rxn_id::String)::Dict{String, Float64}

Return the reaction equation of reaction with id `rxn_id` in model. The reaction
equation maps metabolite ids to their stoichiometric coefficients.
"""
function reaction_stoichiometry(m::CoreModel, rxn_id::String)::Dict{String,Float64}
    Dict(
        m.mets[k] => v for
        (k, v) in zip(findnz(m.S[:, first(indexin([rxn_id], m.rxns))])...)
    )
end

"""
    reaction_stoichiometry(model::CoreModel, rxn_id::String)::Dict{String, Float64}

Return the reaction equation of reaction with id `rxn_ind` in model. The reaction
equation maps metabolite ids to their stoichiometric coefficients.
"""
function reaction_stoichiometry(m::CoreModel, rxn_ind::Int)::Dict{String,Float64}
    Dict(m.mets[k] => v for (k, v) in zip(findnz(m.S[:, rxn_ind])...))
end

"""
    Base.convert(::Type{CoreModel}, m::M) where {M <: MetabolicModel}

Make a `CoreModel` out of any compatible model type.
"""
function Base.convert(::Type{CoreModel}, m::M) where {M<:MetabolicModel}
    if typeof(m) == CoreModel
        return m
    end

    (xl, xu) = bounds(m)
    CoreModel(
        stoichiometry(m),
        balance(m),
        objective(m),
        xl,
        xu,
        reactions(m),
        metabolites(m),
    )
end


text:
function part1and2()
	# Input file is small enough that you can just enter it manually instead of parsing the file
	xmin = 257
	xmax = 286
	ymin = -101
	ymax = -57

	highestheight = 0
	validvelocities = 0

	for xvel in 1:xmax
		for yvel in ymin:100
			x, y = 0, 0
			xvelcur = xvel
			yvelcur = yvel
			highesty = 0
			success = true
			while !(xmin <= x <= xmax && ymin <= y <= ymax)
				x += xvelcur
				y += yvelcur
				if xvelcur > 0
					xvelcur -= 1
				elseif xvelcur < 0
					xvelcur += 1
				end
				yvelcur -= 1

				if y > highesty
					highesty = y
				end

				if x > xmax || y < ymin
					success = false
					break
				end
			end

			if success
				if highesty > highestheight
					highestheight = highesty
				end
				validvelocities += 1
			end
		end
	end

	return (highestheight, validvelocities)
end

println(part1and2())


text:
import Oceananigans.Grids: short_show

#####
##### BoundaryCondition
#####

Base.show(io::IO, bc::BC{C, T}) where {C, T} =
    print(io, "BoundaryCondition: type=$C, condition=$(bc.condition)")

#####
##### FieldBoundaryConditions
#####

bctype_str(::FBC)  = "Flux"
bctype_str(::PBC)  = "Periodic"
bctype_str(::NPBC) = "NoPenetration"
bctype_str(::VBC)  = "Value"
bctype_str(::GBC)  = "Gradient"
bctype_str(::NFBC) = "NoFlux"

short_show(fbcs::FieldBoundaryConditions) =
    string("x=(west=$(bctype_str(fbcs.x.left)), east=$(bctype_str(fbcs.x.right))), ",
           "y=(south=$(bctype_str(fbcs.y.left)), north=$(bctype_str(fbcs.y.right))), ",
           "z=(bottom=$(bctype_str(fbcs.z.left)), top=$(bctype_str(fbcs.z.right)))")

show_field_boundary_conditions(bcs::FieldBoundaryConditions, padding="") =
    string("Oceananigans.FieldBoundaryConditions (NamedTuple{(:x, :y, :z)}), with boundary conditions", '\n',
           padding, "├── x: ", typeof(bcs.x), '\n',
           padding, "├── y: ", typeof(bcs.y), '\n',
           padding, "└── z: ", typeof(bcs.z))

Base.show(io::IO, fieldbcs::FieldBoundaryConditions) = print(io, show_field_boundary_conditions(fieldbcs))


text:
# This file was generated, do not modify it. # hide
using Gnuplot, Colors, ColorSchemes
x = LinRange(-2π,2π,200)
#colorbar tricks... # hide
cbwt = 0.02 # hide
rightmargin = 0.875 # hide
cboxp="set colorbox user origin graph 1.01, graph 0 size $cbwt, graph 1" # hide
addmargin="set rmargin at screen $rightmargin" # hide
goffset="set offsets graph .05, graph .05, graph .05, graph .05" # hide
#custom palette, colormap # hide
function custom_palette(colormap=:viridis) # hide
cmap = get(colorschemes[colormap], LinRange(0,1,256)) # hide
ctmp = "0 '#$(hex(cmap[1]))'," # hide
for i in 2:256; ctmp = ctmp*"$(i-1) '#$(hex(cmap[i]))'," end; # hide
"set palette defined("*ctmp[1:end-1]*")" # hide
end # hide
@gp(x, -0.65sin.(3x), x,  "w l lw 3 dt 1 lc palette",
    "set key off", "set auto fix")
    #goffset, cboxp, addmargin) # hide
save(term="pngcairo size 600,400", output="plt_ex8.png") # hide


text:
using Sockets

const UPPER_BOUND = 5_000_000
const PREFIX = 32_338

mutable struct Node
    children::Dict
    terminal::Bool
end
Node(children::Dict) = Node(children, false)
Node() = Node(Dict(), false)

mutable struct Sieve
    limit::Int64
    prime::Vector{Bool}

    function Sieve(limit::Int64)
        return new(limit, falses(limit + 1))
    end
end

function to_list(s::Sieve)::Vector{Int64}
    result = [2, 3]
    @inbounds for p = 1:s.limit+1
        if s.prime[p] == true
            push!(result, p)
        end
    end
    return result
end

function omit_squares!(s::Sieve)
    r = 5
    while r^2 < s.limit
        @inbounds if s.prime[r]
            i = r^2
            while i < s.limit
                s.prime[i] = false
                i += r^2
            end
        end
        r += 1
    end
end

function step1!(s::Sieve, x::Int64, y::Int64)
    n = 4x^2 + y^2
    if n <= s.limit && (mod(n, 12) == 1 || mod(n, 12) == 5)
        @inbounds s.prime[n] = !s.prime[n]
    end
end

function step2!(s::Sieve, x::Int64, y::Int64)
    n = 3x^2 + y^2
    if n <= s.limit && mod(n, 12) == 7
        @inbounds s.prime[n] = !s.prime[n]
    end
end

function step3!(s::Sieve, x::Int64, y::Int64)
    n = 3x^2 - y^2
    if x > y && n <= s.limit && mod(n, 12) == 11
        @inbounds s.prime[n] = !s.prime[n]
    end
end

function loop_y!(s::Sieve, x::Int64)
    y = 1
    while y^2 < s.limit
        step1!(s, x, y)
        step2!(s, x, y)
        step3!(s, x, y)
        y += 1
    end
end

function loop_x!(s::Sieve)
    x = 1
    while x^2 < s.limit
        loop_y!(s, x)
        x += 1
    end
end

function calc(s::Sieve)
    loop_x!(s)
    omit_squares!(s)
    return s
end

function generate_trie(prime_list::Vector{Int64})::Node
    root = Node()
    for p in prime_list
        head = root
        for ch in string(p)
            if !(ch in keys(head.children))
                head.children[ch] = Node()
            end
            head = head.children[ch]
        end
        head.terminal = true
    end
    return root
end

function find(upper_bound::Int64, search_prefix::Int64)::Vector
    primes = calc(Sieve(upper_bound))
    str_prefix = string(search_prefix)

    head = generate_trie(to_list(primes))

    for char in str_prefix
        head = get(head.children, char, nothing)
        head == nothing && return nothing
    end

    queue = Vector{Tuple}([(head, str_prefix)])
    result = Vector{Int64}()

    while !isempty(queue)
        top, prefix = pop!(queue)

        if top.terminal == true
            push!(result, parse(Int64, prefix))
        end

        for (char, v) in top.children
            insert!(queue, 1, (v, prefix * char))
        end
    end

    return sort(result)
end

function notify(msg)
    try
        socket = connect("localhost", 9001)
        write(socket, msg)
        close(socket)
    catch
        # standalone usage
    end
end

function verify()
    left = [2, 23, 29]
    right = find(100, 2)

    if left != right
        println("$left != $right")
        exit()
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    verify()

    notify("Julia\t$(getpid())")
    results = find(UPPER_BOUND, PREFIX)
    notify("stop")

    println(results)
end

max_stars_repo_name: mkitto/benchmarks

max_issues_repo_name: mkitto/benchmarks

max_forks_repo_name: mkitto/benchmarks