text:
## Main functions of the algorithm

"""
    check_dimensionality(c::CrystalNet)

Check that the dimensionality of the net (i.e. the number of independent axes along which
it is periodic) is equal to `D`, or throw a DimensionMismatch otherwise.
"""
function check_dimensionality(c::CrystalNet{D}) where {D}
    s, (d, o) = e
    edgs = [c.pge.pos[d] .+ o .- c.pge.pos[s] for e in edges(c.pge.g)]
    sort!(edgs)
    unique!(edgs)
    mat = reduce(hcat, edgs)
    if D == 3
        isrank3(mat) || throw(DimensionMismatch("Internal error: the input net does not have expected dimensionality 3."))
    elseif D == 2
        isrank2(mat) || throw(DimensionMismatch("Internal error: the input net does not have expected dimensionality 2."))
    elseif D == 1
        isrank1(mat) || throw(DimensionMismatch("Internal error: the input net does not have expected dimensionality 1."))
    else
        throw(AssertionError("1 ≤ D ≤ 3"))
    end
    nothing
end

"""
    possible_translations(c::CrystalNet)

Return a list of tuples `(nz, i_max_den, max_den, t)` where
- `t` is a translation mapping at the origin vertex to another one in the unit cell.
- `max_den` is the maximum denominator in the `D` coefficients of `t`.
- `i_max_den` is the index.
- `nz` is the number of zeros in `t`.

The list is guaranteed to contain all the possible valid translations but may contain some
invalid translations.

See also: [`find_all_valid_translations`](@ref), `PeriodicGraphEmbeddings.check_valid_symmetry`
"""
function possible_translations(c::CrystalNet{D,T}) where {D,T}
    ts = Tuple{Int, Int, Int, SVector{D,T}}[]
    sortedpos = copy(c.pge.pos)
    origin = popfirst!(sortedpos)
    @toggleassert iszero(origin)
    sort!(sortedpos, by=norm)
    for t in sortedpos
        @toggleassert t == back_to_unit.(t)
        max_den, i_max_den = findmax(denominator.(t))
        numerator(t[i_max_den]) == 1 || continue
        nz = count(iszero, t)
        push!(ts, (nz, i_max_den, max_den, t))
    end
    return sort!(ts; by=(x->(x[1], x[2], x[3])))
end

"""
    find_all_valid_translations(c::CrystalNet{D}, collisions) where D

Return a `D`-tuple of list of tuples `(i_max_den, max_den, t)` (see
[`possible_translations`](@ref) for interpretation) where the `n`-th list contains all
valid translations of the net having exactly `n-1` zeros.

A translation is valid if it maps exactly each vertex to a vertex and each edge to an edge.

See also: [`possible_translations`](@ref), `PeriodicGraphEmbeddings.check_valid_symmetry`
"""
function find_all_valid_translations(shrunk_net::CrystalNet{D,T}, collisions) where {D,T}
    ret = NTuple{D, Vector{Tuple{Int, Int, SVector{D,T}}}}(ntuple(_->[], Val(D)))
    check_symmetry = check_symmetry_with_collisions(collisions)
    for (nz, i_max_den, max_den, t) in possible_translations(shrunk_net)
        vmap = check_symmetry(shrunk_net.pge, t, nothing, shrunk_net.types)
        if vmap isa Vector{PeriodicVertex{D}}
            push!(ret[nz+1], (i_max_den, max_den, t))
        end
    end
    return ret
end

"""
    minimal_volume_matrix(translations::NTuple{D}) where {D}

Given the output of [`find_all_valid_translations`](@ref), compute the transformation that
allows reducing the net to its minimal cell.
"""
function minimal_volume_matrix end

function minimal_volume_matrix(translations::Tuple{Vector{Tuple{Int,Int,SVector{1,T}}}}) where T
    nz0 = translations[1]
    denmax = 1
    imax = 0
    for j in 1:length(nz0)
        _, den, _ = nz0[j]
        if den > denmax
            imax = j
            denmax = den
        end
    end
    if imax == 0
        push!(nz0, (1, 1, [1]))
        imax = length(nz0)
    end
    _nz0 = [nz0[imax]]
    empty!(nz0)
    append!(nz0, _nz0)

    n = length(nz0)
    detmin = one(T)
    best = 0
    @inbounds for i in 1:n
        d = nz0[i][3][1]
        d == 0 && continue
        if d < detmin
            detmin = d
            best = i
        end
    end
    @toggleassert !iszero(best)
    ret = hcat(nz0[best][3])
    if ret[1] < 0
        ret = hcat(.-nz0[best][3])
    end
    return ret
end

function minimal_volume_matrix(translations::NTuple{2, Vector{Tuple{Int,Int,SVector{2,T}}}}) where T
    nz0, nz1 = translations
    denmax = [1, 1]
    imax = [0, 0]
    for j in 1:length(nz1)
        i, den, _ = nz1[j]
        if den > denmax[i]
            imax[i] = j
            denmax[i] = den
        end
    end
    for j in 1:2
        if imax[j] == 0
            push!(nz1, (j, 1, [j==1, j==2]))
            imax[j] = length(nz1)
        end
    end
    _nz1 = [nz1[i] for i in imax]
    empty!(nz1)
    append!(nz1, _nz1)

    # TODO optimize
    all = vcat(nz0, nz1)
    n = length(all)
    detmin = one(T)
    best = (0, 0)
    @inbounds for i in 1:n-1
        for j in i+1:n
            d = abs(det(hcat(all[i][3], all[j][3])))
            d == 0 && continue
            if d < detmin
                detmin = d
                best = (i, j)
            end
        end
    end
    @toggleassert !any(iszero.(best))
    i, j = best
    ret = hcat(all[i][3], all[j][3])
    if det(ret) < 0
        ret = hcat(all[i][3], .-all[j][3])
    end
    return ret
end

function minimal_volume_matrix(translations::NTuple{3, Vector{Tuple{Int,Int,SVector{3,T}}}}) where T
    nz0, nz1, nz2 = translations

    denmax = [1, 1, 1]
    imax = [0, 0, 0]
    for j in 1:length(nz2)
        i, den, _ = nz2[j]
        if den > denmax[i]
            imax[i] = j
            denmax[i] = den
        end
    end
    for j in 1:3
        if imax[j] == 0
            push!(nz2, (j, 1, [j==1, j==2, j==3]))
            imax[j] = length(nz2)
        end
    end
    _nz2 = [nz2[i] for i in imax]
    empty!(nz2)
    append!(nz2, _nz2)

    # TODO optimize
    all = vcat(nz0, nz1, nz2)
    n = length(all)
    detmin = one(T)
    best = (0, 0, 0)
    @inbounds for i in 1:n-2
        for j in i+1:n-1
            for k in j+1:n
                d = abs(det(hcat(all[i][3], all[j][3], all[k][3])))
                d == 0 && continue
                if d < detmin
                    detmin = d
                    best = (i, j, k)
                end
            end
        end
    end
    @toggleassert !any(iszero.(best))
    i, j, k = best
    ret = hcat(all[i][3], all[j][3], all[k][3])
    if det(ret) < 0
        ret = hcat(all[i][3], all[j][3], .-all[k][3])
    end
    return ret
end

"""
    reduce_with_matrix(c::CrystalNet, mat, shrunk_net, collisions)

Given the net and the output of `minimal_volume_matrix` computed on the valid translations
of the net, return the new net representing the initial net in the computed unit cell.
"""
function reduce_with_matrix(c::CrystalNet{D,Rational{T}}, mat, collisions) where {D,T}
    if D == 3
        cell = Cell(c.pge.cell, c.pge.cell.mat * mat)
    else
        _mat = SizedMatrix{3,3,BigFloat}(LinearAlgebra.I)
        _mat[1:D,1:D] .= mat
        cell = Cell(c.pge.cell, c.pge.cell.mat * _mat)
    end

    imat = T.(inv(mat)) # The inverse should only have integer coefficients
    poscol = (imat,) .* c.pge.pos
    n = length(poscol)
    offset = Vector{SVector{D,Int}}(undef, n)
    for (i, pos) in enumerate(poscol)
        ofs = floor.(Int, pos)
        offset[i] = ofs
        poscol[i] = pos .- ofs
    end
    I_sort = sort(1:n; by=i->(poscol[i], hash_position(offset[i])))
    _i = popfirst!(I_sort)
    @toggleassert iszero(offset[_i])
    I_kept = Int[_i]
    last_sortedcol = poscol[_i]

    for i in I_sort
        x = poscol[i]
        if x != last_sortedcol
            push!(I_kept, i)
            last_sortedcol = x
        end
    end

    kept_collisions = Int[]
    if !isempty(collisions)
        m = n - length(collisions)
        reorder = Vector{Int}(undef, length(I_kept))
        idx = 0
        for (j, i) in enumerate(I_kept)
            if i > m
                reorder[end-idx] = j
                idx += 1
                push!(kept_collisions, i-m)
            else
                reorder[j-idx] = j
            end
        end
        I_kept = I_kept[reorder]
    end

    sortedcol = SVector{D,Rational{T}}[SVector{D,Rational{T}}(poscol[i]) for i in I_kept]

    vmap = Vector{Int}(undef, n)
    for (i, pos) in enumerate(poscol)
        j = searchsortedfirst(sortedcol, pos)
        @toggleassert j <= length(sortedcol) && sortedcol[j] == pos
        vmap[i] = j
    end

    edges = PeriodicEdge{D}[]
    for i in 1:length(I_kept)
        ofs_i = offset[I_kept[i]]
        for neigh in neighbors(c.pge.g, I_kept[i])
            ofs_x = offset[neigh.v]
            push!(edges, (i, vmap[neigh.v], ofs_x - ofs_i .+ imat*neigh.ofs))
        end
    end

    newcollisions = [CollisionNode(collisions[i], vmap) for i in kept_collisions]
    for newnode in newcollisions
        if !allunique(newnode.neighs)
            newcollisions = nothing
            break
        end
    end

    graph = PeriodicGraph{D}(edges)
    return CrystalNet{D,Rational{T}}(cell, c.types[I_kept], sortedcol, graph, c.options), newcollisions
end

"""
    minimize(net::CrystalNet, [collisions::Vector{CollisionNode}])

Return a CrystalNet representing the same net as the input, but in a unit cell.
If `collisions` is given, also return the corresponding collisions after minimization.

The computed unit cell may depend on the representation of the input, i.e. it is not
topologicallly invariant.
"""
function minimize(net::CrystalNet, collisions::Vector{CollisionNode})
    translations = find_all_valid_translations(net, collisions)
    all(isempty.(translations)) && return net, collisions
    mat = minimal_volume_matrix(translations)
    _net, collisions = reduce_with_matrix(net, mat, collisions)
    @toggleassert all(isempty.(find_all_valid_translations(_net, collisions)))
    return _net, collisions
end
minimize(net::CrystalNet) = minimize(net, CollisionNode[])[1]

function findfirstbasis(offsets)
    newbasis = normal_basis_rational(offsets)
    invbasis = inv(newbasis)
    intcoords = [Int.(invbasis * x) for x in offsets]
    return newbasis, intcoords
end

function findbasis(edges::Vector{Tuple{Int,Int,SVector{D,T}}}) where {D,T}
    m = length(edges)
    positivetrans = SVector{D,T}[]
    map_to_ofs = Vector{Int}(undef, m)
    tmp_map = Int[]
    for i in 1:m
        trans = last(edges[i])
        c = cmp(trans, zero(SVector{D,Int}))
        if iszero(c)
            map_to_ofs[i] = 0
        elseif c < 0
            push!(tmp_map, -i)
            push!(positivetrans, .-trans)
        else
            push!(tmp_map, i)
            push!(positivetrans, trans)
        end
    end

    I_sort = sortperm(positivetrans)
    uniques = SVector{D,T}[]
    last_trans = zero(SVector{D,T})

    for j in I_sort
        i = tmp_map[j]
        trans = positivetrans[j]
        if trans != last_trans
            push!(uniques, trans)
            last_trans = trans
        end
        map_to_ofs[abs(i)] = sign(i) * length(uniques)
    end

    basis, unique_coords = findfirstbasis(uniques)

    newedges = Vector{PeriodicEdge{D}}(undef, m)
    for j in 1:m
        i = map_to_ofs[j]
        src, dst, _ = edges[j]
        if i == 0
            newedges[j] = PeriodicEdge{D}(src, dst, zero(SVector{D,Int}))
        else
            newedges[j] = PeriodicEdge{D}(src, dst, sign(i) .* unique_coords[abs(i)])
        end
    end

    @toggleassert all(z-> begin x, y = z; x == (y.src, y.dst.v, basis*y.dst.ofs) end, zip(edges, newedges))

    return basis, newedges
end

"""
    candidate_key(net::CrystalNet, u, basis, minimal_edgs)

Given the net, a candidate `u => basis` where `u` is the origin and `basis` the triplet of
axes, and `minimal_edgs` the last minimal key (for the pseudo-lexicographical order used),
extract the key corresponding to the current candidate.

The key is the lexicographically ordered list of edges of the graph when its vertices are
numbered according to the candidate. The ordering of keys first compares the list of edges
disregarding the offsets, and then only compares the offsets if the rest is identical.

If the key is larger or equal to `minimal_edgs`, early stop and return two empty lists.
Otherwise, the extracted key is the current best: return the vmap between the initial
vertices and their ordered image in the candidate, as well as the key.

See also: [`find_candidates`](@ref)
"""
function candidate_key(net::CrystalNet{D,T}, u, basis, minimal_edgs) where {D,T}
    n = nv(net.pge.g)
    h = 2 # next node to assign
    origin = net.pge.pos[u]
    newpos = Vector{SVector{D,T}}(undef, n) # positions of the kept representatives
    newpos[1] = zero(SVector{D,T})
    offsets = Vector{SVector{D,Int32}}(undef, n) # offsets of the new representatives with
    # respect to the original one, in the original basis
    offsets[1] = zero(SVector{D,Int32})
    vmap = Vector{Int}(undef, n) # bijection from the old to the new node number
    vmap[1] = u
    rev_vmap = zeros(Int, n) # inverse of vmap
    rev_vmap[u] = 1
    flag_bestedgs = false # marks whether the current list of edges is lexicographically
    # below the best known one. If not by the end of the algorithm, return false
    edgs = Tuple{Int,Int,SVector{D,T}}[]
    bigbasis = T == Rational{BigInt} ? basis : widen(T).(basis)
    mat = T == Rational{BigInt} ? inv(bigbasis) : T.(inv(bigbasis))
    for t in 1:n # t is the node being processed
        neighs = neighbors(net.pge.g, vmap[t])
        ofst = offsets[t]
        pairs = Vector{Tuple{SVector{D,T},Int}}(undef, length(neighs))
        for (i,x) in enumerate(neighs)
            pairs[i] = ((mat*(net.pge.pos[x.v] .+ x.ofs .- origin .+ ofst)), x.v)
        end
        # (x,i) ∈ pairs means that vertex i (in the old numerotation) has position x in the new basis
        order = unique(last.(sort(pairs)))
        inv_order = Vector{Int}(undef, n)
        for i in 1:length(order)
            inv_order[order[i]] = i
        end
        sort!(pairs, lt = ((x, a), (y, b)) -> begin inv_order[a] < inv_order[b] ||
                                                    (inv_order[a] == inv_order[b] && x < y) end)
        # pairs is sorted such that different nodes first appear in increasing order of their position
        # but different representatives of the same node are contiguous and also sorted by position.
        for (coordinate, v) in pairs
            idx = rev_vmap[v]
            if idx == 0 # New node to which h is assigned
                @toggleassert t < h
                vmap[h] = v
                rev_vmap[v] = h
                newpos[h] = coordinate
                offsets[h] = SVector{D,Int32}(bigbasis * coordinate .+ origin .- net.pge.pos[v])
                push!(edgs, (t, h, zero(SVector{D,T})))
                h += 1
            else
                realofs = coordinate .- newpos[idx]
                # offset between this representative of the node and that which was first encountered
                push!(edgs, (t, idx, realofs))
            end
            if !flag_bestedgs
                j = length(edgs)
                c = cmp(minimal_edgs[j], edgs[j])
                c < 0 && return (Int[], Tuple{Int,Int,SVector{D,T}}[]) # early stop
                c > 0 && (flag_bestedgs = true)
            end
        end
    end
    @toggleassert allunique(edgs)
    if !flag_bestedgs # the current list of edges is equal to minimal_edgs
        @toggleassert minimal_edgs == edgs
        return (Int[], Tuple{Int,Int,SVector{D,T}}[])
    end
    return vmap, edgs
end

"""
    partition_by_coordination_sequence(graph, symmetries::AbstractSymmetryGroup=NoSymmetryGroup(graph))

Partition the vertices of the graph into disjoint categories, one for each coordination
sequence. The partition is then sorted by order of coordination sequence.
This partition does not depend on the representation of the graph.
The optional argument `vmaps` is a set of permutations of the vertices that leave the graph
unchanged. In other words, `vmaps` is a set of symmetry operations of the graph.

Return the categories and a list of unique representative for each symmetry class.
"""
function partition_by_coordination_sequence(graph, symmetries::AbstractSymmetryGroup=NoSymmetryGroup(graph))
    n = nv(graph)
    uniques::Vector{Int} = unique(symmetries)
    unique_reprs = [Int[repr] for repr in uniques]
    csequences = [coordination_sequence(graph, repr, 10) for repr in uniques]
    I = sortperm(csequences)
    @toggleassert csequences[I[1]][1] >= 2 # vertices of degree <= 1 should have been removed at input creation

    if symmetries isa NoSymmetryGroup
        categories = [Int[i] for i in 1:n]
    else
        rev_uniques = Vector{Int}(undef, n)
        for (i, repr) in enumerate(uniques)
            rev_uniques[repr] = i
        end
        categories = [Int[] for _ in 1:length(uniques)]
        for i in 1:n
            push!(categories[rev_uniques[symmetries(i)]], i)
        end
    end

    todelete = falses(length(categories))
    last_i = I[1]
    for j in 2:length(I)
        i = I[j]
        if csequences[i] == csequences[last_i]
            todelete[i] = true
            append!(categories[last_i], categories[i])
            push!(unique_reprs[last_i], unique_reprs[i][1])
            # We asserted that unique_reprs[i] only had one element before
            # Note that unique_reprs[last_i] has more than one element now
        else
            last_i = i
        end
    end
    deleteat!(categories, todelete)
    deleteat!(unique_reprs, todelete)
    deleteat!(csequences, todelete)

    num = length(categories)
    numoutgoingedges = Vector{Tuple{Int,Vector{Int}}}(undef, num)
    for i in 1:num
        seq = csequences[i]
        numoutgoingedges[i] = (length(categories[i]) * seq[1], seq)
    end
    sortorder = sortperm(numoutgoingedges)
    # categories are sorted by the total number of outgoing directed edges from each category
    # and by the coordination sequence in case of ex-aequo.
    # categories are thus uniquely determined and ordered independently of the representation of the net

    @toggleassert allunique(csequences)
    # for i in 1:num
    #     @toggleassert all(coordination_sequence(graph, x, 10) == csequences[i] for x in categories[i])
    # end # if enabled, these assertions are somewhat costly (up to ~10% total execution time)
    return categories[sortorder], unique_reprs[sortorder]
end

function check_symmetry_with_collisions(collisions)
    function check_symmetry(pge::PeriodicGraphEmbedding{D,T}, t::SVector{D,T}, r, vtypes) where {D,T}
        vmap = check_valid_symmetry(pge, t, r, vtypes, true)
        vmap isa Nothing && return nothing
        n = length(pge.pos)
        m = n - length(collisions)
        for (i, node) in enumerate(collisions)
            j = vmap[m+i].v
            m + i == j && continue
            j ≤ m && return nothing # a non-collision node is mapped to a collision node
            collisions[j-m] == CollisionNode(node, [x.v for x in vmap]) || return nothing
        end
        return vmap
    end
    check_symmetry
end

"""
    find_candidates(net::CrystalNet{D}, collisions::Vector{CollisionNode}) where D

Return a non-empty set of candidates `u => basis` where `u` is a vertex and `basis` is
matrix whose columns are `D` linearly independent euclidean embeddings of edges.
The returned set is independent of the representation of the graph used in `net`.

Also return a `category_map` linking each vertex to its category number, as defined by
[`partition_by_coordination_sequence`](@ref)

See also: [`candidate_key`](@ref)
"""
function find_candidates(net::CrystalNet{D,T}, collisions::Vector{CollisionNode}) where {D,T}
    L = D*D
    if D == 3
        check_symmetry = check_symmetry_with_collisions(collisions)
        symmetries = find_symmetries(net.pge, net.types, check_symmetry)
        categories, unique_reprs = partition_by_coordination_sequence(net.pge.g, symmetries)
    else
        categories, unique_reprs = partition_by_coordination_sequence(net.pge.g)
    end

    category_map = Vector{Int}(undef, nv(net.pge.g))
    for (i, cat) in enumerate(categories)
        for j in cat
            category_map[j] = i
        end
    end
    # @toggleassert sort.(categories) == sort.(first(partition_by_coordination_sequence(net.pge.g)))
    candidates = Dict{Int,Vector{SMatrix{D,D,T,L}}}()
    for reprs in unique_reprs
        # First, we try to look for triplet of edges all starting from the same vertex within a category
        degree(net.pge.g, first(reprs)) <= D && continue
        candidates = find_candidates_onlyneighbors(net, reprs, category_map)
        !isempty(candidates) && break
    end
    if D >= 3 && isempty(candidates)
        # If we arrive at this point, it means that all vertices only have coplanar neighbors
        # and D >= 3
        # Then we look for all triplets of edges two of which start from a vertex
        # and one does not, stopping at the first pair of categories for which this
        # set of candidates is non-empty.
        if isempty(candidates)
            for reprs in unique_reprs
                candidates = find_candidates_fallback(net, reprs, categories, category_map)
                !isempty(candidates) && break
            end
        end
    end
    if isempty(candidates)
        check_dimensionality(net)
        error("Internal error: no candidate found.")
    end
    if D == 3
        return extract_through_symmetry(candidates, symmetries), category_map
    else
        flattened_candidates = Pair{Int,SMatrix{D,D,T,L}}[]
        for (i, mats) in candidates
            for mat in mats
                push!(flattened_candidates, i => SMatrix{D,D,T,L}(mat))
            end
        end
        return flattened_candidates, category_map
    end
end

"""
    extract_through_symmetry(candidates::Dict{Int,Vector{SMatrix{3,3,T,9}}}, symmetries::AbstractSymmetryGroup) where T

Given the candidates and the list of symmetries of the net, return the flattened list of
candidates after removing candidates that are symmetric images of the kept ones.
"""
function extract_through_symmetry(candidates::Dict{Int,Vector{SMatrix{3,3,T,9}}}, symmetries::AbstractSymmetryGroup) where T
    unique_candidates = Pair{Int,SMatrix{3,3,T,9}}[]
    for (i, mats) in candidates
        @toggleassert i == symmetries(i)
        symms = [symm for symm in symmetries if symm[i] == i] # vmaps for which i is a fixpoint
        min_mats = Set{SVector{9,T}}()
        for mat in mats
            # At this point, `i => mat` is a potential candidate, and so are all its
            # symmetric images. Among those, only the ones of the form `i => mat2` (i.e.
            # with the same origin `i`) have not been eliminated earlier by symmetry.
            # Let's enumerate all such candidates and only keep one of them, for example
            # the lexicographically smallest, since they are all equivalent to `i => mat`.
            min_mat = SVector{9,T}(mat) # flattened to allow lexicographic comparison
            for symm in symms
                new_mat = SVector{9,T}(symm(mat))
                if new_mat < min_mat
                    min_mat = new_mat
                end
            end
            push!(min_mats, min_mat) # min_mats is a Set to eliminate duplicates
        end
        for x in min_mats
            push!(unique_candidates, i => SMatrix{3,3,T,9}(x))
        end
    end
    return unique_candidates
end

"""
    find_initial_candidates(net::CrystalNet{D}, candidates_v, category_map) where D

Given the net, a list of vertices in a given category and the `category_map`, return a
list of pairs `u => (basis, cats)` where `u ∈ candidates_v`, `basis` is a `D`-rank matrix
made by juxtaposing the euclidean embeddings of outgoing edges from `u`, and `cats` are the
categories of the respective neighbors of `u`.

If the `basis` corresponding to vertex `u` is not of rank `D`, it is not included in the
returned list (for instance, if all outgoing edges of a vertex are coplanar with `D == 3`).
"""
function find_initial_candidates(net::CrystalNet{D,T}, candidates_v, category_map) where {D,T}
    @toggleassert 1 ≤ D ≤ 3
    deg = degree(net.pge.g, first(candidates_v)) # The degree is the same for all vertices of the same category
    n = length(candidates_v)
    _initial_candidates = Vector{Pair{Int,Tuple{Matrix{T},Vector{Int}}}}(undef, n)
    valid_initial_candidates = falses(n)

    @threads for i in 1:n
        v = candidates_v[i]
        a = Matrix{T}(undef, D, deg)
        cats = Vector{Int}(undef, deg)
        posi = net.pge.pos[v]
        for (j, x) in enumerate(neighbors(net.pge.g, v))
            a[:,j] .= net.pge.pos[x.v] .+ x.ofs .- posi
            cats[j] = category_map[x.v]
        end
        if (D == 3 && isrank3(a)) || (D == 2 && isrank2(a)) || (D == 1 && ((@toggleassert a[:,1] != 0); true))
            _initial_candidates[i] = (v => (a, cats))
            valid_initial_candidates[i] = true
        end
    end
    return [_initial_candidates[i] for i in 1:n if valid_initial_candidates[i]]
end

"""
    find_candidates_onlyneighbors(net::CrystalNet{D}, candidates_v, category_map) where D

Given the net, a list of vertices in a given category and the `category_map`, return a Dict
whose pairs `u => matv` are such that `u ∈ candidates_v` and `matv` is a list of unique
invertible matrices of size `D` whose columns are euclidean embeddings of outgoing edges
from `u`.
Each such matrix has a category, defined by the `D`-uplet of categories of each
corresponding outneighbor of `u`: the returned Dict is such that all the matrices belonging
to all `matv` share the same category.

The returned Dict is empty iff `find_initial_candidates(net, candidates_v, category_map)`
is empty.
"""
function find_candidates_onlyneighbors end

function find_candidates_onlyneighbors(net::CrystalNet3D{T}, candidates_v, category_map) where T
    initial_candidates = find_initial_candidates(net, candidates_v, category_map)
    candidates = Dict{Int,Vector{SMatrix{3,3,T,9}}}()
    isempty(initial_candidates) && return candidates
    mincats::SVector{3,Int} = SVector{3,Int}((length(category_map), length(category_map), length(category_map)))
    ordertype::Int = 1
    # ordertype designates the kind of ordering of the category of the three edges:
    # ordertype == 1 means c1 == c2 == c3 where ci is the category of edge i
    # ordertype == 2 means c1  < c2 == c3
    # ordertype == 3 means c1 == c2  < c3
    # ordertype == 4 means c1  < c2  < c3

    for (v, (mat, cats)) in initial_candidates
        _, n = size(mat)
        matv = SMatrix{3,3,T,9}[] # the list of bases making a candidate with origin v
        for _i in 1:(n-2), _j in (_i+1):(n-1), _k in (_j+1):n
            m = SMatrix{3,3,T,9}(mat[:,[_i,_j,_k]])
            issingular(m) && continue
            orders = SVector{3,Int}[]
            subcats = cats[SVector{3,Int}(_i, _j, _k)]
            reorder = SVector{3,Int}(begin
                cmp23 = subcats[2] >= subcats[3]
                cmp13 = subcats[1] >= subcats[3]
                if subcats[1] >= subcats[2]
                    cmp23 ? (3,2,1) : (cmp13 ? (2,3,1) : (2,1,3))
                else
                    cmp23 ? (cmp13 ? (3,1,2) : (1,3,2)) : (1,2,3)
                end
            end)
            subcats = subcats[reorder]
            # subcats are the ordered categories of the three edges
            if subcats[1] == subcats[3] # c1 == c2 == c3
                (ordertype != 1 || mincats[1] < subcats[1]) && continue
                if mincats[1] > subcats[1]
                    empty!(matv)
                    empty!(candidates)
                    mincats = subcats
                end
                orders = SVector{3,Int}[(1,2,3), (1,3,2), (2,1,3), (2,3,1), (3,1,2), (3,2,1)]
            elseif subcats[2] == subcats[3] # c1 < c2 == c3
                (ordertype > 2 || (ordertype == 2 && mincats < subcats)) && continue
                if ordertype == 1 || mincats > subcats
                    empty!(matv)
                    empty!(candidates)
                    ordertype = 2
                    mincats = subcats
                end
                orders = SVector{3,Int}[reorder, (reorder[1], reorder[3], reorder[2])]
            elseif subcats[1] == subcats[2] # c1 == c2 < c3
                (ordertype == 4 || (ordertype == 3 && mincats < subcats)) && continue
                if ordertype <= 2 || mincats > subcats
                    empty!(matv)
                    empty!(candidates)
                    ordertype = 3
                    mincats = subcats
                end
                orders = SVector{3,Int}[reorder, (reorder[2], reorder[1], reorder[3])]
            else # c1 < c2 < c3
                ordertype == 4 && mincats < subcats && continue
                if ordertype != 4 || mincats > subcats
                    empty!(matv)
                    empty!(candidates)
                    ordertype = 4
                    mincats = subcats
                end
                orders = SVector{3,Int}[reorder]
            end
            append!(matv, m[:,o] for o in orders)
        end
        candidates[v] = matv
    end

    return candidates
end

function find_candidates_onlyneighbors(net::CrystalNet2D{T}, candidates_v, category_map) where T
    initial_candidates = find_initial_candidates(net, candidates_v, category_map)

    candidates = Dict{Int,Vector{SMatrix{2,2,T,4}}}()
    isempty(initial_candidates) && return candidates
    mincats::SVector{2,Int} = SVector{2,Int}((length(category_map), length(category_map)))
    ordertype::Int = 1

    for (v, (mat, cats)) in initial_candidates
        _, n = size(mat)
        matv = SMatrix{2,2,T,4}[]
        for _i in 1:(n-1), _j in (_i+1):n
            m = SMatrix{2,2,T,4}(mat[:,[_i,_j]])
            issingular(m) && continue
            orders = SVector{2,Int}[]
            subcats = cats[SVector{2,Int}(_i, _j)]
            reorder = SVector{2,Int}(subcats[1] >= subcats[2] ? (2,1) : (1,2))
            subcats = subcats[reorder]
            if subcats[1] == subcats[2]
                (ordertype != 1 || mincats[1] < subcats[1]) && continue
                if mincats[1] > subcats[1]
                    empty!(matv)
                    empty!(candidates)
                    mincats = subcats
                end
                orders = SVector{2,Int}[(1,2), (2,1)]
            else
                ordertype == 2 && mincats < subcats && continue
                if ordertype != 2 || mincats > subcats
                    empty!(matv)
                    empty!(candidates)
                    ordertype = 2
                    mincats = subcats
                end
                orders = SVector{2,Int}[reorder]
            end
            append!(matv, m[:,o] for o in orders)
        end
        candidates[v] = matv
    end

    return candidates
end

function find_candidates_onlyneighbors(net::CrystalNet1D{T}, candidates_v, category_map) where T
    initial_candidates = find_initial_candidates(net, candidates_v, category_map)

    candidates = Dict{Int,Vector{SMatrix{1,1,T,1}}}()
    isempty(initial_candidates) && return candidates
    mincat::Int = length(category_map)

    for (v, (mat, cats)) in initial_candidates
        n = length(mat)
        matv = SMatrix{1,1,T,1}[]
        for i in 1:n
            cat = cats[i]
            mincat < cat && continue
            if mincat > cat
                empty!(matv)
                empty!(candidates)
                mincat = cat
            end
            push!(matv, mat[[i]])
        end
        candidates[v] = matv
    end

    return candidates
end

"""
    find_candidates_fallback(net::CrystalNet3D, reprs, othercats, category_map)

Return candidates in the same form as [`find_candidates_onlyneighbors`](@ref) except that
only two edges start from `u` and one does not.
"""
function find_candidates_fallback(net::CrystalNet3D{T}, reprs, othercats, category_map) where T
    candidates = Dict{Int,Vector{SMatrix{3,3,T,9}}}(u => [] for u in reprs)
    n = length(category_map)
    mincats = SizedVector{3,Int}(fill(n, 3))
    current_cats = SizedVector{3,Int}(fill(0, 3))
    for cat in othercats
        for u in reprs
            mats = SMatrix{3,3,T,9}[]
            posu = net.pge.pos[u]
            neighu = neighbors(net.pge.g, u)
            for (j1, _x1) in enumerate(neighu)
                category_map[_x1.v] > mincats[1] && category_map[_x1.v] > mincats[2] && continue
                for j2 in j1+1:length(neighu)
                    x2 = neighu[j2]
                    x1 = _x1
                    if category_map[x1.v] > category_map[x2.v]
                        x1, x2 = x2, x1
                    end
                    current_cats[1] = category_map[x1.v]
                    current_cats[1] > mincats[1] && continue
                    current_cats[2] = category_map[x2.v]
                    current_cats[1] == mincats[1] && current_cats[2] > mincats[2] && continue
                    vec1 = net.pge.pos[x1.v] .+ x1.ofs .- posu
                    vec2 = net.pge.pos[x2.v] .+ x2.ofs .- posu
                    j = iszero(vec1[1]) ? iszero(vec1[2]) ? 3 : 2 : 1
                    vec1[j] * vec2 == vec2[j] * vec1 && continue # colinear
                    for v in cat
                        posv = net.pge.pos[v]
                        for x3 in neighbors(net.pge.g, v)
                            vec3 = net.pge.pos[x3.v] .+ x3.ofs .- posv
                            mat = SMatrix{3,3,T,9}(hcat(vec1, vec2, vec3))
                            issingular(mat) && continue
                            current_cats[3] = category_map[x3.v]
                            current_cats > mincats && continue
                            if current_cats < mincats
                                mincats .= current_cats
                                foreach(empty!, values(candidates))
                                empty!(mats)
                            end
                            # if d < 0
                            #     mat = hcat(vec2, vec1, vec3)
                            # end
                            push!(mats, mat)
                            if current_cats[1] == current_cats[2]
                                push!(mats, hcat(vec2, vec1, vec3))
                            end
                        end
                    end
                end
            end
            append!(candidates[u], mats)
        end
        # If at this point candidates is not empty, then we have found the first
        # other category cat such that there is a candidate with two edges from
        # maincat and one from cat. This is uniquely determined independently of
        # the representation of the net, so we can simply return the found candidates
        # without needing to check for the remaining categories.
        if !all(isempty, values(candidates))
            for (u, mats) in candidates
                if isempty(mats)
                    delete!(candidates, u)
                end
            end
            return candidates
         end
    end
    @toggleassert all(isempty, values(candidates))
    return Dict{Int,Vector{SMatrix{3,3,T,9}}}()
end

"""
    topological_key(net::CrystalNet)

Return a unique topological key for the net, which is a topological invariant of the net
(i.e. it does not depend on its initial representation).
"""
function topological_key(net::CrystalNet{D}) where D
    isempty(net.pge.pos) && return PeriodicGraph{D}()
    newnet, collisions = collision_nodes(net)
    if collisions isa Nothing
        net.pge.g.width[] = -2 # internal error code
        return net.pge.g
    end
    return topological_key(newnet, collisions::Vector{CollisionNode})
end

function topological_key(net::CrystalNet{D,T}, collisions::Vector{CollisionNode}) where {D,T}
    candidates, category_map = find_candidates(net, collisions)
    v, minimal_basis = popfirst!(candidates)
    n = nv(net.pge.g)
    _edgs = [(n + 1, 0, zero(SVector{D,T}))]
    minimal_vmap, minimal_edgs = candidate_key(net, v, minimal_basis, _edgs)
    for (v, basis) in candidates
        vmap, edgs = candidate_key(net, v, basis, minimal_edgs)
        isempty(vmap) && continue
        @toggleassert edgs < minimal_edgs
        if edgs < minimal_edgs
            minimal_edgs = edgs
            minimal_vmap = vmap
            minimal_basis = basis
        end
    end

    newbasis, edges = findbasis(minimal_edgs)
    graph = PeriodicGraph{D}(n, edges)

    if !isempty(collisions)
        graph = expand_collisions(collisions, graph, minimal_vmap)
    end

    # finalbasis = minimal_basis * newbasis
    # return Int.(finalbasis), minimal_vmap, graph
    return graph
end


text:
@testset "LightObservable" begin

    @testset "General" begin

        # constructor
        @test typeof(LightObservable(Float64)) == LightObservable{Float64, 10}
        @test typeof(LightObservable(ComplexF64)) == LightObservable{ComplexF64, 10}
        @test typeof(LightObservable(zeros(2,2))) == LightObservable{Matrix{Float64}, 10}
        @test typeof(LightObservable(zeros(ComplexF64, 2,2))) == LightObservable{Matrix{ComplexF64}, 10}

        @test typeof(LightObservable(Int64)) == LightObservable{Float64, 10}
        @test typeof(LightObservable(zeros(Int64, 2,2))) == LightObservable{Matrix{Float64}, 10}

        # size etc.
        @test size(LightObservable(rand(10))) == () # TODO
        @test size(LightObservable([rand(2,2) for _ in 1:3])) == (2,2)
        @test ndims(LightObservable(rand(10))) == 0
        @test ndims(LightObservable([rand(2,2) for _ in 1:3])) == 2

        # name
        obs = LightObservable(Float64, name="myobs")
        @test name(obs) == "myobs"
        @test name(LightObservable(ComplexF64, name="julia")) == "julia"

        # adding and reading
        obs = LightObservable(Float64)
        @test inmemory(obs)
        @test isempty(obs)
        @test push!(obs, 1.0) == nothing
        @test length(obs) == 1
        @test append!(obs, 2.0:4.0) == nothing
        @test length(obs) == 4
        @test append!(obs, 5.0:9.0) == nothing
        @test length(obs) == 9
        @test push!(obs, 10.0) == nothing
        @test length(obs) == 10
        @test !isempty(obs)

        # more than alloc test
        @test typeof(LightObservable(Float64, alloc=1023)) == LightObservable{Float64,10}
        @test typeof(LightObservable(Float64, alloc=2^2 - 1)) == LightObservable{Float64,2}
        @test typeof(LightObservable(Float64, alloc=2^5 - 1)) == LightObservable{Float64,5}

        # adding matrix observables
        ots = [rand(ComplexF64, 2,2) for _ in 1:3]
        obs = LightObservable(zero(ots[1]))
        @test push!(obs, ots[1]) == nothing
        @test append!(obs, ots[2:3]) == nothing
        @test_throws MethodError push!(obs, rand(["a", "b"]))
        @test_throws DimensionMismatch push!(obs, rand(ComplexF64, 2,2,3,4))

        # adding vector observables
        ots = [rand(ComplexF64, 3) for _ in 1:3]
        obs = LightObservable(zero(ots[1]))
        @test push!(obs, ots[1]) == nothing
        @test append!(obs, ots[2:3]) == nothing
        @test_throws MethodError push!(obs, rand(["a", "b"]))
        @test_throws DimensionMismatch push!(obs, rand(ComplexF64, 4))
    end

    @testset "Statistics" begin
        @testset "Real Observables" begin
            ots = [0.00124803, 0.643089, 0.183268, 0.799899, 0.0857666, 0.955348, 0.165763, 0.765998, 0.63942, 0.308818]
            obs = LightObservable(ots)
            @test isapprox(mean(ots), mean(obs))
            @test isapprox(std(ots), std(obs))
            @test isapprox(var(ots), var(obs))

            @test_deprecated error(obs)
            @test std_error(obs) == 0.1083461975750141
            @test tau(obs) == 0.0
        end

        @testset "Complex Observables" begin
            ots = Complex{Float64}[0.458585+0.676913im, 0.41603+0.0800011im, 0.439703+0.472044im, 0.86602+0.756838im, 0.615955+0.312498im, 0.916813+0.150829im, 0.434218+0.839293im, 0.888952+0.648892im, 0.799521+0.734382im, 0.678336+0.810805im]
            obs = LightObservable(ots)
            @test isapprox(mean(ots), mean(obs))
            @test isapprox(std(ots), std(obs))
            @test isapprox(var(ots), var(obs))

            @test_deprecated error(obs)
            @test std_error(obs) == 0.10946759231383459
            @test tau(obs) == 0.0
        end

        @testset "Matrix Observables" begin
            ots = Array{Float64,2}[[0.127479 0.144452; 0.0934332 0.465612], [0.716647 0.576685; 0.44389 0.256331], [0.811945 0.457262; 0.634971 0.188656]]
            obs = LightObservable(ots)
            @test isapprox(mean(ots), mean(obs))
            @test isapprox(std(ots), std(obs))
            @test isapprox(var(ots), var(obs))

            @test_deprecated error(obs)
            @test all(isapprox.(std_error(obs), [0.214048 0.128871; 0.158569 0.083361], atol=1e-6))
            @test all(isapprox.(tau(obs), [0.0 0.0; 0.0 0.0], atol=1e-6))
        end
    end

    @testset "IO" begin
        mktempdir() do d
            cd(d) do
                obs = LightObservable(rand(10))

                @test export_result(obs, "myresults.jld", "myobservables") == nothing
                or = load_result("myresults.jld", "myobservables")
                @test typeof(or) == ObservableResult{Float64,Float64}
                rm("myresults.jld")

                # import BinningAnalysis # necessary! can we avoid this?
                # saveobs(obs, "myobs.jld", "myobservables/obs")
                # saveobs(obs, "myobs.jld", "myobservables/obs_again") # test writing to already existing file
                # x = loadobs("myobs.jld", "myobservables/obs")
                # @test mean(x) == mean(obs)
                # @test std_error(x) == std_error(obs)
                # @test "obs" in listobs("myobs.jld", "myobservables/")
                # rmobs("myobs.jld", "obs", "myobservables/")
                # @test !("obs" in listobs("myobs.jld", "myobservables/"))
                # rm("myobs.jld")
            end
        end
    end

end

ext: jl


text:
function simulated_annealing(tuning_run::Run,
                             channel::RemoteChannel;
                             temperature::Function = log_temperature)
    iteration = 1

    function iterate(tuning_run::Run)
        p          = temperature(iteration)
        iteration += 1
        return probabilistic_improvement(tuning_run, threshold = p)
    end

    technique(tuning_run,
              channel,
              iterate,
              name = "Simulated Annealing")
end


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

commit_sha = "11bfbf509088f168622b8268f49c0a59ee81758b"

# Collection of sources required to build Cuba
name = "Cuba"
version = v"4.2a"
sources = [
    "https://github.com/giordano/cuba/archive/$commit_sha.tar.gz" =>
    "9cbb3a9c6ea541d7f3b0efcba6a865082e70536aded6a347bcec4873df3f3cc4",

]

# Bash recipe for building across all platforms
script = """
cd \$WORKSPACE/srcdir
cd cuba-$commit_sha/
./configure --prefix=\$prefix --host=\$target
make shared
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:i686, :glibc),
    Linux(:x86_64, :glibc),
    Linux(:aarch64, :glibc),
    Linux(:armv7l, :glibc, :eabihf),
    Linux(:powerpc64le, :glibc),
    Linux(:i686, :musl),
    Linux(:x86_64, :musl),
    Linux(:aarch64, :musl),
    Linux(:armv7l, :musl, :eabihf),
    MacOS(:x86_64),
    Windows(:i686),
    Windows(:x86_64),
    FreeBSD(:x86_64)
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libcuba", :libcuba)
]

# Dependencies that must be installed before this package can be built
dependencies = []

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


text:
function BurstSelector(arg0::jint, arg1::jdouble, arg2::jdouble, arg3::TimeScale)
    return BurstSelector((jint, jdouble, jdouble, TimeScale), arg0, arg1, arg2, arg3)
end

function select_dates(obj::BurstSelector, arg0::AbsoluteDate, arg1::AbsoluteDate)
    return jcall(obj, "selectDates", List, (AbsoluteDate, AbsoluteDate), arg0, arg1)
end


text:
module RunSimulation

using FFTW
using Random
using Serialization
using Parameters

using ..Models: AbstractModel
using ..SimulationParams: SimulationParameters
using ..GreensFunctions: EstimateGreensFunction, update!
using ..MuFinder: MuTuner, update_μ!
using ..FourierAcceleration: FourierAccelerator
using ..LangevinDynamics: evolve!, Dynamics, EulerDynamics, RungeKuttaDynamics, HeunsDynamics
using ..HMC: HybridMonteCarlo
import ..HMC
using ..SpecialUpdates: SpecialUpdate, NullUpdate, ReflectionUpdate, SwapUpdate, special_update!
using ..Measurements: initialize_measurements_container, initialize_measurement_folders!
using ..Measurements: make_measurements!, process_measurements!, write_measurements!, reset_measurements!

export run_simulation!

"""
Run Langevin simulation.
"""
function run_simulation!(model::AbstractModel, Gr::EstimateGreensFunction, μ_tuner::MuTuner, sim_params::SimulationParameters,
                         simulation_dynamics::Dynamics, burnin_dynamics::Dynamics,
                         sim_reflect_update::SpecialUpdate, burnin_reflect_update::SpecialUpdate,
                         sim_swap_update::SpecialUpdate, burnin_swap_update::SpecialUpdate,
                         fa::FourierAccelerator, container::NamedTuple, preconditioner, sim_stats::Dict,
                         burnin_start::Int=1, sim_start::Int=1,)::Dict

    ###############################################################
    ## PRE-ALLOCATING ARRAYS AND VARIABLES NEEDED FOR SIMULATION ##
    ###############################################################

    # previous epoch time
    t_prev = 0.0

    # current epoch time
    t_new = time()

    ##############################################
    ## RUNNING SIMULATION: THERMALIZATION STEPS ##
    ##############################################

    # frequency with which to update μ if tuning the denisty
    μ_update_freq = max(sim_params.meas_freq,1)

    # iterate over thermalization timesteps
    for t in burnin_start:sim_params.burnin

        # check if checkpoint need to be written
        t_new = time()
        if (t_new-t_prev) > sim_params.chckpnt_freq
            t_prev = t_new
            chkpnt = (model=model, μ_tuner=μ_tuner, container=container,
                      burnin_start=t, sim_start=1, sim_stats=sim_stats)
            sim_stats["write_time"] += @elapsed serialize(joinpath(sim_params.datafolder,"checkpoint.jls"),chkpnt)
        end

        # update phonon fields
        sim_stats["simulation_time"] += @elapsed sim_stats["iters"] += evolve!(model, burnin_dynamics, fa, preconditioner)

        # update chemical potential
        if μ_tuner.active && t%μ_update_freq==0
            sim_stats["simulation_time"] += @elapsed update!(Gr,model,preconditioner)
            sim_stats["simulation_time"] += @elapsed update_μ!(model, μ_tuner, Gr)
        end
    end

    ###########################################
    ## RUNNING SIMULATION: MEASUREMENT STEPS ##
    ###########################################

    # iterate over simulation time steps
    for t in sim_start:sim_params.nsteps

        # check if checkpoint needs to be written
        t_new = time()
        if (t_new-t_prev) > sim_params.chckpnt_freq
            t_prev = t_new
            chkpnt = (model=model, μ_tuner=μ_tuner, container=container,
                      burnin_start=sim_params.burnin+1, sim_start=t, sim_stats=sim_stats)
            sim_stats["write_time"] += @elapsed serialize(joinpath(sim_params.datafolder,"checkpoint.jls"),chkpnt)
        end

        # udpate phonon fields
        sim_stats["simulation_time"] += @elapsed sim_stats["iters"] += evolve!(model, simulation_dynamics, fa, preconditioner)

        # if time to make measurements
        if t%sim_params.meas_freq==0

            # make measurements
            nmeas = div(t,sim_params.meas_freq)
            sim_stats["measurement_time"] += @elapsed make_measurements!(container,model,Gr,nmeas,preconditioner)

            # update chemical potential
            if μ_tuner.active
                sim_stats["simulation_time"] += @elapsed update_μ!(model, μ_tuner, Gr)
            end

            # get measurement number
            nmeas = div(t,sim_params.meas_freq)

            # if bin of measurements full
            if nmeas%sim_params.bin_size==0

                # get bin number
                bin = div(nmeas,sim_params.bin_size)

                # process measurements
                sim_stats["measurement_time"] += @elapsed process_measurements!(container,sim_params,model)

                # write measurements to file
                sim_stats["write_time"] += @elapsed write_measurements!(container,model,bin)

                # reset measurements container
                sim_stats["measurement_time"] += @elapsed reset_measurements!(container,model)

                # write checkpoint
                t_new = time()
                t_prev = t_new
                chkpnt = (model=model, μ_tuner=μ_tuner, container=container,
                          burnin_start=sim_params.burnin+1, sim_start=t+1, sim_stats=sim_stats)
                sim_stats["write_time"] += @elapsed serialize(joinpath(sim_params.datafolder,"checkpoint.jls"),chkpnt)
            end
        end
    end

    # calculating the average number of iterations needed to solve linear system
    sim_stats["iters"] /= (sim_params.nsteps+sim_params.burnin)

    # report timings in units of minutes
    sim_stats["simulation_time"]  /= 60.0
    sim_stats["measurement_time"] /= 60.0
    sim_stats["write_time"]       /= 60.0
    sim_stats["acceptance_rate"]   = 1.0

    chkpnt = (model=model, μ_tuner=μ_tuner, burnin_start=sim_params.burnin+1,
              sim_start=sim_params.nsteps+1, sim_stats=sim_stats)
    sim_stats["write_time"] += @elapsed serialize(joinpath(sim_params.datafolder,"checkpoint.jls"),chkpnt)

    return sim_stats
end

"""
Run Hybrid Monte Carlo simulation.
"""
function run_simulation!(model::AbstractModel, Gr::EstimateGreensFunction, μ_tuner::MuTuner, sim_params::SimulationParameters,
                         simulation_hmc::HybridMonteCarlo, burnin_hmc::HybridMonteCarlo,
                         sim_reflect_update::SpecialUpdate, burnin_reflect_update::SpecialUpdate,
                         sim_swap_update::SpecialUpdate, burnin_swap_update::SpecialUpdate,
                         fa::FourierAccelerator, container::NamedTuple, preconditioner, sim_stats::Dict,
                         burnin_start::Int=1, sim_start::Int=1)::Dict

    ###############################################################
    ## PRE-ALLOCATING ARRAYS AND VARIABLES NEEDED FOR SIMULATION ##
    ###############################################################

    # previous epoch time
    t_prev = 0.0

    # current epoch time
    t_new = time()

    ##############################################
    ## RUNNING SIMULATION: THERMALIZATION STEPS ##
    ##############################################

    # thermalizing system
    for n in burnin_start:sim_params.burnin

        # check if checkpoint needs to be written
        t_new = time()
        if (t_new-t_prev) > sim_params.chckpnt_freq
            t_prev = t_new
            chkpnt = (model=model, μ_tuner=μ_tuner, container=container,
                      burnin_start=n, sim_start=1, sim_stats=sim_stats)
            sim_stats["write_time"] += @elapsed serialize(joinpath(sim_params.datafolder,"checkpoint.jls"),chkpnt)
        end

        # perform hmc update
        sim_stats["simulation_time"] += @elapsed accepted, niters = HMC.update!(model,burnin_hmc,fa,preconditioner)
        sim_stats["iters"]           += niters
        sim_stats["acceptance_rate"] += accepted

        # do reflection update
        if burnin_reflect_update.active && mod(n,burnin_reflect_update.freq)==0
            sim_stats["simulation_time"] += @elapsed accepted = special_update!(model,burnin_hmc,burnin_reflect_update,preconditioner)
            sim_stats["reflect_acceptance_rate"] += accepted
        end

        # do swap update
        if burnin_swap_update.active && mod(n,burnin_swap_update.freq)==0
            sim_stats["simulation_time"] += @elapsed accepted = special_update!(model,burnin_hmc,burnin_swap_update,preconditioner)
            sim_stats["swap_acceptance_rate"] += accepted
        end

        # update chemical potential
        if μ_tuner.active
            sim_stats["simulation_time"] += @elapsed update!(Gr,model,preconditioner)
            sim_stats["simulation_time"] += @elapsed update_μ!(model, μ_tuner, Gr)
        end
    end

    # close log file
    close(burnin_hmc.logfile)

    ###########################################
    ## RUNNING SIMULATION: MEASUREMENT STEPS ##
    ###########################################

    # iterate over hmc updates
    for n in sim_start:sim_params.nsteps

        # check if checkpoint needs to be written
        t_new = time()
        if (t_new-t_prev) > sim_params.chckpnt_freq
            t_prev = t_new
            chkpnt = (model=model, μ_tuner=μ_tuner, container=container,
                      burnin_start=sim_params.burnin+1, sim_start=n, sim_stats=sim_stats)
            sim_stats["write_time"] += @elapsed serialize(joinpath(sim_params.datafolder,"checkpoint.jls"),chkpnt)
        end

        # perform hmc update
        sim_stats["simulation_time"] += @elapsed accepted, niters = HMC.update!(model,simulation_hmc,fa,preconditioner)
        sim_stats["iters"]           += niters
        sim_stats["acceptance_rate"] += accepted

        # do reflection update
        if burnin_reflect_update.active && mod(n,sim_reflect_update.freq)==0
            sim_stats["simulation_time"] += @elapsed accepted = special_update!(model,burnin_hmc,sim_reflect_update,preconditioner)
            sim_stats["reflect_acceptance_rate"] += accepted
        end

        # do swap update
        if burnin_swap_update.active && mod(n,sim_swap_update.freq)==0
            sim_stats["simulation_time"] += @elapsed accepted = special_update!(model,burnin_hmc,sim_swap_update,preconditioner)
            sim_stats["swap_acceptance_rate"] += accepted
        end

        # if time to perform measurements (almost always yes)
        if n%sim_params.meas_freq==0

            # perform measurements
            nmeas = div(n,sim_params.meas_freq)
            sim_stats["measurement_time"] += @elapsed make_measurements!(container,model,Gr,nmeas,preconditioner)

            # update chemical potential
            if μ_tuner.active
                sim_stats["simulation_time"] += @elapsed update_μ!(model, μ_tuner, Gr)
            end

            # get measurement number
            nmeas = div(n,sim_params.meas_freq)

            # if bin of measurements is full
            if nmeas%sim_params.bin_size==0

                # calculate bin number
                bin = div(nmeas,sim_params.bin_size)

                # process measurements
                sim_stats["measurement_time"] += @elapsed process_measurements!(container,sim_params,model)

                # write measurements to file
                sim_stats["write_time"] += @elapsed write_measurements!(container,model,bin)

                # reset measurements container
                sim_stats["measurement_time"] += @elapsed reset_measurements!(container,model)

                # write checkpoint
                t_new = time()
                t_prev = t_new
                chkpnt = (model=model, μ_tuner=μ_tuner, container=container,
                          burnin_start=sim_params.burnin+1, sim_start=n+1, sim_stats=sim_stats)
                sim_stats["write_time"] += @elapsed serialize(joinpath(sim_params.datafolder,"checkpoint.jls"),chkpnt)
            end
        end
    end

    # calculating the average number of iterations needed to solve linear system
    sim_stats["iters"] /= (sim_params.nsteps+sim_params.burnin)

    # calculating the acceptance acceptance rate
    sim_stats["acceptance_rate"] /= (sim_params.nsteps + sim_params.burnin)

    # calculating special acceptance rate
    burnin_freq = burnin_reflect_update.freq
    sim_freq    = sim_reflect_update.freq
    sim_stats["reflect_acceptance_rate"] /= ( div(sim_params.nsteps,sim_freq) + div(sim_params.burnin,burnin_freq) )

    # calculating special acceptance rate
    burnin_freq = burnin_swap_update.freq
    sim_freq    = sim_swap_update.freq
    sim_stats["swap_acceptance_rate"] /= ( div(sim_params.nsteps,sim_freq) + div(sim_params.burnin,burnin_freq) )

    # report timings in units of minutes
    sim_stats["simulation_time"]  /= 60.0
    sim_stats["measurement_time"] /= 60.0
    sim_stats["write_time"]       /= 60.0

    # final checkpoint
    chkpnt = (model=model, μ_tuner=μ_tuner, burnin_start=sim_params.burnin+1,
              sim_start=sim_params.nsteps+1, sim_stats=sim_stats)
    serialize(joinpath(sim_params.datafolder,"checkpoint.jls"),chkpnt)

    # close log files
    close(simulation_hmc.logfile)

    return sim_stats
end

end

max_stars_repo_name: cohensbw/Langevin

max_issues_repo_name: cohensbw/ElPhDynamics

max_forks_repo_name: cohensbw/ElPhDynamics


text:
function pack{T1<:Number,T2<:Integer}(X::SparseMatrixCSC{T1,T2}; symmetric::Bool=false,
                                                                 hermitian::Bool=false,
                                                                 posdef::Bool=false)
    # Packs a sparse array
    _nnz = nnz(X)
    dims = (X.m, X.n)
    rowind = Vector{UInt32}(_nnz)
    colind = Vector{UInt32}(_nnz)
    nzind = Vector{UInt32}(_nnz)
    rowind_wblock = Vector{UInt32}(128)
    colind_wblock = Vector{UInt32}(128)
    nzind_wblock = Vector{UInt32}(128)
    unique_nzval = unique(X.nzval)
    # Fill row and column index
    ind = 1
    for colptr_ind = 1:length(X.colptr)-1
        @inbounds for i = X.colptr[colptr_ind]:X.colptr[colptr_ind+1]-1
            rowind[ind] = X.rowval[ind]
            colind[ind] = colptr_ind
            ind += 1
        end
    end
    # Fill non-zero index
    unique_nzvals_dict = Dict{Float64,UInt32}(zip(unique_nzval, 1:length(unique_nzval)))
    for i = 1:_nnz
        @inbounds nzind[i] = unique_nzvals_dict[X.nzval[i]]
    end
    # Compress integers
    rowind = pack(rowind)
    colind = pack(colind)
    nzind = pack(nzind)
    return TinySparseMat{T1}(f, fc, X.m, X.n, _nnz, dims, rowind, colind, nzind, unique_nzval, rowind_wblock, colind_wblock, nzind_wblock, symmetric, hermitian, posdef)
end

# Multiplication with vector
function f{T<:Number}(y::Vector{T}, A::TinySparseMat{T}, x::Vector{T})
    fill!(y, 0.0)
    nchunks = cld(A._nnz, 128)
    for chunk = 1:nchunks - 1
        unpack!(A.rowind, A.rowind_wblock, chunk)
        unpack!(A.colind, A.colind_wblock, chunk)
        unpack!(A.nzind, A.nzind_wblock, chunk)
        for i::Int64 = 1:128
            @inbounds y[A.rowind_wblock[i]] += x[A.colind_wblock[i]]*A.nzval[A.nzind_wblock[i]]
        end
    end
    unpack!(A.rowind, A.rowind_wblock, nchunks)
    unpack!(A.colind, A.colind_wblock, nchunks)
    unpack!(A.nzind, A.nzind_wblock, nchunks)
    for i = 1:A._nnz-(nchunks-1)*128
        @inbounds y[A.rowind_wblock[i]] += x[A.colind_wblock[i]]*A.nzval[A.nzind_wblock[i]]
    end
end

# Conjugate transpose multiplication with vector
function fc{T<:Number}(y::Vector{T}, A::TinySparseMat{T}, x::Vector{T})
    fill!(y, 0.0)
    nchunks = cld(A._nnz, 128)
    for chunk = 1:nchunks - 1
        unpack!(A.rowind, A.rowind_wblock, chunk)
        unpack!(A.colind, A.colind_wblock, chunk)
        unpack!(A.nzind, A.nzind_wblock, chunk)
        for i = 1:128
            @inbounds y[A.colind_wblock[i]] += x[A.rowind_wblock[i]]*A.nzval[A.nzind_wblock[i]]'
        end
    end
    unpack!(A.rowind, A.rowind_wblock, nchunks)
    unpack!(A.colind, A.colind_wblock, nchunks)
    unpack!(A.nzind, A.nzind_wblock, nchunks)
    for i = 1:A._nnz-(nchunks-1)*128
        @inbounds y[A.colind_wblock[i]] += x[A.rowind_wblock[i]]*A.nzval[A.nzind_wblock[i]]'
    end
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule YASM_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("YASM")
JLLWrappers.@generate_main_file("YASM", UUID("997772c2-56d0-5ccd-9329-3f55f14e5768"))
end  # module YASM_jll


text:
function call(::Types.Who, args::Vector, storage, conn::IO)
    ip_addr = string(first(getpeername(conn)))

    user = storage.active_clients[ip_addr]
    Network.fancy_write(storage, conn, "$(user.name) $(user.ip_addr)\n")
end


text:
################################################################################
#  Copyright 2020, Tom Van Acker                                               #
################################################################################
# MultiStateSystems.jl                                                         #
# A Julia package to solve multi-state system models.                          #
# See http://github.com/timmyfaraday/MultiStateSystems.jl                      #
################################################################################

# load pkgs
using Unitful
using MultiStateSystems

# initialize the state-transition diagram corresponding to the front-end ac/dc
# converter
stdᵃᶜ = STD()

# add the states to the std
add_states!(stdᵃᶜ, name  = ["available", "unavailable"],
                   power = [1.0u"MW", 0.0u"MW"],
                   init  = [1.0, 0.0])

# add the transitions to the std
add_transitions!(stdᵃᶜ, rate = [0.0u"1/hr"  1.14e-4u"1/hr"
                                12.0u"1/hr" 0.0u"1/hr"])

# return the std
return stdᵃᶜ


text:
# TODO: maybe we should get rid of the LearningRate abstraction, make it a number,
#   and then allow sub-learner(s) in the GradientLearner to update the learning rate?

immutable FixedLR <: LearningRate
    lr::Float64
end
value(lr::FixedLR) = lr.lr
update!(lr::FixedLR, err) = lr

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

# "Adapts learning rate based on relative variance of the changes in the test error"
# @with_kw type AdaptiveLR <: LearningRate
#     lr::Float64 = 1e-2
#     ε::Diff{Float64} = Diff()
#     lookback::Int = 20
#     σ²::Variance{BoundedEqualWeight} = Variance(BoundedEqualWeight(lookback))
#     adjpct::Float64 = 1e-2
#     cutoff::Float64 = 1e-1
# end
#
#
# # if the error is decreasing at a large rate relative to the variance, increase the learning rate (speed it up)
# function update!(lr::AdaptiveLR, err)
#     fit!(lr.ε, err)
#     fit!(lr.σ², diff(lr.ε))
#     μ = mean(lr.σ²)
#     σ = std(lr.σ²)
#     if σ > 0
#         pct = lr.adjpct * (μ / σ < -lr.cutoff ? 1.0 : -1.0)
#         lr.lr *= (1.0 + pct)
#     end
#     lr
# end

max_stars_repo_name: JuliaML/StochasticOptimization

max_issues_repo_name: JuliaML/StochasticOptimization

max_forks_repo_name: JuliaML/StochasticOptimization


text:
using DiffEqSensitivity, OrdinaryDiffEq, Zygote
using RecursiveArrayTools: DiffEqArray
using Test
import Tracker

function fiip(du,u,p,t)
  du[1] = dx = p[1]*u[1] - p[2]*u[1]*u[2]
  du[2] = dy = -p[3]*u[2] + p[4]*u[1]*u[2]
end
function foop(u,p,t)
  dx = p[1]*u[1] - p[2]*u[1]*u[2]
  dy = -p[3]*u[2] + p[4]*u[1]*u[2]
  [dx,dy]
end
function foop(u::Tracker.TrackedArray,p,t)
  dx = p[1]*u[1] - p[2]*u[1]*u[2]
  dy = -p[3]*u[2] + p[4]*u[1]*u[2]
  Tracker.collect([dx,dy])
end

p = [1.5,1.0,3.0,1.0]; u0 = [1.0;1.0]
prob = ODEProblem(fiip,u0,(0.0,10.0),p)
proboop = ODEProblem(foop,u0,(0.0,10.0),p)

sol = concrete_solve(prob,Tsit5(),abstol=1e-14,reltol=1e-14)
@test sol isa DiffEqArray
sumsol = sum(sol)
@test sum(concrete_solve(prob,Tsit5(),u0,p,abstol=1e-14,reltol=1e-14)) == sumsol
@test sum(concrete_solve(prob,Tsit5(),u0,p,abstol=1e-14,reltol=1e-14,sensealg=ForwardDiffSensitivity())) == sumsol
@test sum(concrete_solve(prob,Tsit5(),u0,p,abstol=1e-14,reltol=1e-14,sensealg=BacksolveAdjoint())) == sumsol

###
### adjoint
###

_sol = solve(prob,Tsit5(),abstol=1e-14,reltol=1e-14)
ū0,adj = adjoint_sensitivities(_sol,Tsit5(),((out,u,p,t,i) -> out .= -1),0.0:0.1:10,abstol=1e-14,
                                  reltol=1e-14,iabstol=1e-14,ireltol=1e-12)
du01,dp1 = Zygote.gradient((u0,p)->sum(concrete_solve(prob,Tsit5(),u0,p,abstol=1e-14,reltol=1e-14,saveat=0.1,sensealg=QuadratureAdjoint())),u0,p)
du02,dp2 = Zygote.gradient((u0,p)->sum(concrete_solve(prob,Tsit5(),u0,p,abstol=1e-14,reltol=1e-14,saveat=0.1,sensealg=InterpolatingAdjoint())),u0,p)
du03,dp3 = Zygote.gradient((u0,p)->sum(concrete_solve(prob,Tsit5(),u0,p,abstol=1e-14,reltol=1e-14,saveat=0.1,sensealg=BacksolveAdjoint())),u0,p)
du04,dp4 = Zygote.gradient((u0,p)->sum(concrete_solve(prob,Tsit5(),u0,p,abstol=1e-14,reltol=1e-14,saveat=0.1,sensealg=TrackerAdjoint())),u0,p)
@test_broken Zygote.gradient((u0,p)->sum(concrete_solve(prob,Tsit5(),u0,p,abstol=1e-14,reltol=1e-14,saveat=0.1,sensealg=ZygoteAdjoint())),u0,p) isa Tuple
du05,dp5 = Zygote.gradient((u0,p)->sum(concrete_solve(prob,Tsit5(),u0,p,abstol=1e-14,reltol=1e-14,saveat=0.1,save_idxs=1:1,sensealg=InterpolatingAdjoint())),u0,p)
du06,dp6 = Zygote.gradient((u0,p)->sum(concrete_solve(prob,Tsit5(),u0,p,abstol=1e-14,reltol=1e-14,saveat=0.0:0.1:10.0,save_idxs=1:1,sensealg=QuadratureAdjoint())),u0,p)
du07,dp7 = Zygote.gradient((u0,p)->sum(concrete_solve(prob,Tsit5(),u0,p,abstol=1e-14,reltol=1e-14,saveat=0.1,save_idxs=1,sensealg=InterpolatingAdjoint())),u0,p)

@test ū0 ≈ du01 rtol=1e-12
@test ū0 == du02
@test ū0 ≈ du03 rtol=1e-12
@test ū0 ≈ du04 rtol=1e-12
@test ū0 ≈ du05 rtol=1e-12
@test ū0 ≈ du06 rtol=1e-12
@test ū0 ≈ du07 rtol=1e-12
@test adj ≈ dp1' rtol=1e-12
@test adj == dp2'
@test adj ≈ dp3' rtol=1e-12
@test adj ≈ dp4' rtol=1e-12
@test adj ≈ dp5' rtol=1e-12
@test adj ≈ dp6' rtol=1e-12
@test adj ≈ dp7' rtol=1e-12

###
### forward
###

du06,dp6 = Zygote.gradient((u0,p)->sum(concrete_solve(prob,Tsit5(),u0,p,abstol=1e-14,reltol=1e-14,saveat=0.1,sensealg=ForwardSensitivity())),u0,p)
du07,dp7 = Zygote.gradient((u0,p)->sum(concrete_solve(prob,Tsit5(),u0,p,abstol=1e-14,reltol=1e-14,saveat=0.1,sensealg=ForwardDiffSensitivity())),u0,p)
@test_broken du08,dp8 = Zygote.gradient((u0,p)->sum(concrete_solve(prob,Tsit5(),u0,p,abstol=1e-14,reltol=1e-14,saveat=0.1,save_idxs = 1:1,sensealg=ForwardSensitivity())),u0,p)
@test_broken du09,dp9 = Zygote.gradient((u0,p)->sum(concrete_solve(prob,Tsit5(),u0,p,abstol=1e-14,reltol=1e-14,saveat=0.1,save_idxs = 1:1,sensealg=ForwardDiffSensitivity())),u0,p)

@test du06 === du07 === nothing
@test adj ≈ dp6' rtol=1e-12
@test adj ≈ dp7' rtol=1e-12

###
### OOPs
###

###
### adjoint
###

du01,dp1 = Zygote.gradient((u0,p)->sum(concrete_solve(proboop,Tsit5(),u0,p,abstol=1e-14,reltol=1e-14,saveat=0.1,sensealg=QuadratureAdjoint())),u0,p)
du02,dp2 = Zygote.gradient((u0,p)->sum(concrete_solve(proboop,Tsit5(),u0,p,abstol=1e-14,reltol=1e-14,saveat=0.1,sensealg=InterpolatingAdjoint())),u0,p)
du03,dp3 = Zygote.gradient((u0,p)->sum(concrete_solve(proboop,Tsit5(),u0,p,abstol=1e-14,reltol=1e-14,saveat=0.1,sensealg=BacksolveAdjoint())),u0,p)
du04,dp4 = Zygote.gradient((u0,p)->sum(concrete_solve(proboop,Tsit5(),u0,p,abstol=1e-14,reltol=1e-14,saveat=0.1,sensealg=TrackerAdjoint())),u0,p)
@test_broken Zygote.gradient((u0,p)->sum(concrete_solve(proboop,Tsit5(),u0,p,abstol=1e-14,reltol=1e-14,saveat=0.1,sensealg=ZygoteAdjoint())),u0,p) isa Tuple

@test ū0 ≈ du01 rtol=1e-12
@test ū0 ≈ du02 rtol=1e-12
@test ū0 ≈ du03 rtol=1e-12
@test ū0 ≈ du04 rtol=1e-12
@test adj ≈ dp1' rtol=1e-12
@test adj ≈ dp2' rtol=1e-12
@test adj ≈ dp3' rtol=1e-12
@test adj ≈ dp4' rtol=1e-12

###
### forward
###

@test_broken Zygote.gradient((u0,p)->sum(concrete_solve(proboop,Tsit5(),u0,p,abstol=1e-14,reltol=1e-14,saveat=0.1,sensealg=ForwardSensitivity())),u0,p) isa Tuple
du07,dp7 = Zygote.gradient((u0,p)->sum(concrete_solve(proboop,Tsit5(),u0,p,abstol=1e-14,reltol=1e-14,saveat=0.1,sensealg=ForwardDiffSensitivity())),u0,p)

@test du06 === du07 === nothing
@test adj ≈ dp6' rtol=1e-12
@test adj ≈ dp7' rtol=1e-12


text:
using SBMLBioModelsRepository
using Pkg, Test
using SBML, SBMLToolkit
using ModelingToolkit, OrdinaryDiffEq, CSV, DataFrames, BenchmarkTools, Sundials
using Base.Threads, Glob, Dates

mkpath(logdir)

@testset "SBMLBioModelsRepository.jl" begin
    @testset "test_suite" begin include("test_suite.jl") end
    @testset "biomd" begin include("biomd.jl") end
end


text:
using LinearAlgebra

function vandermonde_matrix(x::Vector{T}) where {T}
    local n = length(x)
	local V = Array{T,2}(undef, n, n)

	for i in 1:n
		V[:,i] .= x.^(i-1)
	end

	V
end

function vandermonde_matrix_inverse(x::Vector{T}) where {T}
    local n = length(x)

    local L::Matrix{T} = zeros(n,n)
    local U::Matrix{T} = Matrix{T}(I, n, n)
    local V::Matrix{T}

    L[1,1] = 1
    for i in 2:n
        for j in 1:i
            p = 1
            for k in 1:i
                if k ≠ j
                    p *= (x[j] - x[k])
                end
            end
            L[i,j] = 1/p
        end
    end

    i = 1
    for j in i+1:n
        U[i,j] = - U[i,j-1] * x[j-1]
    end

    for i in 2:n
        for j in i+1:n
            U[i,j] = U[i-1,j-1] - U[i,j-1] * x[j-1]
        end
    end

    V = *(U,L)
end


text:
import Base.rand

"""
    rand(N::ProbabilisticNetwork)

Converts a probabilistic network into a deterministic one by performing random
draws. All interactions are treated as independent Bernoulli events. Note that
this network is *not* check for degeneracy, *i.e.* species can end up with no
interactions.

#### References

Poisot, T., Cirtwill, A.R., Cazelles, K., Gravel, D., Fortin, M.-J., Stouffer, D.B., 2016. The structure of probabilistic networks. Methods in Ecology and Evolution 7, 303–312. https://doi.org/10.1111/2041-210X.12468
"""
function rand(N::ProbabilisticNetwork)
    # Get the correct network type
    newtype = typeof(N) <: AbstractUnipartiteNetwork ? UnipartiteNetwork : BipartiteNetwork
    return newtype(rand(Float64, size(N)).<=N.A, species_objects(N)...)
end

"""
    rand(N::ProbabilisticNetwork, n::T) where {T<:Integer}

Generates a number of random deterministic networks based on a probabilistic
network.

#### References

Poisot, T., Cirtwill, A.R., Cazelles, K., Gravel, D., Fortin, M.-J., Stouffer, D.B., 2016. The structure of probabilistic networks. Methods in Ecology and Evolution 7, 303–312. https://doi.org/10.1111/2041-210X.12468
"""
function rand(N::ProbabilisticNetwork, n::T) where {T<:Integer}
    @assert n > 0
    return map(x -> rand(N), 1:n)
end

"""
    rand(N::ProbabilisticNetwork, S::Tuple{Int64,Int64})

Generates a number of random deterministic networks based on a probabilistic
network, and returns them as a matrix.

#### References

Poisot, T., Cirtwill, A.R., Cazelles, K., Gravel, D., Fortin, M.-J., Stouffer, D.B., 2016. The structure of probabilistic networks. Methods in Ecology and Evolution 7, 303–312. https://doi.org/10.1111/2041-210X.12468
"""
function rand(N::ProbabilisticNetwork, S::Tuple{Int64,Int64})
    @assert minimum(S) > 0
    return reshape(map(x -> rand(N), 1:prod(S)), S)
end


text:
name: FluxBot

on:
  issue_comment:
    types: [created, edited]

jobs:
  build:
    if: contains(github.event.comment.body, '@ModelZookeeper')
    runs-on: ${{ matrix.os }}
    strategy:
      matrix:
        julia-version: [1.5.0]
        julia-arch: [x86]
        os: [ubuntu-latest]
    steps:
      - uses: actions/checkout@af513c7a016048ae468971c52ed77d9562c7c819 # v1.0.0
      - uses: julia-actions/setup-julia@v1
        with:
          version: ${{ matrix.julia-version }}
      - name: Install dependencies
        run: julia --project=.fluxbot/ -e 'using Pkg; Pkg.instantiate()'
      - name: FluxBot.respond
        env:
          FLUXBOT_GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
          BOT_SECRET: ${{ secrets.BOT_SECRET }}
          MODELZOO_TRIGGER_TOKEN: ${{ secrets.MODELZOO_TRIGGER_TOKEN }}
        run: julia --project=.fluxbot -e 'using FluxBot; FluxBot.trial()'

max_stars_repo_name: shawn-kern/model-zoo

max_issues_repo_name: shawn-kern/model-zoo

max_forks_repo_name: shawn-kern/model-zoo


text:
using Test

using MathOptInterface
const MOI = MathOptInterface
const MOIT = MOI.Test
const MOIU = MOI.Utilities
const MOIB = MOI.Bridges

import ECOS
const optimizer = ECOS.Optimizer()
MOI.set(optimizer, MOI.Silent(), true)

@testset "SolverName" begin
    @test MOI.get(optimizer, MOI.SolverName()) == "ECOS"
end

@testset "supports_allocate_load" begin
    @test MOIU.supports_allocate_load(optimizer, false)
    @test !MOIU.supports_allocate_load(optimizer, true)
end

# UniversalFallback is needed for starting values, even if they are ignored by ECOS
const cache = MOIU.UniversalFallback(MOIU.Model{Float64}())
const cached = MOIU.CachingOptimizer(cache, optimizer)

const bridged = MOIB.full_bridge_optimizer(cached, Float64)

# SOC2 requires 1e-4
const config = MOIT.TestConfig(atol=1e-4, rtol=1e-4)

@testset "Unit" begin
    MOIT.unittest(bridged,
                  config,
                  [
        # `NumberOfThreads` not supported.
        "number_threads",
        # `TimeLimitSec` not supported.
        "time_limit_sec",
        # Need https://github.com/JuliaOpt/MathOptInterface.jl/issues/529
        "solve_qp_edge_cases",
        # Integer and ZeroOne sets are not supported
        "solve_integer_edge_cases", "solve_objbound_edge_cases",
        "solve_zero_one_with_bounds_1",
        "solve_zero_one_with_bounds_2",
        "solve_zero_one_with_bounds_3"])
end

@testset "Continuous linear problems" begin
    MOIT.contlineartest(bridged, config)
end

@testset "Continuous quadratic problems" begin
    MOIT.qcptest(bridged, config)
end

@testset "Continuous conic problems" begin
    exclude = ["dualexp", "pow", "dualpow", "sdp", "rootdet", "logdet"]
    MOIT.contconictest(bridged, config, exclude)
end


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

name = "BWA"
version = v"0.7.17"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/lh3/bwa/releases/download/v$(version)/bwa-$(version).tar.bz2",
                  "de1b4d4e745c0b7fc3e107b5155a51ac063011d33a5d82696331ecf4bed8d0fd"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/bwa*/
atomic_patch -p1 ../patches/makefile.patch
atomic_patch -p1 ../patches/0001-Fix-building-against-GCC-10.patch
make -j${nproc}
mkdir -p "${libdir}" "${bindir}"
cp "bwa${exeext}" "${bindir}/bwa${exeext}"
cp "libbwa.${dlext}" "${libdir}/libbwa.${dlext}"
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line.
platforms = supported_platforms(; experimental=true)
# The package uses Intel intrinsics, so we can build only for Intel platforms.
filter!(p -> BinaryBuilder.proc_family(p) == "intel", platforms)
# ...and uses some Unix header files not available on Windows.
filter!(!Sys.iswindows, platforms)

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

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("Zlib_jll"),
]

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

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
using Statistics
using StatsBase:sample, sample!
using EvoTrees
using BenchmarkTools
using CUDA

# prepare a dataset
features = rand(Int(1.25e6), 100)
# features = rand(100, 10)
X = features
Y = rand(size(X, 1))
𝑖 = collect(1:size(X, 1))

# train-eval split
𝑖_sample = sample(𝑖, size(𝑖, 1), replace=false)
train_size = 0.8
𝑖_train = 𝑖_sample[1:floor(Int, train_size * size(𝑖, 1))]
𝑖_eval = 𝑖_sample[floor(Int, train_size * size(𝑖, 1)) + 1:end]

X_train, X_eval = X[𝑖_train, :], X[𝑖_eval, :]
Y_train, Y_eval = Y[𝑖_train], Y[𝑖_eval]

###########################
# Tree CPU
###########################
params_c = EvoTreeRegressor(T=Float32,
    loss=:linear, metric=:none,
    nrounds=100,
    λ=1.0, γ=0.1, η=0.1,
    max_depth=6, min_weight=1.0,
    rowsample=0.5, colsample=0.5, nbins=64);

# params_c = EvoTreeGaussian(T=Float32,
#     loss=:gaussian, metric=:none,
#     nrounds=100,
#     λ=1.0, γ=0.1, η=0.1,
#     max_depth=6, min_weight=1.0,
#     rowsample=0.5, colsample=0.5, nbins=64);

model_c, cache_c = EvoTrees.init_evotree(params_c, X_train, Y_train);

# initialize from cache
params_c = model_c.params
X_size = size(cache_c.X_bin)

# select random rows and cols
sample!(params_c.rng, cache_c.𝑖_, cache_c.nodes[1].𝑖, replace=false, ordered=true);
sample!(params_c.rng, cache_c.𝑗_, cache_c.𝑗, replace=false, ordered=true);
# @btime sample!(params_c.rng, cache_c.𝑖_, cache_c.nodes[1].𝑖, replace=false, ordered=true);
# @btime sample!(params_c.rng, cache_c.𝑗_, cache_c.𝑗, replace=false, ordered=true);

𝑖 = cache_c.nodes[1].𝑖
𝑗 = cache_c.𝑗

# build a new tree
# 897.800 μs (6 allocations: 736 bytes)
@time EvoTrees.update_grads!(params_c.loss, cache_c.δ𝑤, cache_c.pred, cache_c.Y, params_c.α)
# @btime EvoTrees.update_grads!($params_c.loss, $cache_c.δ𝑤, $cache_c.pred_cpu, $cache_c.Y_cpu, $params_c.α)
# ∑ = vec(sum(cache_c.δ[𝑖,:], dims=1))
# gain = EvoTrees.get_gain(params_c.loss, ∑, params_c.λ)
# assign a root and grow tree
# train_nodes[1] = EvoTrees.TrainNode(UInt32(0), UInt32(1), ∑, gain)

# 62.530 ms (7229 allocations: 17.43 MiB)
# 1.25e5: 9.187 ms (7358 allocations: 2.46 MiB)
tree = EvoTrees.Tree(params_c.max_depth, model_c.K, zero(typeof(params_c.λ)))
@time EvoTrees.grow_tree!(tree, cache_c.nodes, params_c, cache_c.δ𝑤, cache_c.edges, cache_c.𝑗, cache_c.left, cache_c.left, cache_c.right, cache_c.X_bin, cache_c.K)
@btime EvoTrees.grow_tree!($EvoTrees.Tree(params_c.max_depth, model_c.K, zero(typeof(params_c.λ))), $cache_c.nodes, $params_c, $cache_c.δ𝑤, $cache_c.edges, $cache_c.𝑗, $cache_c.left, $cache_c.left, $cache_c.right, $cache_c.X_bin, $cache_c.K)

@time EvoTrees.grow_tree!(EvoTrees.Tree(params_c.max_depth, model_c.K, params_c.λ), params_c, cache_c.δ, cache_c.hist, cache_c.histL, cache_c.histR, cache_c.gains, cache_c.edges, 𝑖, 𝑗, 𝑛, cache_c.X_bin);
@btime EvoTrees.grow_tree!(EvoTrees.Tree($params_c.max_depth, $model_c.K, $params_c.λ), $params_c, $cache_c.δ, $cache_c.hist, $cache_c.histL, $cache_c.histR, $cache_c.gains, $cache_c.edges, $𝑖, $𝑗, $𝑛, $cache_c.X_bin);
@code_warntype EvoTrees.grow_tree!(EvoTrees.Tree(params_c.max_depth, model_c.K, params_c.λ), params_c, cache_c.δ, cache_c.hist, cache_c.histL, cache_c.histR, cache_c.gains, cache_c.edges, 𝑖, 𝑗, 𝑛, cache_c.X_bin);

# push!(model_c.trees, tree)
# 1.883 ms (83 allocations: 13.77 KiB)
@btime EvoTrees.predict!(model_c.params.loss, cache_c.pred_cpu, tree, cache_c.X, model_c.K)

δ𝑤, K, edges, X_bin, nodes, out, left, right = cache_c.δ𝑤, cache_c.K, cache_c.edges, cache_c.X_bin, cache_c.nodes, cache_c.out, cache_c.left, cache_c.right;

# 9.613 ms (81 allocations: 13.55 KiB)
# 1.25e5: 899.200 μs (81 allocations: 8.22 KiB)
@time EvoTrees.update_hist!(params_c.loss, nodes[1].h, δ𝑤, X_bin, 𝑖, 𝑗, K)
@btime EvoTrees.update_hist!($params_c.loss, $nodes[1].h, $δ𝑤, $X_bin, $𝑖, $𝑗, $K)
@btime EvoTrees.update_hist!($nodes[1].h, $δ𝑤, $X_bin, $nodes[1].𝑖, $𝑗)
@code_warntype EvoTrees.update_hist!(hist, δ, X_bin, 𝑖, 𝑗, 𝑛)

j = 1
# 8.399 μs (80 allocations: 13.42 KiB)
n = 1
nodes[1].∑ .= vec(sum(δ𝑤[:, 𝑖], dims=2))
EvoTrees.update_gains!(params_c.loss, nodes[n], 𝑗, params_c, K)
nodes[1].gains
# findmax(nodes[1].gains) #1.25e5: 36.500 μs (81 allocations: 8.22 KiB)
@btime EvoTrees.update_gains!($params_c.loss, $nodes[n], $𝑗, $params_c, $K)
@code_warntype EvoTrees.update_gains!(params_c.loss, nodes[n], 𝑗, params_c, K)

#1.25e5: 14.100 μs (1 allocation: 32 bytes)
best = findmax(nodes[n].gains)
@btime best = findmax(nodes[n].gains)
@btime best = findmax(view(nodes[n].gains, :, 𝑗))

tree.cond_bin[n] = best[2][1]
tree.feat[n] = best[2][2]

Int.(tree.cond_bin[n])
# tree.cond_bin[n] = 32

# 204.900 μs (1 allocation: 96 bytes)
offset = 0
@time EvoTrees.split_set!(left, right, 𝑖, X_bin, tree.feat[n], tree.cond_bin[n], offset)
@btime EvoTrees.split_set!($left, $right, $𝑖, $X_bin, $tree.feat[n], $tree.cond_bin[n], $offset)
@code_warntype EvoTrees.split_set!(left, right, 𝑖, X_bin, tree.feat[n], tree.cond_bin[n])

# 1.25e5: 227.200 μs (22 allocations: 1.44 KiB)
@time EvoTrees.split_set_threads!(out, left, right, 𝑖, X_bin, tree.feat[n], tree.cond_bin[n], offset)
@btime EvoTrees.split_set_threads!($out, $left, $right, $𝑖, $X_bin, $tree.feat[n], $tree.cond_bin[n], $offset, Int(2e15))

###################################################
# GPU
###################################################
params_g = EvoTreeRegressor(T=Float32,
    loss=:linear, metric=:none,
    nrounds=100,
    λ=1.0, γ=0.1, η=0.1,
    max_depth=6, min_weight=1.0,
    rowsample=0.5, colsample=0.5, nbins=64);

model_g, cache_g = EvoTrees.init_evotree_gpu(params_g, X_train, Y_train);

params_g = model_g.params;
X_size = size(cache_g.X_bin);

# select random rows and cols
𝑖c = cache_g.𝑖_[sample(params_g.rng, cache_g.𝑖_, ceil(Int, params_g.rowsample * X_size[1]), replace=false, ordered=true)]
𝑖 = CuVector(𝑖c)
𝑗c = cache_g.𝑗_[sample(params_g.rng, cache_g.𝑗_, ceil(Int, params_g.colsample * X_size[2]), replace=false, ordered=true)]
𝑗 = CuVector(𝑗c)
sum(𝑖)
cache_g.nodes[1].𝑖 = 𝑖
cache_g.nodes[1].𝑖 .= 𝑖
sum(𝑗)
cache_g.𝑗 .= 𝑗
sum(cache_g.𝑗)

# build a new tree
# 144.600 μs (23 allocations: 896 bytes) - 5-6 X time faster on GPU
@time CUDA.@sync EvoTrees.update_grads_gpu!(params_g.loss, cache_g.δ𝑤, cache_g.pred, cache_g.Y)
# @btime CUDA.@sync EvoTrees.update_grads_gpu!($params_g.loss, $cache_g.δ𝑤, $cache_g.pred_gpu, $cache_g.Y_gpu)
# sum Gradients of each of the K parameters and bring to CPU

# 33.447 ms (6813 allocations: 307.27 KiB)
tree = EvoTrees.TreeGPU(params_g.max_depth, model_g.K, params_g.λ)
sum(cache_g.δ𝑤[:, cache_g.nodes[1].𝑖], dims=2)
CUDA.@time EvoTrees.grow_tree_gpu!(tree, cache_g.nodes, params_g, cache_g.δ𝑤, cache_g.edges, CuVector(cache_g.𝑗), cache_g.out, cache_g.left, cache_g.right, cache_g.X_bin, cache_g.K)
@btime EvoTrees.grow_tree_gpu!(EvoTrees.TreeGPU(params_g.max_depth, model_g.K, params_g.λ), cache_g.nodes, params_g, $cache_g.δ𝑤, $cache_g.edges, $𝑗, $cache_g.out, $cache_g.left, $cache_g.right, $cache_g.X_bin, $cache_g.K);
@code_warntype EvoTrees.grow_tree_gpu!(EvoTrees.TreeGPU(params_g.max_depth, model_g.K, params_g.λ), params_g, cache_g.δ, cache_g.hist, cache_g.histL, cache_g.histR, cache_g.gains, cache_g.edges, 𝑖, 𝑗, 𝑛, cache_g.X_bin);

push!(model_g.trees, tree);
# 2.736 ms (93 allocations: 13.98 KiB)
@time CUDA.@sync EvoTrees.predict!(cache_g.pred_gpu, tree, cache_g.X_bin)
@btime CUDA.@sync EvoTrees.predict!($cache_g.pred_gpu, $tree, $cache_g.X_bin)

###########################
# Tree GPU
###########################
δ𝑤, K, edges, X_bin, nodes, out, left, right = cache_g.δ𝑤, cache_g.K, cache_g.edges, cache_g.X_bin, cache_g.nodes, cache_g.out, cache_g.left, cache_g.right;

# 2.571 ms (1408 allocations: 22.11 KiB)
# 𝑗2 = CuArray(sample(UInt32.(1:100), 50, replace=false, ordered=true))
@time EvoTrees.update_hist_gpu!(params_g.loss, nodes[1].h, δ𝑤, X_bin, 𝑖, 𝑗, K)
@btime EvoTrees.update_hist_gpu!($params_g.loss, $nodes[1].h, $δ𝑤, $X_bin, $𝑖, $𝑗, $K)
@btime EvoTrees.update_hist_gpu!($nodes[1].h, $δ𝑤, $X_bin, $nodes[1].𝑖, $𝑗)
@code_warntype EvoTrees.update_hist_gpu!(hist, δ, X_bin, 𝑖, 𝑗, 𝑛)

# 72.100 μs (186 allocations: 6.00 KiB)
n = 1
nodes[1].∑ .= vec(sum(δ𝑤[:, 𝑖], dims=2))
CUDA.@time EvoTrees.update_gains_gpu!(params_g.loss, nodes[n], 𝑗, params_g, K)
@btime EvoTrees.update_gains_gpu!($params_g.loss, $nodes[n], $𝑗, $params_g, $K)

tree = EvoTrees.TreeGPU(params_g.max_depth, model_g.K, params_g.λ)
best = findmax(nodes[n].gains)
if best[2][1] != params_g.nbins && best[1] > nodes[n].gain + params_g.γ
    tree.gain[n] = best[1]
    tree.cond_bin[n] = best[2][1]
    tree.feat[n] = best[2][2]
    tree.cond_float[n] = edges[tree.feat[n]][tree.cond_bin[n]]
end

tree.split[n] = tree.cond_bin[n] != 0
tree.feat[n]
Int(tree.cond_bin[n])

# 673.900 μs (600 allocations: 29.39 KiB)
offset = 0
_left, _right = EvoTrees.split_set_threads_gpu!(out, left, right, 𝑖, X_bin, tree.feat[n], tree.cond_bin[n], offset)
@time EvoTrees.split_set_threads_gpu!(out, left, right, 𝑖, X_bin, tree.feat[n], tree.cond_bin[n], offset)
@btime EvoTrees.split_set_threads_gpu!($out, $left, $right, $𝑖, $X_bin, $tree.feat[n], $tree.cond_bin[n], $offset)


text:
#!/usr/bin/env julia

using Luxor

using Test

using Random
Random.seed!(42)

function testbezierstroke(fname)
    Random.seed!(3)
    currentwidth = 1200
    currentheight = 1200
    Drawing(currentwidth, currentheight, fname)
    origin()
    background("pink")
    setopacity(0.5)
    path1 = [Point(-50, 50), Point(60, -60)]
    path4 = ngon(O, 40, 2, vertices=true)
    path2 = ngon(O, 40, 3, vertices=true)
    path3 = ngon(O, 40, 4, vertices=true)
    path5 = ngon(O, 40, 5, vertices=true)
    path6 = ngon(O, 40, 6, vertices=true)

    paths = [path1, path2, path3, path4, path5, path6]
    paths = vcat(paths, paths, paths)

    t = Table(3, 6, 120, 150, Point(0, 300))
    for (i, path) in enumerate(paths)
        @layer begin
        translate(t[i])
        poly(path, :path)
        bpatharray = pathtobezierpaths()
        for (n, bpath) in enumerate(bpatharray)
            randomhue()
            drawbezierpath(bpath, :stroke)
            for bps in bpath
                randomhue()
                drawbezierpath(setbezierhandles(bps, angles=[2rand(), 2rand()], handles=[3rand(), 2rand()]), :fill)
            end
            # close needs doing manually
            if length(bpath) >= 1
                closingsegment = BezierPathSegment(bpath[end].cp2, bpath[end].p2, bpath[1].p1, bpath[1].cp1)
                seg = setbezierhandles(closingsegment, angles=[2rand(), -2rand()], handles=[3rand(), 2rand()])
                drawbezierpath(seg, :fill)
            end
        end
        end
    end

    # Julia!
    fontsize(500)
    fontface("Times-Bold")
    translate(-450, 100)

    textpath("julia")
    bp = pathtobezierpaths()
    newpath() # clear current path
    for i in 1:5
    for b in bp
        for bps in b
            randomhue()
            drawbezierpath(setbezierhandles(bps,
                angles=[3rand(), 3rand()], handles=[4rand(), 4rand()]),
                :fill)
        end
        # close needs doing manually
        if length(b) >= 1
            closingsegment = BezierPathSegment(b[end].cp2, b[end].p2, b[1].p1, b[1].cp1)
            seg = setbezierhandles(closingsegment, angles=[4rand(), -4rand()], handles=[4rand(), 4rand()])
            drawbezierpath(seg, :fill)
        end
    end

    setmode("overlay")
    textpath("julia")
    sethue("blue")
    fillpath()

end
    @test finish() == true
end

fname = "bezierstroke.png"
testbezierstroke(fname)
println("...finished bezier stroke tests: output in $(fname)")


text:
__precompile__()

module HigherPrecision

    import Base: +, -, *, /, ^, <, ==, <=

    include("basics.jl")
    include("computemode.jl")
    include("double.jl")
    # include("quad.jl")

end # module


text:
# HSA runtime and ROCm External Libraries
## copied from CUDAdrv/src/CUDAdrv.jl
if !parse(Bool, get(ENV, "JULIA_AMDGPU_DISABLE_ARTIFACTS", "false"))
    using hsa_rocr_jll
    using hsakmt_roct_jll
end

try
    include("ext.jl")
catch err
    if !isfile(joinpath(@__DIR__, "ext.jl"))
        @warn "Didn't find $ext, please build AMDGPU.jl"
        @eval const hsa_configured = false
        @eval const ext_libs_configured = false
    else
        rethrow(err)
    end
end

# Default (non-functional) values for critical variables,
# making it possible to load the package at all times.
if !hsa_configured
    const libhsaruntime_version = v"0.0"
    const libhsaruntime_vendor = "none"
    const libhsaruntime_path = nothing
end
if !ext_libs_configured
    const librocblas = nothing
    const librocsparse = nothing
    const librocalution = nothing
    const librocfft = nothing
    const librocrand = nothing
    const libmiopen = nothing
    const libhip = nothing
end

# ROCm-Device-Libs
device_libs_deps = joinpath(@__DIR__, "deps.jl")
isfile(device_libs_deps) && include(device_libs_deps)
const device_libs_path = joinpath(@__DIR__, "usr", "lib")

const configured = hsa_configured


text:
module NirsGadgets

using DSP
using Loess
using StatsBase
using Statistics
using Parameters
using MultivariateStats

include("./check_data.jl")
include("./outliers.jl")
include("./types.jl")
include("./rm_glb.jl")
include("./rm_spk.jl")
include("./rm_trd.jl")
include("./re_smp.jl")
include("./fl_pas.jl")
include("./ts_hbt.jl")
include("./cc_avg.jl")
include("./cc_ana.jl")
include("./cc_add.jl")
include("./rs_bnd.jl")

export outliers
export to_1ds
export rm_glb
export rm_spk
export rm_trd
export re_smp
export fl_pas
export an_hbt
export cc_avg
export cc_mod
export cc_arg
export cc_rea
export cc_add
export rs_bnd

export Neuro1DRealSignal
export Neuro1DAnalyticSignal

end # module


text:
using MittagLefflerFunctions
using SpecialFunctions: erfcx
using PyPlot

T = BigFloat

function g(y)
    x = (1+y) / (1-y)
    return ( 2*erfcx(x) / (1-y) - 1 ) / (1+y)
end

function erfcx_approx(x)
    y = (x-1) / (x+1)
    return (1-y)*(1+(1+y)*g_approx(y)) / 2
end

nmax = 30
M = nmax + 6
a = chebyshev_coefs(T, g, nmax, M)

y = range(-one(T), one(T), length=201)
g_approx(y) = chebyshev_sum(a, y)

figure(1)
plot(y, g.(y))
grid(true)
xlabel(L"$y$")
ylabel(L"$g(y)$")

figure(2)
plot(y, g.(y)-g_approx.(y))
grid(true)
xlabel(L"$y$")
ylabel(L"error in $g$")

figure(3)
semilogy(abs.(a))
grid(true)
xlabel(L"$n$")
ylabel(L"$|a_n|$")

figure(4)
x = range(zero(T), big(10.0), length=201)
plot(x, erfcx_approx.(x) - erfcx.(x))
grid(true)
xlabel(L"$x$")
ylabel("error in erfcx")


text:
@testset "Examples" begin
  include("Examples/MultiQuad.jl")
end


text:
function mw_s_cov_ha(cov_h0, param, rnd_delta; big=true)#, rmatrix_string)
  # s_cov_ha - Generates nsim HA covariance matrices for S-test
  # param = local parameter
  # nsim = number of simulations
  # rmatrix_string = name of matrix with standard random variables (on disk)
  nsim = size(rnd_delta, 2)
  q = size(cov_h0, 1)
  param = 5.0
  qq =Int( q*(q+1)/2 )
  delta = rnd_delta
  delta = cumsum(delta[1:q, 1:nsim])/sqrt(q)
  delta = delta*(param/sqrt(q))
  gam = exp.(delta)
  cov_ha = NaN*zeros(qq, nsim)
  big_cov_ha = NaN*zeros(q, q, nsim)
  for i ∈ 1:nsim
    tmp = (cov_h0.*gam[:, i]).*gam[:, i]'
    cov_ha[:, i] = vec(tril(tmp))[vec(tril(tmp).!=0.0)]
    big_cov_ha[:, :, i] = Hermitian(tmp, :L)
  end
  (big==true) && return big_cov_ha
  return cov_ha
end
#
function random_gen_delta(t1, t2; seed=89871924)
  srand(seed)
  rnd_delta = rand(Normal(), t1, t2)
  return rnd_delta
end


text:
"""
TODO: document
"""
struct MvEpd{T <: Real, Cov<:AbstractPDMat, Mean<:AbstractVector, shape<:Real} <: ContinuousMultivariateDistribution
    μ::Mean
    Σ::Cov
    p::shape
end

### Construction
function MvEpd(μ::AbstractVector{T}, Σ::AbstractPDMat{T}, p::T) where {T <: Real}
    dim(Σ) == length(μ) || throw(DimensionMismatch("The dimensions of mu and Sigma are inconsistent."))
    p > 0 || throw(DomainError("p must be positive"))
    MvEpd{T, typeof(Σ), typeof(μ), typeof(p)}(μ, Σ, p)
end

function MvEpd(μ::AbstractVector{<:Real}, Σ::AbstractPDMat, p::Real)
    R = Base.promote_eltype(μ, Σ)
    MvEpd(convert(AbstractArray{R}, μ), convert(AbstractArray{R}, Σ), p)
end

# constructor with general covariance matrix
MvEpd(μ::AbstractVector{<:Real}, Σ::AbstractMatrix{<:Real}, p::Real) = MvEpd(μ, PDMat(Σ), p)
MvEpd(μ::AbstractVector{<:Real}, Σ::Diagonal{<:Real}, p::Real) = MvEpd(μ, PDiagMat(diag(Σ)), p)
MvEpd(μ::AbstractVector{<:Real}, Σ::UniformScaling{<:Real}, p::Real) =
    MvEpd(μ, ScalMat(length(μ), Σ.λ), p)

# constructor with vector of standard deviations
MvEpd(μ::AbstractVector{<:Real}, σ::AbstractVector{<:Real}, p::Real) = MvEpd(μ, PDiagMat(abs2.(σ)), p)

# constructor with scalar standard deviation
MvEpd(μ::AbstractVector{<:Real}, σ::Real, p::Real) = MvEpd(μ, ScalMat(length(μ), abs2(σ)), p)

# constructor without mean vector
MvEpd(Σ::AbstractVecOrMat{<:Real}, p::Real) = MvEpd(Zeros{eltype(Σ)}(size(Σ, 1)), Σ, p)

# special constructor
MvEpd(d::Int, σ::Real, p::Real) = MvEpd(Zeros{typeof(σ)}(d), σ, p)

### Basic statistics
length(d::MvEpd) = length(d.μ)
location(d::MvEpd) = d.μ
params(d::MvEpd) = (d.μ, d.Σ, d.p)
@inline partype(d::MvEpd{T}) where {T<:Real} = T
Base.eltype(::Type{<:MvEpd{T}}) where {T} = T

cov(d::MvEpd) = d.Σ * 2^(1/d.p) *
    gamma((length(d.μ)+2)/(2*d.p)) / (length(d.μ) * gamma(length(d.μ)/(2*d.p))) |> Matrix

invcov(d::MvEpd) = Matrix(inv(cov(d)))

### Evaluation
insupport(d::MvEpd, x::AbstractVector{T}) where {T<:Real} =
    length(d) == length(x) && all(isfinite, x)

sqmahal(d::MvEpd, x::AbstractVector{<:Real}) = invquad(d.Σ, x - d.μ)

function sqmahal!(r::AbstractArray, d::MvEpd, x::AbstractMatrix{<:Real})
    invquad!(r, d.Σ, x .- d.μ)
end

sqmahal(d::MvEpd, x::AbstractMatrix{T}) where {T<:Real} = sqmahal!(Vector{T}(undef, size(x, 2)), d, x)

function mvepd_consts(d::MvEpd)
    n = length(d.μ)
    log(n) + loggamma(n/2) - n/2 * log(π) - loggamma(1+n/(2*d.p)) - (1+n/(2*d.p))*log(2) - 0.5*logdet(d.Σ)
end

_logpdf(d::MvEpd, x::AbstractVector{T}) where {T<:Real} = mvepd_consts(d) - 0.5 * sqmahal(d, x)^d.p

function _logpdf!(r::AbstractArray, d::MvEpd, x::AbstractMatrix)
    sqmahal!(r, d, x)
    println(r)
    k = mvepd_consts(d)
    for i = 1:size(x, 2)
        @inbounds r[i] = k - 0.5 * r[i]^d.p
    end
    r
end

_pdf!(r::AbstractArray, d::MvEpd, x::AbstractVector{T}) where {T<:Real} = exp.(r, _logpdf!(r, d, x))

function runifsphere(rng::AbstractRNG, n::Integer, k::Integer)
    k >= 2 || throw(DomainError(k, "k must be [2, ∞)"))
    X = reshape(rand(rng, Normal(), n*k), (n,k))
    rownorms = .√sum(X.^2, dims=2)
    for i in 1:size(X, 2)
        X[:,i] = X[:,i] ./ rownorms
    end
    X
end

function _rand!(rng::AbstractRNG, d::MvEpd, x::AbstractMatrix{T}) where T<:Real
    μ, Σ, p = params(d)
    k = size(Σ, 1)
    n = size(x, 2)
    val, ev = eigen(Σ)
    Σsqrt = ev *  diagm(.√val) * ev'
    radius = (rand(rng,Gamma(k/(2*p), 1/2), n)).^(1/(2*p))
    un = reshape(rand(rng, Normal(), n*k), (n,k))
    rownorms = .√sum(un.^2, dims=2)
    for i in 1:size(un, 2)
        un[:,i] = un[:,i] ./ rownorms
    end
    x = μ .+ reshape(radius .* un * Σsqrt, (k,n))
    x
end


text:
addIntegrator_! =
Dict{String,Function}(
    "Euler"=>addIntegratorEuler_!,
    "Heun"=>addIntegratorHeun_!,
)

"""
    function updateVariables_!(abm::Agent,p::Program)

Adds to p.update those variables that should be updated.
"""
function updateVariables_!(abm::Agent,p::Program_)

    #Check updated
    up = symbols_(abm,abm.declaredUpdates["Equation"])
    up = up[(up[:,"assigned"].==true),:]

    var = [Meta.parse(string(i,"̇")) for i in abm.declaredSymbols["Local"]]
    for (i,j) in zip(var,abm.declaredSymbols["Local"])
        if i in up.Symbol
            push!(p.update,j)
        end
    end

    return
end

max_stars_repo_name: gatocor/embryogenesisJulia

max_issues_repo_name: gatocor/embryogenesisJulia

max_forks_repo_name: gatocor/embryogenesisJulia


text:
# Autogenerated wrapper script for AsmJit_jll for i686-w64-mingw32
export libasmjit

JLLWrappers.@generate_wrapper_header("AsmJit")
JLLWrappers.@declare_library_product(libasmjit, "libasmjit.dll")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libasmjit,
        "bin\\libasmjit.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""
Collection of miscellaneous utilities

"""
module SLIMutils

# Vector{Vector} nested blocks
include("nested_blocks.jl")

end # module


text:
using IRTracker.GraphViz

function pushstmt!(stmts, ref, stmt::Constant)
    value = escape_string(sprint(showstmt, ref, stmt))
    push!(stmts, GraphViz.Node(string(ref.number), label=value, shape="rectangle"))
end
function pushstmt!(stmts, ref, stmt::Call)
    argstring = join(stmt.args, ", ")
    label = escape_string(sprint(showstmt, ref, stmt))
    push!(stmts, GraphViz.Node(string(ref.number), label=label, shape="rectangle"))
end
function pushstmt!(stmts, ref, stmt::Union{Assumption, Observation})
    label = escape_string(sprint(showstmt, ref, stmt))
    push!(stmts, GraphViz.Node(string(ref.number), label=label, shape="circle"))
end

function Base.convert(::Type{GraphViz.Graph}, graph::Graph)
    stmts = Vector{GraphViz.Statement}()

    for (ref, stmt) in graph
        pushstmt!(stmts, ref, stmt)

        for dep in dependencies(stmt)
            from = GraphViz.NodeID(string(ref.number))
            to = GraphViz.NodeID(string(dep.number))
            push!(stmts, GraphViz.Edge([from, to]))
        end

        if ref isa NamedReference
            from = GraphViz.NodeID(string(ref.number))
            @show ref
            for index in DynamicPPL.getindexing(ref.vn), ix in index
                if ix isa Reference
                    to = GraphViz.NodeID(string(ix.number))
                    @show ix
                    push!(stmts, GraphViz.Edge([from, to]))
                end
            end
        end
    end

    graph_attrs = Dict(:ordering => "in",    # edges sorted by incoming
                       :rankdir => "BT",     # order nodes from right to left
                       )
    edge_attrs = Dict{Symbol, String}()
    node_attrs = Dict(:shape => "plaintext")

    return GraphViz.DiGraph(stmts,
                            graph_attrs = graph_attrs,
                            edge_attrs = edge_attrs,
                            node_attrs = node_attrs)
end

function savedot(fn::AbstractString, graph::Graph)
    open(fn, "w") do fp
        GraphViz.pprint(fp, Base.convert(GraphViz.Graph, graph))
    end
end

# dot /tmp/graph.dot -Tpdf -Nfontname="DejaVu Sans Mono" -Efontname="DejaVu Sans Mono" > /tmp/graph.pdf