text:
"""
Use krylovkit to perform exact diagonalization
"""
function exact_diagonalization(opp::MPOHamiltonian;sector = first(sectors(oneunit(opp.pspaces[1]))),len::Int = opp.period,num::Int = 1,which::Symbol=:SR,alg::KrylovKit.KrylovAlgorithm = Lanczos())
    left = ℂ[sector=>1];
    right = oneunit(left);

    middle_site = Int(round(len/2));

    Ot = eltype(opp);

    mpst_type = tensormaptype(spacetype(Ot),2,1,eltype(Ot));
    mpsb_type = tensormaptype(spacetype(Ot),1,1,eltype(Ot));
    CLs = Vector{Union{Missing,mpsb_type}}(missing,len+1);
    ALs = Vector{Union{Missing,mpst_type}}(missing,len);
    ARs = Vector{Union{Missing,mpst_type}}(missing,len);
    ACs = Vector{Union{Missing,mpst_type}}(missing,len);

    for i in 1:middle_site-1
        ALs[i] = isomorphism(storagetype(Ot),left*opp.pspaces[i],fuse(left*opp.pspaces[i]));
        left = _lastspace(ALs[i])';
    end
    for i in len:-1:middle_site+1
        ARs[i] = _permute_front(isomorphism(storagetype(Ot),fuse(opp.pspaces[i]'*right'),opp.pspaces[i]'*right));
        right = _firstspace(ARs[i]);
    end
    ACs[middle_site] = TensorMap(rand,ComplexF64,left*opp.pspaces[middle_site],right);
    norm(ACs[middle_site]) == 0 && throw(ArgumentError("invalid sector"));
    normalize!(ACs[middle_site]);

    #construct the largest possible finite mps of that length
    state = FiniteMPS{mpst_type,mpsb_type}(ALs,ARs,ACs,CLs);
    envs = environments(state,opp);

    #optimize the middle site. Because there is no truncation, this single site captures the entire possible hilbert space
    (vals,vecs,convhist) = eigsolve(state.AC[middle_site],num,which,alg) do x
        ac_prime(x,middle_site,state,envs)
    end

    state_vecs = map(vecs) do v
        cs = copy(state);
        cs.AC[middle_site] = v;
        cs
    end

    return vals,state_vecs,convhist
end


text:
module DocExample

export summa, calcQuadRoots,
        showText, showFirst

"""
    calcQuadRoots(a,b,c)

Returns double the number `x` plus `1`.
"""

include("mathFoo.jl")
include("textFoo.jl")

end # module

ext: jl


text:
# This file was generated, do not modify it. # hide
allknntime = @elapsed knns, dists = allknn(G, k; parallel=true)
nothing # hide

max_stars_repo_name: sadit/SimilaritySearchDemos

max_issues_repo_name: sadit/SimilaritySearchDemos

max_forks_repo_name: sadit/SimilaritySearchDemos


text:
import SparseArrays: SparseMatrixCSC, SparseMatrixCSCInterface, is_hermsym
is_hermsym(A::Number, check::Function) = A == check(A)

function is_hermsym3(A::SparseMatrixCSCInterface{Tv,Ti}, check::Function) where {Tv,Ti}
    m, n = size(A)
    m == n || return false
    rowval = rowvals(A)
    nzval = nonzeros(A)
    start = Vector{Int}(undef, n)
    toend = A isa SparseMatrixCSC || last(A.indices[1]) == size(parent(A), 1)
    stop = toend ? parent(A).colptr : Vector{Ti}(undef, n+1)
    for col = 1:n
        nzr = nzrange(A, col)
        if !toend
            stop[col+1] = last(nzr) + 1
        end
        start[col] = stop[col+1]
        for j in nzr
            row = rowval[j]
            row == 0 && continue
            if row >= col
                if row == col
                    is_hermsym(nzval[j], check) || return false
                    start[col] = j + 1
                 else
                    start[col] = j
                end
                break
            else # row < col
                nz = nzval[j]
                found = false
                for k = start[row]:stop[row+1]-1
                    mcol = rowval[k]
                    if mcol < col
                        iszero(nzval[k]) || return false
                    elseif mcol == col
                        nzval[k] == check(nz) || return false
                        start[row] = k + 1
                        found = true
                        break
                    else
                        start[row] = k
                        break
                    end
                end
                found || iszero(nz) || return false
            end
        end
    end
    # return all(iszero(nzval[k]) for row = 1:n, k = start[row]:stop[row+1]-1)
    for row = 1:n, k = start[row]:stop[row+1]-1
        iszero(nzval[k]) || return false
    end
    return true
end

function is_hermsym2(A::SubArray{<:Any,2,<:SparseMatrixCSC,<:Tuple{<:StepRange{<:Integer,<:Integer},<:AbstractVector{<:Integer}},false}, check::Function)

    st = step(A.indices[1])
    if st == 1
        ind = (first(A.indices[1]):last(A.indices[1]), A.indices[2])
    elseif st == -1
        ind = (last(A.indices[1]):first(A.indices[1]), reverse(A.indices[2]))
    else
        throw(ArgumentError("view step of first index is $st not in {1, -1}"))
    end
    is_hermsym2(view(parent(A), ind...), check)
end

function is_hermsym2(A::SparseMatrixCSCInterface, check::Function)
    m, n = size(A)
    if m != n; return false; end

    rowval = rowvals(A)
    nzval = nonzeros(A)
    tracker, stop = tracker_stop(A)
    for col = 1:n
        # `tracker` is updated such that, for symmetric matrices,
        # the loop below starts from an element at or below the
        # diagonal element of column `col`"
        for p = tracker[col]:stop[col+1]-1
            val = nzval[p]
            row = rowval[p]

            # Ignore stored zeros
            if iszero(val)
                continue
            end

            # If the matrix was symmetric we should have updated
            # the tracker to start at the diagonal or below. Here
            # we are above the diagonal so the matrix can't be symmetric.
            if row < col
                return false
            end

            # Diagonal element
            if row == col
                if !is_hermsym(val, check)
                    return false
                end
            else
                offset = tracker[row]

                # If the matrix is unsymmetric, there might not exist
                # a rowval[offset]
                if offset > length(rowval)
                    return false
                end

                row2 = rowval[offset]

                # row2 can be less than col if the tracker didn't
                # get updated due to stored zeros in previous elements.
                # We therefore "catch up" here while making sure that
                # the elements are actually zero.
                while row2 < col
                    if !iszero(nzval[offset])
                        return false
                    end
                    offset += 1
                    row2 = rowval[offset]
                    tracker[row] += 1
                end

                # Non zero A[i,j] exists but A[j,i] does not exist
                if row2 > col
                    return false
                end

                # A[i,j] and A[j,i] exists
                if row2 == col
                    if val != check(nzval[offset])
                        return false
                    end
                    tracker[row] += 1
                end
            end
        end
    end
    return true
end

function tracker_stop(A::SparseMatrixCSC)
    copy(A.colptr), A.colptr
end
function tracker_stop(A::SparseMatrixCSCInterface{Tv,Ti}) where {Tv,Ti}
    n = size(A, 2)
    stop = Vector{Ti}(undef, n+1)
    stop[1] = 1
    tracker = Vector{Ti}(undef, n)
    for col = 1:n
        nzr = nzrange(A, col)
        tracker[col] = first(nzr)
        stop[col+1] = last(nzr) + 1
    end
    tracker, stop
end


text:
## header to provide surface and context
using Cairo
c = CairoRGBSurface(256,256);
cr = CairoContext(c);

save(cr);
set_source_rgb(cr,0.8,0.8,0.8);    # light gray
rectangle(cr,0.0,0.0,256.0,256.0); # background
fill(cr);
restore(cr);

save(cr);

set_font_face(cr, "Sans 16")

text(cr,16.0,40.0,"Hamburgefons")
text(cr,16.0,72.0,"sp⁰¹²³,min⁻²,αΑβΒφϕΦγΓ")
text(cr,16.0,104.0,"TextBoldItalicsuper-2",markup=true)

text(cr,40.0,224.0,"Es geht aufwärts !",markup=true,angle=30.0)

#using textwidth and height

set_font_face(cr, "Sans 12")

a = "A"
aheight = textheight(cr,a)
awidth = textwidth(cr,a)
atext = @sprintf("%s wd=%2.2f,ht=%2.2f",a,awidth,aheight)
text(cr,16.0,240.0,atext,markup=true)

## mark picture with current date
restore(cr);
move_to(cr,0.0,12.0);
set_source_rgb(cr, 0,0,0);
show_text(cr,Libc.strftime(time()));
write_to_png(c,"sample_pango_text.png");


text:
# AbstractGridFunction

export AbstractGridFunction
"""
    abstract type AbstractGridFunction{S,T,D} end

An abstract grid function with domain `S` and codomain `T` in dimension `D`.

    AbstractGridFunction
    AbstractGridHierarchy
    AbstractGridMesh
    AbstractGridComponent
"""
abstract type AbstractGridFunction{S,T,D} end

Base.ndims(::AbstractGridFunction{S,T,D} where {S,T}) where {D} = D

Base.keytype(::AbstractGridFunction{S,T,D}) where {S,T,D} = SVector{D,S}

Base.eltype(::AbstractGridFunction{S,T}) where {S,T} = T

export domain
domain(::AbstractGridFunction) = error("undefined method")

Base.zero(::AbstractGridFunction) = error("undefined method")
Base.:+(::AbstractGridFunction) = error("undefined method")
Base.:-(::AbstractGridFunction) = error("undefined method")
Base.:+(::AbstractGridFunction{S,T,D} where {T}, ::AbstractGridFunction{S,T,D} where {T}) where {S,D} = error("undefined method")
Base.:-(::AbstractGridFunction{S,T,D} where {T}, ::AbstractGridFunction{S,T,D} where {T}) where {S,D} = error("undefined method")
Base.:*(::Number, ::AbstractGridFunction) = error("undefined method")
Base.:*(::AbstractGridFunction, ::Number) = error("undefined method")
Base.:\(::Number, ::AbstractGridFunction) = error("undefined method")
Base.:/(::AbstractGridFunction, ::Number) = error("undefined method")

export evaluate
evaluate(::AbstractGridFunction{S,T,D}, ::SVector{D,S}) where {S,T,D} = error("undefined method")

export test_AbstractGridFunction
function test_AbstractGridFunction(gf::AbstractGridFunction)
    @testset "$(typeof(gf)) <: AbstractGridFunction" begin
        D = ndims(gf)
        @test D isa Int
        @test D ≥ 0

        @test keytype(gf) <: AbstractArray
        S = eltype(keytype(gf))
        @test S isa Type

        T = eltype(gf)
        @test T isa Type

        gfdomain = domain(gf)
        @test gfdomain isa Domain
        @test ndims(gfdomain) == D
        @test eltype(gfdomain) == S

        len = length(gf)
        @test len isa Int
        @test len ≥ 0

        rng = MersenneTwister(0)

        for val in [1, 1.0, SVector(1.0f0, 1.0f0)]
            T = typeof(val)

            function test_props(gf)
                @test gf isa AbstractGridFunction
                @test eltype(keytype(gf)) == S
                @test eltype(gf) == T
                @test ndims(gf) == D
            end

            e = map(_ -> val, gf)
            test_props(e)

            @test reduce(+, e; init=zero(val)) == len * val
            @test mapreduce(sum, +, e; init=sum(zero(val))) == length(val) * len

            c = copy(e)
            s = similar(e)
            test_props(c)
            test_props(s)
            fill!(s, val)
            @test s == c
            @test s ≈ e

            n = zero(e)
            if T <: AbstractFloat
                x = map(v -> T(v[1]), coordinates(e))
            else
                x = rand(rng, e)
            end
            y = rand(rng, e)
            z = rand(rng, e)
            a = rand(rng, eltype(T))
            b = rand(rng, eltype(T))
            if T <: Integer
                # Prevent integer overflow
                map!(v -> v % 1000, x, x)
                map!(v -> v % 1000, y, y)
                map!(v -> v % 1000, z, z)
                a %= 1000
                b %= 1000
            end

            test_props(n)
            test_props(x)
            test_props(y)
            test_props(z)

            @test map(identity, y) == y
            @test map(a -> a .+ 1, map(a -> 2a, y)) ≈ map((a -> a .+ 1) ∘ (a -> 2a), y)

            @test +x == x
            @test n + x == x
            @test x + n == x
            @test (x + y) + z ≈ x + (y + z)

            @test x + (-x) == n
            @test x - y == x + (-y)

            @test zero(a) * x == n
            @test one(a) * x == x

            @test a * n == n

            @test a * x == x * a

            @test (a * b) * x ≈ a * (b * x)
            @test (x * a) * b ≈ x * (a * b)

            @test a * (x + y) ≈ a * x + a * y
            @test (a + b) * x ≈ a * x + b * x

            if a ≠ 0
                @test inv(a) * (a * x) ≈ x
                @test a \ x ≈ inv(a) * x
                @test x / a ≈ x * inv(a)
            end
        end

        T = eltype(gf)
        for n in 1:10
            x = map((a, b) -> (a + b) / 2 + (b - a) / 2 * (2 * rand(S) - 1), min(gfdomain), max(gfdomain))
            @test x isa SVector{D,S}
            fx = evaluate(gf, x)
            @test fx isa T
        end
    end
end


text:
using Stheno: ConstantMean, ITMean, ITKernel, map, pairwise

@testset "input_transform" begin

    # Test ITMean.
    let
        rng, N, D = MersenneTwister(123456), 10, 2
        μ, f, X = ConstantMean(randn(rng)), x->sum(abs2, x), ColsAreObs(randn(rng, D, N))
        μf = ITMean(μ, f)

        @test ITMean(μ, identity) == μ
        @test length(μf) == Inf
        @test μf(X[1]) == (μ ∘ f)(X[1])
        mean_function_tests(μf, X)
    end

    # Test ITKernel.
    let
        rng, N, N′, D = MersenneTwister(123456), 10, 11, 2
        k, f = EQ(), x->sum(abs2, x)
        X0, X1, X2 = ColsAreObs(randn(rng, D, N)), ColsAreObs(randn(rng, D, N)), ColsAreObs(randn(rng, D, N′))
        kf = ITKernel(k, f)

        @test ITKernel(k, identity) == k
        @test size(kf, 1) == Inf && size(kf, 2) == Inf

        @test kf(X0[1], X1[1]) == k(f(X0[1]), f(X1[1]))
        kernel_tests(kf, X0, X1, X2, 1e-6)
    end

    # Test convenience code.
    let
        rng, N, N′, D = MersenneTwister(123456), 10, 11, 2
        m, k = ConstantMean(randn(rng)), EQ()
        X0, X1, X2 = ColsAreObs(randn(rng, D, N)), ColsAreObs(randn(rng, D, N)), ColsAreObs(randn(rng, D, N′))
        m1, k1 = pick_dims(m, 1), pick_dims(k, 1)

        @test m1(X0[1]) == m(X0[1][1])
        @test k1(X0[1], X0[2]) == k(X0[1][1], X0[2][1])
        mean_function_tests(m1, X0)
        kernel_tests(k1, X0, X1, X2)

        x0, x1, x2 = randn(rng, N), randn(rng, N), randn(rng, N′)
        mean_function_tests(periodic(m, 0.1), x0)
        kernel_tests(periodic(k, 0.1), x0, x1, x2)
    end
end


text:
function solve_modelST_TR(nlp_stop, X :: PDataST, δ:: T; cgtol :: T = 0.1) where T
# cas particulier Steihaug-Toint
    ϵ = sqrt(eps(T)) # * 100.0
    n = length(X.g)
    cgtol = max(ϵ, min(cgtol, 9 * cgtol / 10, 0.01 * norm(X.g)^(1.0 + X.ζ)))

    (d, cg_stats) = cg(X.H, -X.g, atol = cgtol, rtol = ϵ,
                     radius = δ,  itmax = max(2 * n, 50),
                     verbose = false)

    λ = 0.0  #  dummy for this variant

    return d, NaN * rand(length(d)), λ#, cgtol
end


text:
using Documenter, DLM

makedocs(;
    modules = [DLM],
    format = Documenter.HTML(),
    pages = [
        "Home" => "index.md",
        "Example" => "example.md",
        "Library" => "library.md"
    ],
    sitename = "DLM.jl",
    authors = "Taylor McDonnell <taylormcd@byu.edu>",
)

deploydocs(
    repo="github.com/byuflowlab/DLM.jl.git",
)


text:
using AerosolMTK
using Test

@testset "AerosolMTK.jl" begin
    # Write your tests here.
end


text:
module Percolation

export dfspercolate, dfspercolate!, colorpercolate, hkpercolate,
    clustgyration, hklabel, hkmaxclust!, genclust, clustfractal

using DataStructures

"""
Find out whether or not percolation is possible in a 2D grid using depth-first search

Percolation is assumed to occur from top to bottom.

# Argument

`grid::Matrix{Bool}`: Represents the grid cells; True if porous (part of the path;
    can pass through) and False if blocked (can't pass through; no path passes through)

# Returns

boolean indicating whether percolation is possible
"""
function dfspercolate(grid::BitMatrix)
    height, width = size(grid)

    undiscovered = BitMatrix(undef, height + 2, width + 2)
    # "pad" the array to avoid bound checking for neigboring
    undiscovered[:, 1] .= false
    undiscovered[1, :] .= false
    undiscovered[end, :] .= false
    undiscovered[:, end] .= false
    # insert grid into undiscovered array
    undiscovered[2:end-1, 2:end-1] = grid

    stack = Stack{Tuple{Integer, Integer}}()
    for root in 2:width+1
        if undiscovered[2, root] # top row
            # search this path
            push!(stack, (2, root))
            if dfs!(height + 1, stack, undiscovered)
                return true
            end
        end
    end
    return false
end

"""
Depth-first search algorithm adopted for the dfspercolate function
"""
@inline function dfs!(destination::Integer, stack::Stack, undiscovered::BitMatrix)
    while !isempty(stack)
        vertex = pop!(stack)
        undiscovered[vertex...] = false

        # check if destination is reached
        if vertex[1] == destination
            return true
        end

        # check neighboring vertices
        for step in ((-1, 0), (0, -1), (0, 1), (1, 0))
            neighbor = vertex .+ step
            if undiscovered[neighbor...]
                push!(stack, neighbor)
            end
        end
    end
    return false
end

"""
Find out whether or not percolation is possible in a 2D grid using depth-first search

Percolation is assumed to occur from top to bottom. this version marks points visited by
the dfs in the visited

# Arguments

- `grid::Matrix{Bool}`: Represents the grid cells; True if porous (part of the path;
    can pass through) and False if blocked (can't pass through; no path passes through)
- `visited::BitMatrix`: whether each point is visited by the dfs

# Returns

boolean indicating whether percolation is possible
"""
function dfspercolate!(grid::BitMatrix, visited::BitMatrix)
    height, width = size(grid)

    undiscovered = BitMatrix(undef, height + 2, width + 2)
    # "pad" the array to avoid bound checking for neigboring
    undiscovered[:, 1] .= false
    undiscovered[1, :] .= false
    undiscovered[end, :] .= false
    undiscovered[:, end] .= false
    # insert grid into undiscovered array
    undiscovered[2:end-1, 2:end-1] = grid

    stack = Stack{Tuple{Integer, Integer}}()
    for root in 2:width+1
        if undiscovered[2, root] # top row
            # search this path
            push!(stack, (2, root))
            if dfs!(height + 1, stack, undiscovered, visited)
                return true
            end
        end
    end
    return false
end

"""
Depth-first search algorithm adopted for the dfspercolate function

this version modifies the `visited` matrix
"""
@inline function dfs!(destination::Integer, stack::Stack, undiscovered::BitMatrix,
        visited::BitMatrix)
    while !isempty(stack)
        vertex = pop!(stack)
        undiscovered[vertex...] = false
        visited[(vertex .- (1, 1))...] = true

        # check if destination is reached
        if vertex[1] == destination
            return true
        end

        # check neighboring vertices
        for step in ((-1, 0), (0, -1), (0, 1), (1, 0))
            neighbor = vertex .+ step
            if undiscovered[neighbor...]
                push!(stack, neighbor)
            end
        end
    end
    return false
end

"""
labels the clusters in a 2D grid for solving the percolation problem

Uses and inefficient O((length * width)^2) coloring algorithm.
Percolation is assumed to occur from top to bottom.

# Argument

`grid::Matrix{Bool}`: Represents the grid cells; True if porous (part of the path;
    can pass through) and False if blocked (can't pass through; no path passes through)

# Returns

a boolean indicating whether percolation is possible and a matrix "colored" with integers.
Points with value 1 are connected to the starting row.
"""
function colorpercolate(grid::BitMatrix)
    height, width = size(grid)
    colors = zeros(Integer, height + 1, width)
    colors[1, :] .= 1
    counter = 2

    for row in 2:height+1
        for col in 1:width
            if grid[row - 1, col]
                neighbors = neighborhood(colors, row, col)
                nonzero = filter(n -> colors[n...] != 0, neighbors)
                if length(nonzero) == 0
                    colors[row, col] = counter
                    counter += 1
                elseif length(nonzero) == 1
                    colors[row, col] = colors[nonzero[1]...]
                else
                    minneighbor = minimum([colors[n...] for n in nonzero])
                    colors[row, col] = minneighbor
                    for neighbor in nonzero
                        if colors[neighbor...] != minneighbor
                            colors[colors .== colors[neighbor...]] .= minneighbor
                        end
                    end
                end
            end
        end
    end
    return any(x -> x == 1, colors[end, :]), colors[2:end, :]
end

"""
Retruns a vector of tuples of the neighbors of the point at [row, col]
"""
@inline function neighborhood(grid::Matrix, row::Integer, col::Integer)
    neighbors = [(row + 1, col), (row - 1, col), (row, col + 1), (row, col - 1)]
    return filter(n -> checkbounds(Bool, grid, n...), neighbors)
end

"""
Uses the Hoshen-Kopelman algorithm to label the clusters in a 2D grid for percolation

Percolation is assumed to occur from top to bottom.

# Arguments

- `grid::Matrix{Bool}`: Represents the grid cells; True if porous (part of the path;
    can pass through) and False if blocked (can't pass through; no path passes through)
- `color::Bool` (default = `false`): whether to color the grid after finishing the algorithm

# Returns

[a boolean indicating whether percolation is possible or, if qinfty is true,
the sum of the sizes of the open connecting clusters (integer)],
an integer matrix where clusters, are marked by unique integers, and the cluster lengths.
"""
function hkpercolate(grid::BitMatrix; color::Bool=false, qinfty::Bool=false)
    height, width = size(grid)
    labels = zeros(Integer, height, width)
    list = Vector{Integer}()
    clusters = Vector{Integer}()

    # first row, special case because of top boundary
    hkfirstrow!(grid, labels, list, clusters)

    for row in 2:height
        # first column, special case because of left boundary
        if grid[row, 1]
            hkfirstcol!(labels, list, clusters, row)
        end

        # "middle" points, with appropriate neighbors
        for col in 2:width
            if grid[row, col]
                hkiter!(labels, list, clusters, row, col)
            end
        end
    end

    if color
        hkcolor!(labels, list)
    end
    if qinfty
        return hkresult(labels, list, clusters), labels, clusters
    else
        return hkbool(labels, list), labels, clusters
    end
end

"""
Finds the radius of gyration for the largest closed cluster in the percolation grid

Percolation is assumed to occur from top to bottom. Uses the Hoshen-Kopelman algorithm.

# Argument

`grid::Matrix{Bool}`: Represents the grid cells; True if porous (part of the path;
    can pass through) and False if blocked (can't pass through; no path passes through)

# Returns

The gyration radius of the largest closed cluster (the largest cluster that doesn't connect
the top and bottom of the grid)
"""
function clustgyration(grid::BitMatrix)
    labels, list, clusters = hklabel(grid)
    if length(clusters) == 0
        return 0
    end

    maxlabel = hkmaxclust!(labels, list, clusters)
    if maxlabel == 0
        return 0
    end

    maxclust = Tuple.(findall(x -> x != 0 && hkfind(list, x) == maxlabel, labels))
    mass = length(maxclust)
    if mass == 0
        return 0
    end

    centerofmass = reduce(.+, maxclust) ./ mass
    return √(sum(pos -> reduce(.+, (pos .- centerofmass).^2),
        maxclust) / mass)
end

"""
Uses the Hoshen-Kopelman algorithm to label the grid

returns the labels, label mapping list, and cluster lengths
"""
@inline function hklabel(grid::BitMatrix)
    height, width = size(grid)
    labels = zeros(Integer, height, width)
    list = Vector{Integer}()
    clusters = Vector{Integer}()

    # first row, special case because of top boundary
    hkfirstrow!(grid, labels, list, clusters)

    for row in 2:height
        # first column, special case because of left boundary
        if grid[row, 1]
            hkfirstcol!(labels, list, clusters, row)
        end

        # "middle" points, with appropriate neighbors
        for col in 2:width
            if grid[row, col]
                hkiter!(labels, list, clusters, row, col)
            end
        end
    end

    return labels, list, clusters
end

"""
Finds the label of the largest closed cluster (a.k.a. finite cluster).

Also, sets the size of open clusters (a.k.a. infinite cluster) to zero, hence the !
"""
@inline function hkmaxclust!(labels::Matrix{Integer}, list::Vector{Integer},
        clusters::Vector{Integer})
    for root in labels[1, :]
        if root > 0
            for dest in labels[end, :]
                if dest > 0
                    toplabel = hkfind(list, root)
                    bottomlabel = hkfind(list, dest)
                    if toplabel == bottomlabel
                        clusters[toplabel] = 0
                    end
                end
            end
        end
    end

    if maximum(clusters) > 0
        return argmax(clusters)
    else
        return 0
    end
end

"""
Uses a depth-first search algorithm to generate a cluster with the given site probability

# Arguments

- `probability`: the probability a site is generated when growing the cluster
- `maxsize::Integer` (optional, defaults to 10000): the maximum size of the cluster
    (breaks the loop after the cluster is grown to this size)

# Returns

an integer matrix (8-bit) with 0 where no site exists, 1 where a cell is generated, and 2
where the growth is stopped
"""
function genclust(probability::Float64; maxsize::Integer=10000, crop::Bool=true)
    cluster = zeros(Int8, maxsize, maxsize)
    # create "fences" to stop the cluster form reaching the bounds of the matrix
    cluster[1, :] .= 2
    cluster[:, 1] .= 2
    cluster[end, :] .= 2
    cluster[:, end] .= 2

    stack = Stack{Tuple{Integer, Integer}}()
    seed = (maxsize÷2, maxsize÷2)
    push!(stack, seed)
    cluster[seed...] = 1

    clustdfs!(stack, cluster, probability, maxsize)

    if crop
        return cropclust!(cluster)
    else
        return cluster
    end
end

"""
The depth-first search algorithm used for `genclust`
"""
@inline function clustdfs!(stack::Stack, cluster::Matrix{<:Integer},
        probability::Float64, maxsize::Integer)
    size = 1
    while !isempty(stack) && size < maxsize
        site = pop!(stack)
        for step in ((-1, 0), (0, -1), (0, 1), (1, 0))
            neighbor = site .+ step
            if cluster[neighbor...] == 0
                if rand() <= probability
                    cluster[neighbor...] = 1
                    push!(stack, neighbor)
                    size += 1
                else
                    cluster[neighbor...] = 2
                end
            end
        end
    end
end

"""
finds the mass and radius of gyration of the cluster from the output of `genclust`
"""
function clustfractal(cluster::Matrix{Int8})
    clust = Tuple.(findall(x -> x == 1, cluster))
    mass = length(clust)
    if mass == 0
        return 0
    end

    centerofmass = reduce(.+, clust) ./ mass
    return mass, √(sum(pos -> reduce(.+, (pos .- centerofmass).^2), clust) / mass)
end

"""
crop the `genclust` function output to remove excess zeros
"""
@inline function cropclust!(cluster::Matrix{<:Integer})
    height, width = size(cluster)
    top = bottom = height ÷ 2
    left = right = width ÷ 2

    while 1 in cluster[top + 1, 2:end-1] || 2 in cluster[top + 1, 2:end-1]
        top += 1
    end
    while 1 in cluster[bottom - 1, 2:end-1] || 2 in cluster[bottom - 1, 2:end-1]
        bottom -= 1
    end
    while 1 in cluster[2:end-1, left - 1] || 2 in cluster[2:end-1, left - 1]
        left -= 1
    end
    while 1 in cluster[2:end-1, right + 1] || 2 in cluster[2:end-1, right + 1]
        right += 1
    end

    return cluster[bottom:top, left:right]
end

"""
Helper function for `hkpercolate`; Adds new label.
"""
@inline function hknew!(labels::Matrix{<:Integer}, list::Vector{<:Integer},
        clusters::Vector{<:Integer}, row::Integer, col::Integer)
    labels[row, col] = length(list) + 1
    push!(list, labels[row, col])
    push!(clusters, 0)
end

"""
Helper function for `hkpercolate`; Applies the Hoshen-Kopelman algorithm to the first row.
"""
@inline function hkfirstrow!(grid::BitMatrix, labels::Matrix{<:Integer},
        list::Vector{<:Integer}, clusters::Vector{<:Integer})
    width = size(labels)[2]
    if grid[1, 1]
        hknew!(labels, list, clusters, 1, 1)
        clusters[labels[1, 1]] += 1
    end
    for col in 2:width
        if grid[1, col]
            if labels[1, col - 1] != 0
                labels[1, col] = hkfind!(list, labels[1, col - 1])
            else
                hknew!(labels, list, clusters, 1, col)
            end
            clusters[labels[1, col]] += 1
        end
    end
end

"""
Helper function for `hkpercolate`; Applies the Hoshen-Kopelman algorithm to the first column.
"""
@inline function hkfirstcol!(labels::Matrix{<:Integer}, list::Vector{<:Integer},
        clusters::Vector{<:Integer}, row::Integer)
    if labels[row - 1, 1] != 0
        labels[row, 1] = hkfind!(list, labels[row - 1, 1])
    else
        hknew!(labels, list, clusters, row, 1)
    end
    clusters[labels[row, 1]] += 1
end

"""
Helper function for `hkpercolate`; Applies one iteration of the Hoshen-Kopelman algorithm.
"""
@inline function hkiter!(labels::Matrix{<:Integer}, list::Vector{<:Integer},
        clusters::Vector{<:Integer}, row::Integer, col::Integer)
    above, left = labels[row - 1, col], labels[row, col - 1]
    if above == 0
        if left == 0
            hknew!(labels, list, clusters, row, col)
        else
            labels[row, col] = hkfind!(list, left)
        end
    elseif left == 0
        labels[row, col] = hkfind!(list, above)

    else
        labels[row, col] = hkunion!(list, clusters, above, left)
    end
    clusters[labels[row, col]] += 1
end

"""
Union function of the union-find algorithm used by the Hoshen-Kopelman algorithm
"""
@inline function hkunion!(list::Vector{<:Integer}, clusters::Vector{<:Integer},
        x::Integer, y::Integer)
    root = hkfind!(list, x)
    dest = hkfind!(list, y)
    if dest != root
        list[dest] = root
        clusters[root] += clusters[dest]
        clusters[dest] = 0
    end
    return root
end

"""
Find function of the union-find algorithm used by the Hoshen-Kopelman algorithm
"""
@inline function hkfind!(list::Vector{<:Integer}, x::Integer)
    y = x;
    while list[y] != y
        y = list[y];
    end
    while list[x] != x
        z = list[x];
        list[x] = y;
        x = z;
    end
    return y;
end

"""
Same as hkfind!, but doesn't update labels
"""
@inline function hkfind(list::Vector{<:Integer}, x::Integer)
    while x != list[x]
        x = list[x];
    end
    return x
end

"""
Color resultant grid of the Hoshen-Kopelman algorithm appropriately
"""
@inline function hkcolor!(labels::Matrix{<:Integer}, list::Vector{<:Integer})
    width, height = size(labels)
    for row in 1:height
        for col in 1:width
            if labels[row, col] != 0
                labels[row, col] = hkfind(list, labels[row, col])
            end
        end
    end
end

"""
Find the sum of the sizes of the open connecting clusters (integer)
"""
@inline function hkresult(labels::Matrix{<:Integer}, list::Vector{<:Integer},
        clusters::Vector{<:Integer})
    sumsize = 0
    # to avoid over-counting the clusters, save the label of each cluster that is counted
    counted = []

    for root in labels[1, :]
        if root > 0
            for dest in labels[end, :]
                if dest > 0
                    toplabel = hkfind(list, root)
                    bottomlabel = hkfind(list, dest)
                    if toplabel == bottomlabel && !(toplabel in counted)
                        sumsize += clusters[toplabel]
                        append!(counted, toplabel)
                    end
                end
            end
        end
    end
    return sumsize
end

"""
Find out whether percolation is possible or not in hkpercolate output
"""
@inline function hkbool(labels::Matrix{<:Integer}, list::Vector{<:Integer})
    for root in labels[1, :]
        if root > 0
            for dest in labels[end, :]
                if dest > 0 && hkfind(list, dest) == hkfind(list, root)
                    return true
                end
            end
        end
    end
    return false
end

end

max_stars_repo_name: slhshamloo/comp-phys

max_issues_repo_name: slhshamloo/comp-phys

max_forks_repo_name: slhshamloo/comp-phys


text:
@sde_model OrnsteinUhlenbeck dx = θ*(μ-x)*dt + σ*dW

function _ou_transition_distr(m, s, x0)
  @unpack θ, μ, σ = m
  m = x0 * exp(-θ * s.Δt) + μ * (1.0 - exp(-θ * s.Δt))
  d = sqrt(1.0 - exp(-2θ * s.Δt)) * σ / sqrt(2θ)
  Normal(m, d)
end

function sample(model::OrnsteinUhlenbeck, scheme::Exact, t0, x0)
  d = _ou_transition_distr(model, scheme, x0)
  rand(d)
end

function transition(model::OrnsteinUhlenbeck, scheme::Exact, t0, x0, x1)
  d = _ou_transition_distr(model, scheme, x0)
  pdf(d, x1)
end

function logtransition(model::OrnsteinUhlenbeck, scheme::Exact, t0, x0, x1)
  d = _ou_transition_distr(model, scheme, x0)
  logpdf(d, x1)
end


text:
module TestStoreApi

using ..MyPetStore
using Swagger
using Test
using Dates

function test(uri)
    @info("StoreApi")
    client = Swagger.Client(uri)
    api = StoreApi(client)

    @info("StoreApi - getInventory")
    inventory = getInventory(api)
    @test isa(inventory, Dict{String,Int32})
    @test !isempty(inventory)

    @info("StoreApi - placeOrder")
    @test_throws Swagger.ValidationException Order(; id=10, petId=10, quantity=2, shipDate=DateTime(2017, 03, 12), status="invalid_status", complete=false)
    order = Order(; id=10, petId=10, quantity=2, shipDate=DateTime(2017, 03, 12), status="placed", complete=false)
    neworder = placeOrder(api, order)
    @test neworder.id == 10

    @info("StoreApi - getOrderById")
    @test_throws Swagger.ValidationException getOrderById(api, 0)
    order = getOrderById(api, 10)
    @test isa(order, Order)
    @test order.id == 10

    @info("StoreApi - getOrderById (async)")
    response_channel = Channel{Order}(1)
    @test_throws Swagger.ValidationException getOrderById(api, response_channel, 0)
    @sync begin
        @async begin
            resp = getOrderById(api, response_channel, 10)
            @test (200 <= resp.status <= 206)
        end
        @async begin
            order = take!(response_channel)
            @test isa(order, Order)
            @test order.id == 10
        end
    end

    # a closed channel is equivalent of cancellation of the call,
    # no error should be thrown, but response can be nothing if call was interrupted immediately
    @test !isopen(response_channel)
    resp = getOrderById(api, response_channel, 10)
    @test (resp === nothing) || (200 <= resp.status <= 206)

    @info("StoreApi - deleteOrder")
    @test deleteOrder(api, 10) === nothing

    nothing
end

end # module TestStoreApi


text:
@compat abstract type abc_plan_type end

mutable struct abc_pmc_plan_saveable_type <: abc_plan_type
   gen_data::String
   calc_summary_stats::String
   calc_dist::String
   prior::Distribution
   make_proposal_dist::String
   normalize::String
   is_valid::String
   num_part::Int64
   num_max_attempt::Int64
   num_max_attempt_valid::Int64
   num_max_times::Int64
   num_dist_per_obs::Int64
   epsilon_init::Float64
   init_epsilon_quantile::Float64
   epsilon_reduction_factor::Float64
   target_epsilon::Float64
   tau_factor::Float64
   frac_rejects_save::Float64   # Not implemented yet
   param_active::Vector{Int64}   # Not fully tested yet
   adaptive_quantiles::Bool      # Not fully tested yet
   stop_on_decreasing_efficiency::Bool # Not implemented yet
   save_params::Bool   # Not implemented yet
   save_summary_stats::Bool   # Not implemented yet
   save_distances::Bool   # Not implemented yet
   in_parallel::Bool
end

mutable struct abc_pmc_plan_type <: abc_plan_type
   gen_data::Function
   calc_summary_stats::Function
   calc_dist::Function
   prior::Distribution
   make_proposal_dist::Function
   normalize::Function
   is_valid::Function
   num_part::Int64
   num_max_attempt::Int64
   num_max_attempt_valid::Int64
   num_max_times::Int64
   num_dist_per_obs::Int64
   epsilon_init::Float64
   init_epsilon_quantile::Float64
   epsilon_reduction_factor::Float64
   target_epsilon::Float64
   tau_factor::Float64
   frac_rejects_save::Float64   # Not implemented yet
   param_active::Vector{Int64}   # Not fully tested yet
   adaptive_quantiles::Bool      # Not fully tested yet
   stop_on_decreasing_efficiency::Bool # Not implemented yet
   save_params::Bool   # Not implemented yet
   save_summary_stats::Bool   # Not implemented yet
   save_distances::Bool   # Not implemented yet
   in_parallel::Bool

"""
   abc_pmc_plan_type(gd::Function,css::Function,cd::Function,p::Distribution; <keyword arguements>)

Create a plan for an ABC simulations.

# Arguements
- `gd`:  Generates the data
- `css`: Calculates the summary statistics from a simulated dataset
- `cd`:  Computes the distance between two sets of summary statistiscs
- `p`:   Prior distribution for model parameters
"""
function abc_pmc_plan_type(gd::Function,css::Function,cd::Function,p::Distribution;
     make_proposal_dist::Function = make_proposal_dist_gaussian_full_covar, param_active::Vector{Int64} = collect(1:length(Distributions.rand(p))),
     normalize::Function = noop, is_valid::Function = noop,
     num_part::Integer = 10*length(Distributions.rand(p))^2, num_max_attempt::Integer = 1000, num_max_attempt_valid::Integer = 50, num_max_times::Integer = 100, num_dist_per_obs::Int64 = 1,
     epsilon_init::Float64 = 1.0, init_epsilon_quantile::Float64 = 0.75, epsilon_reduction_factor::Float64 = 0.9,
     target_epsilon::Float64 = 0.01, tau_factor::Float64 = 2.0,
     adaptive_quantiles::Bool = false, stop_on_decreasing_efficiency::Bool = false,
     save_params::Bool = true, save_summary_stats::Bool = false, save_distances::Bool = true, frac_rejects_save::Float64 = 0.0, in_parallel::Bool = false)
     @assert(num_part>=length(Distributions.rand(p)))
     @assert(num_max_attempt>=1)
     @assert(num_max_attempt_valid>=1)
     @assert(num_max_times>0)
     @assert(1<=num_dist_per_obs<=10000)
     @assert(epsilon_init>0.01)
     @assert(0<init_epsilon_quantile<=1.0)
     @assert(0.5<epsilon_reduction_factor<1.0)
     @assert(target_epsilon>0.0)
     @assert(1.0<=tau_factor<=4.0)
     @assert(0.0<=frac_rejects_save<=1.0)
     new(gd,css,cd,p,make_proposal_dist,normalize,is_valid, num_part,num_max_attempt,num_max_attempt_valid,num_max_times,num_dist_per_obs,epsilon_init, init_epsilon_quantile,epsilon_reduction_factor,target_epsilon,tau_factor,frac_rejects_save,param_active,adaptive_quantiles,stop_on_decreasing_efficiency,save_params,save_summary_stats,save_distances,in_parallel)
   end
end

function abc_pmc_plan_saveable_type(plan::abc_pmc_plan_type)
  abc_pmc_plan_saveable_type(string(plan.gen_data),string(plan.calc_summary_stats),string(plan.calc_dist),plan.prior,string(plan.make_proposal_dist),string(plan.normalize),string(plan.is_valid),plan.num_part,plan.num_max_attempt,plan.num_max_attempt_valid,plan.num_max_times,plan.num_dist_per_obs,plan.epsilon_init,plan.init_epsilon_quantile,plan.epsilon_reduction_factor,plan.target_epsilon,plan.tau_factor,plan.frac_rejects_save,plan.param_active,plan.adaptive_quantiles,plan.stop_on_decreasing_efficiency,plan.save_params,plan.save_summary_stats,plan.save_distances,plan.in_parallel)
end

# WARNING: This doesn't really work yet.  There seems to be some issue with the scope of the saved function and what's currently imported
function abc_pmc_plan_type(plan::abc_pmc_plan_saveable_type)
  abc_pmc_plan_saveable_type(eval(Symbol(plan.gen_data)),eval(Symbol(plan.calc_summary_stats)),eval(Symbol(plan.calc_dist)),plan.prior,eval(Symbol(plan.make_proposal_dist)),eval(Symbol(plan.normalize)),eval(Symbol(plan.is_valid)),plan.num_part,plan.num_max_attempt,plan.num_max_attempt_valid,plan.num_max_times,plan.num_dist_per_obs,plan.epsilon_init,plan.init_epsilon_quantile,plan.epsilon_reduction_factor,plan.target_epsilon,plan.tau_factor,plan.frac_rejects_save,plan.param_active,plan.adaptive_quantiles,plan.stop_on_decreasing_efficiency,plan.save_params,plan.save_summary_stats,plan.save_distances,plan.in_parallel)
end

mutable struct abc_log_type   # Not implemented yet
   theta::Array{Array{Float64,1},1}
   ss::Array{Any,1}
   dist::Array{Float64,1}
   epsilon::Array{Float64,1}
   generation_starts_at::Array{Int64,1}

   function abc_log_type(t::Array{Array{Float64,1},1}, ss::Array{Any,1}, d::Array{Float64,1}, eps::Array{Float64,1}, gsa::Array{Int64,1} = ones(Int64,1))
      @assert length(t) == length(ss) == length(d) == length(eps)
      new( t, ss, d, eps, gsa )
   end
end

function abc_log_type()
   abc_log_type( Array{Array{Float64,1}}(undef,0), Array{Any}(undef,0), Array{Float64}(undef,0), Array{Float64}(undef,0), Array{Int64}(undef,0) )
end

mutable struct abc_population_type_old # Deprecated, but kept for now, in case need to read in data from a file
   theta::Array{Float64,2}
   weights::Array{Float64,1}
   dist::Array{Float64,1}
   log::abc_log_type
   repeats::Array{Int64,1}

   function abc_population_type_old(t::Array{Float64,2}, w::Array{Float64,1}, d::Array{Float64,1}, l::abc_log_type = abc_log_type(), repeats::Array{Int64,1} = zeros(Int64,length(w)) )
      @assert(length(w)==length(d)==size(t,2)==length(repeats))
      new(t,w,d,l,repeats)
   end
end

mutable struct abc_population_type
   theta::Array{Float64,2}
   weights::Array{Float64,1}
   dist::Array{Float64,1}
   epsilon::Float64
   logpdf::Array{Float64,1}
   accept_log::abc_log_type
   reject_log::abc_log_type
   repeats::Array{Int64,1}

   function abc_population_type(t::Array{Float64,2}, w::Array{Float64,1}, d::Array{Float64,1}, eps::Float64, lp::Array{Float64,1} = ones(Float64,length(w))/length(w), la::abc_log_type = abc_log_type(), lr::abc_log_type = abc_log_type(), repeats::Array{Int64,1} = zeros(Int64,length(w)) )
      @assert(length(w)==length(d)==length(lp)==size(t,2)==length(repeats))
      new(t,w,d,eps,lp,la,lr,repeats)
   end
end

"""
   abc_population_type(num_param::Integer, num_particles::Integer; accept_log::abc_log_type = abc_log_type(), reject_log::abc_log_type = abc_log_type(), repeats::Array{Int64,1} = zeros(Int64,num_particles) )

num_param:  Number of model parameters for generating simulated data
num_particles: Number of particles for sequential importance sampler

Optional parameters:
accept_log: Log of accepted parameters/summary statistics/distances
reject_log: Log of rejected parameters/summary statistics/distances
repeats:    Array indicating which particles have been repeated from previous generation

"""
function abc_population_type(num_param::Integer, num_particles::Integer; accept_log::abc_log_type = abc_log_type(), reject_log::abc_log_type = abc_log_type(), repeats::Array{Int64,1} = zeros(Int64,num_particles) )
   abc_population_type(zeros(num_param,num_particles), zeros(num_particles), zeros(num_particles), 0.0, ones(Float64,num_particles)/num_particles, accept_log, reject_log, repeats)
end


text:
module TestLibrary
include("preamble.jl")

@testset "Cat" begin
    # Inner transducer is stateful:
    @testset for xs in iterator_variants([1:5, 1:2])
        xf = DropWhile(x -> x < 3)
        @test collect(Cat() |> xf, xs) == collect(xf, vcat(xs...))
    end

    # Inner transducer is stateful and requires complete:
    @testset for xs in iterator_variants([1:5, 1:2])
        xf = PartitionBy(isequal(3)) |> Map(copy)
        @test collect(Cat() |> xf, xs) == collect(xf, vcat(xs...))
    end
end

@testset "MapCat" begin
    # https://clojuredocs.org/clojure.core/mapcat
    # (transduce (mapcat reverse) conj [[3 2 1 0] [6 5 4] [9 8 7]])
    @testset for xs in iterator_variants(
            [[3, 2, 1, 0], [6, 5, 4], [9, 8, 7]]
            )
        @test collect(MapCat(reverse), xs) == 0:9
    end
end

@testset "PartitionBy" begin
    # https://clojuredocs.org/clojure.core/partition-by#example-542692c7c026201cdc3269da
    xf = PartitionBy(isequal(3)) |> Map(copy)
    @testset for xs in iterator_variants(1:5)
        @test collect(xf, xs) == [[1, 2], [3], [4, 5]]
    end
end

@testset "Scan" begin
    @testset for xs in iterator_variants(1:10)
        @test collect(Scan(+), 1:10) == cumsum(1:10)
    end

    ed = eduction(Scan(+), 1:10)
    @test finaltype(ed.rf) === Int

    xs0 = [0, -1, 3, -2, 1]
    @testset for xs in [xs0, collect(xs0)]
    # @testset for xs in iterator_variants([0, -1, 3, -2, 1])  # TODO: fix
        @test collect(Scan(max), xs) == [0, 0, 3, 3, 3]
        @test collect(Scan(min), xs) == [0, -1, -1, -2, -2]
    end
end

@testset "TeeZip" begin
    @testset for xs in iterator_variants(1:5)
        @test collect(TeeZip(Filter(isodd) |> Map(inc)), xs) ==
            collect(zip(1:2:5, 2:2:6))
    end

    ed = eduction(TeeZip(Filter(isodd) |> Map(x -> x + 1)), 1:5)
    @test finaltype(ed.rf) === Tuple{Int,Int}

    xf = Map(inc) |> TeeZip(Filter(isodd)) |> Map(first)
    @testset for xs in iterator_variants(1:6)
        @test collect(xf, xs) == 3:2:7
    end
end

@testset "Replace" begin
    @testset for xs in iterator_variants(1:3)
        @test collect(Replace(Dict(1 => 10, 2 => 20)), xs) == [10, 20, 3]
    end
end

@testset "TakeWhile" begin
    @testset for xs in iterator_variants(1:5)
        @test collect(TakeWhile(x -> x < 3), xs) == 1:2
    end
end

@testset "TakeNth" begin
    @testset for xs in iterator_variants(1:10)
        @test collect(TakeNth(3), xs) == 1:3:10
    end
end

@testset "Drop" begin
    @testset for xs in iterator_variants(1:10)
        @test collect(Drop(3), xs) == 4:10
    end
end

@testset "DropLast" begin
    @testset for xs in iterator_variants(1:10)
        @test collect(DropLast(3), xs) == 1:7
    end
end

@testset "DropWhile" begin
    @testset for xs in iterator_variants(1:10)
        @test collect(DropWhile(x -> x % 3 != 0), xs) == 3:10
    end
end

# https://clojuredocs.org/clojure.core/keep
@testset "Keep" begin
    xf = Keep() do x
        if x % 3 == 0
            x
        end
    end
    @testset for xs in iterator_variants(1:10)
        @test collect(xf, xs) == 3:3:10
    end
end

# https://clojuredocs.org/clojure.core/distinct
@testset "Distinct" begin
    @testset for xs in iterator_variants([1, 1, 2, 1, 3, 2])
        @test collect(Distinct(), xs) == [1, 2, 3]
    end
end

# https://clojuredocs.org/clojure.core/interpose
@testset "Interpose" begin
    @testset for xs in iterator_variants(1:3)
        @test collect(Interpose(0), xs) == [1, 0, 2, 0, 3]
    end
end

# https://clojuredocs.org/clojure.core/dedupe
@testset "Dedupe" begin
    @testset for xs in iterator_variants([1, 1, 1, 2, 3, 3, 1])
        @test collect(Dedupe(), xs) == [1:3; 1]
    end
end

# https://clojuredocs.org/clojure.core/partition-all
@testset "Window" begin
    @testset for xs in iterator_variants(1:10)
        @test collect(Window(3) |> Map(copy), xs) == [[1:3;], [4:6;], [7:9;]]
    end
    @testset for xs in iterator_variants(1:6)
        @test collect(Window(3, 1) |> Map(copy), xs) ==
            [[1:3;], [2:4;], [3:5;], [4:6;]]
    end
    @testset "w=3, 7:$i" for i in 7:9
        @testset for xs in iterator_variants(1:i)
            @test collect(Window(3, flush=true) |> Map(copy), xs) ==
                [[1:3;], [4:6;], [7:i;]]
        end
    end
    @testset "w=4, 5:$i" for i in 5:8
        @testset for xs in iterator_variants(1:i)
            @test collect(Window(4, flush=true) |> Map(copy), xs) ==
                [[1:4;], [5:i;]]
        end
    end
    @testset "w=4, 9:$i" for i in 9:12
        @testset for xs in iterator_variants(1:i)
            @test collect(Window(4, flush=true) |> Map(copy), xs) ==
                [[1:4;], [5:8;], [9:i;]]
        end
    end
end

end  # module


text:
# Part 1
answers = split.(split(read("day06_input.txt", String), "\n\n"))
any_count(group) = length(∪(Set.(group)...))
println("Part 1: ∑(counts) = $(sum(any_count.(answers)))")

# Part 2
all_count(group) = length(∩(Set.(group)...))
println("Part 2: ∑(counts) = $(sum(all_count.(answers)))")

max_stars_repo_name: ericgreveson/adventofcode

max_issues_repo_name: ericgreveson/adventofcode

max_forks_repo_name: ericgreveson/adventofcode


text:
"""
exports @printf_code and printf_code, which returns
the code generated for printf.
"""
module Meta
export @printf_code, printf_code
import Base: @pure
import ..GenPrintf: printf_impl

"""
returns the code that printf will generate.
e.g. `@printf_code(fmt"abc")` will return
```
quote
    print("abc")
end
```

"""
macro printf_code(fmt, args...)
  quote
    flatten_body(
      remove_line_nodes(
        printf_impl(
          typeof($(esc(fmt))),
          $(length(args))
        )
      )
    )
  end
end

"""
returns the code that printf will generate.
e.g. `printf_code(fmt"abc")` will return
```
quote
    print("abc")
end
```
"""
@pure function printf_code(fmt, args...)
  flatten_body(
    remove_line_nodes(
      printf_impl(
        typeof(fmt),
        length(args)
      )
    )
  )
end

@pure remove_line_nodes(x) = x

@pure function remove_line_nodes(x::Expr)
  Expr(
    x.head,
    map(
      remove_line_nodes,
      filter(!x->x isa LineNumberNode, x.args)
    )...
  )
end

@pure flatten_body(x) = x

@pure function flatten_body(x::Expr)
  result_args = []

  for arg in x.args
    if arg isa Expr && arg.head == :block
      append!(result_args, flatten_body(arg).args)
     else
      push!(result_args, arg)
    end
  end

  Expr(:block, result_args...)
end

end


text:
module F1Method

#======================================================================
This package (the code below) implements the F-1 method as described
in the work of Pasquier et al. (2019). The numbers in parentheses
refer to the Equation numbers in the above manuscript. A bibtex
citation file is available in the GitHub repository.
======================================================================#

using LinearAlgebra, DualNumbers, HyperDualNumbers, DiffEqBase

"""
    Mem

Memory cache to store reusable objects.
Contains
- `s`   the steady-state solution, 𝒔(𝒑)
- `A`   the factors of 𝐀 = ∇ₓ𝑭(𝒔,𝒑)
- `∇s`  the derivative ∇𝒔(𝒑)
- `∇ₓf` the derivative ∇ₓ𝑓(𝒔,𝒑)
- `p`   the parameters 𝒑
The `Mem`-type object should be initialized with `initialize_mem`.
"""
mutable struct Mem
    s     # 𝒔(𝒑)
    A     # factors of 𝐀 = ∇ₓ𝑭(𝒔,𝒑)
    ∇s    # ∇𝒔(𝒑)
    ∇ₓf   # ∇ₓ𝑓(𝒔,𝒑)
    p     # 𝒑
end

function update_mem!(f, F, ∇ₓf, ∇ₓF, mem, p, alg; options...)
    if p ≠ mem.p                      # only update mem if 𝒑 has changed
        update_solution!(F, ∇ₓF, mem, p, alg; options...)
        s, m = mem.s.u, length(p)
        ∇ₚF = reduce(hcat, [𝔇(F(s, p + ε * e(j,m))) for j in 1:m]) # (2.7)
        mem.A = factorize(∇ₓF(s,p))   # update factors of ∇ₓ𝑭(𝒔,𝒑)
        mem.∇s .= mem.A \ -∇ₚF        # update ∇𝒔 (2.2)
        mem.∇ₓf .= ∇ₓf(s,p)           # update ∇ₓ𝑓(𝒔,𝒑)
        mem.p = p                     # update 𝒑
    end
end

function update_solution!(F, ∇ₓF, mem, p, alg; options...)
    if ~(mem.s isa SteadyStateSolution) || p ≠ mem.s.prob.p
        mem.s isa SteadyStateSolution ? x = mem.s.u : x = mem.s
        prob = SteadyStateProblem(F, ∇ₓF, x, p)       # define problem
        mem.s = solve(prob, alg; options...)          # update 𝒔
    end
end

"""
    objective(f, F, ∇ₓF, mem, p, alg; options...)

Returns `f(x,p)` such that `F(x,p)=0` using the F-1 method.

Specifically, `objective(f, F, ∇ₓF, mem, p, alg; options...)`
evaluates the objective function defined by `f̂(p) = f(s(p),p)`, where
`s(p)`, which is the steady-state solution (i.e., such that `F(s(p),p)=0`)
is computed by the iterative Newton-type solver `alg`.
The Jacobian, `∇ₓF`, and the memory cache `mem` must be supplied.
"""
function objective(f, F, ∇ₓF, mem, p, alg; options...)
    update_solution!(F, ∇ₓF, mem, p, alg; options...)
    return f(mem.s,p)
end

"""
    gradient(f, F, ∇ₓf, ∇ₓF, mem, p, alg; options...)

Returns the gradient of the `objective` function using the F-1 method.
"""
function gradient(f, F, ∇ₓf, ∇ₓF, mem, p, alg; options...)
    update_mem!(f, F, ∇ₓf, ∇ₓF, mem, p, alg; options...)
    s, ∇s, m = mem.s, mem.∇s, length(p)
    ∇ₚf = [𝔇(f(s,p + ε * e(j,m))) for j in 1:m]'    # (2.6)
    return mem.∇ₓf * ∇s + ∇ₚf                       # (2.1)
end

"""
    hessian(f, F, ∇ₓf, ∇ₓF, mem, p, alg; options...)

Returns the Hessian of the `objective` function using the F-1 method.
"""
function hessian(f, F, ∇ₓf, ∇ₓF, mem, p, alg; options...)
    update_mem!(f, F, ∇ₓf, ∇ₓF, mem, p, alg; options...)
    s, A, ∇s, m = mem.s, mem.A, mem.∇s, length(p)
    A⁻ᵀ∇ₓfᵀ = vec(A' \ mem.∇ₓf') # independent of (𝑗,𝑘)
    H, xⱼₖ = zeros(m,m), Vector{Hyper{Float64}}(undef, length(s))
    for j in 1:m, k in j:m       # loop upper triangle (symmetry)
        pⱼₖ = p + ε₁ * e(j,m) + ε₂ * e(k,m)              # hyperdual 𝒑
        @views xⱼₖ .= s + ε₁ * ∇s[:,j] + ε₂ * ∇s[:,k]    # hyperdual 𝒙
        H[j,k] = ℌ(f(xⱼₖ,pⱼₖ)) - ℌ(F(xⱼₖ,pⱼₖ))' * A⁻ᵀ∇ₓfᵀ    # (2.8)
        j ≠ k ? H[k,j] = H[j,k] : nothing # Hessian symmetry
    end
    return H
end

"""
    initialize_mem(x, p)

Initializes the memory cache for the F-1 method.
"""
function initialize_mem(x, p)
    n, m = length(x), length(p)
    return Mem(copy(x), nothing, zeros(n,m), zeros(1,n), nothing)
end

e(j, m) = [i == j for i in 1:m]      # 𝑗th basis vector of ℝᵐ
𝔇(x) = DualNumbers.dualpart.(x)      # dual part
ℌ(x) = HyperDualNumbers.ε₁ε₂part.(x) # hyperdual part

end


text:
module SpeckleTest

using Speckles
using Test

################################################################################
function TestSpeckles()
    @testset "Speckles.jl" begin
        # Write your tests here.
    end
end

export TestSpeckles

################################################################################
function TestLightSource()
    @testset "LightSource.jl" begin
        include("testLightSource.jl")
    end
end

export TestLightSource

################################################################################
function TestDetector()
    @testset "Detector.jl" begin
        include("testDetector.jl")
    end
end

export TestDetector

################################################################################
function TestAll()
    TestSpeckles()
    TestLightSource()
    TestDetector()
end

export TestAll

end

import .SpeckleTest

SpeckleTest.TestDetector()


text:
using Documenter
using KrylovKit

makedocs(modules=[KrylovKit],
            format=:html,
            sitename="KrylovKit.jl",
            pages = [
                "Home" => "index.md",
                "Manual" => [
                    "man/intro.md",
                    "man/linear.md",
                    "man/eig.md",
                    "man/svd.md",
                    "man/matfun.md",
                    "man/algorithms.md",
                    "man/implementation.md"
                ]
            ])

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    deps = nothing,
    make = nothing,
    target = "build",
    repo = "github.com/Jutho/KrylovKit.jl.git",
    julia = "1.0"
)


text:
using LinearAlgebra
using Random
using Test
using TrackedDistributions
using Tracker

include("distributions.jl")


text:
# This file is part of RLinearAlgebra.jl

module ProceduralTestLinSysSampler

using Test, RLinearAlgebra, Random

@testset "Linear System Sampler Abstractions" begin

    @test LinSysSketch == LinSysSampler #Check alias
    @test LinSysSelect == LinSysSampler #Check alias

    @test supertype(LinSysVecRowSampler) == LinSysSampler #Check parent type
    @test LinSysVecRowSketch == LinSysVecRowSampler #Check alias
    @test LinSysVecRowSelect == LinSysVecRowSampler #Check alias

    @test supertype(LinSysVecColSampler) == LinSysSampler #Check parent type
    @test LinSysVecColSketch == LinSysVecColSampler #Check alias
    @test LinSysVecColSelect == LinSysVecColSampler #Check alias

    @test supertype(LinSysBlkRowSampler) == LinSysSampler #Check parent type
    @test LinSysBlkRowSketch == LinSysBlkRowSampler #Check alias
    @test LinSysBlkRowSelect == LinSysBlkRowSampler #Check alias

    @test supertype(LinSysBlkColSampler) == LinSysSampler #Check parent type
    @test LinSysBlkColSketch == LinSysBlkColSampler #Check alias
    @test LinSysBlkColSelect == LinSysBlkColSampler #Check alias
end

end


text:
using GraphCircuits
if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end
using LightGraphs

@testset "Cycles" begin
    hadwick_james_ex1 = DiGraph(5)
    add_edge!(hadwick_james_ex1, 1, 2)
    add_edge!(hadwick_james_ex1, 2, 3)
    add_edge!(hadwick_james_ex1, 3, 4)
    add_edge!(hadwick_james_ex1, 3, 5)
    add_edge!(hadwick_james_ex1, 4, 5)
    add_edge!(hadwick_james_ex1, 5, 2)

    # johnson_circuits = find_circuits(hadwick_james_ex1, Val{:johnson})
    hadwick_circuits = find_circuits(hadwick_james_ex1, Val{:hadwick_james})

    expected_circuits = Vector{Int}[
        [2, 3, 4, 5],
        [2, 3, 5]
    ]
    @test issubset(expected_circuits, hadwick_circuits)
    @test issubset(hadwick_circuits, expected_circuits)

    add_edge!(hadwick_james_ex1, 1, 1)
    add_edge!(hadwick_james_ex1, 3, 3)

    hadwick_circuits_self = find_circuits(hadwick_james_ex1, Val{:hadwick_james})

    @test issubset(expected_circuits, hadwick_circuits_self)
    @test [1] ∈ hadwick_circuits_self && [3] ∈ hadwick_circuits_self

    ex2_size = 10
    hadwick_james_ex2 = PathDiGraph(ex2_size)
    @test isempty(find_circuits(hadwick_james_ex2, Val{:hadwick_james}))

    ex3_size = 5
    hadwick_james_ex3 = CompleteDiGraph(ex3_size)
    @test length(find_circuits(hadwick_james_ex3, Val{:hadwick_james})) == length(unique(find_circuits(hadwick_james_ex3, Val{:hadwick_james})))

    # Almost fully connected DiGraph
    ex4 = DiGraph(9)
    for (src, dest) in [(1,2),(1,5),(1,7),(1,8),(2,9),(3,4),(3,6),
                        (4,5),(4,7),(5,6),(6,7),(6,8),(7,9),(8,9)]
        add_edge!(ex4, src, dest)
        add_edge!(ex4, dest, src)
    end
    ex4_output = simplecycles_hadwick_james(ex4)
    @test [1, 2] ∈ ex4_output && [8, 9] ∈ ex4_output
end


text:
function  nll_data(m::AbstractModel,df)
    # -sum(v*Log_LikeliHood(m, k...) for (k, v) in cm)
    -sum(log.(Likelihood(m,df)))
end

function check_params(m::AbstractModel)
    (all(params(m) .> 1e-5) && (all([(v < l) for (v,l) in zip(params(m), m.limits)])))
end

function Distributions.fit(m::WaltonForaging.AbstractModel,df::DataFrames.AbstractDataFrame)
    res = optimize(params(m), NelderMead(), Optim.Options(iterations = 5000)) do param
        (all(param .> 1e-5) && (all([(v < l) for (v,l) in zip(param, m.limits)]))) || return 1e10
        updateparams!(m,param)
        return nll_data(m, df)
    end
    updateparams!(m,Optim.minimizer(res))
    return m, res
end

function Distributions.fit(::Type{T},df::DataFrames.AbstractDataFrame) where T <:WaltonForaging.AbstractModel
    m = init(T)
    Distributions.fit(m,df)
end


text:
using GeCKOReader
using Base.Test
using Bio.Seq
using Query
using HypothesisTests
using DataFrames

@testset "Hamming" begin

    # 1 error
    seq_a = dna"ATCGT"
    seq_b = dna"ATAGT"
    @test GeCKOReader.hamming(seq_a, seq_b) == 1

    # 0 error
    seq_a = dna"ATCGT"
    seq_b = dna"ATCGT"
    @test GeCKOReader.hamming(seq_a, seq_b) == 0

    # All error
    seq_a = dna"AAAAA"
    seq_b = dna"TTTTT"
    @test GeCKOReader.hamming(seq_a, seq_b) == length(seq_a)

    # Empty sequence test
    seq_a = dna""
    seq_b = dna""
    @test GeCKOReader.hamming(seq_a, seq_b) == 0

    # Comparing DNA sequences of different lengths is undefined
    seq_a = dna"CGA"
    seq_b = dna"CGAA"
    @test_throws DimensionMismatch GeCKOReader.hamming(seq_a, seq_b) == 0

end

@testset "Frequencies" begin
    # generate fake sample
    testsample = ScreenSample("test", 51, 41)
    testsample.counts["HGLibA_64614"] = 20
    testsample.counts["HGLibA_64466"] = 30
    testsample.counts["HGLibA_65230"] = 25
    testsample.counts["HGLibA_54150"] = 50

    callable_bc = readtable("callable_sgrnas.csv")

    freqs = GeCKOReader.get_rel_freqs(testsample, callable_bc)

    trimmed_freqs = @from i in freqs begin
        @where i.freqs > 0
        @select i
        @collect DataFrame
    end

    @test sum(trimmed_freqs[:freqs]) == 1

    @test all(isapprox.(Array(trimmed_freqs[:rel_freqs]), [2.0, 1.2, 0.8, 1.0]))
end

@testset "Log2 FC" begin
    df = DataFrame(gene=["A", "B", "A", "A", "B", "A", "B", "C", "C"],
                   log2fc=[5.0, 0.375, 3.0, 6.0, 1.5, Inf, -Inf, Inf, Inf])
    negcontrols = [2.58155,2.29046,1.66488,-0.149982,1.67709]

    result = GeCKOReader.compute_gene_stats(df, negcontrols)

    @test size(result) == (2, 5)

    @test result[1, :med_log2fc] == 5.0
    @test result[2, :med_log2fc] == 0.9375

    pval = -log10(pvalue(MannWhitneyUTest([5.0, 3.0, 6.0], negcontrols)))

    @test result[1, :pval] == pval
end

@testset "Line fit" begin
    t0 = DataFrame(id=["1", "2", "3", "4", "5"],
                   gene=["A", "A", "A", "B", "B"],
                   rel_freqs=[0.5, 0.25, 1, 0.1, 0.25],
                   counts=[10, 5, 20, 2, 5])
    t1 = DataFrame(id=["1", "2", "3", "4", "5"],
                   gene=["A", "A", "A", "B", "B"],
                   rel_freqs=[0.625, 0.875, 2.125, 0.25, 0.375],
                   counts=[5, 7, 17, 2, 3])
    t2 = DataFrame(id=["1", "2", "3", "4", "5"],
                   gene=["A", "A", "A", "B", "B"],
                   rel_freqs=[0.04, 0.24, 0.8, 0.04, 0.4],
                   counts=[1, 6, 20, 1, 10])

    output = GeCKOReader.fit_lines(t0, t1, t2)

    @test all(isapprox.(output[:intercept], [0.618333, 0.46, 1.40833, 0.16, 0.26667], atol=0.0001))
    @test all(isapprox.(output[:slope], [-0.23, -0.005, -0.1, -0.03, 0.075], atol=0.0001))
end


text:
abstract type Excitation end

struct VoltageSource <: Excitation
    w::Wire
    m::Int      # segment number
    v::Complex{Float64}
end

struct LinearWave <: Excitation
    n_thetas::Int       # num theta angles for incident plane wave
    n_phis::Int         # num phi angles for incident plane wave

    theta::Float64      # theta in degrees
    phi::Float64        # phi in degrees
    eta::Float64        # eta in degrees

    d_theta::Float64    # theta angle stepping increment in degrees
    d_phi::Float64      # phi angle stepping increment in degrees
end


text:
using CUDA.CUSPARSE, SparseArrays

@testset "sparse" begin
    n, m = 4, 4
    I = [1,2,3] |> cu
    J = [2,3,4] |> cu
    V = Float32[1,2,3] |> cu

    dense = rand(3,3) |> cu

    # check defaults
    @test sparse(I, J, V) isa CuSparseMatrixCSC
    @test sparse(dense) isa CuSparseMatrixCSC

    for (fmt, T) in  [(:coo, CuSparseMatrixCOO),
                      (:csc, CuSparseMatrixCSC),
                      (:csr, CuSparseMatrixCSR),
                      (:bsr, CuSparseMatrixBSR)
                     ]
        if fmt != :bsr # bsr not supported
            x = sparse(I, J, V; fmt=fmt)
            @test x isa T{Float32}
            @test size(x) == (3, 4)

            x = sparse(I, J, V, m, n; fmt=fmt)
            @test x isa T{Float32}
            @test  size(x) == (4, 4)
        end

        if fmt != :coo # dense to COO not implemented
            x = sparse(dense; fmt=fmt)
            @test x isa T{Float32}
            @test collect(x) == collect(dense)
        end
    end
end


text:
using CpuId, LLVM

features = split(unsafe_string(LLVM.API.LLVMGetHostCPUFeatures()), ',')
avx512f = any(isequal("+avx512f"), features)

register_size = avx512f ? 64 : 32
register_count = avx512f ? 32 : 16
cache_size = CpuId.cachesize()
num_cores = CpuId.cpucores()

avx2 = any(isequal("+avx2"), features)

# Should I just add all the flags in features?
cpu_info_string = """
const REGISTER_SIZE = $register_size
const REGISTER_COUNT = $register_count
const FP256 = $(cpufeature(CpuId.FP256)) # Is AVX2 fast?
const CACHELINE_SIZE = $(cachelinesize())
const CACHE_SIZE = $cache_size
const NUM_CORES = $num_cores
const FMA3 = $(any(isequal("+fma"), features))
const AVX2 = $(avx2)
const AVX512F = $(avx512f)
const AVX512ER = $(any(isequal("+avx512er"), features))
const AVX512PF = $(any(isequal("+avx512pf"), features))
const AVX512VL = $(any(isequal("+avx512vl"), features))
const AVX512BW = $(any(isequal("+avx512bw"), features))
const AVX512DQ = $(any(isequal("+avx512dq"), features))
const AVX512CD = $(any(isequal("+avx512cd"), features))
const SIMD_NATIVE_INTEGERS = $(avx2)

"""


text:
@testset "FixedTimeZone" begin
    utc = TimeZone("UTC")
    fixed = FixedTimeZone("UTC-06:00")
    dt = DateTime(2016, 8, 11, 2, 30)

    for t in (utc, fixed)
        @test ZonedDateTime(LaxZonedDateTime(dt, t)) == ZonedDateTime(dt, t)
        for p in (Hour(1), Day(1))
            @test LaxZonedDateTime(dt, t) + p == LaxZonedDateTime(dt + p, t)
            @test LaxZonedDateTime(dt, t) - p == LaxZonedDateTime(dt - p, t)
        end
    end
end

@testset "VariableTimeZone" begin
    t = VariableTimeZone("Testing", [
        Transition(DateTime(1800,1,1), FixedTimeZone("TST",0,0)),
        Transition(DateTime(1960,4,1,2), FixedTimeZone("TDT",0,7200)),   # 1960-04-01 02:00
        Transition(DateTime(1960,4,1,3), FixedTimeZone("TET",0,10800)),  # 1960-04-01 05:00
        Transition(DateTime(1960,4,1,10), FixedTimeZone("TDT",0,7200)),  # 1960-04-02 12:00
    ])

    @test LaxZonedDateTime(DateTime(1960,3,31,3), t) + Day(1) == LaxZonedDateTime(DateTime(1960,4,1,3), t, NonExistent())

    valid_a = LaxZonedDateTime(DateTime(1960,4,1,1), t)
    valid_b = LaxZonedDateTime(DateTime(1960,4,1,6), t)
    non_existent_a = LaxZonedDateTime(DateTime(1960,4,1,2), t, NonExistent())
    non_existent_b = LaxZonedDateTime(DateTime(1960,4,1,3), t, NonExistent())
    null = LaxZonedDateTime()
    ambiguous = LaxZonedDateTime(DateTime(1960,4,1,12),t)

    @test valid_a + Hour(2) == valid_b
    @test isequal(non_existent_a + Hour(1), null)
    @test isequal(non_existent_b + Hour(1), null)

    @test isequal(non_existent_a + Hour(0), null)
    @test isequal(non_existent_b + Hour(0), null)
    @test isequal(ambiguous + Hour(0), null)

    @test isequal(ambiguous - Hour(1), null)
    @test isequal(ambiguous + Hour(1), null)
end

@testset "non-existent" begin
    lzdt = LaxZonedDateTime(ZonedDateTime(2015,3,7,2, winnipeg))
    lzdt += Day(1)
    @test lzdt == LaxZonedDateTime(DateTime(2015,3,8,2), winnipeg, NonExistent())
    lzdt += Day(1)
    @test lzdt == LaxZonedDateTime(ZonedDateTime(2015,3,9,2,winnipeg))

    non_existent = LaxZonedDateTime(DateTime(2015,3,8,2), winnipeg, NonExistent())
    @test isequal(non_existent + Hour(1), LaxZonedDateTime())
    @test isequal(non_existent - Hour(1), LaxZonedDateTime())
    @test isequal(non_existent + Hour(0), LaxZonedDateTime())
    @test isequal(non_existent - Hour(0), LaxZonedDateTime())

    @test isequal((non_existent + Hour(1)) + Hour(1), LaxZonedDateTime())
    @test isequal((non_existent + Hour(1)) + Day(1), LaxZonedDateTime())
end

@testset "compare" begin

    @testset "basic" begin
        a = LaxZonedDateTime(ZonedDateTime(2014,winnipeg))
        b = ZonedDateTime(2015,winnipeg)
        @test a < b

        a = ZonedDateTime(2016, 11, 6, 1, 30, winnipeg, 1)
        b = ZonedDateTime(2016, 11, 6, 1, winnipeg, 2)

        @test LaxZonedDateTime(a) < b
    end

    @testset "ambiguous" begin
        dt = DateTime(2016, 11, 6, 1, 45)
        amb1 = LaxZonedDateTime(ZonedDateTime(dt, tz"America/Winnipeg", 1))
        amb2 = LaxZonedDateTime(ZonedDateTime(dt, tz"America/Winnipeg", 2))
        amb = LaxZonedDateTime(dt, tz"America/Winnipeg")

        # Note: the fuzzy matching allows ranges to work correctly. Possibly we should
        # use a different operator for this.
        @test amb1 != amb2
        @test amb1 != amb
        @test amb2 != amb

        @test isless(amb1, amb2)
        @test !isless(amb1, amb)
        @test !isless(amb2, amb)

        @test amb1 < amb2
        @test !(amb1 < amb)
        @test !(amb2 < amb)

        @test amb1 <= amb2
        @test amb1 <= amb  # fuzzy
        @test amb2 <= amb  # fuzzy

        @test !(amb1 > amb2)
        @test !(amb1 > amb)
        @test !(amb2 > amb)

        @test !(amb1 >= amb2)
        @test amb1 >= amb  # fuzzy
        @test amb2 >= amb  # fuzzy
    end

    @testset "null equality" begin
        # Unrepresentable LZDTs are treated like NaNs
        null = LaxZonedDateTime()
        @test null != null
        @test isequal(null, null)

        # Test that unrepresentable LZDTs initialized with different values hash the same
        utc = TimeZone("UTC")
        null_2 = LaxZonedDateTime(DateTime(2013, 2, 13, 0, 30), utc, utc, false)
        @test isequal(null, null_2)
        @test hash(null) == hash(null_2)
    end

    @testset "ZonedDateTime" begin
        zdt = ZonedDateTime(DateTime(2016, 11, 10, 1, 45), winnipeg)
        lzdt = LaxZonedDateTime(DateTime(2016, 11, 10, 1, 45), winnipeg)
        dne = LaxZonedDateTime(DateTime(2015, 3, 8, 2), winnipeg)
        nonamb = LaxZonedDateTime(ZonedDateTime(2017, 11, 6, 1, 45, winnipeg, 1))
        amb1 = LaxZonedDateTime(DateTime(2016, 11, 6, 1, 45), winnipeg)
        amb2 = LaxZonedDateTime(DateTime(2017, 11, 5, 1, 45), winnipeg)

        @testset "isequal" begin
            @test isequal(zdt, lzdt)
            @test isequal(lzdt, zdt)
            @test !isequal(zdt, dne)
            @test isequal(dne, dne)
            @test !isequal(amb1, amb2)
            @test isequal(amb1 + Hour(1), amb2 + Hour(1))
            @test !isequal(amb1, nonamb)
        end

        @testset "hash" begin
            @test hash(zdt) == hash(lzdt)
            @test hash(zdt) != hash(dne)
            @test hash(dne) == hash(dne)
            @test hash(amb1) != hash(amb2)
            @test hash(amb1 + Hour(1)) == hash(amb2 + Hour(1))
        end
    end
end

@testset "convert to ZonedDateTime" begin

    @testset "basic" begin
        lzdt = LaxZonedDateTime(DateTime(2020, 1, 1, 9, 30), winnipeg)
        expected =  ZonedDateTime(DateTime(2020, 1, 1, 9, 30), winnipeg)
        @test ZonedDateTime(lzdt) == expected
    end

    @testset "Ambiguous" begin
        # Occurs during Fall DST - it's the same hour just in two different timezones
        amb = LaxZonedDateTime(DateTime(2016, 11, 6, 1, 45), winnipeg)

        zdt = ZonedDateTime(amb, :first)
        exp_first = ZonedDateTime(DateTime(2016, 11, 6, 1, 45), winnipeg, 1)
        @test zdt == exp_first

        zdt = ZonedDateTime(amb, :last)
        exp_last = ZonedDateTime(DateTime(2016, 11, 6, 1, 45), winnipeg, 2)
        @test zdt == exp_last

        @test_throws AmbiguousTimeError ZonedDateTime(amb, :throw)
        @test_throws AmbiguousTimeError ZonedDateTime(amb)

    end

    @testset "NonExistent" begin
        # Occurs during Spring DST
        dne = LaxZonedDateTime(DateTime(2016, 3, 13, 2, 45), winnipeg)

        # :first gets the time just before the transition occurred
        zdt = ZonedDateTime(dne, :first)
        exp_first = ZonedDateTime(2016, 3, 13, 1, 59, 59, 999, winnipeg)
        @test zdt == exp_first

        # :last gets the hour after the transition occurred
        zdt = ZonedDateTime(dne, :last)
        exp_last = ZonedDateTime(2016, 3, 13, 3, 0, winnipeg)
        @test zdt == exp_last

        @test_throws NonExistentTimeError ZonedDateTime(dne, :throw)
        @test_throws NonExistentTimeError ZonedDateTime(dne)
    end

end


text:
using DrWatson
@quickactivate "GCN_HM_GRN-Integration"
using Turing, PyCall
using DataFrames
using GraphRecipes
using Plots
using LightGraphs
using Statistics
using Distributions
using Zygote
theme(:default)
gr(fmt=:svg)

Turing.setadbackend(:tracker)

Zygote.@adjoint function pycall(f, x...; kw...)
    x = map(py, x)
    y = pycall(f, x...; kw...)
    y.detach().numpy(), function (ȳ)
        y.backward(gradient = py(ȳ))
        (nothing, map(x->x.grad.numpy(), x)...)
    end
end

abstract type BayesianInferenceModel end

function load_pyg(path::String, file::String)
    pushfirst!(PyVector(pyimport("sys")."path"), path);
    train_model = pyimport(file).train_model
    run_model = pyimport(file).run_model
    extract_subgraph = pyimport(file).extract_subgraph
    pyimport("sys")."stdout" = PyTextIO(stdout)
    pyimport("sys")."stderr" = PyTextIO(stderr)
    return train_model, run_model, extract_subgraph
end

struct PyGInferenceModel <: BayesianInferenceModel
    C::Int32
    path::String
    file::String
    node::Int32
    k::Int32
    ei::Array{Int32, 2}
    base_pred::Array{Float32, 1}
    model
    x
    y
end

function PyGInferenceModel(path::String, file::String, node::Integer, k::Integer, α::Float64 = 2.0, β::Float64 = 20.0, stdev::Float64 = 0.0001)
    node = convert(Int32, node)
    k = convert(Int32, k)
    train_model, run_model, extract_subgraph = load_pyg(path, file)
    model, x, y, edge_index = train_model();
    edge_index, node_idx = extract_subgraph(node, k, edge_index)
    edge_mask = repeat([1], get(edge_index."shape", 1));
    @model graph_mask(y) = begin
        N = get(edge_index."shape", 1)
        edge_mask = Vector{Int8}(undef, N)
        edge_prob ~ filldist(Beta(α, β), N)
        edge_mask .~ Bernoulli.(edge_prob)
        pred = run_model(edge_mask, edge_index, model, node_idx)
        y = y/sum(y)
        y ~ MvNormal(pred/sum(pred), repeat([stdev], length(y)))
    end
    base = run_model(repeat([1], get(edge_index."shape", 1)), edge_index, model, node_idx)
    return PyGInferenceModel(get(edge_index."shape", 1), path, file, node, k,
            edge_index."detach"()."cpu"()."numpy"(), base, graph_mask, x."detach"()."cpu"()."numpy"(), y."detach"()."cpu"()."numpy"())
end

function sample(m::T where T <: BayesianInferenceModel, iters::Integer)
    return Turing.sample(m.model(m.base_pred), HMC(3.5, 5, :edge_prob), iters)
end

function final_summary(s::Chains)
    return DataFrame(describe(s)[1])
end

function plot_result(model::T where T <: BayesianInferenceModel, result::Chains, k::Number, dlpy::Bool = true)
    ei = model.ei;
    vals = DataFrame(describe(result)[1])[size(ei, 2)+1:end, 2];
    sub_idx = (mean(vals) + k * var(vals)) .< vals;
    vals = vals[sub_idx]
    ei = ei[:, sub_idx]

    weigh = Dict((ei[1, i], ei[2, i]) => vals[i] for i ∈ 1:length(vals));
    uni = unique(ei);
    conv = Dict(uni[i] => i for i ∈ 1:length(uni));
    conv_rev = Dict(i => uni[i] for i ∈ 1:length(uni));

    ei_conv = map(x -> conv[x], ei)
    inp = ei_conv
    classes = model.y[unique(map(x -> conv_rev[x], inp)) .+ 1]

    g = DiGraph(Edge.(zip(inp[1, :], inp[2, :])))

    plt = plot(graphplot(g, names=uni, arrow=true, nodecolor=classes, edgewidth=(s, d, w)->
        2*weigh[(conv_rev[s],conv_rev[d])], fontsize=8), size=(1000, 1000), dpi=300)
    if dlpy
        display(plt)
    end
    return plt
end


text:
include("task_func.jl")

function draw(r)
    putmarker!(r)
    for side = (Nord, Sud)
        listOfShags = []
        k = 1
        k2, k3, lastk2, lastk3 = 0, 0, 0, 0
        while !isborder(r, side)
            push!(listOfShags, goToSideWithShag!(r, side, 1, save=true))

            lastk2 = k2
            k2 = goToSideWithShag!(r, West, k, save = true, checkBorder = true)[1]
            if lastk2!=k2 putmarker!(r) end

            lastk3 = k3
            k3 = goToSideWithShag!(r, Ost, k+k2, save = true, checkBorder = true)[1]
            if lastk3!=k3 putmarker!(r) end

            goToSideWithShag!(r, West, k3-k2, save = false, checkBorder = false)
            k+=1
        end
        goToStart!(r, listOfShags)
    end
end

max_stars_repo_name: Trilonka/MireaAl

max_issues_repo_name: Trilonka/MireaAl

max_forks_repo_name: Trilonka/MireaAl


text:
# refer to benchmarks/array-board/run-benchmarks.jl
using Game2048Core
using Game2048Core: bitboard_to_array, move!

function ok()
    for _ in 1:1_000_000
        bboard = Bitboard(rand(UInt64))
        arr_board = bitboard_to_array(bboard)
        if any(arr_board .== 15)
            break
        end

        new_board = move(bboard, left)

        move!(arr_board, left)

        if  !(bitboard_to_array(new_board) == arr_board)
            return bboard
        end
    end
end

@time a = ok()

move(a, left)
b = bitboard_to_array(a)
move!(b, left)

b


text:
for p in ("ArgParse","JLD","Knet")
    if Pkg.installed(p) == nothing; Pkg.add(p); end
end

module RNNLM
using ArgParse,JLD,Knet
using AutoGrad: getval
using Knet: sigm_dot, tanh_dot
logprint(x)=join(STDERR,[Dates.format(now(),"HH:MM:SS"),x,'\n'],' ')
macro run(i,x) :(if loglevel>=$i; $(esc(x)); end) end
macro msg(i,x) :(if loglevel>=$i; logprint($(esc(x))); end) end
macro log(i,x) :(if loglevel>=$i; logprint($(string(x))); end; $(esc(x))) end

function lstm(weight,bias,hidden,cell,input)            # 2:991  1:992:1617 (id:forw:back)
    gates   = weight * hidden .+ input .+ bias          # 2:312  1:434:499 (43+381+75) (cat+mmul+badd)
    h       = size(hidden,1)                            #
    forget  = sigm_dot(gates[1:h,:])                        # 2:134  1:98:99  (62+37) (index+sigm)
    ingate  = sigm_dot(gates[1+h:2h,:])                     # 2:99   1:73:123 (77+46)
    outgate = sigm_dot(gates[1+2h:3h,:])                    # 2:113  1:66:124 (87+37)
    change  = tanh_dot(gates[1+3h:4h,:])                    # 2:94   1:51:179 (130+49) replace end with 4h?
    cell    = cell .* forget + ingate .* change         # 2:137  1:106:202 (104+93+5) (bmul+bmul+add)
    hidden  = outgate .* tanh_dot(cell)                     # 2:100  1:69:194 (73+121) (tanh+bmul)
    return (hidden,cell)
end

# sequence[t]::Vector{Int} minibatch of tokens
function rnnlm(model, state, sequence; pdrop=0, range=1:(length(sequence)-1), keepstate=nothing, stats=nothing)
    index = vcat(sequence[range]...)
    input = Wm(model)[:,index]                          # 2:15
    for n = 1:nlayers(model)
        input = dropout(input, pdrop)
        input = Wx(model,n) * input                     # 2:26
        w,b,h,c = Wh(model,n),bh(model,n),hdd(state,n),cll(state,n)
        output = []
        j1 = j2 = 0
        for t in range
            j1 = j2 + 1
            j2 = j1 + length(sequence[t]) - 1
            input_t = input[:,j1:j2]                    # 2:35
            (h,c) = lstm(w,b,h,c,input_t)               # 2:991
            push!(output,h)
        end
        if keepstate != nothing
            keepstate[2n-1] = getval(h)
            keepstate[2n] = getval(c)
        end
        input = hcat(output...)                         # 2:39
    end
    input = dropout(input,pdrop)
    input = Wy(model) * input .+ by(model)
    input = logp(input,1)                               # 2:354  1:1067:673
    golds = vcat(sequence[range+1]...)
    index = golds + size(input,1)*(0:(length(golds)-1))
    logp1 = input[index]
    total = sum(logp1)
    count = length(logp1)
    if stats != nothing
        stats[1]=total
        stats[2]=count
    end
    batch = length(sequence[1])
    # return -total / count # per token loss: scale does not depend on sequence length or minibatch
    return -total / batch # per sequence loss: does not depend on minibatch, larger loss for longer seq
    # return -total 	    # total loss: longer sequences and larger minibatches have higher loss
end

rnnlmgrad = grad(rnnlm)

# data[t][b] contains word[(b-1)*T+t]
function bptt(model, data, optim; pdrop=0, slen=20)
    T = length(data)
    B = length(data[1])
    state = initstate(model,B)
    @run 2 begin
        wnorm = zeros(length(model))
        gnorm = zeros(length(model))
        count = 0
    end
    for i = 1:slen:(T-1)
        j = i+slen-1
        if j >= T; break; end
        grads = rnnlmgrad(model, state, data; pdrop=pdrop, range=i:j, keepstate=state)
        update!(model, grads, optim)
        @run 2 begin
            gnorm += map(vecnorm,grads)
            wnorm += map(vecnorm,model)
            count += 1
        end
    end
    @msg 2 string("wnorm=",wnorm./count)
    @msg 2 string("gnorm=",gnorm./count)
end

function loss(model, data; slen=20)
    T = length(data)
    B = length(data[1])
    state = initstate(model,B)
    rat=zeros(2); tot=zeros(2)
    for i = 1:slen:(T-1)
        j = i+slen-1
        if j >= T; break; end
        rnnlm(model, state, data; stats=rat, range=i:j, keepstate=state)
        tot += rat
    end
    return (-tot[1],tot[2])
end

nlayers(model)=div(length(model)-3,3)
Wm(model)=model[1]
Wx(model,n)=model[3n-1]
Wh(model,n)=model[3n]
bh(model,n)=model[3n+1]
Wy(model)=model[end-1]
by(model)=model[end]
hdd(state,n)=state[2n-1]
cll(state,n)=state[2n]

function initmodel(atype, hidden, vocab, embed)
    init(d...)=atype(xavier(Float32,d...))
    bias(d...)=atype(zeros(Float32,d...))
    N = length(hidden)
    model = Array{Any}(3N+3)
    model[1] = init(embed,vocab) # Wm
    X = embed
    for n = 1:N
        H = hidden[n]
        model[3n-1] = init(4H,X) # Wx
        model[3n]   = init(4H,H) # Wh
        model[3n+1] = bias(4H,1) # bh
        model[3n+1][1:H] = 1     # forget gate bias = 1
        X = H
    end
    model[3N+2] = init(vocab,hidden[end]) # Wy
    model[3N+3] = bias(vocab,1)           # by
    return model
end

let blank = nothing; global initstate
function initstate(model, batch)
    N = nlayers(model)
    state = Array{Any}(2N)
    for n = 1:N
        bias = bh(model,n)
        hidden = div(length(bias),4)
        if typeof(blank)!=typeof(bias) || size(blank)!=(hidden,batch)
            blank = fill!(similar(bias, hidden, batch),0)
        end
        state[2n-1] = state[2n] = blank
    end
    return state
end
end

# initoptim creates optimization parameters for each numeric weight
# array in the model.  This should work for a model consisting of any
# combination of tuple/array/dict.
initoptim{T<:Number}(::KnetArray{T},otype)=eval(parse(otype))
initoptim{T<:Number}(::Array{T},otype)=eval(parse(otype))
initoptim(a::Associative,otype)=Dict(k=>initoptim(v,otype) for (k,v) in a)
initoptim(a,otype)=map(x->initoptim(x,otype), a)

# Q: charlm flat minibatch style or sort sentences? => flat
# Q: padding or masking? => no need
# Q: initial state zero or learnt? => zero
# Q: partial batches => no need

function minibatch(data, B)
    T = div(length(data),B)
    batches = Array{Vector{Int32}}(T)
    for t = 1:T
        batch = Array{Int32}(B)
        for b = 1:B
            batch[b] = data[(b-1)*T+t]
        end
        batches[t] = batch
    end
    return batches
end

function initvocab()
    global EOS = Int32(1) # use as both SOS in input and EOS in output.
    Dict{String,Int32}("<s>"=>EOS)
end

function loaddata(file, vocab)
    data = Int32[EOS]; nw = ns = 0
    for l in eachline(file); ns+=1
        for w in split(l); nw+=1
            push!(data, get!(vocab, w, 1+length(vocab)))
        end
        push!(data,EOS)
    end
    @msg 1 "$file: $ns sentences, $nw words, vocab=$(length(vocab)), corpus=$(length(data))"
    return data
end

function mikolovptb()
    files = [ Knet.dir("data","ptb.$x.txt") for x in ("train","valid","test") ]
    if any(!isfile(f) for f in files)
        tgz = Knet.dir("data","simple-examples.tgz")
        if !isfile(tgz)
            url = "http://www.fit.vutbr.cz/~imikolov/rnnlm/simple-examples.tgz"
            download(url,tgz)
        end
        run(`tar --strip-components 3 -C $(Knet.dir("data")) -xzf $tgz ./simple-examples/data/ptb.train.txt ./simple-examples/data/ptb.valid.txt ./simple-examples/data/ptb.test.txt`)
    end
    return files
end

function main(args=ARGS)
    global model, text, data, tok2int, o
    s = ArgParseSettings()
    s.description="rnnlm.jl: LSTM language model\n"
    s.exc_handler=ArgParse.debug_handler
    @add_arg_table s begin
        ("--datafiles"; nargs='+'; help="If provided, use first file for training, second for early stop, others for test. If not provided use mikolovptb files.")
        ("--loadfile"; help="Initialize model from file")
        ("--savefile"; help="Save final model to file")
        ("--bestfile"; help="Save best model to file")
        ("--epochs"; arg_type=Int; default=1; help="Number of epochs for training.")
        ("--hidden"; nargs='+'; arg_type=Int; default=[256]; help="Sizes of one or more LSTM layers.")
        ("--embed"; arg_type=Int; default=128; help="Size of the embedding vector.")
        ("--batchsize"; arg_type=Int; default=64; help="Number of sequences to train on in parallel.")
        ("--bptt"; arg_type=Int; default=20; help="Number of steps to unroll for bptt.")
        ("--optimization"; default="Adam()"; help="Optimization algorithm and parameters.")
        ("--dropout"; arg_type=Float64; default=0.0; help="Dropout probability.")
        ("--gcheck"; arg_type=Int; default=0; help="Check N random gradients.")
        ("--seed"; arg_type=Int; default=-1; help="Random number seed.")
        ("--atype"; default=(gpu()>=0 ? "KnetArray{Float32}" : "Array{Float32}"); help="array type: Array for cpu, KnetArray for gpu")
        ("--fast"; action=:store_true; help="skip loss printing for faster run")
        ("--loglevel"; arg_type=Int; default=1; help="display progress messages")
        # TODO: ("--generate"; arg_type=Int; default=0; help="If non-zero generate given number of tokens.")
    end
    isa(args, AbstractString) && (args=split(args))
    o = parse_args(args, s; as_symbols=true)
    global loglevel = o[:loglevel]
    if o[:seed] > 0; setseed(o[:seed]); end
    if isempty(o[:datafiles]); o[:datafiles] = mikolovptb(); end
    @msg 1 string(s.description,"opts=",[(k,v) for (k,v) in o]...)
    global vocab = initvocab()
    global text = map(f->loaddata(f,vocab), o[:datafiles])
    global data = map(t->minibatch(t, o[:batchsize]), text)
    global model = initmodel(eval(parse(o[:atype])), o[:hidden], length(vocab), o[:embed])
    function report(ep)
        l = [ loss(model,d;slen=o[:bptt]) for d in data ]
        l1 = Float32[ exp(x[1]/x[2]) for x in l ]
        l2 = [ x[2] for x in l ]
        if ep==0; @msg 1 (:epoch,ep,:perp,l1...,:size,l2...)
        else; @msg 1 (:epoch,ep,:perp,l1...); end
        return l1
    end
    if length(data) > 1; devset=2; else devset=1; end
    if !o[:fast]; @log 1 (losses = report(0)); devbest = devlast = losses[devset]; end
    global optim = initoptim(model,o[:optimization])
    Knet.knetgc(); gc() # TODO: fix this otherwise curand cannot initialize no memory left!
    for epoch=1:o[:epochs]
        @log 1 bptt(model, data[1], optim; pdrop=o[:dropout], slen=o[:bptt])
        if o[:fast]; continue; end
        @log 1 (losses = report(epoch))
        if o[:bestfile] != nothing && losses[devset] < devbest
            devbest = losses[devset]
            @log 1 save(o[:bestfile], "model", model, "vocab", vocab)
        end
        # if epoch > 6 # losses[devset] > devlast && isa(optim[1], Sgd)
        #     for p in optim; p.lr /= 1.2; end; @msg 1 "lr=$(optim[1].lr)"
        # end
        devlast = losses[devset]
        if o[:gcheck] > 0
            gradcheck(rnnlm, model, initstate(model,o[:batchsize]), data[1]; gcheck=o[:gcheck], verbose=true, kwargs=[(:range,1:o[:bptt])])
        end
    end
    if o[:savefile] != nothing
        @log 1 save(o[:savefile], "model", model, "vocab", vocab)
    end
    return model
end

# This allows both non-interactive (shell command) and interactive calls like:
# $ julia rnnlm.jl --epochs 10
# julia> RNNLM.main("--epochs 10")
if VERSION >= v"0.5.0-dev+7720"
    if basename(PROGRAM_FILE)==basename(@__FILE__); main(ARGS); end
else
    !isinteractive() && !isdefined(Core.Main,:load_only) && main(ARGS)
end

end  # module


text:
# Draw main plots for execution time & completeness
# Author: Thomas Minier
using Gadfly
using RDatasets
using CalliPlots

custom_theme = Theme(
    key_position = :top,
    line_width=2px,
    major_label_font_size=15px,
    key_title_font_size=18px,
    key_label_font_size=15px,
    major_label_font="Computer Modern",
    major_label_color=colorant"#000000",
    minor_label_color=colorant"#000000",
    key_label_color=colorant"#000000"
)

Gadfly.push_theme(custom_theme)

function colors()
 return Scale.color_discrete_manual(colorant"#990000", colorant"#FF5733", colorant"#FFC300", colorant"#b1d073")
end

# Execution times

time_ref = computeMean(
    :time,
    "scripts/curio/run1/execution_times_ref.csv",
    "scripts/curio/run2/execution_times_ref.csv",
    "scripts/curio/run3/execution_times_ref.csv"
)
time_ref[:approach] = "1 (TPF) "

time_2s = computeMean(
    :time,
    "scripts/curio/run1/execution_times_2servers.csv",
    "scripts/curio/run2/execution_times_2servers.csv",
    "scripts/curio/run3/execution_times_2servers.csv"
)
time_2s[:approach] = "2"

time_3s = computeMean(
    :time,
    "scripts/curio/run1/execution_times_3servers.csv",
    "scripts/curio/run2/execution_times_3servers.csv",
    "scripts/curio/run3/execution_times_3servers.csv"
)
time_3s[:approach] = "3"

time_20s = computeMean(
    :time,
    "scripts/curio/run1/execution_times_20servers.csv",
    "scripts/curio/run2/execution_times_20servers.csv",
    "scripts/curio/run3/execution_times_20servers.csv"
)
time_20s[:approach] = "20"

all = [time_ref;time_2s;time_3s;time_20s]

layer_ref = layer(time_ref, y=:time, x=:clients, size=[5px], Geom.point, Geom.line, Theme(point_shapes=Function[Gadfly.diamond], default_color=colorant"#990000"))
layer_2s = layer(time_2s, y=:time, x=:clients, size=[5px], Geom.point, Geom.line, Theme(point_shapes=Function[Gadfly.square], line_style=:dash, default_color=colorant"#FF5733"))
layer_3s = layer(time_3s, y=:time, x=:clients, size=[5px], Geom.point, Geom.line, Theme(point_shapes=Function[Gadfly.cross], default_color=colorant"#FFC300"))
layer_20s = layer(time_20s, y=:time, x=:clients, size=[5px], Geom.point, Geom.line, Theme(point_shapes=Function[Gadfly.dtriangle], line_style=:dash, default_color=colorant"#b1d073"))

color_guide = Guide.manual_color_key("Number of servers",
                            ["1 (TPF) ", "2", "3", "4"],
                            [colorant"#990000", colorant"#FF5733", colorant"#FFC300", colorant"#b1d073"])

p = plot(layer_ref,layer_2s, layer_3s, layer_20s, Guide.xlabel("Number of concurrent clients executing the query"), Guide.ylabel("avg. execution time (s)", orientation=:vertical), Guide.xticks(ticks=[1, 10, 20, 50, 100]), color_guide)

# RAW version
# p = plot([time_ref;time_mixed;time_ulysses], y=:time, x=:clients, color=:approach, Geom.line, Guide.xlabel("Number of clients"), Guide.ylabel("Execution time (s)"), Scale.x_discrete)
# Smooth version
# p = plot(all, y=:time, x=:clients, color=:approach, Geom.point, Geom.line(), #method=:loess, smoothing=0.75
# )

draw(PDF("scripts/curio/execution_time_with_load.pdf", 7inch, 3inch), p)
draw(PNG("scripts/curio/execution_time_with_load.png", 7inch, 3inch), p)

max_stars_repo_name: Callidon/ulysses-tpf

max_issues_repo_name: Callidon/ulysses-tpf

max_forks_repo_name: Callidon/ulysses-tpf


text:
using StatsBase
using YaoArrayRegister
export UnitaryChannel

# TODO: replace this with uweight
# when StatsBase 0.33 is out
_uweights(n) = weights(ones(n))

function check_nqubits(operators)
    N = nqubits(first(operators))
    for each in operators
        N == nqubits(each) || throw(AssertionError("number of qubits mismatch"))
    end
    return N
end

# NOTE: this can have a better support in YaoIR directly with StatsBase
# as an compiled instruction, but store the weights is the best solution
# here

"""
    UnitaryChannel(operators[, weights])

Create a unitary channel, optionally weighted from an list of weights.
The unitary channel is defined as below in Kraus representation

```math
ϕ(ρ) = \\sum_i U_i ρ U_i^†
```

!!! note
    Unitary channel will only normalize the weights when calculating the matrix form,
    thus you should be careful when you need this condition for other purpose.

!!! note
    when applying a `UnitaryChannel` on the register, a unitary will be sampled
    uniformly or optionally from given weights, then this unitary will be applied
    to the register.

# Example

```jldoctest; setup=:(using YaoBlocks, YaoArrayRegister)
julia> UnitaryChannel([X, Y, Z])
nqubits: 1
unitary_channel
├─ [1.0] X
├─ [1.0] Y
└─ [1.0] Z
```

Or with weights

```jldoctest; setup=:(using YaoBlocks, YaoArrayRegister)
julia> UnitaryChannel([X, Y, Z], [0.1, 0.2, 0.7])
nqubits: 1
unitary_channel
├─ [0.1] X
├─ [0.2] Y
└─ [0.7] Z
```
"""
struct UnitaryChannel{N,W<:AbstractWeights} <: CompositeBlock{N}
    operators::Vector{AbstractBlock{N}}
    weights::W

    function UnitaryChannel(operators::Vector)
        w = _uweights(length(operators))
        N = check_nqubits(operators)
        new{N,typeof(w)}(operators, w)
    end

    function UnitaryChannel(operators::Vector, w::AbstractWeights)
        N = check_nqubits(operators)
        new{N,typeof(w)}(operators, w)
    end

    function UnitaryChannel(operators::Vector, w::AbstractVector)
        w = weights(w)
        N = check_nqubits(operators)
        new{N,typeof(w)}(operators, w)
    end
end

UnitaryChannel(it, weights) = UnitaryChannel(collect(it), weights)
UnitaryChannel(it) = UnitaryChannel(collect(it))

function apply!(r::AbstractRegister, x::UnitaryChannel)
    apply!(r, sample(x.operators, x.weights))
end

function mat(::Type{T}, x::UnitaryChannel) where {T}
    U = sample(x.operators, x.weights)
    return mat(T, U)
end

subblocks(x::UnitaryChannel) = x.operators
chsubblocks(x::UnitaryChannel{N}, it) where {N} = UnitaryChannel{N}(collect(it), x.weights)
occupied_locs(x::UnitaryChannel) = union(occupied_locs.(x.operators)...)

function cache_key(x::UnitaryChannel)
    key = hash(x.weights)
    for each in x.operators
        key = hash(each, key)
    end
    return key
end

function Base.:(==)(lhs::UnitaryChannel{N}, rhs::UnitaryChannel{N}) where {N}
    return (lhs.weights == rhs.weights) && (lhs.operators == rhs.operators)
end

Base.adjoint(x::UnitaryChannel) = UnitaryChannel(adjoint.(x.operators), x.weights)


text:
function gen_qp_functions(solver::Solver)
    n,m,N = get_sizes(solver)
    p,pI,pE = get_num_constraints(solver)
    p_N,pI_N,pE_N = get_num_terminal_constraints(solver)

    c_function!, c_jacobian!, c_labels, cI!, cE! = generate_constraint_functions(solver.obj)
    f = solver.fd
    s_x = reshape(1:(N*n),n,N)
    s_u = reshape(N*n .+ (1:(m*(N-1))),m,N-1)
    s_ν = (N*n + (N-1)*m) .+ (1:N*n)
    sze = (X=s_x, U=s_u, ν=s_ν)

    function newton_cost(V::Vector)
        X,U = unpackV(V)
        J = cost(solver,X,U)
        return J
    end

    function cI(V)
        X,U = unpackV(V)
        C = zeros(eltype(V),pI,N-1)
        for k = 1:N-1
            cI!(view(C,1:pI,k),X[:,k],U[:,k])
        end
        CN = zeros(eltype(V),pI_N)
        cI!(CN,X[:,N])
        CI = [vec(C); CN]
        return CI
    end

    function cE(V)
        X,U = unpackV(V)
        C = zeros(eltype(V),pE,N-1)
        for k = 1:N-1
            cE!(view(C,1:pE,k),X[:,k],U[:,k])
        end
        CN = zeros(eltype(V),pE_N)
        cE!(CN,X[:,N])
        CE = [vec(C); CN]
        return CE
    end

    function d(V)
        X,U = unpackV(V)
        D = zeros(eltype(V),n,N)
        D[:,1] = X[:,1] - solver.obj.x0
        for k = 2:N
            f(view(D,1:n,k),X[:,k-1],U[:,k-1])
            D[:,k] -= X[:,k]
        end
        return vec(D)
    end

    function active_set(V,eps=1e-2)
        X,U,λ,μ,ν,s = unpackV(V)
        a = trues(length(V))
        ci = cI(V)
        c_inds = ci .>= -eps
        a[sze.μ] = c_inds
        a[sze.s] = c_inds
        return a
    end

    function packV(X,U)
        V = [vec(X);vec(U)]
    end

    function unpackV(V)
        X = V[sze.X]
        U = V[sze.U]
        return X,U,λ,μ,ν,s
    end
    return newton_cost, packV, unpackV, cI, cE, d, active_set, sze
end

model,obj = Dynamics.dubinscar
solver = Solver(model,obj,N=31)
n,m,N = get_sizes(solver)
U0 = rand(m,N-1)
res,stats = solve(solver,U0)
# λ_update_default!(res,solver);

newton_cost, packV, unpackV, cI, cE, d, active_set = gen_qp_functions(solver)
X = to_array(res.X)
U = to_array(res.U)
ν = vec(to_array(res.s))
V = packV(X,U)
Nz = length(V)

using JuMP, NLopt
using Gurobi
H = ForwardDiff.hessian(newton_cost,V)
g = ForwardDiff.gradient(newton_cost,V)
c = newton_cost(V)
function quad_cost(δV)
    newton_cost(V) + g'δV + 0.5*δV'H*δV
end

∇d = ForwardDiff.jacobian(d,V)
e = d(V)
d(V)
function quad_dynamics(δV)
    d(V) + ∇d*δV
end

function qp_step(V)
    H = ForwardDiff.hessian(newton_cost,V)
    g = ForwardDiff.gradient(newton_cost,V)
    c = newton_cost(V)

    ∇d = ForwardDiff.jacobian(d,V)
    e = d(V)

    mod = JuMP.Model(solver=GurobiSolver())
    @variable(mod, Z[1:Nz])
    obj = 0.5*Z'H*Z + g'Z + c
    @objective(mod, Min, obj)
    dyn = e + ∇d*Z;
    @constraint(mod, dyn .== 0)
    status = JuMP.solve(mod)
    return getvalue(Z)
end

V = packV(X,U)
newton_cost(V)
norm(d(V),Inf)

max_d = Float64[]

for i = 1:20
    V_ = qp_step(V) + V
    newton_cost(V_)
    push!(max_d,norm(d(V_),Inf))
    copyto!(V,V_)
end

plot(max_d)


text:
type BGZF
    box::Int32
    cache_size::Int32
    block_length::Int32
    block_offset::Int32
    block_address::Int64
    uncompressed_address::Int64
    uncompressed_block::Ptr{Void}
    compressed_block::Ptr{Void}
    cache::Ptr{Void}
    fp::Ptr{Void}
    mt::Ptr{Void}
    idx::Ptr{Void}
    idx_build_otf::Int32
    gz_stream::Ptr{Void}
end

type Bam_hdr_t
    n_targets::Int32
    ignore_sam_err::Int32
    l_text::UInt32
    target_len::UInt32
    cigar_tab::Ptr{Int8}
    target_name::Ptr{Ptr{Char}}
    text::Ptr{Char}
    sdict::Ptr{Void}
end

type Bam1_core_t
    tid::Int32
    pos::Int32
    box1::UInt32
    box2::UInt32
    l_qseq::Int32
    mtid::Int32
    mpos::Int32
    isize::Int32
end

type Bam1_t
    core::Bam1_core_t
    l_data::Int32
    m_data::Int32
    data::Ptr{UInt8}
    id::UInt64
end

type KStr
    l::Csize_t
    m::Csize_t
    s::Ptr{Cchar}
end

function get_ptr_of_null_kstr(num::Int64)
    ccall((:get_ptr_of_kstr, LIBKSTR),Ptr{KStr}, (Int64,), num)
end

function kstrToASCII(kstr::KStr)
    strptr(kstr.s, 1, kstr.l)
end

function kstrToASCII(pkstr::Ptr{KStr})
    kstr = unsafe_load(pkstr)
    kstrToASCII(kstr)
end
type Hts_itr_t
end
type HtsFile
end
type BGZF
end
typealias HtsFile SamFile

type Bam_pileup1_t
end

type Bam_plp_t
end
type Bam_mplp_t
end


text:
using SafeTestsets

@safetestset "Chain Environment Tests" begin include("chain.jl") end