text:
module Day09

import ..data_dir # from parent module
input = read(joinpath(data_dir, "day09"), String)

export part1, part2

"""
Horrendous logic for converting this string of integers to a matrix.
There must be a better way, but I can't find it.
"""
function string_block_to_matrix(string_block::AbstractString)
    as_arrays = split.(split(string_block, "\n"), "")
    as_int_arrays = [parse.(Int64, array) for array in as_arrays]
    # Transposing changes the type!
    Matrix(transpose(hcat(as_int_arrays...)))
end

function out_of_bounds(heatmap::Matrix{Int64}, row::Int64, column::Int64)
    heatmap_height, heatmap_width = size(heatmap)
    row < 1 || row > heatmap_height || column < 1 || column > heatmap_width
end

function identify_neighbours(heatmap::Matrix{Int64}, row::Int64, column::Int64)
    candidates = [
        (row - 1, column), (row + 1, column),
        (row, column - 1), (row, column + 1)
    ]
    filter(candidate -> !out_of_bounds(heatmap, candidate...), candidates)
end

function risk_level(heatmap::Matrix{Int64}, row::Int64, column::Int64)
    value = heatmap[row, column] # matrices are indexed row then column
    if value == 9
        return (0)
    end
    neighbours = identify_neighbours(heatmap, row, column)
    is_minimum = all(value < heatmap[neighbour...] for neighbour in neighbours)
    is_minimum * (1 + value) # non-minima don't get a risk value
end

function add_to_basin(
    heatmap::Matrix{Int64},
    unclaimed::Vector{Tuple{Int64,Int64}},
    row::Int64,
    column::Int64
)
    # can't pop by value, so I need to get the index
    index_in_unclaimed = findfirst(x -> x == (row, column), unclaimed)
    if isnothing(index_in_unclaimed) # already claimed by another path
        return Set{Tuple{Int64,Int64}}()
    end
    in_basin = Set([popat!(unclaimed, index_in_unclaimed)])
    neighbours = [neighbour for neighbour in identify_neighbours(heatmap, row, column)
                  if neighbour in unclaimed]
    if length(neighbours) == 0
        return (in_basin)
    end
    for neighbour in neighbours
        in_basin = in_basin ∪ add_to_basin(heatmap, unclaimed, neighbour...)
    end
    in_basin
end

function part1(input = input)
    heatmap = string_block_to_matrix(input)
    n_rows, n_cols = size(heatmap)
    total_risk_level = 0
    for row = 1:n_rows, column = 1:n_cols
        total_risk_level += risk_level(heatmap, row, column)
    end
    total_risk_level
end

"""
Two things make this second part manageable:

* I didn't need to consider the actual values of the heatmap,
    as long as I knew which ones were 9 and which ones aren't
* Every non-9 element of the heatmap is in exactly one basin
"""
function part2(input = input)
    heatmap = string_block_to_matrix(input)
    n_rows, n_cols = size(heatmap)
    unclaimed = [(row, column) for row in 1:n_rows, column in 1:n_cols if heatmap[row, column] != 9]
    basins = Vector{Set{Tuple{Int64,Int64}}}()
    while length(unclaimed) > 0
        push!(basins, add_to_basin(heatmap, unclaimed, unclaimed[1]...))
    end
    basin_sizes = map(x -> length(x), basins)
    last(sort(basin_sizes), 3) |> prod
end

end #module


text:
export CDS_MVp

"""
compute single-thread matrix vector product with vector x, output is vector y: y=A*x
MVP is in the compressed diagonal format.
R is a tall matrix N by ndiagonals, corresponding to a square matrix A
offsets indicate offset of diagonal compared to the main diagonal in A (which is 0)
"""
function CDS_MVp(
                        N      ::Integer,
                        ndiags ::Integer,
                        R      ::Array{TF,2},
                        offset ::Vector{TI},
                        x      ::Vector{TF},
                        y      ::Vector{TF}) where {TF<:Real,TI<:Integer}

  for i = 1 : ndiags
      d  = offset[i]
      r0 = max(1, 1-d)
      r1 = min(N, N-d)
      c0 = max(1, 1+d)
      for r = r0 : r1
        c = r - r0 + c0 #original
        @inbounds y[r] = y[r] + R[r,i] * x[c]#original
      end
  end
  return y
end

# #the versions below is insignificantly faster than the original version above
# function CDS_MVp(
#   N::Integer,
#   ndiags::Integer,
#   R::Array{TF,2},
#   offset::Vector{TI},
#   x::Vector{TF},
#   y::Vector{TF}) where {TF<:Real,TI<:Integer}

#   for i = 1 : ndiags
#     d = offset[i]
#     r0 = max(1, 1-d)
#     r1 = min(N, N-d)
#     c0 = max(1, 1+d)
#     #r0_plus_c0 = r0 + c0
#     c = deepcopy(c0)
#     for r = r0 : r1
#       @inbounds y[r] = y[r] + R[r,i] * x[c]#original
#       c += 1
#     end
#   end
# return y
# end

# function CDS_MVp2(
#   N      ::Integer,
#   ndiags ::Integer,
#   R      ::Array{TF,2},
#   offset ::Vector{TI},
#   x      ::Vector{TF},
#   y      ::Vector{TF}) where {TF<:Real,TI<:Integer}

#   for i = 1 : ndiags
#     d  = offset[i]
#     r0 = max(1, 1-d)
#     r1 = min(N, N-d)
#     c0 = max(1, 1+d)
#     # for r = r0 : r1
#     #   c = r - r0 + c0 #original
#     #   @inbounds y[r] = y[r] + R[r,i] * x[c]#original
#     # end
#     @inbounds y[r0 : r1] .= y[r0 : r1] .+ R[r0:r1,i].*x[r0 - r0 + c0:r1 - r0 + c0]
#   end

#   return y
# end

ext: jl

max_stars_repo_name: slimgroup/SetIntersectionProjection

max_issues_repo_name: slimgroup/SetIntersectionProjection

max_forks_repo_name: slimgroup/SetIntersectionProjection


text:
#    Copyright (c) 2021 J.A. Duffek
#    Copyright (c) 2000 D.M. Spink
#
#    This program is free software; you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation; either version 2 of the License, or
#    (at your option) any later version.
#
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
#    You should have received a copy of the GNU General Public License
#    along with this program; if not, see <http://www.gnu.org/licenses/>.

"""
    nrbtestsrf()::NURBS2D

Construct a simple test surface.

# Output:
- `srf`: NURBS data structure of a surface

# Examples:
```julia
julia> srf = nrbtestsrf()
```
"""
function nrbtestsrf()::NURBS2D
# allocate multi-dimensional array of control points
pnts = zeros(typeof(1.0),(3,5,5));

# define a grid of control points
# in this case a regular grid of u,v points
# pnts[1:3,u,v]

pnts[:,:,1] = [ 0.0  3.0  5.0  8.0 10.0;  # w*x
                0.0  0.0  0.0  0.0  0.0;  # w*y
                2.0  2.0  7.0  7.0  8.0]; # w*z

pnts[:,:,2] = [ 0.0  3.0  5.0  8.0 10.0;
                3.0  3.0  3.0  3.0  3.0;
                0.0  0.0  5.0  5.0  7.0];

pnts[:,:,3] = [ 0.0  3.0  5.0  8.0 10.0;
                5.0  5.0  5.0  5.0  5.0;
                0.0  0.0  5.0  5.0  7.0];

pnts[:,:,4] = [ 0.0  3.0  5.0  8.0 10.0;
                8.0  8.0  8.0  8.0  8.0;
                5.0  5.0  8.0  8.0 10.0];

pnts[:,:,5] = [ 0.0  3.0  5.0  8.0 10.0;
               10.0 10.0 10.0 10.0 10.0;
                5.0  5.0  8.0  8.0 10.0];

# knots
knots = Vector{Vector{typeof(1.0)}}(undef,2);
knots[1] = vec([0 0 0 1/3 2/3 1 1 1]); # knots along u
knots[2] = vec([0 0 0 1/3 2/3 1 1 1]); # knots along v

# make nurbs surface
return nrbmak(pnts,knots);
end # nrbtestsrf


text:
module StringInterpolation

export @interpolate, interp_parse

"""
Interpolates an escaped string, e.g.

julia> x = "World"; @interpolate \"Hello \\ \b\$ \bx\"

"Hello World"
"""
macro interpolate(args...)
    return interp_parse(args...)
end

"""
String interpolation parsing

Resurrected from Julia base:
https://github.com/JuliaLang/julia/blob/deab8eabd7089e2699a8f3a9598177b62cbb1733/base/string.jl
"""
function interp_parse(s::AbstractString, unescape::Function, printer::Function)
    sx = []
    i = j = start(s)
    while !done(s,j)
        c, k = next(s,j)
        if c == '$'
            if !isempty(s[i:j-1])
                push!(sx, unescape(s[i:j-1]))
            end
            ex, j = parse(s,k,greedy=false)
            if isa(ex,Expr) && is(ex.head,:continue)
                throw(ParseError("Incomplete expression"))
            end
            push!(sx, esc(ex))
            i = j
        elseif c == '\\' && !done(s,k)
            if s[k] == '$'
                if !isempty(s[i:j-1])
                    push!(sx, unescape(s[i:j-1]))
                end
                i = k
            end
            c, j = next(s,k)
        else
            j = k
        end
    end
    if !isempty(s[i:end])
        push!(sx, unescape(s[i:j-1]))
    end
    length(sx) == 1 && isa(sx[1],ByteString) ? sx[1] :
    Expr(:call, :sprint, printer, sx...)
end

interp_parse(s::AbstractString, u::Function) = interp_parse(s, u, print)
interp_parse(s::AbstractString) = interp_parse(s, x -> isvalid(UTF8String,unescape_string(x)) ? unescape_string(x) : error("Invalid UTF-8 sequence"))

end


text:
using DataFrames
using Statistics

import Dates
import LibGit2
import Pkg
import Pkg.TOML

parse_toml_config(config::T where T <: AbstractString) =
    parse_toml_config(TOML.parsefile(config))

function get_default_config()
    return deepcopy(Dict{String, Any}(
        "Data" => Dict{String, Any}(
            "x-column" => "x",
            "y-column" => "y",
            "z-column" => "z",
            "gene-column" => "gene",
            "min-molecules-per-gene" => 1,
            "min-molecules-per-cell" => 3,
            "estimate-scale-from-centers" => true,
            "scale" => nothing,
            "scale-std" => "25%",
            "min-molecules-per-segment" => nothing
        ),
        "Sampling" => Dict{String, Any}(
            "new-component-weight" => 0.2,
            "new-component-fraction" => 0.3
        ),
        "Plotting" => Dict{String, Any}(
            "gene-composition-neigborhood" => nothing,
            "min-pixels-per-cell" => 15
        )
    ))
end

function parse_toml_config(config::Dict{AS, Any}) where AS <: AbstractString
    res_config = get_default_config()
    for (k,v) in config
        if !(k in keys(res_config))
            error("Unexpected value in the config: '$k'")
        end

        cur_def = res_config[k]

        for (k2,v2) in v
            if !(k2 in keys(cur_def))
                error("Unexpected value in the config: '$k' -> '$k2'")
            end

            cur_def[k2] = v2
        end
    end

    return res_config
end

function extend_params_with_config!(params::Dict, config::Dict)
    for sub_cfg in values(config)
        for (k, v) in sub_cfg
            if !(k in keys(params)) || params[k] === nothing
                params[k] = v
            end
        end
    end
end

function load_df(args::Dict; kwargs...)
    exc_genes = (args["exclude-genes"] === nothing) ? String[] : String.(strip.(Base.split(args["exclude-genes"], ",")))

    df_spatial, gene_names = load_df(args["coordinates"]; x_col=args["x-column"], y_col=args["y-column"], z_col=args["z-column"], gene_col=args["gene-column"],
        min_molecules_per_gene=args["min-molecules-per-gene"], exclude_genes=exc_genes, kwargs...)

    if :z in propertynames(df_spatial)
        if ("force-2d" in keys(args)) && args["force-2d"] || (length(unique(df_spatial.z)) < 2)
            select!(df_spatial, Not(:z))
        end
    elseif args["z-column"] != :z
        error("z-column $(args["z-column"]) not found in the data")
    end

    return df_spatial, gene_names
end

append_suffix(output::String, suffix) = "$(splitext(output)[1])_$suffix"

function get_baysor_run_str()::String
    pkg_str = "v$pkg_version"
    try
        pkg_info = Pkg.dependencies()[Base.UUID("cc9f9468-1fbe-11e9-0acf-e9460511877c")]
        repo = LibGit2.GitRepo(pkg_info.source) # In case it's a development version with a commit instead of a release
        hash_str = "$(LibGit2.GitShortHash(LibGit2.peel(LibGit2.GitCommit, LibGit2.head(repo))))"
        pkg_str = "$pkg_str [$hash_str]"
    catch err
    end

    return "($(Dates.Date(Dates.now()))) Run Baysor $pkg_str"
end

run_cli(args::String) = run_cli(String.(Base.split(args)))

function run_cli(args::Vector{String}=ARGS)::Cint
    help_message = "Usage: baysor <command> [options]\n\nCommands:\n\trun\t\trun segmentation of the dataset\n\tpreview\t\tgenerate preview diagnostics of the dataset\n"

    debug = false
    if "--debug" in args
        args = args[args .!= "--debug"]
        debug = true
    end

    try
        if (length(args) == 0) || (length(args) == 1) && (args[1] == "-h" || args[1] == "--help")
            println(help_message)
            return 0
        end

        if args[1] == "run"
            return run_cli_main(args[2:end])
        end

        if args[1] == "preview"
            return run_cli_preview(args[2:end])
        end

        @error "Can't parse argument $(args[1])"
        println(help_message)
        return 1
    catch err
        if debug
            rethrow()
        else
            @error("$err\n\n" * join(["$s" for s in stacktrace(catch_backtrace())], "\n"))
        end
    end

    return 2
end

julia_main()::Cint = run_cli()

max_stars_repo_name: repushko/Baysor

max_issues_repo_name: repushko/Baysor

max_forks_repo_name: repushko/Baysor


text:
# TODO: Make tests pass for
# - T ∈ [Float32, Float64]
# - Setup ∈ [PlateSetup, CylinderSetup, SphereSetup, NeuronSetup]
# - Gradient ∈ [GeneralGradient, ScalarGradient] (all time profiles)

# Floating point type for simulations
T = Float64

@testset "Setup recipes" begin
    @testset "PlateSetup" begin
        setup = get_setup(PlateSetup{T})
        coeffs = get_coeffs(setup)
        femesh, = create_geometry(setup; recreate = true)
        model = Model(; mesh = femesh, coeffs...)
    end

    @testset "CylinderSetup" begin
        setup = get_setup(CylinderSetup{T})
        coeffs = get_coeffs(setup)
        femesh, = create_geometry(setup; recreate = true)
        model = Model(; mesh = femesh, coeffs...)
    end

    @testset "SphereSetup" begin
        setup = get_setup(SphereSetup{T})
        coeffs = get_coeffs(setup)
        @test_broken femesh, = create_geometry(setup; recreate = true)
        @test_broken model = Model(; mesh = femesh, coeffs...)
    end

    @testset "NeuronSetup" begin
        setup = get_setup(NeuronSetup{T})
        coeffs = get_coeffs(setup)
        femesh, = create_geometry(setup; recreate = true)
        model = Model(; mesh = femesh, coeffs...)
    end
end

# Geometrical setup
setup = get_setup(CylinderSetup{T})
coeffs = get_coeffs(setup)

# Get compartimentalized coefficient vectors
femesh, = create_geometry(setup; recreate = true)
model = Model(; mesh = femesh, coeffs...)
volumes = get_cmpt_volumes(model.mesh)
D_avg = 1 / 3 * tr.(model.D)' * volumes / sum(volumes)
ncompartment = length(model.mesh.points)
matrices = assemble_matrices(model);

## Gradients
dir = [1.0, 1.0, 1.0]
profile = PGSE(2000.0, 6000.0)
b = 1000
g = √(b / int_F²(profile)) / coeffs.γ
pgse_gradient = ScalarGradient(dir, profile, g)

dir = [1.0, 1.0, 1.0]
profile = CosOGSE(5000.0, 5000.0, 2)
b = 1000
g = √(b / int_F²(profile)) / coeffs.γ
ogse_gradient = ScalarGradient(dir, profile, g)

TE = 5000
φ = -π / 6
R = [cos(φ) sin(φ) 0; -sin(φ) cos(φ) 0; 0 0 1]
g⃗(t) = 1.0 * R * [sin(2π * t / TE), sin(20π * t / TE) / 5, cos(2π * t / TE)]
general_gradient = GeneralGradient{T,typeof(g⃗)}(; g⃗, TE)

@testset "Signal" begin
    ρ = initial_conditions(model)
    @test compute_signal(matrices.M, ρ) isa Complex{T}
    @test compute_signal.(matrices.M_cmpts, split_field(model.mesh, ρ)) isa
          Vector{Complex{T}}
end

@testset "ADC STA" begin
    @test compute_adc_sta(model, pgse_gradient) isa Vector{T}
    @test compute_adc_sta(model, ogse_gradient) isa Vector{T}
    @test_throws MethodError compute_adc_sta(model, general_gradient)
end

@testset "HADC" begin
    hadc = HADC(;
        model,
        matrices,
        odesolver = QNDF(autodiff = false),
        reltol = 1e-4,
        abstol = 1e-6,
    )
    @test_throws MethodError solve(hadc, general_gradient)
    @test solve(hadc, pgse_gradient) isa Vector{T}
    @test solve(hadc, ogse_gradient) isa Vector{T}
end

@testset "BTPDE" begin
    btpde = GeneralBTPDE(;
        model,
        matrices,
        reltol = 1e-4,
        abstol = 1e-6,
        odesolver = QNDF(autodiff = false),
    )
    @test solve(btpde, general_gradient) isa Vector{Complex{T}}
    @test solve(btpde, ogse_gradient) isa Vector{Complex{T}}

    btpde = IntervalConstanBTPDE{T}(; model, matrices, θ = 0.5, timestep = 5)
    @test_throws MethodError solve(btpde, general_gradient) isa Vector{Complex{T}}
    @test_throws ErrorException solve(btpde, ogse_gradient) isa Vector{Complex{T}}
    @test solve(btpde, pgse_gradient) isa Vector{Complex{T}}
end

@testset "Karger" begin
    # Compute HADC and fit difftensors
    directions = unitsphere(10)
    gradients = [
        ScalarGradient(collect(d), pgse_gradient.profile, pgse_gradient.amplitude) for
        d ∈ eachcol(directions)
    ]
    hadc = HADC(;
        model,
        matrices,
        odesolver = QNDF(autodiff = false),
        reltol = 1e-4,
        abstol = 1e-6,
    )
    adcs, = solve_multigrad(hadc, gradients)
    difftensors = fit_tensors(directions, adcs)

    # Solve Karger
    karger = Karger(; model, difftensors, odesolver = MagnusGL6(), timestep = 5.0)
    signal = solve(karger, pgse_gradient)
end

@testset "Matrix formalism" begin
    # Perform Laplace eigendecomposition
    laplace = Laplace{T}(; model, matrices, neig_max = 400)
    lap_eig = solve(laplace)
    length_scales = eig2length.(lap_eig.values, D_avg)

    # Truncate basis at minimum length scale
    length_scale = 3
    λ_max = length2eig(length_scale, D_avg)
    lap_eig = limit_lengthscale(lap_eig, λ_max)
    @test all(λ -> eig2length(λ, D_avg) ≥ length_scale, lap_eig.values)

    mf = MatrixFormalism(; model, matrices, lap_eig, ninterval = 500)
    @test solve(mf, general_gradient) isa Vector{Complex{T}}
    @test solve(mf, pgse_gradient) isa Vector{Complex{T}}
    @test solve(mf, ogse_gradient) isa Vector{Complex{T}}
end

@testset "Analytical" begin
    length_scale = 1.0
    eigstep = 1e-8
    eiglim = length2eig(length_scale, D_avg)
    analytical_coeffs = analytical_coefficients(setup, coeffs)
    analytical_laplace = AnalyticalLaplace(; analytical_coeffs..., eiglim, eigstep)
    lap_mat = solve(analytical_laplace)
    analytical_mf = AnalyticalMatrixFormalism(; analytical_laplace, lap_mat, volumes)
    @test_throws MethodError solve(analytical_mf, general_gradient)
    @test_throws ErrorException solve(analytical_mf, ogse_gradient)
    @test solve(analytical_mf, pgse_gradient) isa Complex{T}
end


text:
using Test
using Qecsim.App
using Qecsim.BasicModels:FiveQubitCode
using Qecsim.GenericModels:BitFlipErrorModel, DepolarizingErrorModel, NaiveDecoder
using Qecsim.Model:Model, ErrorModel, Decoder, DecodeResult, logical_xs, logical_zs
using Qecsim.PauliTools:to_bsf
using JSON
using Random:MersenneTwister

# TODO: profiling / type-stability checks

# test stub error model that generates a fixed given error
struct _FixedErrorModel <: ErrorModel
    error::BitVector
end
Model.label(::_FixedErrorModel) = "fixed"
Model.generate(error_model::_FixedErrorModel, x...) = error_model.error
# test stub decoder that decodes to a fixed given decode result
struct _FixedDecoder <: Decoder
    result::DecodeResult
end
Model.label(::_FixedDecoder) = "fixed"
Model.decode(decoder::_FixedDecoder, x...; kwargs...) = decoder.result
# test stub decoder that decodes to a cycle of given decode results
struct _CycleDecoder <: Decoder
    result_iter  # Iterator{DecodeResult}
    _CycleDecoder(results) = new(Iterators.Stateful(Iterators.Cycle(results)))
end
Model.label(::_CycleDecoder) = "cycle"
Model.decode(decoder::_CycleDecoder, x...; kwargs...) = popfirst!(decoder.result_iter)

# test stub code using duck-typing
struct _DuckCode end
Model.label(::_DuckCode) = "duck"
Model.stabilizers(::_DuckCode) = to_bsf(["XZZXI", "IXZZX", "XIXZZ", "ZXIXZ"])
Model.logical_xs(::_DuckCode) = to_bsf(["XXXXX"])
Model.logical_zs(::_DuckCode) = to_bsf(["ZZZZZ"])
Model.nkd(::_DuckCode) = (5, 1, 3)
# test stub error model using duck-typing
struct _DuckErrorModel end
Model.label(::_DuckErrorModel) = "duck"
Model.generate(::_DuckErrorModel, x...) = to_bsf("IIIII")
# test stub decoder using duck-typing
struct _DuckDecoder end
Model.label(::_DuckDecoder) = "duck"
Model.decode(::_DuckDecoder, x...; kwargs...) = DecodeResult(recovery=to_bsf("IIIII"))

@testset "RunResult" begin
    r0 = RunResult(false, 3, BitVector([0, 1]), nothing)  # reference
    r1 = RunResult(false, 3, BitVector([0, 1]), nothing)  # same
    r2 = RunResult(true, 3, BitVector([0, 1]), nothing)   # different success
    r3 = RunResult(false, 3, BitVector([0, 0]), nothing)  # different logical_commutations
    r4 = RunResult(false, 4, BitVector([0, 1]), nothing)  # different error_weight
    @test r1 == r0 && hash(r1) == hash(r0)
    @test r2 != r0 && hash(r2) != hash(r0)
    @test r3 != r0 && hash(r3) != hash(r0)
    @test r4 != r0 && hash(r4) != hash(r0)
    @test isequal(r1, r0) && hash(r1) == hash(r0)
    @test !isequal(r2, r0) && hash(r2) != hash(r0)
    @test !isequal(r3, r0) && hash(r3) != hash(r0)
    @test !isequal(r4, r0) && hash(r4) != hash(r0)
end

@testset "qec_run_once" begin
    # simulation parameters
    code = FiveQubitCode()
    error_model = DepolarizingErrorModel()
    decoder = NaiveDecoder()
    # float p
    data = qec_run_once(code, error_model, decoder, 0.1)
    @test !data.success || !any(data.logical_commutations)  # success -> zero commutator
    # rational p
    data = qec_run_once(code, error_model, decoder, 1 // 10)
    @test !data.success || !any(data.logical_commutations)  # success -> zero commutator
    # seeded run
    data1 = qec_run_once(code, error_model, decoder, 0.1, MersenneTwister(13))
    data2 = qec_run_once(code, error_model, decoder, 0.1, MersenneTwister(13))
    @test data1 == data2
    # warn: RECOVERY DOES NOT RETURN TO CODESPACE
    error_model = _FixedErrorModel(to_bsf("IIIII"))
    decoder = _FixedDecoder(DecodeResult(recovery=to_bsf("IIIIX")))
    @test_logs (:warn, "RECOVERY DOES NOT RETURN TO CODESPACE") #=
        =# data = qec_run_once(code, error_model, decoder, 0.1)
    @test !data.success
end

@testset "qec_run_once-duck-typing" begin
    # full duck-type simulation
    data = qec_run_once(_DuckCode(), _DuckErrorModel(), _DuckDecoder(), 0.1)
    @test data.success
    @test !any(data.logical_commutations)
    @test data.error_weight == 0
    @test isnothing(data.custom_values)
    # mixed duck-type simulation (no exceptions thrown)
    qec_run_once(_DuckCode(), _DuckErrorModel(), NaiveDecoder(), 0.1)  # no exception
    qec_run_once(_DuckCode(), BitFlipErrorModel(), NaiveDecoder(), 0.1)  # no exception
end

@testset "qec_run_once-override" begin
    # common simulation parameters
    code = FiveQubitCode()
    identity = to_bsf("IIIII")
    error_model = _FixedErrorModel(identity)
    p = 0.1
    # tests with and without overrides
    for (decode_result, expected) in [
        # identity recovery
        (DecodeResult(recovery=identity),
            RunResult(true, 0, BitVector([0, 0]), nothing)),
        # logical_x recovery
        (DecodeResult(recovery=logical_xs(code)[1,:]),
            RunResult(false, 0, BitVector([0, 1]), nothing)),
        # logical_z recovery
        (DecodeResult(recovery=logical_zs(code)[1,:]),
            RunResult(false, 0, BitVector([1, 0]), nothing)),
        # identity but override success=false
        (DecodeResult(success=false, recovery=identity),
            RunResult(false, 0, BitVector([0, 0]), nothing)),
        # identity but override logical_commutations=[1, 1]
        (DecodeResult(recovery=identity, logical_commutations=[1, 1]),
            RunResult(true, 0, BitVector([1, 1]), nothing)),
        # identity but override success=false and logical_commutations=[1, 1]
        (DecodeResult(success=false, recovery=identity, logical_commutations=[1, 1]),
            RunResult(false, 0, BitVector([1, 1]), nothing)),
        # no-recovery but override success=false
        (DecodeResult(success=false),
            RunResult(false, 0, nothing, nothing)),
        # no-recovery but override success=false and logical_commutations=[1, 1]
        (DecodeResult(success=false, logical_commutations=[1, 1]),
            RunResult(false, 0, BitVector([1, 1]), nothing)),
        # no-recovery but override success=false and custom_values=[1, 2]
        (DecodeResult(success=false, custom_values=[1, 2]),
            RunResult(false, 0, nothing, [1, 2])),
        # no-recovery but override success=false and custom_values=[1., 2.]
        (DecodeResult(success=false, custom_values=[1., 2.]),
            RunResult(false, 0, nothing, [1., 2.])),
        # no-recovery but override success=false and custom_values=[1//2, 1//3]
        (DecodeResult(success=false, custom_values=[1 // 2, 1 // 3]),
            RunResult(false, 0, nothing, [1 // 2, 1 // 3])),
    ]
        decoder = _FixedDecoder(decode_result)
        data = qec_run_once(code, error_model, decoder, p)
        @test data == expected
    end
end

@testset "qec_run" begin
    # simulation parameters
    code = FiveQubitCode()
    error_model = DepolarizingErrorModel()
    decoder = NaiveDecoder()
    max_runs = 1000
    for p in [0.25, 1 // 4]  # test with float and rational p
        data = qec_run(code, error_model, decoder, p; max_runs=max_runs)
        # JSON.print(data, 4)
        expected_keys = Set([:code, :n_k_d, :time_steps, :error_model, :decoder,
            :error_probability, :measurement_error_probability, :n_run, :n_success, :n_fail,
            :n_logical_commutations, :custom_totals, :error_weight_total,
            :error_weight_pvar, :logical_failure_rate, :physical_error_rate, :wall_time])
        @test keys(data) == expected_keys
        @test data[:error_probability] == p
        @test typeof(data[:error_probability]) == typeof(p)
        @test data[:n_run] == max_runs
        @test data[:n_success] + data[:n_fail] == data[:n_run]
        @test data[:n_success] >= 0 && data[:n_fail] >= 0
        @test data[:n_fail] <= sum(data[:n_logical_commutations])
        @test data[:logical_failure_rate] == data[:n_fail] / data[:n_run]
        # physical_error_rate
        p_rate = data[:physical_error_rate]
        p_rate_std = sqrt(data[:error_weight_pvar] / (data[:n_k_d][1]^2))
        @test p_rate - p_rate_std < p < p_rate + p_rate_std
        # run count
        data = qec_run(code, error_model, decoder, p)
        @test data[:n_run] == 1
        data = qec_run(code, error_model, decoder, p; max_runs=10)
        @test data[:n_run] == 10
        data = qec_run(code, error_model, decoder, p; max_failures=2)
        @test data[:n_fail] == 2
        data = qec_run(code, error_model, decoder, p; max_runs=10, max_failures=3)
        @test ((data[:n_run] == 10 && data[:n_fail] <= 3)
            || (data[:n_run] <= 10 && data[:n_fail] == 3))
    end
    # seeded run
    data1 = qec_run(code, error_model, decoder, 0.1, 13; max_runs=max_runs)
    data2 = qec_run(code, error_model, decoder, 0.1, 13; max_runs=max_runs)
    delete!(data1, :wall_time)
    delete!(data2, :wall_time)
    @test data1 == data2
end

@testset "qec_run-duck-typing" begin
    # full duck-type simulation
    data = qec_run(_DuckCode(), _DuckErrorModel(), _DuckDecoder(), 0.1; max_runs=10)
    @test data[:n_success] == 10
    @test all(data[:n_logical_commutations] .== 0)
    @test data[:error_weight_total] == 0
    @test isnothing(data[:custom_totals])
    # mixed duck-type simulation (no exceptions thrown)
    qec_run(_DuckCode(), _DuckErrorModel(), NaiveDecoder(), 0.1; max_runs=10)
    qec_run(_DuckCode(), BitFlipErrorModel(), NaiveDecoder(), 0.1; max_runs=10)
end

@testset "qec_run-override" begin
    # common simulation parameters
    code = FiveQubitCode()
    identity = to_bsf("IIIII")
    error_model = _FixedErrorModel(identity)
    p = 0.1
    # tests with and without overrides
    for (decoder, max_runs, expected) in [
        # identity recovery
        (_FixedDecoder(DecodeResult(recovery=identity)), 1,
            Dict(:n_success => 1, :n_fail => 0, :n_logical_commutations => [0, 0])),
        # logical_x recovery
        (_FixedDecoder(DecodeResult(recovery=logical_xs(code)[1,:])), 2,
            Dict(:n_success => 0, :n_fail => 2, :n_logical_commutations => [0, 2])),
        # logical_z recovery
        (_FixedDecoder(DecodeResult(recovery=logical_zs(code)[1,:])), 3,
            Dict(:n_success => 0, :n_fail => 3, :n_logical_commutations => [3, 0])),
        # identity but override success=false
        (_FixedDecoder(DecodeResult(success=false, recovery=identity)), 4,
            Dict(:n_success => 0, :n_fail => 4, :n_logical_commutations => [0, 0])),
        # identity but override logical_commutations=[1, 1]
        (_FixedDecoder(DecodeResult(recovery=identity, logical_commutations=[1, 1])), 5,
            Dict(:n_success => 5, :n_fail => 0, :n_logical_commutations => [5, 5])),
        # identity but override success=false and logical_commutations=[1, 1]
        (_FixedDecoder(DecodeResult(success=false, recovery=identity,
            logical_commutations=[1, 1])), 6,
            Dict(:n_success => 0, :n_fail => 6, :n_logical_commutations => [6, 6])),
        # no-recovery but override success=false
        (_FixedDecoder(DecodeResult(success=false)), 7,
            Dict(:n_success => 0, :n_fail => 7, :n_logical_commutations => nothing)),
        # no-recovery but override success=false and logical_commutations=[1, 1]
        (_FixedDecoder(DecodeResult(success=false, logical_commutations=[1, 1])), 8,
            Dict(:n_success => 0, :n_fail => 8, :n_logical_commutations => [8, 8])),
        # no-recovery but override success=false and custom_values=[1, 1]
        (_FixedDecoder(DecodeResult(success=false, custom_values=[1, 1])), 9,
            Dict(:n_success => 0, :n_fail => 9, :custom_totals => [9, 9])),
        # no-recovery but override success=false and custom_values=[1., 1.]
        (_FixedDecoder(DecodeResult(success=false, custom_values=[1., 1.])), 10,
            Dict(:n_success => 0, :n_fail => 10, :custom_totals => [10., 10.])),
        # no-recovery but override success=false and custom_values=[1//2, 1//3]
        (_FixedDecoder(DecodeResult(success=false, custom_values=[1 // 2, 1 // 3])), 11,
            Dict(:n_success => 0, :n_fail => 11, :custom_totals => [11 // 2, 11 // 3])),
        # no-recovery but override success=false/true and custom_values=[1, 1]/[2, 3]
        (_CycleDecoder([DecodeResult(success=false, custom_values=[1, 1]),
            DecodeResult(success=true, custom_values=[2, 3])]), 12,
            Dict(:n_success => 6, :n_fail => 6, :custom_totals => [18, 24])),
        # no-recovery but override success=false/true and custom_values=[1, [1.]]/[2, [3.]]
        (_CycleDecoder([DecodeResult(success=false, custom_values=[1, [1.]]),
            DecodeResult(success=true, custom_values=[2, [3.]])]), 2,
            Dict(:n_success => 1, :n_fail => 1, :custom_totals => [3, [1., 3.]])),
        # no-recovery but override success=false/true and custom_values=[1,[1 2]]/[2,[3 4]]
        (_CycleDecoder([DecodeResult(success=false, custom_values=[1, [1 2]]),
            DecodeResult(success=true, custom_values=[2, [3 4]])]), 2,
            Dict(:n_success => 1, :n_fail => 1, :custom_totals => [3, [1 2; 3 4]])),
    ]
        data = qec_run(code, error_model, decoder, p; max_runs=max_runs)
        @test issubset(expected, data)
    end
    # test that inconsistent overrides fail
    for (decode_results, expected) in [
        ([DecodeResult(success=true, logical_commutations=nothing),
            DecodeResult(success=true, logical_commutations=[1, 0])], MethodError),
        ([DecodeResult(success=true, logical_commutations=[1, 1, 0]),
            DecodeResult(success=true, logical_commutations=[1, 0])], DimensionMismatch),
        ([DecodeResult(success=true, custom_values=[1, 0]),
            DecodeResult(success=true, custom_values=[1.1, 0.1])], InexactError),
    ]
        decoder = _CycleDecoder(decode_results)
        @test_throws expected qec_run(code, error_model, decoder, p; max_runs=5)
    end
end

@testset "qec_merge-similar" begin
    # similar simulation parameters
    code = FiveQubitCode()
    error_model = DepolarizingErrorModel()
    decoder = NaiveDecoder()
    p = 0.20
    max_runs = 10
    data1 = qec_run(code, error_model, decoder, p; max_runs=max_runs)
    # JSON.print(data1, 4)
    data2 = qec_run(code, error_model, decoder, p; max_runs=max_runs)
    # JSON.print(data2, 4)
    # data1 and data2 have same keys
    @test keys(data1) == keys(data2)
    merged_data_list = qec_merge(data1, data2)
    # JSON.print(merged_data_list, 4)
    # merged_data_list is a vector with one entry
    @test length(merged_data_list) == 1
    merged_data = merged_data_list[1]
    # merged_data drops the error_weight_pvar key
    @test keys(merged_data) ⊆ keys(data1)
    @test setdiff(keys(data1), keys(merged_data)) == Set([:error_weight_pvar])
    # merged_data preserves p and sums max_runs
    @test merged_data[:error_probability] == p
    @test merged_data[:n_run] == max_runs * 2
    # all value types match
    @test all(typeof(merged_data[k]) == typeof(data1[k]) for k in keys(merged_data))
end

@testset "qec_merge-distinct" begin
    # distinct simulation parameters
    code = FiveQubitCode()
    error_model = DepolarizingErrorModel()
    decoder = NaiveDecoder()
    p1, p2 = 0.1, 0.2
    max_runs = 10
    data1 = qec_run(code, error_model, decoder, p1; max_runs=max_runs)
    data2 = qec_run(code, error_model, decoder, p2; max_runs=max_runs)
    # data1 and data2 have same keys
    @test keys(data1) == keys(data2)
    merged_data_list = qec_merge(data1, data2)
    # merged_data_list is a vector with two entries
    @test length(merged_data_list) == 2
    merged_data1 = merged_data_list[1]
    merged_data2 = merged_data_list[2]
    # merged_data preserves all values except :error_weight_pvar
    delete!(data1, :error_weight_pvar)
    @test merged_data1 == data1
    delete!(data2, :error_weight_pvar)
    @test merged_data2 == data2
end

@testset "qec_merge-zero" begin
    # no data to merge
    merged_data_list = qec_merge()
    # merged_data_list is a vector with no entries
    @test length(merged_data_list) == 0
end

@testset "qec_merge-one" begin
    # single simulation
    data = qec_run(FiveQubitCode(), BitFlipErrorModel(), NaiveDecoder(), 0.1; max_runs=10)
    merged_data_list = qec_merge(data)
    # merged_data_list is a vector with one entry
    @test length(merged_data_list) == 1
    merged_data = merged_data_list[1]
    # merged_data preserves all values except :error_weight_pvar
    delete!(data, :error_weight_pvar)
    @test merged_data == data
end

@testset "qec_merge-custom_totals" begin
    # similar simulation parameters
    code = FiveQubitCode()
    error_model = DepolarizingErrorModel()
    decoder = _CycleDecoder([
        DecodeResult(success=false, custom_values=[1, [1]]),
        DecodeResult(success=true, custom_values=[2, [3]])
    ])
    p = 0.20
    max_runs = 3
    data1 = qec_run(code, error_model, decoder, p; max_runs=max_runs)
    # JSON.print(data1, 4)
    data2 = qec_run(code, error_model, decoder, p; max_runs=max_runs)
    # JSON.print(data2, 4)
    merged_data_list = qec_merge(data1, data2)
    # JSON.print(merged_data_list, 4)
    # merged_data_list is a vector with one entry
    @test length(merged_data_list) == 1
    merged_data = merged_data_list[1]
    # check custom_totals as expected
    @test merged_data[:custom_totals] == [9, [1, 3, 1, 3, 1, 3]]
end

@testset "qec_write_read" begin
    # multiple simulation parameters
    code = FiveQubitCode()
    error_model = DepolarizingErrorModel()
    decoder = NaiveDecoder()
    p1, p2 = 0.08, 0.10
    max_runs = 10
    data = []
    push!(data, qec_run(code, error_model, decoder, p1; max_runs=max_runs))
    push!(data, qec_run(code, error_model, decoder, p2; max_runs=max_runs))
    # JSON.print(data, 4)
    # write data
    filename = tempname()
    qec_write(filename, data...)
    # read data
    data_read = qec_read(filename)
    # JSON.print(data_read, 4)
    @test data_read == data
end

@testset "qec_write-existing-file" begin
    # single simulation
    data = qec_run(FiveQubitCode(), BitFlipErrorModel(), NaiveDecoder(), 0.1; max_runs=10)
    # write data
    filename = tempname()
    touch(filename)
    # refuse to overwrite existing file
    @test_throws ErrorException qec_write(filename, data...)
end

@testset "qec_write_read-empty-data" begin
    # write empty data
    filename = tempname()
    qec_write(filename)
    # read empty data
    data_read = qec_read(filename)
    @test data_read == Dict{Symbol,Any}[]
end

@testset "qec_read-invalid-data" begin
    # write invalid data
    filename = tempname()
    touch(filename)
    # fail read invalid data
    @test_throws Exception qec_read(filename)
    # write invalid data
    filename = tempname()
    open(io -> JSON.print(io, Dict(1 => '1')), filename, "w")
    # fail read invalid data
    @test_throws Exception qec_read(filename)
end

@testset "qec_read-qecsim-py-data" begin
    # read data from Python version of Qecsim
    data_py = qec_read(joinpath(@__DIR__, "./App_files/qecsim-py-data.json"))
    # create local data
    data = qec_run(FiveQubitCode(), BitFlipErrorModel(), NaiveDecoder(), 0.1; max_runs=10)
    for d in data_py
        # keys match
        @test keys(d) == keys(data)
        # all value types match
        @test all(typeof(d[k]) == typeof(data[k]) for k in keys(data))
    end
end


text:
"""
`module Constraints`

TODO: write documentation
"""
module Constraints

import JuLIP: Dofs, AbstractConstraint,
         dofs, project!, set_positions!, positions,
         mat, vecs, JVecs,
         AbstractAtoms

export FixedCell

function zeros_free{T}(n::Integer, x::Vector{T}, free::Vector{Int})
   z = zeros(T, n)
   z[free] = x
   return z
end

function insert_free!{T}(p::Array{T}, x::Vector{T}, free::Vector{Int})
   p[free] = x
   return p
end

"""
`FixedCell`: no constraints are placed on the motion of atoms, but the
cell shape is fixed

Constructor:
```julia
FixedCell(at::AbstractAtoms; free=..., clamp=..., mask=...)
```
Set at most one of the kwargs:
* no kwarg: all atoms are free
* `free` : list of free atom indices (not dof indices)
* `clamp` : list of clamped atom indices (not dof indices)
* `mask` : TODO
"""
type FixedCell <: AbstractConstraint
   ifree::Vector{Int}
end

function FixedCell(at::AbstractAtoms; free=nothing, clamp=nothing, mask=nothing)
   if length(find((free != nothing, clamp != nothing, mask != nothing))) > 1
      error("FixedCell: only one of `free`, `clamp`, `mask` may be provided")
   elseif all( (free == nothing, clamp == nothing, mask == nothing) )
      # in this case (default) all atoms are free
      return FixedCell(collect(1:3*length(at)))
   end

   # determine free dof indices
   Nat = length(at)
   if clamp != nothing
      # revert to setting free
      free = setdiff(1:Nat, clamp)
   end
   if free != nothing
      # revert to setting mask
      mask = Matrix{Bool}(3, Nat)
      fill!(mask, false)
      mask[:, free] = true
   end
   return FixedCell(find( mask[:] ))
end

dofs{T}( at::AbstractAtoms, cons::FixedCell, v::JVecs{T}) = mat(v)[cons.ifree]
#    !!!!!!!! this may end up being a problem !!!!!!!
# dofs{T}( at::AbstractAtoms, cons::FixedCell, p::JPts{T}) = mat(p)[cons.ifree]

vecs(cons::FixedCell, at::AbstractAtoms, dofs::Dofs) =
      zeros_free(length(at), dofs, cons.ifree) |> vecs

positions(cons::FixedCell, at::AbstractAtoms, dofs::Dofs) =
      insert_free!(positions(at) |> mat, dofs, cons.ifree) |> vecs

project!(cons::FixedCell, at::AbstractAtoms) = at

# TODO: this is a temporaruy hack, and I think we need to
#       figure out how to do this for more general constraints
project!(cons::FixedCell, A::SparseMatrixCSC) = A[cons.ifree, cons.ifree]

end


text:
extent(grid::AbstractGrid) = (grid.Lx, grid.Ly, grid.Lz)

function get_grid_params(datapath::String)
    file = jldopen(datapath, "r")

    Nx = file["grid/Nx"]
    Ny = file["grid/Ny"]
    Nz = file["grid/Nz"]

    Lx = file["grid/Lx"]
    Ly = file["grid/Ly"]
    Lz = file["grid/Lz"]

    close(file)
    return Nx, Ny, Nz, Lx, Ly, Lz
end

get_grid_params(grid::AbstractGrid) = tuple(size(grid)..., extent(grid)...)


text:
module CompressHashDisplace

using MurmurHash3
export FrozenUnsafeDict, FrozenDict

mmhash(str::String) = last(mmhash128_a(sizeof(str), pointer(str), 0%UInt32))
mmhash(str::String, d::UInt32) = last(mmhash128_a(sizeof(str), pointer(str), d))
mmhash(key) = hash(key)
mmhash(key, d) = hash(key, UInt64(d))

struct FrozenUnsafeDict{V}
    G::Vector{Int}
    values::Vector{V}
    sz::UInt64
end

optimal_size(dict) = 1 << (leading_zeros(0) - leading_zeros(length(dict)))

"""
    FrozenUnsafeDict(dict)

Computes a perfect hash table using given `dict`. `FrozenUnsafeDict` do not hold key values,
and wouldn't do any validations during lookup procedure.
"""
function FrozenUnsafeDict(dict::Dict{K, V}) where {K, V}
    sz = optimal_size(dict)
    szmask = (sz - 1) % UInt64

    # Step 1: Place all of the keys into buckets
    buckets = [K[] for _ in 1:sz]
    G = zeros(Int, sz)
    # values = Vector{Union{Int, Nothing}}(undef, sz)
    values = Vector{V}(undef, sz)
    keyset = BitArray(undef, sz)
    for i in 1:sz
        keyset[i] = false
    end

    for key in keys(dict)
        idx = mmhash(key, 0 % UInt32) & szmask + 1
        push!(buckets[idx], key)
    end

    # Step 2: Sort the buckets and process the ones with the most items first.
    sort!(buckets, by = length, rev = true)
    slots = Vector{Int}()
    sizehint!(slots, length(first(buckets)))
    b = 1
    while b <= sz
        bucket = buckets[b]
        length(bucket) <= 1 && break

        d = 1 % UInt32
        item = 1

        # Repeatedly try different values of d until we find a hash function
        # that places all items in the bucket into free slots
        while item <= length(bucket)
            # slot = mmhash(bucket[item], d) & szmask + 1
            slot = (mmhash(bucket[item], 0%UInt32) >> d) & szmask + 1
            if keyset[slot] || (slot in slots)
                d += 1 % UInt32
                item = 1
                empty!(slots)
            else
                push!(slots, slot)
                item += 1
            end
        end

        # We already know this value, why do we calculate it again??
        G[mmhash(first(bucket), 0 % UInt32) & szmask + 1] = d
        for i in 1:length(bucket)
            values[slots[i]] = dict[bucket[i]]
            keyset[slots[i]] = true
        end
        empty!(slots)
        b += 1
    end

    idx = 1
    for b2 in b:sz
        bucket = buckets[b2]
        isempty(bucket) && break
        bucket = first(bucket)
        while keyset[idx]
            idx += 1
        end

        # why we calculate it second time?? It makes no sense at all
        G[mmhash(bucket, 0 % UInt32) & szmask + 1] = -idx
        values[idx] = dict[bucket]
        idx += 1
    end

    return FrozenUnsafeDict{V}(G, values, sz)
end

# Look up a value in the hash table, defined by G and V.
function Base.:getindex(FD::FrozenUnsafeDict{V}, key) where V
    szmask = (FD.sz - 1) % UInt64
    idx = mmhash(key, 0%UInt32)
    @inbounds d = FD.G[idx & szmask + 1%UInt64]
    if d < 0
        @inbounds return FD.values[-d]
    else
        # @inbounds return FD.values[mmhash(key, d % UInt32) & szmask + 1%UInt64]
        @inbounds return FD.values[(idx >> d) & szmask + 1%UInt64]
    end
end

struct FrozenDict{K, V} <: AbstractDict{K, V}
    ks::Vector{K}
    G::Vector{Int}
    values::Vector{V}
    mask::BitArray{1}
    sz::UInt64
    length::Int
end

"""
    FrozenDict(dict)

Computes a perfect hash table using given `dict`. It returns read-only dictionary.
"""
function FrozenDict(dict::Dict{K, V}) where {K, V}
    sz = optimal_size(dict)
    szmask = (sz - 1) % UInt64
    ks = Vector{K}(undef, sz)

    # Step 1: Place all of the keys into buckets
    buckets = [K[] for _ in 1:sz]
    G = zeros(Int, sz)
    # values = Vector{Union{Int, Nothing}}(undef, sz)
    values = Vector{V}(undef, sz)
    keyset = BitArray(undef, sz)
    for i in 1:sz
        keyset[i] = false
    end

    for key in keys(dict)
        idx = (mmhash(key, 0 % UInt32) & szmask) % Int + 1
        push!(buckets[idx], key)
    end

    # Step 2: Sort the buckets and process the ones with the most items first.
    sort!(buckets, by = length, rev = true)
    slots = Vector{Int}()
    sizehint!(slots, length(first(buckets)))
    b = 1
    while b <= sz
        bucket = buckets[b]
        length(bucket) <= 1 && break

        resize!(slots, length(bucket))
        d = 1 % UInt32
        item = 1

        # Repeatedly try different values of d until we find a hash function
        # that places all items in the bucket into free slots
        while item <= length(bucket)
            slot = (mmhash(bucket[item], 0%UInt32) >> d) & szmask + 1
            slot == 0 && throw(KeyError(bucket[item]))
            if keyset[slot] || (slot in slots)
                d += 1 % UInt32
                item = 1
            else
               slots[item] = slot
                item += 1
            end
        end

        # We already know this value, why do we calculate it again??
        G[mmhash(first(bucket), 0 % UInt32) & szmask + 1] = d
        for i in 1:length(bucket)
            values[slots[i]] = dict[bucket[i]]
            keyset[slots[i]] = true
            ks[slots[i]] = bucket[i]
        end
        b += 1
    end

    idx = 1
    for b2 in b:sz
        bucket = buckets[b2]
        isempty(bucket) && break
        key = first(bucket)
        while keyset[idx]
            idx += 1
        end

        # why we calculate it second time?? It makes no sense at all
        G[mmhash(key, 0 % UInt32) & szmask + 1] = -idx
        values[idx] = dict[key]
        ks[idx] = key
        keyset[idx] = true
        idx += 1
    end

    return FrozenDict{K, V}(ks, G, values, keyset, sz, length(dict))
end

# Look up a value in the hash table, defined by G and V.
function Base.:getindex(FD::FrozenDict{K, V}, key) where {K, V}
    szmask = UInt64(FD.sz - 1)
    idx = mmhash(key, 0%UInt32)
    @inbounds d = FD.G[idx & szmask + 1%UInt32]
    idx = d < 0 ? -d%UInt64 : (idx >> d) & szmask + 1%UInt64

    @inbounds if FD.ks[idx] === key || isequal(FD.ks[idx], key)
        return FD.values[idx]
    else
        throw(KeyError(key))
    end
end

Base.length(FD::FrozenDict) = FD.length

function Base.iterate(FD::FrozenDict, state=0)
    for i in (state + 1):FD.sz
        if FD.mask[i]
            return (FD.ks[i], FD.values[i]), i
        end
    end
    return nothing
end

end # module


text:
# NOTE setup too for `cthulhu_ast`
include("sandbox.jl")
using .CthulhuTestSandbox
Revise.track(CthulhuTestSandbox, normpath(@__DIR__, "sandbox.jl"))

@testset "printer test" begin
    _, src, infos, mi, rt, slottypes = process(testf_revise);
    tf = (true, false)

    @testset "codeview: $codeview" for codeview in Cthulhu.CODEVIEWS
        if !@isdefined(Revise)
            codeview == Cthulhu.cthulhu_ast && continue
        end
        @testset "optimize: $optimize" for optimize in tf
            @testset "debuginfo: $debuginfo" for debuginfo in instances(Cthulhu.DInfo.DebugInfo)
                params = Cthulhu.current_params()
                config = Cthulhu.CONFIG

                io = IOBuffer()
                codeview(io, mi, optimize, debuginfo, params, config)
                @test !isempty(String(take!(io))) # just check it works
            end
        end
    end

    @testset "debuginfo: $debuginfo" for debuginfo in instances(Cthulhu.DInfo.DebugInfo)
        @testset "iswarn: $iswarn" for iswarn in tf
            @testset "hide_type_stable: $hide_type_stable" for hide_type_stable in tf
                @testset "inline_cost: $inline_cost" for inline_cost in tf
                    @testset "type_annotations: $type_annotations" for type_annotations in tf
                        io = IOBuffer()
                        Cthulhu.cthulhu_typed(io, debuginfo,
                            src, rt, mi;
                            iswarn, hide_type_stable, inline_cost, type_annotations)
                        @test !isempty(String(take!(io))) # just check it works
                    end
                end
            end
        end
    end
end

using PerformanceTestTools
PerformanceTestTools.@include("irshow.jl")

@testset "hide type-stable statements" begin
    let # optimize code
        _, src, infos, mi, rt, slottypes = @eval Module() begin
            const globalvar = Ref(42)
            $process() do
                a = sin(globalvar[])
                b = sin(undefvar)
                return (a, b)
            end
        end
        function prints(; kwargs...)
            io = IOBuffer()
            Cthulhu.cthulhu_typed(io, :none, src, rt, mi; kwargs...)
            return String(take!(io))
        end

        let # by default, should print every statement
            s = prints()
            @test occursin("globalvar", s)
            @test occursin("undefvar", s)
        end
        let # should omit type stable statements
            s = prints(; hide_type_stable=true)
            @test !occursin("globalvar", s)
            @test occursin("undefvar", s)
        end
    end

    let # unoptimize code
        _, src, infos, mi, rt, slottypes = @eval Module() begin
            const globalvar = Ref(42)
            $process(; optimize=false) do
                a = sin(globalvar[])
                b = sin(undefvar)
                return (a, b)
            end
        end
        function prints(; kwargs...)
            io = IOBuffer()
            Cthulhu.cthulhu_typed(io, :none, src, rt, mi; kwargs...)
            return String(take!(io))
        end

        let # by default, should print every statement
            s = prints()
            @test occursin("globalvar", s)
            @test occursin("undefvar", s)
        end
        let # should omit type stable statements
            s = prints(; hide_type_stable=true)
            @test !occursin("globalvar", s)
            @test occursin("undefvar", s)
        end

        # should work for warn mode
        let
            s = prints(; iswarn=true)
            @test occursin("globalvar", s)
            @test occursin("undefvar", s)
        end
        let
            s = prints(; iswarn=true, hide_type_stable=true)
            @test !occursin("globalvar", s)
            @test occursin("undefvar", s)
        end
    end
end


text:
export scaleby, inverse, from_to, gettransform!

function Transform(;scale=1, rotation=rad(0), translation=(0, 0))
    Transform(scale, rotation, translation)
end

Base.Broadcast.broadcastable(t::Transform) = Ref(t)

function rotmat(ang::Angle)
    SMatrix{2, 2}(cos(ang.rad), sin(ang.rad), -sin(ang.rad), cos(ang.rad))
end

function Base.:*(t1::Transform, t2::Transform)
    rmat = rotmat(t1.rotation)
    scale = t1.scale * t2.scale
    rotation = t2.rotation + t1.rotation
    translation = t1.scale * (rmat * t2.translation) + t1.translation
    Transform(scale, rotation, translation)
end

function Base.:*(t::Transform, r::Real)
    t.scale * r
end

function inverse(t::Transform)
    scale = 1 / t.scale
    rotation = -1 * t.rotation
    rmat = rotmat(rotation)

    translation = -scale * rmat * t.translation
    Transform(scale, rotation, translation)
end

function from_to(t1::Transform, t2::Transform)
    t2 * inverse(t1)
end

function scaleby(t::Transform, s::Real)
    Transform(t.scale * s, t.rotation, t.translation)
end

function Base.:*(t::Transform, p::Point)
    rmat = rotmat(t.rotation)
    Point(rmat * t.scale * p.xy + t.translation)
end

function Base.:*(t::Transform, l::Line)
    Line(t * l.from, t * l.to)
end

function Base.:*(t::Transform, c::Circle)
    Circle(t * c.center, t.scale * c.radius)
end

function Base.:*(t::Transform, r::Rect)
    Rect(t * r.center, t.scale * r.width, t.scale * r.height, r.angle + t.rotation)
end

function Base.:*(t::Transform, b::Bezier)
    Bezier((t .* (b.from, b.c1, b.c2, b.to))...)
end

function Base.:*(t::Transform, x::T) where T <: Union{Move, Lineto}
    T(t * x.p)
end

function Base.:*(t::Transform, x::Close)
    x
end

function Base.:*(t::Transform, c::CurveTo)
    CurveTo(t * c.c1, t * c.c2, t * c.p)
end

function Base.:*(t::Transform, bp::Path)
    Path([t * c for c in bp.commands])
end

function Base.:*(t::Transform, p::Polygon)
    Polygon(t .* p.points)
end

function Base.:*(tr::Transform, t::Txt)
    Txt(
        tr * t.pos,
        t.text,
        tr.scale * t.size,
        t.halign,
        t.valign,
        tr.rotation + t.angle,
        t.font,
        tr.scale * t.extent
    )
end

function Base.:*(t::Transform, a::Arc)
    Arc(
        t * a.center,
        t.scale * a.radius,
        t.rotation + a.start_angle,
        t.rotation + a.end_angle,
    )
end


text:
"Distance between two values"
function dist end

# zt: these choices seem a little adhoc.  Should I specialise by type

@inline dist(x::Real, y::Real) = (xy = (x - y); xy * xy)
# @inline dist(x::Vector{<:Real}, y::Vector{<:Real}) = norm(x - y)
@inline dist(x::AbstractVector{<:Real}, y::AbstractVector{<:Real}) = sum(dist.(x, y))
@inline dist(x::NTuple{N,<: Real}, y::NTuple{N,<:Real}) where N = sum(dist.(x, y))
@inline dist(x::AbstractArray{<:Real}, y::AbstractArray{<:Real}) = norm(x[:] - y[:])

"Distance from x to [a, b]"
function bound_loss(x, a, b)
  # @pre b >= a
  if x < a
    a - x
  elseif x > b
    x - b
  else
    zero(x)
  end
end

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


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

name = "Capstone"
version = v"4.0.2"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/aquynh/capstone/archive/4.0.2.tar.gz", "7c81d798022f81e7507f1a60d6817f63aa76e489aa4e7055255f21a22f5e526a")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd capstone-4.0.2/
mkdir build
cd build/
cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} -DCMAKE_BUILD_TYPE=Release ..
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

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

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

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

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
using Video4Linux
using Colors, ColorTypes, FixedPointNumbers
using Test

@testset "Video4Linux" begin

    #conversions
    testu8 = zeros(UInt8,5)
    for i = 1:5
        testu8[i] = 0x01 << (i-1)
    end
    @test convertY10BtoU16(testu8) == [0x0004, 0x0020, 0x0102, 0x0010]

    # Reference
    A = [0.0  0.0 32.0;
       128.0 64.0 96.0;
        64.0  0.0 32.0;
       192.0 64.0 96.0]

    testu8 = zeros(UInt8,8)
    y = 0x00
    for i = 2:2:8
        testu8[i] = UInt8(y)
        y += 0x40
    end
    y = 0x00
    for i = 1:2:8
        testu8[i] = UInt8(y)
        y += 0x20
    end

    im = convertUYVYtoYCbCr(testu8, 2, 2)

    @test im == reshape(YCbCr.(A[:,1], A[:,2], A[:,3]),(2,2))

    #using convert
    frame = Video4Linux.RawFrame(testu8, 2, 2, :UYVY)
    @test Video4Linux.convert(YCbCr, frame) == reshape(YCbCr.(A[:,1], A[:,2], A[:,3]),(2,2))
    @test Video4Linux.convert(Gray, frame)  == reinterpret(Gray{N0f8},UInt8.(reshape((A[:,1]),(2,2))))

    testu8 = zeros(UInt8,8)
    y = 0x00
    for i = 1:2:8
        testu8[i] = UInt8(y)
        y += 0x40
    end
    y = 0x00
    for i = 2:2:8
        testu8[i] = UInt8(y)
        y += 0x20
    end
    im = convertYUYVtoYCbCr(testu8, 2, 2)
    @test im == reshape(YCbCr.(A[:,1], A[:,2], A[:,3]),(2,2))

    im = Video4Linux.convertYUYVtoArray(testu8, 2, 2)
    @test im == reshape(A,(2,2,3))

    #using convert
    frame = Video4Linux.RawFrame(testu8, 2, 2, :YUYV)
    @test Video4Linux.convert(YCbCr, frame) == reshape(YCbCr.(A[:,1], A[:,2], A[:,3]),(2,2))
    @test Video4Linux.convert(Gray, frame)  == reinterpret(Gray{N0f8},UInt8.(reshape((A[:,1]),(2,2))))

    #wrapper
    @test set_io_method() == nothing
    @test set_io_method(Int32(0)) == nothing
    @test set_io_method(Video4Linux.IO_METHOD_MMAP) == nothing

    #for now these should just not break julia
    fid = open_device("/dev/video0")
    init_device(Int32(0))
    start_capturing(Int32(0))
    stop_capturing(Int32(0))
    uninit_device(Int32(0))
    close_device(Int32(0))
    temptestvalue = 0
    try
        # this should through error to protect julia from crashing
        mainloop( Int32(-1), 1 )
    catch
        temptestvalue = 1
    end
    @test 1 == temptestvalue

    # some more api
    d1 = Video4Linux.UYVY(640,480)
    d2 = Video4Linux.UYVYonlyY(640,480)
    d3 = Video4Linux.YUYV(640,480)
    d4 = Video4Linux.YUYVonlyY(640,480)
    d5 = Video4Linux.Y10B(640,480)

    #also just test not te break julia
    vidchan = Channel((c::Channel) -> videoproducer(c, d1))
    vidchan = Channel((c::Channel) -> videoproducer(c, d2, devicename = "/dev/video0"))
    vidchan = Channel((c::Channel) -> videoproducer(c, d3, devicename = "/dev/video0", iomethod = Video4Linux.IO_METHOD_READ))

end


text:
mutable struct SchellingAgent <: AbstractAgent
    id::Int # The identifier number of the agent
    pos::Tuple{Int,Int} # The x, y location of the agent
    mood::Bool # whether the agent is happy in its node. (true = happy)
    group::Int # The group of the agent,
    # determines mood as it interacts with neighbors
end

function schelling(; numagents = 320, griddims = (20, 20), min_to_be_happy = 3)
    space = Space(griddims, moore = true)
    properties = Dict(:min_to_be_happy => min_to_be_happy)
    schelling =
        ABM(SchellingAgent, space; properties = properties, scheduler = random_activation)
    for n in 1:numagents
        agent = SchellingAgent(n, (1, 1), false, n < numagents / 2 ? 1 : 2)
        add_agent_single!(agent, schelling)
    end
    return schelling
end

@testset "Agent distribution" begin
    model = schelling()
    groupcolor(a) = a.group == 1 ? :blue : :orange
    groupmarker(a) = a.group == 1 ? :circle : :square
    p = plotabm(model; ac = groupcolor, am = groupmarker)
    @test typeof(p) <: AgentsPlots.AbstractPlot
end


text:
# This test case is copied from http://hakank.org/julia/constraints/monks_and_doors.jl
# which is a website by Håkan Kjellerstrand
# please check it out :)

@testset "Monks & Doors" begin
    all_solutions = true
    model = Model(optimizer_with_attributes(CS.Optimizer,
        "all_solutions"=> all_solutions,
        # "all_optimal_solutions"=>all_solutions,
        "logging"=>[],

        "traverse_strategy"=>:BFS,
        "branch_split"=>:InHalf,

        "branch_strategy" => :IMPS, # default
    ))

    num_doors = 4
    num_monks = 8
    @variable(model, doors[1:num_doors], Bin)
    da,db,dc,dd = doors
    door_names = ["A","B","C","D"]

    @variable(model, monks[1:num_monks], Bin)
    m1,m2,m3,m4,m5,m6,m7,m8 = monks

    # Monk 1: Door A is the exit.
    # M1 #= A (Picat constraint)
    @constraint(model, m1 == da)

    #  Monk 2: At least one of the doors B and C is the exit.
    # M2 #= 1 #<=> (B #\/ C)
    @constraint(model, m2 := { db == 1 || dc == 1})

    #  Monk 3: Monk 1 and Monk 2 are telling the truth.
    # M3 #= 1 #<=> (M1 #/\ M2)
    @constraint(model, m3 := { m1 == 1 && m2 == 1})

    #  Monk 4: Doors A and B are both exits.
    # M4 #= 1 #<=> (A #/\ B)
    @constraint(model, m4 := { da == 1 && db == 1})

    #  Monk 5: Doors A and C are both exits.
    # M5 #= 1 #<=> (A #/\ C)
    @constraint(model, m5 := { da == 1 && dc == 1})

    #  Monk 6: Either Monk 4 or Monk 5 is telling the truth.
    # M6 #= 1 #<=> (M4 #\/ M5)
    @constraint(model, m6 := { m4 == 1|| m5 == 1})

    #  Monk 7: If Monk 3 is telling the truth, so is Monk 6.
    # M7 #= 1 #<=> (M3 #=> M6)
    @constraint(model, m7 := { m3 => {m6 == 1}})

    #  Monk 8: If Monk 7 and Monk 8 are telling the truth, so is Monk 1.
    # M8 #= 1 #<=> ((M7 #/\ M8) #=> (M1))
    b1 = @variable(model, binary=true)
    @constraint(model, b1 := {m7 == 1 && m8 == 1})
    @constraint(model, m8 := {b1 => {m1 == 1}})

    # Exactly one door is an exit.
    # (A + B + C + D) #= 1
    @constraint(model, da + db + dc + dd == 1)

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    # println("status:$status")
    num_sols = 0
    @test status == MOI.OPTIMAL
    num_sols = MOI.get(model, MOI.ResultCount())
    @test num_sols == 1
    doors_val = convert.(Integer,JuMP.value.(doors))
    monks_val = convert.(Integer,JuMP.value.(monks))
    # exit door is A
    @test doors_val[1] == 1
    @test sum(doors_val) == 1
    @test sum(monks_val) == 3
    # monks 1,7,8 are the only ones who tell the truth
    @test monks_val[1] == monks_val[7] == monks_val[8] == 1
end