text:
struct CSTShapeFunction{N,T <: AbstractFloat}
    bino::NTuple{N,Int}
    CSTShapeFunction(N::Int, ::Type{T}) where T = new{N,T}(tuple((binomial(N, i) for i in 1:N)...))
end

struct CSTClassFunction{N1,N2,T}
    CSTClassFunction(N1::T, N2::T) where T = new{N1,N2,T}()
end

shape_f(s::CSTShapeFunction{N,T}, i::Int, x::T) where {N,T} =
if i == 0
    (1 - x)^N
elseif i == N
    x^N
else
    s.bino[i] * x^i * (1 - x)^(N - i)
end

shape_d(s::CSTShapeFunction{N,T}, i::Int, x::T) where {N,T} =
if i == 0
    -N * (1 - x)^(N - 1)
elseif i == N
    N * x^(N - 1)
else
    s.bino[i] * (i * x^(i - 1) * (1 - x)^(N - i) - (N - i) * x^i * (1 - x)^(N - i - 1))
end

class_f(::CSTClassFunction{N1,N2,T},x::T) where {N1,N2,T} = x^N1 * (1 - x)^N2

class_d(::CSTClassFunction{N1,N2,T},x::T) where {N1,N2,T} = N1 * x^(N1 - 1) * (1 - x)^N2 - N2 * x^N1 * (1 - x)^(N2 - 1)

function cst_preprocess(xx::AbstractVector{T}, yy::AbstractVector{T}) where T
    xte = (xx[1] + xx[end]) / 2
    if yy[end-1] > yy[2]
        reverse!(xx)
        reverse!(yy)
    end
    imin = argmin(xx)
    x0 = xx[imin]
    y0 = yy[imin]
    cc = xte - x0
    xu_ = @view(xx[imin:-1:1])
    yu_ = @view(yy[imin:-1:1])
    xl_ = @view(xx[imin:end] )
    yl_ = @view(yy[imin:end] )
    xu = (xu_ .- x0) ./ cc
    yu = (yu_ .- y0) ./ cc
    xl = (xl_ .- x0) ./ cc
    yl = (yl_ .- y0) ./ cc
    return xu, yu, xl, yl
end

function cst_matrix(
    class::CSTClassFunction{N1,N2,T},
    shape::CSTShapeFunction{N,T},
    xu::V, yu::V, xl::V, yl::V)　where {T,V <: AbstractVector{T},N1,N2,N}
    NP = length(xu) + length(xl) - 2
    A = zeros(T, NP, 2N + 3)
    b = vcat(yu[2:end], yl[2:end])
    ip = 1
    @inbounds for x in @view xu[2:end]
        c = class_f(class, x)
        for i in 1:N
            A[ip, i] = c * shape_f(shape, i, x) # cstShapeF(N, i, x)
        end
        A[ip, 2N + 1] = c * shape_f(shape, 0, x) # a_le
        A[ip, 2N + 2] = x # y_tu
        ip += 1
    end
    @inbounds for x in @view xl[2:end]
        c = class_f(class, x)
        for i in 1:N
            A[ip, i + N] = c * shape_f(shape, i, x)
        end
        A[ip, 2N + 1] = -c * shape_f(shape, 0, x) # -a_le
        A[ip, 2N + 3] = x # y_tl
        ip += 1
    end
    return A, b
end

function cst_fit(
    class::CSTClassFunction{N1,N2},
    shape::CSTShapeFunction{N,T},
    xu::AbstractVector{T}, yu::AbstractVector{T},
    xl::AbstractVector{T}, yl::AbstractVector{T} ) where {T,N,N1,N2}
    A, b = cst_matrix(class, shape, xu, yu, xl, yl)
    U, S, V = svd(A)
    params = V * ((transpose(U) * b) ./ S)
    au = params[1:N]
    al = params[N + 1:2N]
    a0 = params[2N + 1]
    yu = params[2N + 2]
    yl = params[2N + 3]
    return au, al, a0, yu, yl
end


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

name = "OpenSpiel"
version = v"0.1.1"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/findmyway/open_spiel/archive/v0.1.1.tar.gz", "8791ab9902aa46fd8263b79ea789f34d17e3a38dbc9e7ecbd46c5423922f0124"),
    ArchiveSource("https://github.com/findmyway/dds/archive/v0.1.0.tar.gz", "81070b8e96779b5b2303185642753013aa874bffbd58b9cc599204aee064292d"),
    ArchiveSource("https://github.com/findmyway/abseil-cpp/archive/v0.1.0.tar.gz", "7b612c1fed278250b5d1a4e29ddb410145b26a0e7c781c1ca4ac03d092179202"),
    ArchiveSource("https://github.com/findmyway/hanabi-learning-environment/archive/v0.1.0.tar.gz", "6126936fd13a95f8cadeacaa69dfb38a960eaf3bd588aacc8893a6e07e4791a3"),
    ArchiveSource("https://github.com/findmyway/project_acpc_server/archive/v0.1.0.tar.gz", "e29f969dd62ba354b7019cae3f7f1dbfbd9a744687ea4a8f7494c2bb1ee87382"),
    ArchiveSource("https://julialang-s3.julialang.org/bin/linux/armv7l/1.3/julia-1.3.1-linux-armv7l.tar.gz", "965c8fab2214f8ce1b3d449d088561a6de61be42543b48c3bbadaed5b02bf824"; unpack_target="julia-arm-linux-gnueabihf"),
    ArchiveSource("https://julialang-s3.julialang.org/bin/linux/x64/1.3/julia-1.3.1-linux-x86_64.tar.gz", "faa707c8343780a6fe5eaf13490355e8190acf8e2c189b9e7ecbddb0fa2643ad"; unpack_target="julia-x86_64-linux-gnu"),
    ArchiveSource("https://github.com/Gnimuc/JuliaBuilder/releases/download/v1.3.0/julia-1.3.0-x86_64-apple-darwin14.tar.gz", "f2e5359f03314656c06e2a0a28a497f62e78f027dbe7f5155a5710b4914439b1"; unpack_target="julia-x86_64-apple-darwin14"),
    ArchiveSource("https://github.com/Gnimuc/JuliaBuilder/releases/download/v1.3.0/julia-1.3.0-x86_64-w64-mingw32.tar.gz", "c7b2db68156150d0e882e98e39269301d7bf56660f4fc2e38ed2734a7a8d1551"; unpack_target="julia-x86_64-w64-mingw32"),
]

# Bash recipe for building across all platforms
script = raw"""

case "$target" in
	arm-linux-gnueabihf|x86_64-linux-gnu)
        Julia_PREFIX=${WORKSPACE}/srcdir/julia-$target/julia-1.3.1
        ;;
    x86_64-apple-darwin14|x86_64-w64-mingw32)
        Julia_PREFIX=${WORKSPACE}/srcdir/julia-$target/juliabin
        ;;
esac

mv open_spiel-0.1.1 open_spiel
mv abseil-cpp-0.1.0/ open_spiel/open_spiel/abseil-cpp
mv dds-0.1.0 open_spiel/open_spiel/games/bridge/double_dummy_solver
mv hanabi-learning-environment-0.1.0 open_spiel/open_spiel/games/hanabi/hanabi-learning-environment
mv project_acpc_server-0.1.0 open_spiel/open_spiel/games/universal_poker/acpc

mkdir open_spiel/build
cd open_spiel/build
BUILD_WITH_PYTHON=OFF BUILD_WITH_JULIA=ON BUILD_WITH_HANABI=ON BUILD_WITH_ACPC=OFF cmake -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} -DJulia_PREFIX=$Julia_PREFIX -DCMAKE_PREFIX_PATH=$prefix/destdir/lib/cmake/JlCxx -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_FIND_ROOT_PATH=${prefix} ../open_spiel/
make -j${nproc}
make install
install_license ${WORKSPACE}/srcdir/open_spiel/LICENSE
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:x86_64, libc=:glibc, compiler_abi=CompilerABI(;cxxstring_abi=:cxx11)),
]

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

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

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

ext: jl

max_stars_repo_name: jbrea/Yggdrasil

max_issues_repo_name: jbrea/Yggdrasil

max_forks_repo_name: jbrea/Yggdrasil


text:
import Base.Test

using DebugPrintln
nodebug()

type FunctionCountingAstWalker <: AstWalker
	count::Int
end
FunctionCountingAstWalker() = FunctionCountingAstWalker(0)

# We count number of visits to function signatures...
function visit_expr_function_sig(fc::FunctionCountingAstWalker, sig, funcex)
	fc.count += 1
	sig
end

astfunc1 = :(
	function f(a)
		a+1
	end
)

fc = FunctionCountingAstWalker()
visit(fc, astfunc1)
@assert fc.count == 1

astfunc2 = quote
	function f1(a)
		a+1
	end
	function f2(b)
		b+2
	end
end

fc2 = FunctionCountingAstWalker()
visit(fc2, astfunc2)
@assert fc2.count == 2

type LineCountingAstWalker <: AstWalker
	count::Int
end
LineCountingAstWalker() = LineCountingAstWalker(0)

function visit_line_node(lc::LineCountingAstWalker, l::LineNumberNode)
	lc.count += 1
	l
end

lc1 = LineCountingAstWalker()
visit(lc1, astfunc1)
@assert lc1.count == 1 # Only actual lines in the body of the funcs are counted

lc2 = LineCountingAstWalker()
visit(lc2, astfunc2)
@assert lc2.count == 4 # Only actual lines in the body of the funcs are counted


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module Base64

using Base: require_one_based_indexing

export
    Base64EncodePipe,
    base64encode,
    Base64DecodePipe,
    base64decode,
    stringmime

# Base64EncodePipe is a pipe-like IO object, which converts into base64 data
# sent to a stream. (You must close the pipe to complete the encode, separate
# from closing the target stream).  We also have a function base64encode(f,
# args...) which works like sprint except that it produces base64-encoded data,
# along with base64encode(args...)  which is equivalent to base64encode(write,
# args...), to return base64 strings.  A Base64DecodePipe object can be used to
# decode base64-encoded data read from a stream , while function base64decode is
# useful for decoding strings

include("buffer.jl")
include("encode.jl")
include("decode.jl")

"""
    stringmime(mime, x; context=nothing)

Returns an `AbstractString` containing the representation of `x` in the
requested `mime` type. This is similar to [`repr(mime, x)`](@ref) except
that binary data is base64-encoded as an ASCII string.

The optional keyword argument `context` can be set to `:key=>value` pair
or an `IO` or [`IOContext`](@ref) object whose attributes are used for the I/O
stream passed to `show`.
"""
stringmime(m::MIME, x; context=nothing) = istextmime(m) ? Base.Multimedia._textrepr(m, x, context) : _binstringmime(m, x, context)
stringmime(m::AbstractString, x; context=nothing) = stringmime(MIME(m), x; context=context)

_binstringmime(m::MIME, x, context) = Base64.base64encode(show, m, x; context=context)
_binstringmime(m::MIME, x::Vector{UInt8}, context) = Base64.base64encode(write, x; context=context)

end

max_stars_repo_name: EricForgy/julia

max_issues_repo_name: EricForgy/julia

max_forks_repo_name: EricForgy/julia


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct TransportProtocol2 <: SwaggerModel

    function TransportProtocol2(;)
        o = new()
        o
    end
end # type TransportProtocol2

const _property_map_TransportProtocol2 = Dict{Symbol,Symbol}()
const _property_types_TransportProtocol2 = Dict{Symbol,String}()
Base.propertynames(::Type{ TransportProtocol2 }) = collect(keys(_property_map_TransportProtocol2))
Swagger.property_type(::Type{ TransportProtocol2 }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_TransportProtocol2[name]))}
Swagger.field_name(::Type{ TransportProtocol2 }, property_name::Symbol) =  _property_map_TransportProtocol2[property_name]

function check_required(o::TransportProtocol2)
    true
end

function validate_property(::Type{ TransportProtocol2 }, name::Symbol, val)
end


text:
using Test
function test_state(v, acc, s=101.0)
    ego = Vehicle(vehicle_state(s, main_lane(mdp.env), v, mdp.env.roadway), VehicleDef(), EGO_ID)
    veh1 = Vehicle(vehicle_state(70.0, main_lane(mdp.env), 20.0, mdp.env.roadway), VehicleDef(), EGO_ID + 1)
    veh2 = Vehicle(vehicle_state(150.0, main_lane(mdp.env), 20.0, mdp.env.roadway), VehicleDef(), EGO_ID + 2)
    scene = Scene()
    push!.(Ref(scene), [ego, veh1, veh2])

    return AugScene(scene, (acc=acc,))
end

AutoViz.render(s.scene, mdp.env.roadway, [IDOverlay()], cam=CarFollowCamera(1, 10.0))

s = test_state(20.0, 0.0)
@test speed_limit_actions(mdp, s) == actions(mdp)
s = test_state(25.0, 0.0)
@test speed_limit_actions(mdp, s) == [1,2,3,4,7]
s = test_state(26.0, 0.0)
@test speed_limit_actions(mdp, s) == [1]

# test acceleration limit
s = test_state(20.0, 0.0)
@test acceleration_limit_actions(mdp, s) == [2,3,4,5,6,7]
s = test_state(20.0, 2.0)
acceleration_limit_actions(mdp, s) == [2,3,4,7]
s = test_state(20.0, -2.0)
acceleration_limit_actions(mdp, s) == [4,5,6,7]

# test idm_limit_actions
s = test_state(20.0, 0.0)
idm_limit_actions(mdp, s) == [1,2,3,4,5,6,7]
s = test_state(20.0, 0.0, 120.0)
idm_limit_actions(mdp, s) == [1]
s = test_state(20.0, 0.0, 110.0)
idm_limit_actions(mdp, s) == [1,2,3,4,7]

# test joint mask
s = test_state(20.0, 0.0)
safe_actions(mdp, s) == [2,3,4,5,6,7]
s = test_state(25.0, 0.0, 110.0)
safe_actions(mdp, s) == [1]
s = test_state(20.0, 0.0, 110.0)
safe_actions(mdp, s) == [2,3,4,7]


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

MOIU.@model(SDModelData,
            (),
            (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan),
            (MOI.Zeros, MOI.Nonnegatives, MOI.Nonpositives,
             MOI.PositiveSemidefiniteConeTriangle),
            (),
            (MOI.SingleVariable,),
            (MOI.ScalarAffineFunction,),
            (MOI.VectorOfVariables,),
            (MOI.VectorAffineFunction,))
# UniversalFallback is needed for starting values, even if they are ignored by SDPA
const cache = MOIU.UniversalFallback(SDModelData{Float64}())
const optimizer = MOIU.CachingOptimizer(cache,
                                        SDPA.Optimizer())
# test 1e-3 because of rsoc3 test, otherwise, 1e-5 is enough
const config = MOIT.TestConfig(atol=1e-3, rtol=1e-3)

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

@testset "Linear tests" begin
    MOIT.contlineartest(MOIB.SplitInterval{Float64}(optimizer), config,
                        ["linear12"])
end
@testset "Conic tests" begin
    MOIT.contconictest(MOIB.RootDet{Float64}(MOIB.GeoMean{Float64}(MOIB.RSOCtoPSD{Float64}(MOIB.SOCtoPSD{Float64}(optimizer)))),
                       config,
                       ["lin3", "soc3", "psds", "rootdets", "logdet", "exp"])
end


text:
using HDF5
using JLD
using Base.Test

# Define variables of different types
x = 3.7
A = reshape(1:15, 3, 5)
Aarray = Vector{Float64}[[1.2,1.3],[2.2,2.3,2.4]]
str = "Hello"
stringsA = ASCIIString["It", "was", "a", "dark", "and", "stormy", "night"]
stringsU = UTF8String["It", "was", "a", "dark", "and", "stormy", "night"]
strings16 = convert(Array{UTF16String}, stringsA)
strings16_2d = reshape(strings16[1:6], (2,3))
empty_string = ""
empty_string_array = ASCIIString[]
empty_array_of_strings = ASCIIString[""]
tf = true
TF = A .> 10
B = [-1.5 sqrt(2) NaN 6;
     0.0  Inf eps() -Inf]
AB = Any[A, B]
t = (3, "cat")
c = float32(3)+float32(7)im
cint = 1+im  # issue 108
C = reinterpret(Complex128, B, (4,))
emptyA = zeros(0,2)
emptyB = zeros(2,0)
try
    global MyStruct
    type MyStruct
        len::Int
        data::Array{Float64}
        MyStruct(len::Int) = new(len)
        MyStruct(len::Int, data::Array{Float64}) = new(len, data)
    end
catch
end
ms = MyStruct(2, [3.2, -1.7])
msempty = MyStruct(5, Float64[])
sym = :TestSymbol
syms = [:a, :b]
d = Dict([(syms[1],"aardvark"), (syms[2], "banana")])
ex = quote
    function incrementby1(x::Int)
        x+1
    end
end
T = Uint8
char = 'x'
unicode_char = '\U10ffff'
α = 22
β = Any[[1, 2], [3, 4]]  # issue #93
vv = Vector{Int}[[1,2,3]]  # issue #123
typevar = Array{Int}[[1]]
typevar_lb = Vector{TypeVar(:U, Integer)}[[1]]
typevar_ub = Vector{TypeVar(:U, Int, Any)}[[1]]
typevar_lb_ub = Vector{TypeVar(:U, Int, Real)}[[1]]
undef = cell(1)
undefs = cell(2, 2)
ms_undef = MyStruct(0)
# Unexported type:
cpus = Base.Sys.cpu_info()
# Immutable type:
rng = 1:5
# Type with a pointer field (#84)
immutable ObjWithPointer
    a::Ptr{Void}
end
objwithpointer = ObjWithPointer(0)
# Custom BitsType (#99)
bitstype 64 MyBT
bt = reinterpret(MyBT, int64(55))
# Symbol arrays (#100)
sa_asc = [:a, :b]
sa_utf8 = [:α, :β]
# SubArray (to test tuple type params)
subarray = sub([1:5], 1:5)
# Array of empty tuples (to test tuple type params)
arr_empty_tuple = ()[]
immutable EmptyImmutable end
emptyimmutable = EmptyImmutable()
type EmptyType end
emptytype = EmptyType()
immutable EmptyII
    x::EmptyImmutable
end
emptyii = EmptyII(EmptyImmutable())
immutable EmptyIT
    x::EmptyType
end
emptyit = EmptyIT(EmptyType())
type EmptyTI
    x::EmptyImmutable
end
emptyti = EmptyTI(EmptyImmutable())
type EmptyTT
    x::EmptyType
end
emptytt = EmptyTT(EmptyType())
immutable EmptyIIOtherField
    x::EmptyImmutable
    y::Float64
end
emptyiiotherfield = EmptyIIOtherField(EmptyImmutable(), 5.0)

# Unicode type field names (#118)
type MyUnicodeStruct☺{τ}
    α::τ
    ∂ₓα::τ
    MyUnicodeStruct☺(α::τ, ∂ₓα::τ) = new(α, ∂ₓα)
end
unicodestruct☺ = MyUnicodeStruct☺{Float64}(1.0, -1.0)
# Arrays of matrices (#131)
array_of_matrices = Matrix{Int}[[1 2; 3 4], [5 6; 7 8]]
# Tuple of arrays and bitstype
tup = (1, 2, [1, 2], [1 2; 3 4], bt)
# Empty tuple
empty_tup = ()
# Non-pointer-free immutable
immutable MyImmutable{T}
    x::Int
    y::Vector{T}
    z::Bool
end
nonpointerfree_immutable_1 = MyImmutable(1, [1., 2., 3.], false)
nonpointerfree_immutable_2 = MyImmutable(2, Any[3., 4., 5.], true)
immutable MyImmutable2
    x::Vector{Int}
    MyImmutable2() = new()
end
nonpointerfree_immutable_3 = MyImmutable2()
# Immutable with a non-concrete datatype (issue #143)
immutable Vague
    name::ByteString
end
vague = Vague("foo")
# Immutable with a union of BitsTypes
immutable BitsUnion
    x::Union(Int64, Float64)
end
bitsunion = BitsUnion(5.0)
# Immutable with a union of Types
immutable TypeUnionField
    x::Union(Type{Int64}, Type{Float64})
end
typeunionfield = TypeUnionField(Int64)
# Generic union type field
immutable GenericUnionField
    x::Union(Vector{Int},Int)
end
genericunionfield = GenericUnionField(1)
# Array references
arr_contained = [1, 2, 3]
arr_ref = typeof(arr_contained)[]
push!(arr_ref, arr_contained, arr_contained)
# Object references
type ObjRefType
    x::ObjRefType
    y::ObjRefType
    ObjRefType() = new()
    ObjRefType(x, y) = new(x, y)
end
ref1 = ObjRefType()
obj_ref = ObjRefType(ObjRefType(ref1, ref1), ObjRefType(ref1, ref1))
# Immutable that requires padding between elements in array
immutable PaddingTest
    x::Int64
    y::Int8
end
padding_test = PaddingTest[PaddingTest(i, i) for i = 1:8]
# Empty arrays of various types and sizes
empty_arr_1 = Int[]
empty_arr_2 = Array(Int, 56, 0)
empty_arr_3 = Any[]
empty_arr_4 = cell(0, 97)
# Moderately big dataset (which will be mmapped)
bigdata = [1:10000]
# BigFloats and BigInts
bigints = big(3).^(1:100)
bigfloats = big(3.2).^(1:100)
# None
none = Union()
nonearr = Array(Union(), 5)

# some data big enough to ensure that compression is used:
Abig = kron(eye(10), rand(20,20))
Sbig = "A test string "^1000

iseq(x,y) = isequal(x,y)
iseq(x::MyStruct, y::MyStruct) = (x.len == y.len && x.data == y.data)
iseq(x::MyImmutable, y::MyImmutable) = (isequal(x.x, y.x) && isequal(x.y, y.y) && isequal(x.z, y.z))
iseq(x::Union(EmptyTI, EmptyTT), y::Union(EmptyTI, EmptyTT)) = isequal(x.x, y.x)
iseq(c1::Array{Base.Sys.CPUinfo}, c2::Array{Base.Sys.CPUinfo}) = length(c1) == length(c2) && all([iseq(c1[i], c2[i]) for i = 1:length(c1)])
function iseq(c1::Base.Sys.CPUinfo, c2::Base.Sys.CPUinfo)
    for n in Base.Sys.CPUinfo.names
        if getfield(c1, n) != getfield(c2, n)
            return false
        end
    end
    true
end
iseq(x::MyUnicodeStruct☺, y::MyUnicodeStruct☺) = (x.α == y.α && x.∂ₓα == y.∂ₓα)
iseq(x::Array{None}, y::Array{None}) = size(x) == size(y)
macro check(fid, sym)
    ex = quote
        let tmp
            try
                tmp = read($fid, $(string(sym)))
            catch e
                warn("Error reading ", $(string(sym)))
                rethrow(e)
            end
            if !iseq(tmp, $sym)
                written = $sym
                error("For ", $(string(sym)), ", read value $tmp does not agree with written value $written")
            end
            written_type = typeof($sym)
            if typeof(tmp) != written_type
                error("For ", $(string(sym)), ", read type $(typeof(tmp)) does not agree with written type $(written_type)")
            end
        end
    end
    esc(ex)
end

# Test for equality of expressions, skipping line numbers
checkexpr(a, b) = @assert a == b
function checkexpr(a::Expr, b::Expr)
    @assert a.head == b.head
    i = 1
    j = 1
    while i <= length(a.args) && j <= length(b.args)
        if isa(a.args[i], Expr) && a.args[i].head == :line
            i += 1
            continue
        end
        if isa(b.args[j], Expr) && b.args[j].head == :line
            j += 1
            continue
        end
        checkexpr(a.args[i], b.args[j])
        i += 1
        j += 1
    end
    @assert i >= length(a.args) && j >= length(b.args)
end

fn = joinpath(tempdir(),"test.jld")

# Issue #106
module Mod106
bitstype 64 Typ{T}
typ{T}(x::Int64, ::Type{T}) = Base.box(Typ{T}, Base.unbox(Int64,x))
abstract UnexportedT
end

for compress in (true,false)
    fnc = compress ? fn*".c" : fn # workaround #176
    fid = jldopen(fnc, "w", compress=compress)
    @write fid x
    @write fid A
    @write fid Aarray
    @write fid str
    @write fid stringsA
    @write fid stringsU
    @write fid strings16
    @write fid strings16_2d
    @write fid empty_string
    @write fid empty_string_array
    @write fid empty_array_of_strings
    @write fid tf
    @write fid TF
    @write fid AB
    @write fid t
    @write fid c
    @write fid cint
    @write fid C
    @write fid emptyA
    @write fid emptyB
    @write fid ms
    @write fid msempty
    @write fid sym
    @write fid syms
    @write fid d
    if compress # work around #176
        @write fid ex
    end
    @write fid T
    @write fid char
    @write fid unicode_char
    @write fid α
    @write fid β
    @write fid vv
    @write fid cpus
    @write fid rng
    @write fid typevar
    @write fid typevar_lb
    @write fid typevar_ub
    @write fid typevar_lb_ub
    @write fid undef
    @write fid undefs
    @write fid ms_undef
    @test_throws JLD.PointerException @write fid objwithpointer
    @write fid bt
    @write fid sa_asc
    @write fid sa_utf8
    @write fid subarray
    @write fid arr_empty_tuple
    @write fid emptyimmutable
    @write fid emptytype
    @write fid emptyii
    @write fid emptyit
    @write fid emptyti
    @write fid emptytt
    @write fid emptyiiotherfield
    @write fid unicodestruct☺
    @write fid array_of_matrices
    @write fid tup
    @write fid empty_tup
    @write fid nonpointerfree_immutable_1
    @write fid nonpointerfree_immutable_2
    @write fid nonpointerfree_immutable_3
    @write fid vague
    @write fid bitsunion
    @write fid typeunionfield
    @write fid genericunionfield
    @write fid arr_ref
    @write fid obj_ref
    @write fid padding_test
    @write fid empty_arr_1
    @write fid empty_arr_2
    @write fid empty_arr_3
    @write fid empty_arr_4
    @write fid bigdata
    @write fid bigfloats
    @write fid bigints
    @write fid none
    @write fid nonearr
    @write fid Abig
    @write fid Sbig
    # Make sure we can create groups (i.e., use HDF5 features)
    g = g_create(fid, "mygroup")
    i = 7
    @write g i
    write(fid, "group1/x", Any[1])
    write(fid, "group2/x", Any[2])
    close(fid)

    # mmapping currently fails on Windows; re-enable once it can work
    for mmap = (@windows ? false : (false, true))
        fidr = jldopen(fnc, "r", mmaparrays=mmap)
        @check fidr x
        @check fidr A
        @check fidr Aarray
        @check fidr str
        @check fidr stringsA
        @check fidr stringsU
        @check fidr strings16
        @check fidr strings16_2d
        @check fidr empty_string
        @check fidr empty_string_array
        @check fidr empty_array_of_strings
        @check fidr tf
        @check fidr TF
        @check fidr AB
        @check fidr t
        @check fidr c
        @check fidr cint
        @check fidr C
        @check fidr emptyA
        @check fidr emptyB
        @check fidr ms
        @check fidr msempty
        @check fidr sym
        @check fidr syms
        @check fidr d
        if compress # work around #176
            exr = read(fidr, "ex")   # line numbers are stripped, don't expect equality
            checkexpr(ex, exr)
        end
        @check fidr T
        @check fidr char
        @check fidr unicode_char
        @check fidr α
        @check fidr β
        @check fidr vv
        @check fidr cpus
        @check fidr rng
        @check fidr typevar
        @check fidr typevar_lb
        @check fidr typevar_ub
        @check fidr typevar_lb_ub

        # Special cases for reading undefs
        undef = read(fidr, "undef")
        if !isa(undef, Array{Any, 1}) || length(undef) != 1 || isdefined(undef, 1)
            error("For undef, read value does not agree with written value")
        end
        undefs = read(fidr, "undefs")
        if !isa(undefs, Array{Any, 2}) || length(undefs) != 4 || any(map(i->isdefined(undefs, i), 1:4))
            error("For undefs, read value does not agree with written value")
        end
        ms_undef = read(fidr, "ms_undef")
        if !isa(ms_undef, MyStruct) || ms_undef.len != 0 || isdefined(ms_undef, :data)
            error("For ms_undef, read value does not agree with written value")
        end

        @check fidr bt
        @check fidr sa_asc
        @check fidr sa_utf8
        @check fidr subarray
        @check fidr arr_empty_tuple
        @check fidr emptyimmutable
        @check fidr emptytype
        @check fidr emptyii
        @check fidr emptyit
        @check fidr emptyti
        @check fidr emptytt
        @check fidr emptyiiotherfield
        @check fidr unicodestruct☺
        @check fidr array_of_matrices
        @check fidr tup
        @check fidr empty_tup
        @check fidr nonpointerfree_immutable_1
        @check fidr nonpointerfree_immutable_2
        @check fidr nonpointerfree_immutable_3
        vaguer = read(fidr, "vague")
        @test typeof(vaguer) == typeof(vague) && vaguer.name == vague.name
        @check fidr bitsunion
        @check fidr typeunionfield
        @check fidr genericunionfield

        arr = read(fidr, "arr_ref")
        @test arr == arr_ref
        @test arr[1] === arr[2]

        obj = read(fidr, "obj_ref")
        @test obj.x.x === obj.x.y == obj.y.x === obj.y.y
        @test obj.x !== obj.y

        @check fidr padding_test
        @check fidr empty_arr_1
        @check fidr empty_arr_2
        @check fidr empty_arr_3
        @check fidr empty_arr_4
        @check fidr bigdata
        @check fidr bigfloats
        @check fidr bigints
        @check fidr none
        @check fidr nonearr
        @check fidr Abig
        @check fidr Sbig

        x1 = read(fidr, "group1/x")
        @assert x1 == Any[1]
        x2 = read(fidr, "group2/x")
        @assert x2 == Any[2]

        close(fidr)
    end
end # compress in (true,false)

for compress in (false,true)
    # object references in a write session
    x = ObjRefType()
    a = [x, x]
    b = [x, x]
    @save fn a b
    jldopen(fn, "r") do fid
        a = read(fid, "a")
        b = read(fid, "b")
        @test a[1] === a[2] === b[2] === a[1]

        # Let gc get rid of a and b
        a = nothing
        b = nothing
        gc()

        a = read(fid, "a")
        b = read(fid, "b")
        @test typeof(a[1]) == ObjRefType
        @test a[1] === a[2] === b[2] === a[1]
    end

    # do syntax
    jldopen(fn, "w"; compress=compress) do fid
        g_create(fid, "mygroup") do g
            write(g, "x", 3.2)
        end
    end
    fid = jldopen(fn, "r")
    @assert names(fid) == ASCIIString["mygroup"]
    g = fid["mygroup"]
    @assert names(g) == ASCIIString["x"]
    @assert read(g, "x") == 3.2
    close(g)
    close(fid)

    # Function load() and save() syntax
    d = Dict([("x",3.2), ("β",β), ("A",A)])
    save(fn, d, compress=compress)
    d2 = load(fn)
    @assert d == d2
    β2 = load(fn, "β")
    @assert β == β2
    β2, A2 = load(fn, "β", "A")
    @assert β == β2
    @assert A == A2

    save(fn, "x", 3.2, "β", β, "A", A, compress=compress)
    d3 = load(fn)
    @assert d == d3

    # #71
    jldopen(fn, "w", compress=compress) do file
        file["a"] = 1
    end
    jldopen(fn, "r") do file
        @assert read(file, "a") == 1
    end

    # Issue #106
    save(fn, "i106", Mod106.typ(int64(1), Mod106.UnexportedT), compress=compress)
    i106 = load(fn, "i106")
    @assert i106 == Mod106.typ(int64(1), Mod106.UnexportedT)

    # bracket syntax for datasets
    jldopen(fn, "w", compress=compress) do file
        file["a"] = [1:100]
        file["b"] = [x*y for x=1:10,y=1:10]
        file["c"] = Any[1, 2, 3]
        file["d"] = [1//2, 1//4, 1//8]
    end
    jldopen(fn, "r+", compress=compress) do file
        @test(file["a"][1:50] == [1:50])
        file["a"][1:50] = 1:2:100
        @test(file["a"][1:50] == [1:2:100])
        @test(file["b"][5,6][1]==5*6)
        @test(file["c"][1:2] == [1, 2])
        file["c"][2:3] = [5, 7]
        @test(read(file, "c") == [1, 5, 7])
        @test(file["d"][2:3] == [1//4, 1//8])
        file["d"][1:1] = [9]
        @test(read(file, "d") == [9, 1//4, 1//8])
    end

    # bracket syntax when created by HDF5
    println("The following unrecognized JLD file warning is a sign of normal operation.")
    if compress
        h5open(fn, "w") do file
            file["a", "blosc",5] = [1:100]
            file["a"][51:100] = [1:50]
            file["b", "blosc",5] = [x*y for x=1:10,y=1:10]
        end
    else
        h5open(fn, "w") do file
            file["a"] = [1:100]
            file["a"][51:100] = [1:50]
            file["b"] = [x*y for x=1:10,y=1:10]
        end
    end
    jldopen(fn, "r") do file
        @assert(file["a"][1:50] == [1:50])
        @assert(file["a"][:] == [[1:50],[1:50]])
        @assert(file["b"][5,6][1]==5*6)
    end

    # delete!
    jldopen(fn, "w", compress=compress) do file
        file["ms"] = ms
        delete!(file, "ms")
        file["ms"] = β
        g = g_create(file,"g")
        file["g/ms"] = ms
        @test_throws ErrorException delete!(file, "_refs/g/ms")
        delete!(file, "g/ms")
        file["g/ms"] = ms
        delete!(file, "/g/ms")
        g["ms"] = ms
        delete!(g,"ms")
        g["ms"] = ms
        delete!(g["ms"])
        g["ms"] = ms
        delete!(g)
        g = g_create(file,"g")
        g["ms"] = ms
        delete!(g)
    end
    jldopen(fn, "r") do file
        @assert(read(file["ms"]) == β)
        @assert(!exists(file, "g/ms"))
        @assert(!exists(file, "g"))
    end

end # compress in (false,true)

# mismatched and missing types
module JLDTemp1
using JLD
import ..fn, Core.Intrinsics.box

type TestType1
    x::Int
end
type TestType2
    x::Int
end
immutable TestType3
    x::TestType2
end

type TestType4
    x::Int
end
type TestType5
    x::TestType4
end
type TestType6 end
bitstype 8 TestType7
immutable TestType8
    a::TestType4
    b::TestType5
    c::TestType6
    d::TestType7
end

jldopen(fn, "w") do file
    truncate_module_path(file, JLDTemp1)
    write(file, "x1", TestType1(1))
    write(file, "x2", TestType3(TestType2(1)))
    write(file, "x3", TestType4(1))
    write(file, "x4", TestType5(TestType4(2)))
    write(file, "x5", [TestType5(TestType4(i)) for i = 1:5])
    write(file, "x6", TestType6())
    write(file, "x7", reinterpret(TestType7, 0x77))
    write(file, "x8", TestType8(TestType4(2), TestType5(TestType4(3)),
                                TestType6(), reinterpret(TestType7, 0x12)))
end
end

type TestType1
    x::Float64
end
type TestType2
    x::Int
end
immutable TestType3
    x::TestType1
end

import Core.Intrinsics.box, Core.Intrinsics.unbox
jldopen(fn, "r") do file
    @test_throws JLD.TypeMismatchException read(file, "x1")
    @test_throws TypeError read(file, "x2")
    println("The following missing type warnings are a sign of normal operation.")
    @test read(file, "x3").x == 1
    @test read(file, "x4").x.x == 2
    x = read(file, "x5")
    for i = 1:5
        @test x[i].x.x == i
    end
    @test typeof(read(file, "x6")).names == ()
    @test reinterpret(Uint8, read(file, "x7")) == 0x77
    x = read(file, "x8")
    @test x.a.x == 2
    @test x.b.x.x == 3
    @test typeof(x.c).names == ()
    @test reinterpret(Uint8, x.d) == 0x12
end


text:
"""
    is_leap_year(year)

Return `true` if `year` is a leap year in the gregorian calendar.

"""
function is_leap_year(year)
           tst1 = year % 4 == 0
           tst2 = year % 100 != 0
           tst3 = year % 400 == 0
           return tst1 && (tst2 || tst3)
end

max_stars_repo_name: stepinski/julia-in

max_issues_repo_name: stepinski/julia-init

max_forks_repo_name: stepinski/julia-init


text:
"""
    reverse_pure_string(text::String)

Reverse only the non-special characters (alphabet: 'a', to 'z' and 'A' to 'Z') and letting
the rest untouched. For this purpose, `reverse_pure_string` is scanning through the string
as an array, and if the current character is not special, it will reverse the elements by
swapping; vice versa, if the current character is a special element, the boundaries will be
modified.

# Arguments
- `text::String`: String with special character like !,<#

# Examples
```julia-repl
julia> import ClassicAlgorithmsCollections
julia> string = "a!!!b.c79.d,e'f,ghi3###""
julia> ClassicAlgorithmsCollectionsreverse_pure_string(string)
"a!!!b.c79.d,e'f,ghi3###""
```
"""
function reverse_pure_string(text::String)
    low = 1
    high = sizeof(text)

    text_array = collect(text)
    while low < high
        # Check if lower char is special character and change the boundary
        if !isletter(text_array[low])
            low += 1
            # Check if higher char is special character and change the boundary
        elseif !isletter(text_array[high])
            high -= 1
        else
            # Reverse the string by swapping
            text_array[low], text_array[high] = text_array[high], text_array[low]
            low += 1
            high -= 1

        end
    end
    return join(text_array)
end

max_stars_repo_name: Anselmoo/ClassicAlgorithmsCollections

max_issues_repo_name: Anselmoo/ClassicAlgorithmsCollections

max_forks_repo_name: Anselmoo/ClassicAlgorithmsCollections


text:
export ProdAnalyticalRuleAvailable, ProdAnalyticalRuleUnknown
export prod_analytical_rule

import Base: show, prod
import Distributions

struct ProdAnalyticalRuleAvailable end
struct ProdAnalyticalRuleUnknown end

"""
    prod_analytical_rule(::Type, ::Type)

Returns either `ProdAnalyticalRuleAvailable` or `ProdAnalyticalRuleUnknown` for two given distribution types.
Returns `ProdAnalyticalRuleUnknown` by default.

See also: [`prod`](@ref), [`ProdAnalytical`](@ref), [`ProdGeneric`](@ref)
"""
prod_analytical_rule(::Type, ::Type) = ProdAnalyticalRuleUnknown()


text:
struct CartesianDDMRAS{T,C,A} <: DDMPreconditioner{T}
    context::C
    parent::A
end

parent(P::CartesianDDMRAS) = getproperty(P, :parent)
getcontext(P::CartesianDDMRAS) = getproperty(P, :context)

function CartesianDDMRAS(A::CartesianDDMMatrix, init=factorize)
    (; context, parent) = A

    out = map(parent) do el
        init(el)
    end

    CartesianDDMRAS{eltype(A),typeof(context),typeof(out)}(context, out)
end

function size(x::CartesianDDMRAS)
    (; context, parent) = x
    (; dims) = context
    ntuple(i -> mapreduce(getdof, *, dims), 2)
end

"""
    ldiv!(y, P, x)

Computes `P \\ x` in-place of `y`.

"""
function ldiv!(y::CartesianDDMVector, P::CartesianDDMRAS, x::CartesianDDMVector)
    map(parent(y), parent(P), parent(x)) do y_, P_, x_
        y_ .= P_ \ x_
    end
    makecoherent!(y)
    y
end

"""
    ldiv!(P, x)

Computes `P \\ x` in-place of `x`.

"""
function ldiv!(P::CartesianDDMRAS, x::CartesianDDMVector)
    map(parent(P), parent(x)) do P_, x_
        x_ .= P_ \ x_
    end
    makecoherent!(x)
    x
end

"""
    \\(P, x)

Computes `P \\ x`.

"""
function (\)(P::CartesianDDMRAS, x::CartesianDDMVector)
    y = similar(x)
    ldiv!(y, P, x)
end


text:
import Base: copy, *, +, -
export ComputState, rawcoeffs, coeff, BellState

type ComputState{A<:AbstractVector,T,N} <: AbstractQuState{T,N}
    coeffs::A
end

# constructors
ComputState{A<:AbstractVector}(state_vec::A) = ComputState{A,eltype(A),length(state_vec)|>log2|>Int}(state_vec)

coefftype{A,T,N}(state::ComputState{A,T,N}) = A
rawcoeffs(state::ComputState) = state.coeffs
coeff(state::ComputState) = state.coeffs

# Base functions overload
copy(state::ComputState) = ComputState(copy(state.coeffs))
similar_type{Q<:ComputState}(::Type{Q}) = ComputState

# Shortcuts
BellState(n::Int) = ComputState(Complex128[1/sqrt(2^n) for i=1:2^n])

# Reload Operators
(*)(a::Number,B::ComputState) = a*coeff(B)
(*)(B::ComputState,a::Number) = a*B
(*)(A::ComputState,B::ComputState) = coeff(A)*coeff(B)
(+)(A::ComputState,B::ComputState) = coeff(A)+coeff(B)
(-)(A::ComputState,B::ComputState) = coeff(A)-coeff(B)

max_stars_repo_name: Roger-luo/openQC

max_issues_repo_name: Roger-luo/openQC

max_forks_repo_name: Roger-luo/openQC


text:
using Base.Test
using StrainRecon

@testset "randomM" begin

# check that the columns are unique
(m,n,N) = (3,2,100)
M = randomM(m,n,N=N)
@test size(M) == (m,n,N)
uniqcols = zeros(Int,N)
for idx = 1:N
	uniqcols[idx] = size(unique(M[:,:,idx],2),2)
end
@test all(uniqcols .== n)

# check that the column sums of each block are at most 1
(m,n,p) = (24,2,3)
M = randomM(m,n,p=p)
@test size(M) == ((p-1)*m,n)
blockcolsums = sum(reshape(M,(p-1,:)),1)
@test maximum(blockcolsums) <= 1

end


text:
module PyInspector
export pyfuncs, ispyfunc, pydoc, PyArgSpec

using PyCall
@pyimport inspect

#const PY_FUNC_TYPE = pyeval("type(lambda:1)")

"True if the object passed is a Python function"
ispyfunc(o) = isa(o, PyObject) && inspect.isroutine(o) # isfunction, ismethod, ismethoddescriptor

"True if the object passed is a Python type"
ispytype(o) = isa(o, PyObject) && inspect.isclass(o)

"Given a module, return the list of its top-level Python functions"
pyfuncs(m::Module) = filter(ispyfunc, [eval(m, n) for n in names(m)])

"Given a module, return the list of top-level Python types"
pytypes(m::Module) = filter(ispytype, [eval(m, n) for n in names(m)])

"Given a module, return the list of __init__ methods of top-level Python types"
pyinits(m::Module) = filter(t->haskey(t, :__init__) && inspect.ismethod(t[:__init__]), pytypes(m))

"Return the docstring for a Python object"
pydoc(o::PyObject) = o[:__doc__]

immutable PyArgSpec
  args::Array{Any}
  varargs::Union{AbstractString, Void}
  varkw::Union{AbstractString, Void}
  defaults::Union{Tuple, Void}
end

PyArgSpec(f::PyObject) = begin
  args, varargs, varkw, defaults = inspect.getargspec(f)
  PyArgSpec(args, varargs, varkw, (defaults == nothing) ? () : defaults)
end

end


text:
using CSV
using Dates
using DataFrames
using DelimitedFiles
using XLSX
using YAML

using SLiDE  # see src/SLiDE.jl

READ_DIR = joinpath("data", "readfiles")
NROW = 4;

# ******************************************************************************************
files_parse = XLSXInput("generate_yaml.xlsx", "parse", "M1:M180", "parse")
files_parse = write_yaml(READ_DIR, files_parse)

y = read_file(files_parse[1])
files = [[y[k] for k in collect(keys(y))[occursin.("Input", keys(y))]]...;]
file = files[1]

# Find which of these edits are represented in the yaml file of defined edits.
EDITS = ["Rename", "Group", "Match", "Melt", "Add", "Map", "Replace", "Drop", "Operate", "Order"]
KEYS = intersect(EDITS, collect(keys(y)))
if ("Drop" in KEYS && all(get_val.(y["Drop"]) .== "all"))
    KEYS = ["Drop"; setdiff(KEYS, ["Drop"])]
end

# Read data file and shorten to appropriate length.
df = read_file(y["PathIn"], file)
df = df[df[:,:State] .== "WISCONSIN", :][1:NROW,:]

#
desc = "nass"
pathout = joinpath(SLIDE_DIR, "data", "output_example")
CSV.write(joinpath(pathout, desc * "_00.csv"), df)

for ii in 1:length(KEYS)-1
    # editname = lowercase(KEYS[ii])
    editfile = desc * "_" * lpad(ii,2,'0') * "_" * lowercase(KEYS[ii]) * ".csv"

    global df = edit_with(df, y[KEYS[ii]])

    # Ugly, but fine for now.
    # if (desc == "nass" && editname == "operate")
    #     cols = [:yr,:r,:n,:value,:units,:value_0,:factor,:units_0]
    #     global df = df[:,cols]
    # end
    # CSV.write(joinpath(pathout, editfile), df)
    # show(first(df,4))
end

x = y["Operate"]
<<<<<<< HEAD
cols = [setdiff(propertynames(df), unique([x.from; x.to; x.input; x.output]))[1:end-2]; x.output; x.from]
=======
cols = [setdiff(names(df), unique([x.from; x.to; x.input; x.output]))[1:end-2]; x.output; x.from]
>>>>>>> dev
[inp in ensurearray(x.output) ? Symbol(inp, :_0) : inp for inp in x.input]

max_stars_repo_name: matbesancon/SLiDE

max_issues_repo_name: matbesancon/SLiDE

max_forks_repo_name: matbesancon/SLiDE


text:
import Random

@testset "SimpleGraphs" begin
    adjmx1 = [0 1 0; 1 0 1; 0 1 0] # graph
    adjmx2 = [0 1 0; 1 0 1; 1 1 0] # digraph
    # specific concrete generators - no need for loop
    @test @inferred(eltype(SimpleGraph())) == Int
    @test @inferred(eltype(SimpleGraph(adjmx1))) == Int
    @test_throws ArgumentError SimpleGraph(adjmx2)

    @test_throws ErrorException badj(DummySimpleGraph())

    @test @inferred(ne(SimpleGraph(path_digraph(5)))) == 4
    @test @inferred(!is_directed(SimpleGraph))
    @test @inferred(!is_directed(SimpleGraph{Int}))

    @test @inferred(eltype(SimpleDiGraph())) == Int
    @test @inferred(eltype(SimpleDiGraph(adjmx2))) == Int
    @test @inferred(ne(SimpleDiGraph(path_graph(5)))) == 8
    @test @inferred(is_directed(SimpleDiGraph))
    @test @inferred(is_directed(SimpleDiGraph{Int}))

    for gbig in [SimpleGraph(0xff), SimpleDiGraph(0xff)]
        @test @inferred(!add_vertex!(gbig))    # overflow
        @test @inferred(add_vertices!(gbig, 10) == 0)
    end

    gdx = path_digraph(4)
    gx = SimpleGraph()
    for g in testgraphs(gx)
        T = eltype(g)
        @test sprint(show, g) == "{0, 0} undirected simple $T graph"
        @test @inferred(add_vertices!(g, 5) == 5)
        @test sprint(show, g) == "{5, 0} undirected simple $T graph"
    end
    gx = SimpleDiGraph()
    for g in testdigraphs(gx)
        T = eltype(g)
        @test sprint(show, g) == "{0, 0} directed simple $T graph"
        @test @inferred(add_vertices!(g, 5) == 5)
        @test sprint(show, g) == "{5, 0} directed simple $T graph"
    end

    gx = path_graph(4)
    for g in testgraphs(gx)
        @test @inferred(vertices(g)) == 1:4
        @test Edge(2, 3) in edges(g)
        @test @inferred(nv(g)) == 4
        @test @inferred(fadj(g)) == badj(g) == adj(g) == g.fadjlist
        @test @inferred(fadj(g, 2)) == badj(g, 2) == adj(g, 2) == g.fadjlist[2]

        @test @inferred(has_edge(g, 2, 3))
        @test @inferred(!has_edge(g, 20, 3))
        @test @inferred(!has_edge(g, 2, 30))
        @test @inferred(has_edge(g, 3, 2))

        gc = copy(g)
        @test @inferred(add_edge!(gc, 4 => 1)) && gc == cycle_graph(4)
        @test @inferred(has_edge(gc, 4 => 1)) && has_edge(gc, 0x04 => 0x01)
        gc = copy(g)
        @test @inferred(add_edge!(gc, (4, 1))) && gc == cycle_graph(4)
        @test @inferred(has_edge(gc, (4, 1))) && has_edge(gc, (0x04, 0x01))
        gc = copy(g)
        @test add_edge!(gc, 4, 1) && gc == cycle_graph(4)

        @test @inferred(inneighbors(g, 2)) == @inferred(outneighbors(g, 2)) == @inferred(neighbors(g, 2)) == [1, 3]
        @test @inferred(add_vertex!(gc))   # out of order, but we want it for issubset
        @test @inferred(g ⊆ gc)
        @test @inferred(has_vertex(gc, 5))

        @test @inferred(ne(g)) == 3

        @test @inferred(rem_edge!(gc, 1, 2)) && @inferred(!has_edge(gc, 1, 2))
        ga = @inferred(copy(g))
        @test @inferred(rem_vertex!(ga, 2)) && ne(ga) == 1
        @test @inferred(!rem_vertex!(ga, 10))

        @test @inferred(zero(g)) == SimpleGraph{eltype(g)}()

        # concrete tests below

        @test @inferred(eltype(g)) == eltype(fadj(g, 1)) == eltype(nv(g))
        T = @inferred(eltype(g))
        @test @inferred(nv(SimpleGraph{T}(6))) == 6

        @test @inferred(eltype(SimpleGraph(T))) == T
        @test @inferred(eltype(SimpleGraph{T}(adjmx1))) == T

        ga = SimpleGraph(10)
        @test @inferred(eltype(SimpleGraph{T}(ga))) == T

        for gd in testdigraphs(gdx)
            U = eltype(gd)
            @test @inferred(eltype(SimpleGraph(gd))) == U
        end

        @test @inferred(edgetype(g)) == SimpleGraphEdge{T}
        @test @inferred(copy(g)) == g
        @test @inferred(!is_directed(g))

        e = first(edges(g))
        @test @inferred(has_edge(g, e))
    end

    gdx = path_digraph(4)
    for g in testdigraphs(gdx)
        @test @inferred(vertices(g)) == 1:4
        @test Edge(2, 3) in edges(g)
        @test !(Edge(3, 2) in edges(g))
        @test @inferred(nv(g)) == 4
        @test @inferred(fadj(g)[2]) == fadj(g, 2) == [3]
        @test @inferred(badj(g)[2]) == badj(g, 2) == [1]
        @test_throws MethodError adj(g)

        @test @inferred(has_edge(g, 2, 3))
        @test @inferred(!has_edge(g, 3, 2))
        @test @inferred(!has_edge(g, 20, 3))
        @test @inferred(!has_edge(g, 2, 30))

        gc = copy(g)
        @test @inferred(add_edge!(gc, 4 => 1)) && gc == cycle_digraph(4)
        @test @inferred(has_edge(gc, 4 => 1)) && has_edge(gc, 0x04 => 0x01)
        gc = copy(g)
        @test @inferred(add_edge!(gc, (4, 1))) && gc == cycle_digraph(4)
        @test @inferred(has_edge(gc, (4, 1))) && has_edge(gc, (0x04, 0x01))
        gc = @inferred(copy(g))
        @test @inferred(add_edge!(gc, 4, 1)) && gc == cycle_digraph(4)

        @test @inferred(inneighbors(g, 2)) == [1]
        @test @inferred(outneighbors(g, 2)) == @inferred(neighbors(g, 2)) == [3]
        @test @inferred(add_vertex!(gc))   # out of order, but we want it for issubset
        @test @inferred(g ⊆ gc)
        @test @inferred(has_vertex(gc, 5))

        @test @inferred(ne(g)) == 3

        @test @inferred(!rem_edge!(gc, 2, 1))
        @test @inferred(rem_edge!(gc, 1, 2)) && @inferred(!has_edge(gc, 1, 2))
        ga = @inferred(copy(g))
        @test @inferred(rem_vertex!(ga, 2)) && ne(ga) == 1
        @test @inferred(!rem_vertex!(ga, 10))

        @test @inferred(zero(g)) == SimpleDiGraph{eltype(g)}()

        # concrete tests below

        @test @inferred(eltype(g)) == eltype(@inferred(fadj(g, 1))) == eltype(nv(g))
        T = @inferred(eltype(g))
        @test @inferred(nv(SimpleDiGraph{T}(6))) == 6

        @test @inferred(eltype(SimpleDiGraph(T))) == T
        @test @inferred(eltype(SimpleDiGraph{T}(adjmx2))) == T

        ga = SimpleDiGraph(10)
        @test @inferred(eltype(SimpleDiGraph{T}(ga))) == T

        for gu in testgraphs(gx)
            U = @inferred(eltype(gu))
            @test @inferred(eltype(SimpleDiGraph(gu))) == U
        end

        @test @inferred(edgetype(g)) == SimpleDiGraphEdge{T}
        @test @inferred(copy(g)) == g
        @test @inferred(is_directed(g))

        e = first(@inferred(edges(g)))
        @test @inferred(has_edge(g, e))
    end

    gx = complete_graph(4)
    for g in testgraphs(gx)
        h = Graph(g)
        @test g == h
        @test rem_vertex!(g, 2)
        @test nv(g) == 3 && ne(g) == 3
        @test g != h
    end

    gdx = complete_digraph(4)
    for g in testdigraphs(gdx)
        h = DiGraph(g)
        @test g == h
        @test rem_vertex!(g, 2)
        @test nv(g) == 3 && ne(g) == 6
        @test g != h
    end
    # tests for #820
    g = complete_graph(3)
    add_edge!(g, 3, 3)
    rem_vertex!(g, 1)
    @test nv(g) == 2 && ne(g) == 2 && has_edge(g, 1, 1)

    g = path_digraph(3)
    add_edge!(g, 3, 3)
    rem_vertex!(g, 1)
    @test nv(g) == 2 && ne(g) == 2 && has_edge(g, 1, 1)

     # Tests for constructors from iterators of edges
    let
        g_undir = erdos_renyi(200, 100; seed=0)
        add_edge!(g_undir, 200, 1) # ensure that the result uses all vertices
        add_edge!(g_undir, 2, 2) # add a self-loop
        for g in testgraphs(g_undir)
            # We create an edge list, shuffle it and reverse half of its edges
            # using this edge list should result in the same graph
            edge_list = [e for e in edges(g)]
            shuffle!(MersenneTwister(0), edge_list)
            for i in rand(MersenneTwister(0), 1:length(edge_list), length(edge_list) ÷ 2)
                e = edge_list[i]
                Te = typeof(e)
                edge_list[i] = Te(dst(e), src(e))
            end

            edge_iter = (e for e in edge_list)
            edge_set = Set(edge_list)
            edge_set_any = Set{Any}(edge_list)

            g1 = @inferred SimpleGraph(edge_list)
            # we can't infer the return type of SimpleGraphFromIterator at the moment
            g2 = SimpleGraphFromIterator(edge_list)
            g3 = SimpleGraphFromIterator(edge_iter)
            g4 = SimpleGraphFromIterator(edge_set)
            g5 = SimpleGraphFromIterator(edge_set_any)

            @test g == g1
            @test g == g2
            @test g == g3
            @test g == g4
            @test g == g5
            @test edgetype(g) == edgetype(g1)
            @test edgetype(g) == edgetype(g2)
            @test edgetype(g) == edgetype(g3)
            @test edgetype(g) == edgetype(g4)
            @test edgetype(g) == edgetype(g5)
        end
        g_dir = erdos_renyi(200, 100; is_directed=true, seed=0)
        add_edge!(g_dir, 200, 1)
        add_edge!(g_dir, 2, 2)
        for g in testdigraphs(g_dir)
            # We create an edge list and shuffle it
            edge_list = [e for e in edges(g)]
            shuffle!(MersenneTwister(0), edge_list)

            edge_iter = (e for e in edge_list)
            edge_set = Set(edge_list)
            edge_set_any = Set{Any}(edge_list)

            g1 = @inferred SimpleDiGraph(edge_list)
            # we can't infer the return type of SimpleDiGraphFromIterator at the moment
            g2 = SimpleDiGraphFromIterator(edge_list)
            g3 = SimpleDiGraphFromIterator(edge_iter)
            g4 = SimpleDiGraphFromIterator(edge_set)
            g5 = SimpleDiGraphFromIterator(edge_set_any)

            @test g == g1
            @test g == g2
            @test g == g3
            @test g == g4
            @test g == g5
            @test edgetype(g) == edgetype(g1)
            @test edgetype(g) == edgetype(g2)
            @test edgetype(g) == edgetype(g3)
            @test edgetype(g) == edgetype(g4)
            @test edgetype(g) == edgetype(g5)
        end

        # SimpleGraphFromIterator of an empty iterator should result
        # in an empty graph of default edgetype
        empty_iter = (x for x in [])
        @test SimpleGraphFromIterator(empty_iter) == SimpleGraph(0)
        @test SimpleDiGraphFromIterator(empty_iter) == SimpleDiGraph(0)
        @test edgetype(SimpleGraphFromIterator(empty_iter)) == edgetype(SimpleGraph(0))
        @test edgetype(SimpleDiGraphFromIterator(empty_iter)) == edgetype(SimpleDiGraph(0))

        # check if multiple edges && multiple self-loops result in the
        # correct number of edges & vertices
        # edges using integers < 1 should be ignored
        g_undir = SimpleGraph(0)
        for g in testgraphs(g_undir)
            T = edgetype(g)
            edge_list = T.([(4, 4),(1, 2),(4, 4),(1, 2),(4, 4),(2, 1),(0, 1),(1, 0),(0, 0)])
            edge_iter = (e for e in edge_list)
            edge_set = Set(edge_list)
            edge_set_any = Set{Any}(edge_list)

            g1 = @inferred SimpleGraph(edge_list)
            g2 = SimpleGraphFromIterator(edge_list)
            g3 = SimpleGraphFromIterator(edge_iter)
            g4 = SimpleGraphFromIterator(edge_set)
            g5 = SimpleGraphFromIterator(edge_set_any)

            @test nv(g1) == 4
            @test nv(g2) == 4
            @test nv(g3) == 4
            @test nv(g4) == 4
            @test nv(g5) == 4

            @test ne(g1) == 2
            @test ne(g2) == 2
            @test ne(g3) == 2
            @test ne(g4) == 2
            @test ne(g5) == 2
        end
        g_dir = SimpleDiGraph(0)
        for g in testdigraphs(g_dir)
            T = edgetype(g)
            edge_list = T.([(4, 4),(1, 2),(4, 4),(1, 2),(4, 4),(2, 1),(0, 1),(1, 0),(0, 0)])
            edge_iter = (e for e in edge_list)
            edge_set = Set(edge_list)
            edge_set_any = Set{Any}(edge_list)

            g1 = @inferred SimpleDiGraph(edge_list)
            g2 = SimpleDiGraphFromIterator(edge_list)
            g3 = SimpleDiGraphFromIterator(edge_iter)
            g4 = SimpleDiGraphFromIterator(edge_set)
            g5 = SimpleDiGraphFromIterator(edge_set_any)

            @test nv(g1) == 4
            @test nv(g2) == 4
            @test nv(g3) == 4
            @test nv(g4) == 4
            @test nv(g5) == 4

            @test ne(g1) == 3
            @test ne(g2) == 3
            @test ne(g3) == 3
            @test ne(g4) == 3
            @test ne(g5) == 3
        end

        # test for iterators where the type of the elements can only be determined at runtime
        g_undir = SimpleGraph(0)
        for g in testgraphs(g_undir)
            T = edgetype(g)
            edge_list_good = Any[ T.(1, 2), T.(3, 4) ]
            edge_list_bad =  Any[ T.(1, 2), Int64(1) ]

            g1 = SimpleGraphFromIterator(edge_list_good)
            @test edgetype(g1) == T
            @test_throws ArgumentError SimpleGraphFromIterator(edge_list_bad)
        end
        g_dir = SimpleDiGraph(0)
        for g in testdigraphs(g_dir)
            T = edgetype(g)
            edge_list_good = Any[ T.(1, 2), T.(3, 4) ]
            edge_list_bad =  Any[ T.(1, 2), Int64(1) ]

            g1 = SimpleDiGraphFromIterator(edge_list_good)
            @test edgetype(g1) == T
            @test_throws ArgumentError SimpleDiGraphFromIterator(edge_list_bad)
        end

        # If there are edges of multiple types, they should be propagated
        # to a common supertype
        edge_list_1 = Any[Edge{Int8}(1, 2), Edge{Int16}(3, 4)]
        edge_list_2 = Any[Edge{Int16}(1, 2), Edge{Int8}(3, 4)]
        g1_undir = SimpleGraphFromIterator(edge_list_1)
        g2_undir = SimpleGraphFromIterator(edge_list_2)
        g1_dir = SimpleGraphFromIterator(edge_list_1)
        g2_dir = SimpleGraphFromIterator(edge_list_2)

        @test Int8 <: eltype(g1_undir)
        @test Int16 <: eltype(g1_undir)
        @test Int8 <: eltype(g2_undir)
        @test Int16 <: eltype(g2_undir)
        @test Int8 <: eltype(g1_dir)
        @test Int16 <: eltype(g1_dir)
        @test Int8 <: eltype(g2_dir)
        @test Int16 <: eltype(g2_dir)
    end

    # test for rem_vertices!
    let
        g_undir = complete_graph(5)
        g_dir = complete_digraph(5)
        for g in (testgraphs(g_undir) ∪ testdigraphs(g_dir))
            T = eltype(g)

            g5 = copy(g)
            vmap = @inferred rem_vertices!(g5, T[], keep_order=true)
            @test g5 == g
            @test vmap == 1:5
            @test isvalid_simplegraph(g5)

            g4 = copy(g)
            vmap = rem_vertices!(g4, T[3], keep_order=true)
            @test g4 == (is_directed(g) ? complete_digraph(T(4)) : complete_graph(T(4)))
            @test vmap == [1, 2, 4, 5]
            @test isvalid_simplegraph(g4)

            g4 = copy(g)
            add_edge!(g4, 1, 1) # some self_loops
            add_edge!(g4, 2, 2)
            vmap = rem_vertices!(g4, T[1, 1, 1], keep_order=false)
            @test ne(g4) == (is_directed(g) ? 13 : 7)
            @test sort(vmap) == [2, 3, 4, 5]
            @test isvalid_simplegraph(g4)

            g2 = copy(g)
            vmap = rem_vertices!(g2, T[2, 1, 4], keep_order=false)
            @test g2 == (is_directed(g) ? complete_digraph(T(2)) : complete_graph(T(2)))
            @test sort(vmap) == [3, 5]
            @test isvalid_simplegraph(g2)

            g0 = copy(g)
            vmap = rem_vertices!(g0, T[1, 3, 2, 3, 5, 4], keep_order=false)
            @test g0 == (is_directed(g) ? SimpleDiGraph(T(0)) : SimpleGraph(T(0)))
            @test isempty(vmap)
            @test isvalid_simplegraph(g0)
            vmap = rem_vertices!(g0, T[], keep_order=false)
            @test g0 == (is_directed(g) ? SimpleDiGraph(T(0)) : SimpleGraph(T(0)))
            @test isempty(vmap)
            @test isvalid_simplegraph(g0)

            g5 = copy(g)
            @test_throws ArgumentError rem_vertices!(g5, T[2, 6], keep_order=true)
            g5 = copy(g)
            @test_throws ArgumentError rem_vertices!(g5, T[3, 0], keep_order=false)
        end

        g_undir = erdos_renyi(10, 0.5)
        g_dir = erdos_renyi(10, 0.5, is_directed=true)
        for u = 1:2:10
            add_edge!(g_undir, u, u)
            add_edge!(g_dir, u, u)
        end
        a = Random.randsubseq(1:10, 0.5)
        for g in (testgraphs(g_undir) ∪ testdigraphs(g_dir))
            T = eltype(g)

            gt = copy(g)
            gf = copy(g)
            a_converted = convert(Vector{T}, a)
            vmap_t = rem_vertices!(gt, a_converted, keep_order=true)
            vmap_f = rem_vertices!(gf, a_converted, keep_order=false)
            @test issorted(vmap_t)
            @test allunique(vmap_t)
            @test allunique(vmap_f)
            @test sort(vmap_f) == vmap_t
            @test length(vmap_f) == nv(g) - length(a_converted)
            @test vmap_t == setdiff(collect(vertices(g)), a_converted)
            @test LightGraphs.Experimental.has_isomorph(gt, gf)
            gi = g[setdiff(collect(vertices(g)), a_converted)]
            @test LightGraphs.Experimental.has_isomorph(gf, gi)
            @test isvalid_simplegraph(gt)
            @test isvalid_simplegraph(gf)
        end
    end
    # codecov for has_edge(::AbstractSimpleGraph, x, y)
    @test @inferred has_edge(DummySimpleGraph(), 1, 2)
end


text:
import Base: getindex, setindex!, lastindex

Base.IndexStyle(::Type{AFArray}) = Base.IndexCartesian()

export allowslow

const _allowslow = Ref(true)
allowslow(::Type{AFArray}, flag = true) = (_allowslow[] = flag)
function allowslow(f, ::Type{AFArray}, flag = true)
    old, _allowslow[] = _allowslow[], flag
    try
        return f()
    finally
        _allowslow = old
    end
end

assertslow(op) = _allowslow[] || error("$op is disabled")

create_seq(r::UnitRange) = af_seq(r.start-1, r.stop-1, 1)
create_seq(r::StepRange) = af_seq(r.start-1, r.stop-1, r.step)
create_seq(i::Int) = af_seq(i-1, i-1, 1)
create_seq(::Colon) = af_seq(1, 1, 0)

set_indexer!(indexers, i, s::Union{AbstractRange,Int,Colon}) = set_seq_indexer(indexers, create_seq(s), i, true)
set_indexer!(indexers, i, s::AFArray{Bool}) = set_array_indexer(indexers, af_where(s), i)
set_indexer!(indexers, i, s::AFArray) = set_array_indexer(indexers, s-UInt32(1), i)
set_indexer!(indexers, i, s::Vector) = set_indexer!(indexers, i, AFArray(s))

function set_seq_indexer(indexer, idx, dim::dim_t, is_batch::Bool)
    _error(ccall((:af_set_seq_indexer,af_lib),af_err,
                 (Ptr{af_index_t},Ptr{af_seq},dim_t,Bool),
                 indexer, RefValue{af_seq}(idx), dim, is_batch), false)
end

function set_array_indexer(indexer, idx, dim::dim_t)
    _error(ccall((:af_set_array_indexer,af_lib),af_err,
                 (Ptr{af_index_t},af_array,dim_t),
                 indexer,idx.arr,dim), false)
end

function release_indexers(indexers)
    _error(ccall((:af_release_indexers,af_lib),af_err,(Ptr{af_index_t},), indexers), false)
end

function assign_gen!(lhs::AFArray,ndims::dim_t,indices,rhs::AFArray)
    out = RefValue{af_array}(lhs.arr)
    _error(ccall((:af_assign_gen,af_lib),af_err,(Ptr{af_array},af_array,dim_t,Ptr{af_index_t},af_array),out,lhs.arr,ndims,indices,rhs.arr), false)
    if lhs.arr != out[]
        release_array(lhs)
        lhs.arr = out[]
    end
end

function create_indexers(idx)
    indexers = create_indexers()
    for (i, thing) in enumerate(idx)
        set_indexer!(indexers, i-UInt32(1), thing)
    end
    indexers
end

function getindex(a::AFArray{T}, idx::Union{AbstractRange,Colon,AFArray,Vector}, idx1::Int...) where {T}
    assertslow("getindex")
    indexers = create_indexers((idx, idx1...))
    out = index_gen_1(a, length(idx1)+1, indexers)
    release_indexers(indexers)
    out
end

function getindex(a::AFArray{T}, idx0::Union{AbstractRange,Colon,AFArray,Vector,Int},
                  idx::Union{AbstractRange,Colon,AFArray,Vector}, idx1::Int...) where {T}
    assertslow("getindex")
    indexers = create_indexers((idx0, idx, idx1...))
    out = index_gen_2(a, length(idx1)+2, indexers)
    release_indexers(indexers)
    out
end

function getindex(a::AFArray{T}, idx0::Union{AbstractRange,Colon,AFArray,Vector,Int},
                  idx::Union{AbstractRange,Colon,AFArray,Vector}, idx1::Union{AbstractRange,Colon,AFArray,Vector}) where {T}
    assertslow("getindex")
    indexers = create_indexers((idx0, idx, idx1))
    out = index_gen_3(a, 3, indexers)
    release_indexers(indexers)
    out
end

function index_gen_1(_in::AFArray{T,N},ndims::dim_t,indices) where {T,N}
    out = RefValue{af_array}(0)
    _error(ccall((:af_index_gen,af_lib),
                 af_err,(Ptr{af_array},af_array,dim_t,Ptr{af_index_t}),
                 out,_in.arr,ndims,indices), false)
    AFArray{T,1}(out[])
end

function index_gen_2(_in::AFArray{T,N},ndims::dim_t,indices) where {T,N}
    out = RefValue{af_array}(0)
    _error(ccall((:af_index_gen,af_lib),
                 af_err,(Ptr{af_array},af_array,dim_t,Ptr{af_index_t}),
                 out,_in.arr,ndims,indices), false)
    AFArray{T,2}(out[])
end

function index_gen_3(_in::AFArray{T,N},ndims::dim_t,indices) where {T,N}
    out = RefValue{af_array}(0)
    _error(ccall((:af_index_gen,af_lib),
                 af_err,(Ptr{af_array},af_array,dim_t,Ptr{af_index_t}),
                 out,_in.arr,ndims,indices), false)
    AFArray{T,3}(out[])
end

function getindex(a::AFArray{T}, idx::Int...) where {T}
    assertslow("getindex")
    @assert length(idx) <= length(size(a))
    indexers = create_indexers(idx)
    out = index_gen(a, length(idx), indexers)
    release_indexers(indexers)
    Array(out)[1]
end

function setindex!(lhs::AFArray{T}, rhs::AFArray{S}, idx::Union{AbstractRange,Int,Colon,AFArray}...) where {T,S}
    assertslow("setindex!")
    @assert length(idx) <= length(size(lhs))
    indexers = create_indexers(idx)
    if T == S
        assign_gen!(lhs, length(idx), indexers, rhs)
    else
        assign_gen!(lhs, length(idx), indexers, convert(AFArray{T}, rhs))
    end
    release_indexers(indexers)
    rhs
end

function setindex!(lhs::AFArray{T}, val::S, idx::Union{AbstractRange,Int,Colon,AFArray}...) where {T,S}
    assertslow("setindex!")
    sz = get_sizes(idx, lhs)
    rhs = constant(T(val), sz)
    setindex!(lhs, rhs, idx...)
    val
end

function get_sizes(idx::Tuple, lhs::AFArray)
    s = Array{Int}(undef, length(idx))
    for i = 1:length(idx)
        if typeof(idx[i]) <: AbstractRange
            s[i] = length(idx[i])
        elseif idx[i] == Colon()
            s[i] = size(lhs,i)
        elseif typeof(idx[i]) <: Integer
            s[i] = 1
        elseif typeof(idx[i]) <: AFArray{Bool}
            s[i] = sum(idx[i])
        elseif typeof(idx[i]) <: AFArray
            s[i] = length(idx[i])
        end
    end
    (s...,)
end

lastindex(a::AFArray) = length(a)
lastindex(a::AFArray, d) = size(a, d)


text:
using Formatting

# port of the SimHashIndex class

_fe = FormatExpr("{1:x},{2}")

type SimhashIndex
    # from py : list of tuples
    objs::Array{Tuple{String,Simhash}}
    f::Int
    k::Int
    count::Int

    # from py : default dict of sets
    bucket::Dict{String, Set}

    # generate the offsets from the f, k
    offsets::Array{Int}

    function SimhashIndex(objs::Array{Tuple{String,Simhash}}, f::Int, k::Int)
        bucket = Dict{String, Set}(Set())
        offsets = _offsets(f, k)
        for obj in objs
            add_to_index(bucket, obj, offsets)
        end
        new(objs, f, k, length(objs), bucket, offsets)
    end
end

function add_to_index(bucket, obj, offsets)
    # obj = Tuple(string, simhash)
    obj_id = obj[1]
    simhash = obj[2]
    for key in _get_keys(simhash, offsets)
        v = format(_fe, simhash.value, obj_id)
        if !haskey(bucket, key)
            bucket[key] = Set()
        end
        push!(bucket[key], v)
    end
end

function _offsets(f, k)
    # per the orig.: You may optimize this method according to
    #     <http://www.wwwconference.org/www2007/papers/paper215.pdf>
    return [fld(f, (k + 1)) * (i - 1) for i in 1:(k + 1)]
end

function _get_keys(simhash, offsets)
    keys = []
    fe = FormatExpr("{1:x}:{2:x}")
    for (i, offset) in enumerate(offsets)
        if i == length(offsets)
            m = 2 ^ (simhash.f - offset) - 1
        else
            m = 2 ^ (offsets[i+1] - offset) - 1
        end
        c = (simhash.value >> offset) & m
        push!(keys, format(fe, c, i-1))
    end
    return keys
end

function get_near_dups(index, simhash)
    # simhash = simhash obj to check against index
    ans = Set()

    for key in _get_keys(simhash, index.offsets)
        if !haskey(index.bucket, key)
            continue
        end
        dups = index.bucket[key]

        if length(dups) > 200
            println("Warning: big bucket found: ", key)
        end

        for dup in dups
            simhash2, obj_id = split(dup, ',', limit=2)
            simhash2 = Simhash(parse(BigInt, simhash2, 16), simhash.f)
            d = distance(simhash, simhash2)
            if d <= index.k
                push!(ans, obj_id)
            end
        end
    end

    return collect(ans)
end

function delete_from_index(index, obj_id, simhash)
    # string, simhash
    for key in _get_keys(simhash, index.offsets)
        v = format(_fe, simhash.value, obj_id)
        delete!(index.bucket[key], v)
    end
end


text:
# Recipe implementation of Simon Byrne editable table
editablefield(f; kwargs...) = editablefield(widget(f)::AbstractWidget; kwargs...)

function editablefield(w::AbstractWidget; editing = false, format = InteractBase.format)
    editing isa AbstractObservable || (editing = Observable(editing))
    data = [:widget => w, :editing => editing]
    wdg = Widget{:editablefield}(data; output = w)
    @layout! wdg Observables.@map &(:editing) ? :widget : map(format, :widget)
end

function editbutton(save = () -> nothing; editing = false)
    editing isa AbstractObservable || (editing = Observable(editing))
    data = [:edit => button("Edit"), :save => button("Save"), :editing => editing]
    changestate() = (editing[] = !editing[])
    wdg = Widget{:editbutton}(data)
    Observables.on(t -> changestate(), wdg[:edit])
    Observables.on(t -> (save(); changestate()), wdg[:save])
    @layout! wdg Observables.@map &(:editing) ? :save : :edit
end


text:
export Sector

#= Sector (for a base shape) =#

mutable struct Sector <: Shape{2,4}  # M = 2K
    c::SVector{2,Float64}  # center of circle
    r::Float64  # radius of circle
    ϕ₀::Float64  # center angle bisecting sector: -π ≤ ϕ₀ < π  (π excluded)
    ∆ϕ2::Float64  # "radius" in angle dimension: 0 ≤ ∆ϕ2 ≤ π (sector spans from ϕ₀ - ∆ϕ2 to ϕ₀ + ∆ϕ2)

    function Sector(c::AbstractVector{<:Real}, r::Real, ϕ::Real, ∆ϕ::Real)
        r≥0 || throw(ArgumentError("r = $r must be nonnegative."))
        -2π≤∆ϕ≤2π  || throw(ArgumentError("∆ϕ = $∆ϕ must be between -2π and 2π, inclusive."))

        ϕ₀ = rem(ϕ + ∆ϕ/2, 2π, RoundNearest)  # put ϕ₀ in [-π, π)
        ∆ϕ2 = abs(∆ϕ/2)

        return new(c, r, ϕ₀, ∆ϕ2)
    end
end

Base.:(==)(s1::Sector, s2::Sector) = s1.c==s2.c && s1.r==s2.r && s1.ϕ₀==s2.ϕ₀ && s1.∆ϕ2==s2.∆ϕ2
Base.isapprox(s1::Sector, s2::Sector) = s1.c≈s2.c && s1.r≈s2.r && s1.ϕ₀≈s2.ϕ₀ && s1.∆ϕ2≈s2.∆ϕ2
Base.hash(s::Sector, h::UInt) = hash(s.c, hash(s.r, hash(s.ϕ₀, hash(s.∆ϕ2, hash(:Sector, h)))))

distangle(ϕ::Real, ϕ₀:: Real) = rem(ϕ-ϕ₀, 2π, RoundNearest)  # ϕ measured from ϕ₀; result within [-π, π)

function level(x::SVector{2,<:Real}, s::Sector)
    d = x - s.c
    ld = norm(d)
    ϕ = ld==0 ? s.ϕ₀ : atan(d[2], d[1])  # angle to x with respect to c

    return max(ld/s.r, abs(distangle(ϕ, s.ϕ₀)) / s.∆ϕ2) - 1.0
end

function surfpt_nearby(x::SVector{2,<:Real}, s::Sector)
    # Basically mimic the same function for Prism, but proceeds in the (ρ,ϕ) domain.
    d = x - s.c
    ld = norm(d)

    # Calculate the closest point in the ρ dimension and outward normal direciton there.
    r2 = s.r / 2
    ρ = ld - r2  # positive if closer to arc; negative if closer to center
    d̂ = ld ≤ Base.rtoldefault(Float64) * r2  ? SVector(cos(s.ϕ₀),sin(s.ϕ₀)) : normalize(d)

    surfρ = ρ<0 ? 0.0 : s.r  # scalar: closest point to x between center and perimeter point
    noutρ = copysign(1.0,ρ) * d̂  # SVector{2}: outward direction normal at surfρ

    absρ = abs(ρ)
    abs∆ρ = abs(r2 - absρ)  # radial distance between x and either center or perimeter, whichever closer to x

    onbndρ = abs∆ρ ≤ Base.rtoldefault(Float64) * r2  # basically r2 ≈ ρ but faster
    isoutρ = (r2 < absρ) || onbndρ

    # Calculate the closest point in the ϕ dimension and outward normal direciton there.
    ϕ = ld==0 ? 0.0 : distangle(atan(d[2], d[1]), s.ϕ₀)  # positive if closer to end side; negative if closer to start side

    ϕtemp = ϕ<0 ? s.ϕ₀ - s.∆ϕ2 : s.ϕ₀ + s.∆ϕ2
    surfϕ = SVector(cos(ϕtemp), sin(ϕtemp))  # SVector{2}: closest point in ϕ dimension
    cosθ = d̂⋅surfϕ  # cosine of angle between d and surfϕ
    ldcosθ = ld*cosθ  # ldconθ .* surfϕ is surface point
    ldsinθ = ld*sqrt(1-cosθ^2)  # always positive

    ϕtemp += copysign(π/2,ϕ)
    noutϕ = SVector(cos(ϕtemp), sin(ϕtemp))  # SVector{2}: outward direction normal at surfϕ

    absϕ = abs(ϕ)
    abs∆ϕ = abs(s.∆ϕ2 - absϕ)  # angular distance between x and closer side of sector

    onbndϕ = abs∆ϕ ≤ Base.rtoldefault(Float64) * s.∆ϕ2  # basically ∆ϕ2 ≈ ϕ but faster
    isoutϕ = (s.∆ϕ2 < absϕ) || onbndϕ

    # Pick the surface point and outward direction normal depending on the location of x.
    if isoutρ && isoutϕ  # x outside in both ρ and ϕ dimensions
        surf = surfρ .* surfϕ  # one of two end points of arc
        nout = (onbndρ && onbndϕ) ? (noutρ + noutϕ) : (d - surf)
        nout = normalize(nout)
    elseif !isoutρ && isoutϕ  # x inside in ρ dimension, but outside in ϕ dimension
        (surf, nout) = (absϕ < s.∆ϕ2+π/2) ? (ldcosθ .* surfϕ, noutϕ) : (SVector(0.0,0.0), d̂)
    elseif isoutρ && !isoutϕ  # x outside in ρ dimension, but inside in ϕ dimension
        (surf, nout) = (surfρ .* d̂, noutρ)
    else  # !isoutρ && !isoutϕ: x strictly inside sector
        (surf, nout) = (abs∆ρ ≤ ldsinθ) ? (surfρ .* d̂, noutρ) : (ldcosθ .* surfϕ, noutϕ)
    end

    return surf+s.c, nout
end

function bounds(s::Sector)
    # Find the minimum and maximum coordinates among the center and two ends of the arc.
    ϕ = SVector(s.ϕ₀ - s.∆ϕ2, s.ϕ₀ + s.∆ϕ2)  # [start angle, end angle]
    v = s.r .* [SVector(0.0,0.0) [cos.(ϕ) sin.(ϕ)]']  # [center, start point of arc, end point of arc]

    # Consider using the code below once https://github.com/JuliaArrays/StaticArrays.jl/issues/498
    # is resolved:
    # l = minimum(v, dims=Val(2))[:,1]
    # u = maximum(v, dims=Val(2))[:,1]

    xs = v[1,:]
    ys = v[2,:]

    xmin = minimum(xs)
    xmax = maximum(xs)

    ymin = minimum(ys)
    ymax = maximum(ys)

    # Consider the extreme points on the arc in the Cartesian directions.
    abs(distangle(0, s.ϕ₀)) ≤ s.∆ϕ2 && (xmax = max(xmax,s.r))  # sector contains +x-direction from center
    abs(distangle(π/2, s.ϕ₀)) ≤ s.∆ϕ2 && (ymax = max(ymax, s.r))  # sector contains +y-direction from center
    abs(distangle(π, s.ϕ₀)) ≤ s.∆ϕ2  && (xmin = min(xmin, -s.r))  # sector contains -x-direction from center
    abs(distangle(3π/2, s.ϕ₀)) ≤ s.∆ϕ2  && (ymin = min(ymin, -s.r))  # sector contains -y-direction from center

    return (SVector(xmin,ymin)+s.c, SVector(xmax,ymax)+s.c)
end


text:
# NOTE: We use a vector of vectors (of string), to represent sequences with
# each nested vector corresponding to a time point. For example, the pattern
# [["a"], ["b", "c"]] has one item in the first time point and two in the
# second. So it corresponds to "a => b,c" in our original implementation.

import Base.isempty

type Sequence
    sid::Int64                          # sequence ID
    eids::Array{Int, 1}                 # event IDs
    items::Array{Array{String, 1}, 1}   # items at each event (e.g., symptoms)
end

type IDList
    patrn::Array{Array{String, 1}, 1}   # vector of vectors with the elements in `pattern`
    sids::Array{Int, 1}
    eids::Array{Int, 1}
    typ::Symbol                         # pattern type is `:initial`, `:sequence` or `:event`
    supp::Float64
    supp_cnt::Int

    function IDList(patrn, sids, eids, typ, num_sequences)
        res = new(patrn, sids, eids, typ, length(unique(sids))/num_sequences, length(unique(sids)))
        return res
    end
end

type PrefixNode
    pattern::Array{Array{String,1},1}
    # parent::PrefixNode
    seq_ext_children::Array{PrefixNode,1}
    item_ext_children::Array{PrefixNode,1}
    support::Int64

    PrefixNode(pattern) = new(pattern)
end

type SeqRule
    prefix::Array{Array{String,1},1}
    postfix::Array{Array{String,1},1}
    conf::Float64
end

isempty(x::IDList) = isempty(x.sids)

function allempty(x::Array{IDList, 1})
    res = true
    for i = 1:length(x)
        if !isempty(x[i])
            res = false
            break
        end
    end
    return res
end

# This function extracts the suffix
# from the ID list's pattern
suffix(idlist::IDList) = idlist.patrn[end]

# This function extracts the prefix
# from the ID list's pattern
prefix(idlist::IDList) = prefix(idlist.patrn)

function prefix(x::Array{Array{String, 1}, 1})
    res = x[1:end-1]
    if length(x[end]) > 1
        res = [res; [x[end][1:end-1]]]
    end
    res                 # returns array of arrays
end

function suffix(x::Array{Array{String, 1}, 1})
     res = x[end][end]
     res                # returns a string
end

# Given an array of `Sequence` objects, this function
# returns the first (k = 1) id-lists.
function first_idlist(seqs::Array{Sequence, 1}, pattern, num_sequences)
    sids = Array{Int, 1}(0)
    eids = Array{Int, 1}(0)

    for s in seqs
        for j = 1:length(s.eids)
            if pattern ∈ s.items[j]
                push!(sids, s.sid)
                push!(eids, s.eids[j])
            end
        end
    end
    return IDList([[pattern]], sids, eids, :initial, num_sequences)
end

# l1 and l2 are IDList objects
function equality_join(l1, l2, num_sequences)
    sids = Array{Int, 1}(0)
    eids = Array{Int, 1}(0)
    n = length(l1.sids)
    m = length(l2.sids)
    for i = 1:n
        for j = 1:m
            if l1.sids[i] == l2.sids[j] && l1.eids[i] == l2.eids[j]
                push!(sids, l1.sids[i])
                push!(eids, l1.eids[i])
            end
        end
    end
    ## type of `pattern` wasn't being infered by compiler
    # pattern::Array{Array{String, 1}, 1} = [l1.patrn[1:end-1]; [sort([l1.patrn[end]; suffix(l2.patrn)])]]

    ## this method works with type inference
    patrn_head = l1.patrn[1:end-1]
    patrn_tail::Array{String,1} = vcat(l1.patrn[end], suffix(l2.patrn))
    sort!(patrn_tail)
    pattern = vcat(patrn_head, [patrn_tail])

    return IDList(pattern, sids, eids, :event, num_sequences)
end

# This is a helper function used in the various
# join functions. It takes a seqeunce ID and event ID,
# as well as vectors of sequence and event IDs we
# have already seen. It returns a boolean.
function already_seen(sid1, eid2, tm_sids, tm_eids)
    res = false
    n = length(tm_sids)
    for i = 1:n
        if tm_sids[i] == sid1 && tm_eids[i] == eid2
            res = true
            break
        end
    end
    return res
end

# @code_warntype already_seen([1, 1, 1, 2], [10, 15, 20, 10], [1, 2, 3], [10, 20, 15])

# This is the guts of the function below, which performs
# a temporal join for those cases in which both id-lists
# are for sequence patterns.
function temporal_join!(l1, l2, sids1, eids1, sids2, eids2, sids3, eids3, n, m)
    for i = 1:n
        for j = 1:m
            if l1.sids[i] == l2.sids[j]
                if l1.eids[i] < l2.eids[j] && !already_seen(l1.sids[i], l2.eids[j], sids1, eids1)
                    push!(sids1, l1.sids[i])
                    push!(eids1, l2.eids[j])
                elseif l2.eids[j] < l1.eids[i] && !already_seen(l2.sids[j], l1.eids[i], sids2, eids2)
                    push!(sids2, l2.sids[j])
                    push!(eids2, l1.eids[i])
                elseif l1.eids[i] == l2.eids[j] && !already_seen(l1.sids[i], l1.eids[i], sids3, eids3) && suffix(l1.patrn) ≠ suffix(l2.patrn)
                    push!(sids3, l1.sids[i])
                    push!(eids3, l1.eids[i])
                end
            end
        end
    end
end

# This function executes a temporal join for those cases
# in which both id-lists are for sequence patterns.
function temporal_join(l1, l2, ::Type{Val{:sequence}}, ::Type{Val{:sequence}}, num_sequences)
    # initialize 3 pairs of empty `sids` and `eids` arrays
    sids1 = Array{Int,1}(0)
    eids1 = Array{Int,1}(0)
    sids2 = Array{Int,1}(0)
    eids2 = Array{Int,1}(0)
    sids3 = Array{Int,1}(0)
    eids3 = Array{Int,1}(0)

    n = length(l1.sids)
    m = length(l2.sids)

    # this is the workhorse of the function
    temporal_join!(l1, l2, sids1, eids1, sids2, eids2, sids3, eids3, n, m)

    seq_patrn1 = [l1.patrn; [[suffix(l2.patrn)]]]
    patrn_head = l1.patrn[1:end-1]
    patrn_tail::Array{String,1} = vcat(l1.patrn[end], suffix(l2.patrn))
    sort!(patrn_tail)
    event_patrn = vcat(patrn_head, [patrn_tail])

    seq_patrn2 = [l2.patrn; [[suffix(l1.patrn)]]]

    idlist_arr = IDList[IDList(seq_patrn1,
                               sids1,
                               eids1,
                               :sequence,
                               num_sequences),
                        IDList(seq_patrn2,
                               sids2,
                               eids2,
                               :sequence,
                               num_sequences),
                        IDList(event_patrn,
                               sids3,
                               eids3,
                               :event,
                               num_sequences)]
    return idlist_arr
end

# This function executes the temporal join for cases in which
# one id-list is for an event pattern and the other is for a
# sequence pattern.
function temporal_join(l1, l2, ::Type{Val{:event}}, ::Type{Val{:sequence}}, num_sequences)
    sids = Array{Int,1}(0)
    eids = Array{Int,1}(0)
    n = length(l1.sids)
    m = length(l2.sids)

    for i = 1:n
        for j = 1:m
            if l1.sids[i] == l2.sids[j] && l1.eids[i] < l2.eids[j] && !already_seen(l1.sids[i], l2.eids[j], sids, eids)
                push!(sids, l1.sids[i])
                push!(eids, l2.eids[j])
            end
        end
    end
    pattern = [l1.patrn; [[suffix(l2.patrn)]]]

    return IDList[IDList(pattern, sids, eids, :sequence, num_sequences)]
end

temporal_join(l1, l2, ::Type{Val{:sequence}}, ::Type{Val{:event}}, num_sequences) = temporal_join(l2, l1, Val{:event}, Val{:sequence}, num_sequences)

# The first merging operation executes both equality and
# temporal joins on id-lists with atoms of length 1.
function first_merge!(l1::IDList, l2::IDList, eq_sids, eq_eids, tm_sids, tm_eids)
    for i = 1:length(l1.sids)
        for j = 1:length(l2.sids)
            if l1.sids[i] == l2.sids[j]
                # equality join
                if l1.eids[i] == l2.eids[j] && suffix(l1) ≠ suffix(l2)
                    push!(eq_sids, l1.sids[i])
                    push!(eq_eids, l1.eids[i])
                # temporal join
                elseif l1.eids[i] < l2.eids[j] && !already_seen(l1.sids[i], l2.eids[j], tm_sids, tm_eids)
                    push!(tm_sids, l1.sids[i])
                    push!(tm_eids, l2.eids[j])
                end
            end
        end
    end
end

# Given two id-lists, this function executes the first
# merge operation. An array of merged id-list (k = 2)
# is returned.
function first_merge(l1::IDList, l2::IDList, num_sequences, minsupp)
    eq_sids = Array{Int, 1}(0)
    tm_sids = Array{Int, 1}(0)

    eq_eids = Array{Int, 1}(0)
    tm_eids = Array{Int, 1}(0)

    first_merge!(l1, l2, eq_sids, eq_eids, tm_sids, tm_eids)
    patrn_tail::Array{String,1} = vcat(l1.patrn[end], suffix(l2.patrn))
    sort!(patrn_tail)
    event_patrn = vcat(l1.patrn[1:end-1], [patrn_tail])
    seq_patrn = [l1.patrn; [[suffix(l2.patrn)]]]

    event_idlist = IDList(event_patrn, eq_sids, eq_eids, :event, num_sequences)
    seq_idlist = IDList(seq_patrn, tm_sids, tm_eids, :sequence, num_sequences)

    merged_idlists = Array{IDList, 1}(0)

    if !isempty(event_idlist) && event_idlist.supp ≥ minsupp
        push!(merged_idlists, event_idlist)
    end
    if !isempty(seq_idlist) && seq_idlist.supp ≥ minsupp
        push!(merged_idlists, seq_idlist)
    end
    return merged_idlists
end

# This function wraps all our join functions.
function merge_idlists(l1, l2, num_sequences)
    if l1.typ == l2.typ == :event                     # both event patterns
        idlist_arr = IDList[equality_join(l1, l2, num_sequences)]
    else
        idlist_arr = temporal_join(l1, l2, Val{l1.typ}, Val{l2.typ}, num_sequences)
    end
    return idlist_arr          # array of merged ID-lists (often of length 1)
end

# s1 = Sequence(
#     1,
#     [1, 2, 3, 4, 5, 6],
#     [["a", "b", "d"], ["a", "e"], ["a", "b", "e"], ["b", "c", "d"], ["b", "c"], ["b", "d"]])
#
# s2 = Sequence(
#     2,
#     [1, 2, 3, 4, 5],
#     [["a", "c", "d"], ["a"], ["a", "b", "d"], ["a", "b"], ["b", "d"]])
#
#
# seq_arr = [s1, s2]
# alist = first_idlist(seq_arr, "a", 2)
# clist = first_idlist(seq_arr, "c", 2)
# dlist = first_idlist(seq_arr, "d", 2)
#
# @code_warntype first_idlist(seq_arr, "d", 2)
# @time first_idlist(seq_arr, "d", 2);
#
#
# cdlist = first_merge(clist, dlist, 2, 0.1)
# @code_warntype first_merge(clist, dlist, 2, 0.1)
# adlist = first_merge(alist, dlist, 2, 0.1)
#
# @code_warntype merge_idlists(adlist[1], cdlist[2], 0.1)
# @code_warntype temporal_join(cdlist[2], adlist[2], Val{:sequence}, Val{:sequence}, 2)
# @code_warntype temporal_join(cdlist[1], adlist[2], Val{:event}, Val{:sequence}, 2)
# @code_warntype temporal_join(cdlist[2], adlist[1], Val{:sequence}, Val{:event}, 2)
# @code_warntype temporal_join(cdlist[1], adlist[1], Val{:event}, Val{:event}, 2)
#
#
# @time temporal_join(cdlist[1], adlist[1], Val{:sequence}, Val{:sequence}, 2)
#
# @code_warntype equality_join(cdlist[1], adlist[1], 2)

# This function is our workhorse for the spade() function
# below. It is only called starting at `F[3]`. Given `f`,
# which is the vector of IDLists `F[k-1]`, this function
# generates `F[k]` by merging the IDList in `f`. It modifies
# `F` in place.
function spade!(f, F, num_sequences, minsupp)
    n = length(f)
    f_tmp = Array{Array{IDList, 1}, 1}(0)

    for i = 1:n
        for j = i:n
            # If both are event patterns, we will only merge
            # id-lists when the suffixes are not identical.
            if f[i].typ == f[j].typ == :event && suffix(f[i]) == suffix(f[j])
                continue
            elseif prefix(f[i]) == prefix(f[j])
                idlist_arr = merge_idlists(f[i], f[j], num_sequences)
                filter!(x -> x.supp ≥ minsupp, idlist_arr)

                if !allempty(idlist_arr)
                    push!(f_tmp, idlist_arr)
                end
            end
        end
    end
    if !isempty(f_tmp)
        fk = reduce(vcat, f_tmp)
        push!(F, unique(fk))
    end

    return nothing  # fixes flag from @code_warntype
end

"""
    spade(seqs, minsupp, maxlength)

Given a vector of `Sequence` objects, this function executes the SPADE
algorithm (Zaki, 2001).

### Arguments
* `seqs`: a vector of `Sequence` objects
* `minsupp`: minimum level of support for a sequential pattern.
* `maxlength`: maximum number of items in a given pattern

The return value is an array of arrays, `F`, where each element of `F` is
an array of `IDList` objects of length `k`. For example, `F[2]` has an array
of all `IDLists` with patterns of length 2 (e.g., {A,B} or {C},{D})
"""
function spade(seqs::Array{Sequence, 1}, minsupp = 0.1, maxlength = 5)
    F = Array{Array{IDList, 1}, 1}(0)
    f1 = Array{IDList, 1}(0)
    items = Array{String, 1}(0)
    n_seq = length(seqs)

    for i = 1:n_seq
        append!(items, unique_items(seqs[i]))
    end
    uniq_items = unique(items)

    for itm in uniq_items
        push!(f1, first_idlist(seqs, itm, n_seq))
    end
    push!(F, f1)
    push!(F, Array{IDList,1}(0))

    n = length(F[1])

    # first merge is handled differently
    for j = 1:n
        for k = 1:n
            idlist_arr = first_merge(F[1][j], F[1][k], n_seq, minsupp)
            if idlist_arr ≠ nothing
                append!(F[2], idlist_arr)
            end
        end
    end
    F[2] = unique(F[2])
    i = 3

    # We persist until arriving at maxlength or
    # until F[k] was empty on most recent iteration
    while i ≤ maxlength && length(F) == i - 1
        spade!(F[i-1], F, n_seq, minsupp)
        i += 1
    end

    return F
end


text:
##
using Random, MAT, Statistics, Printf, LinearAlgebra
using Plots
doPlot=0

try
  global N_threads = ccall((:openblas_get_num_threads64_, Base.libblas_name), Cint, ())
catch
  global N_threads = 0
end

## Functions
function forward(G,m)
    d= G*(m[:].^(-1))
    return d
end

function log_likelihood(d,d_obs,iCd)
    logL = -0.5*(d-d_obs)'*iCd*(d-d_obs)
    return logL[1]
end

function sample_prior(m0,L,nx,ny)
    r=randn(size(L,1),1);
    m_est = m0[:] + L*r;
    m_est=reshape(m_est,ny,nx);
    return m_est
end

## LOAD SOME DATA
vars = matread("ArrenaesGprTomo_30_60.mat");
#vars = matread("ArrenaesGprTomo_55_115.mat");

G = vars["G"]
d0 = vars["d0"]
Cd = vars["Cd"]
Cd = Cd + I(size(Cd,1))*0.4.^2
m0 = vars["m0"]
Cm = vars["Cm"]
d_obs = vars["d_obs"]
d_std = vars["d_std"]
x = vars["x"]
y = vars["y"]
nx=length(x)
ny=length(y)

# next line for faster log-likelihood
iCd = inv(Cd);
# next line for faster prior sampler
C = cholesky(Cm)
L = C.L

## The extended Metropolis algorithm
N_ite=100000

# store every 'i_save' model
i_save = 10;
n_save = 0;
m_post = zeros(ny,nx,Int(floor(N_ite/10)))
size(m_post)
logL_post=zeros(1,N_ite)

# get an initial current model, and data, and log-likelihood
m_cur=sample_prior(m0,L,nx,ny);
d_cur=forward(G,m_cur);
logL_cur=log_likelihood(d_cur,d_obs,iCd);

t_start=time()
step_final=0.02
i_step_final=10000
global N_acc=0
for i=1:N_ite

    ## adjust step length in the first i_step_final iterations
    if (i<i_step_final)
        stepsize=(1-step_final).*(1-(i/i_step_final))+step_final;
    else
        stepsize=step_final
    end

    ## sample prior
    m_pro=sample_prior(m0,L,nx,ny);
    m_pro = m0 + cos(stepsize*pi/2).*(m_cur-m0) + sin(stepsize*pi/2).*(m_pro-m0);

    ## forward response
    d_pro=forward(G,m_pro);
    #
    ## likelihood
    logL_pro=log_likelihood(d_pro,d_obs,iCd);

    ## Accept move in Markov chain
    P_acc = exp(logL_pro-logL_cur);
    #P_acc = 1.0

    if rand(1)[1] <P_acc[1]
        global d_cur = d_pro;
        global m_cur = m_pro;
        global logL_cur = logL_pro;
        global N_acc=N_acc+1;
    end

    logL_post[i]=logL_cur;

    if (mod(i,20000)==0)
        #@printf("i=%5d/%5d, logL_cur=%5.3f, Nacc=%4d, P_acc=%3.2f\n",i,N_ite, logL_cur, N_acc, P_acc);
        @printf("i=%5d/%5d, logL=%5.3f,%5.3f, Nacc=%4d, P_acc=%3.2f\n",i,N_ite, logL_cur,logL_pro, N_acc, P_acc);
    end

    if (mod(i,10000)==0)&(doPlot>1)
        i1=Int(ceil(i/5))
        p1=heatmap(m_cur)
        p2=plot(i1:i,logL_post[i1:i])
        p3=plot([d_cur[:],d_obs[:]])
        l = @layout [a ; b c]
        plot(p1, p2, p3, layout = l)
    end

    # save model?
    if mod(i,i_save)==0
        global n_save = n_save+1;
        m_post[:,:,n_save]=m_cur;
    end

end
t_stop=time()
t_elapsed  = (t_stop-t_start)

@printf("%24s%6s (%s,%dt): t=%6.2fs, N_ite=%8d, %8d iterations/s\n", "JULIA", VERSION,BLAS.vendor(),N_threads,t_elapsed, N_ite, ceil(N_ite/t_elapsed))

## Burnin
burnin=Int(ceil(10000/i_save))
m_post=m_post[:,:,burnin:end];

m_mean = mean(m_post, dims=3)
m_std = std(m_post, dims=3)

## Plot posterior stats
if (doPlot>0)
    heatmap(m_mean[:,:,1])
    savefig("julia_mean.png")
    heatmap(m_std[:,:,1])
    savefig("julia_std.png")
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/ModelReduction.jl/blob/master/LICENSE

using ModelReduction: block_multiply
using MatrixChainMultiply

function perf_test(ndofs, nvals, nkeep, nblocks)
  X = rand(ndofs, nvals)
  R = rand(nvals, nvals)
  A = rand(ndofs, nkeep)
  @time B = block_multiply(X, R, A; n=nblocks)
  gc()
  @time B = X * R * X' * A
  gc()
  @time B = X * R * (X' * A)
  gc()
  @time B = matrixchainmultiply(X, R, X', A)
end

perf_test(30000, 1, 500, 10)


text:
using Documenter, Metaheuristics
using DocumenterCitations

bib = CitationBibliography(joinpath(@__DIR__, "references.bib"))

makedocs(
         bib,
         format = Documenter.HTML(
                                  prettyurls = get(ENV, "CI", nothing) == "true",
                                  assets = ["assets/favicon.ico", "assets/extra_styles.css"],
                                  analytics = "UA-184071594-1",
                                  collapselevel = 1,
                                 ),
         sitename="Metaheuristics.jl",
         authors = "Jesús Mejía",
         pages = [
                  "Index" => "index.md",
                  "Tutorials" => ["tutorials/simple-tutorial.md", "tutorials/create-metaheuristic.md"],
                  "Examples" =>  "examples.md",
                  "Algorithms" => "algorithms.md",
                  "Problems" => "problems.md",
                  "Performance Indicators" => "indicators.md",
                  "Visualization" => "visualization.md",
                  "API References" => "api.md",
                  "Contributing" => "contributing.md",
                  "References" => "references.md",
                 ]
        )

deploydocs(
           repo = "github.com/jmejia8/Metaheuristics.jl.git",
          )


text:
for Typ in (:Ldiv, :Rdiv)
    @eval begin
        struct $Typ{StyleA, StyleB, AType, BType}
            A::AType
            B::BType
        end

        @inline $Typ{StyleA, StyleB}(A::AType, B::BType) where {StyleA,StyleB,AType,BType} =
            $Typ{StyleA,StyleB,AType,BType}(A,B)

        @inline $Typ(A::AType, B::BType) where {AType,BType} =
            $Typ{typeof(MemoryLayout(AType)),typeof(MemoryLayout(BType)),AType,BType}(A, B)

        @inline BroadcastStyle(::Type{<:$Typ}) = ApplyBroadcastStyle()
        @inline broadcastable(M::$Typ) = M

        similar(A::$Typ, ::Type{T}, axes) where T = similar(Array{T}, axes)
        similar(A::$Typ, ::Type{T}) where T = similar(A, T, axes(A))
        similar(A::$Typ) = similar(A, eltype(A))

        @inline copy(M::$Typ) = copyto!(similar(M), M)
        @inline materialize(M::$Typ) = copy(instantiate(M))
    end
end

@inline _ldivaxes(::Tuple{}, ::Tuple{}) = ()
@inline _ldivaxes(::Tuple{}, Bax::Tuple) = Bax
@inline _ldivaxes(::Tuple{<:Any}, ::Tuple{<:Any}) = ()
@inline _ldivaxes(::Tuple{<:Any}, Bax::Tuple{<:Any,<:Any}) = (OneTo(1),last(Bax))
@inline _ldivaxes(Aax::Tuple{<:Any,<:Any}, ::Tuple{<:Any}) = (last(Aax),)
@inline _ldivaxes(Aax::Tuple{<:Any,<:Any}, Bax::Tuple{<:Any,<:Any}) = (last(Aax),last(Bax))

@inline ldivaxes(A, B) = _ldivaxes(axes(A), axes(B))

@inline axes(L::Ldiv) = ldivaxes(L.A, L.B)
@inline size(L::Ldiv) = map(length, axes(L))
@inline length(L::Ldiv{<:Any,<:Any,<:Any,<:AbstractVector}) =size(L.A, 2)

@inline size(L::Rdiv) = (size(L.A, 1),size(L.B,1))
@inline axes(L::Rdiv) = (axes(L.A, 1),axes(L.B,1))

@inline ndims(L::Ldiv) = ndims(L.B)
@inline ndims(L::Rdiv) = 2
@inline eltype(M::Ldiv) = promote_type(Base.promote_op(inv, eltype(M.A)), eltype(M.B))
@inline eltype(M::Rdiv) = promote_type(eltype(M.A), Base.promote_op(inv, eltype(M.B)))

# Lazy getindex

getindex(L::Ldiv, k...) = _getindex(indextype(L), L, k)
concretize(L::AbstractArray) = convert(Array,L)
concretize(L::Ldiv) = ldiv(concretize(L.A), concretize(L.B))
_getindex(::Type{Tuple{I}}, L::Ldiv, (k,)::Tuple{I}) where I = concretize(L)[k]
_getindex(::Type{Tuple{I,J}}, L::Ldiv, (k,j)::Tuple{Colon,J}) where {I,J} = Ldiv(L.A, L.B[:,j])
_getindex(::Type{Tuple{I,J}}, L::Ldiv, (k,j)::Tuple{I,J}) where {I,J} = L[:,j][k]

check_ldiv_axes(A, B) =
    axes(A,1) == axes(B,1) || throw(DimensionMismatch("First axis of A, $(axes(A,1)), and first axis of B, $(axes(B,1)) must match"))

check_rdiv_axes(A, B) =
    axes(A,2) == axes(B,2) || throw(DimensionMismatch("Second axis of A, $(axes(A,2)), and second axis of B, $(axes(B,2)) must match"))

@inline function instantiate(L::Ldiv)
    check_ldiv_axes(L.A, L.B)
    Ldiv(instantiate(L.A), instantiate(L.B))
end

@inline function instantiate(L::Rdiv)
    check_rdiv_axes(L.A, L.B)
    Rdiv(instantiate(L.A), instantiate(L.B))
end

__ldiv!(::Mat, ::Mat, B) where Mat = error("Overload materialize!(::Ldiv{$(typeof(MemoryLayout(Mat))),$(typeof(MemoryLayout(B)))})")
__ldiv!(::Mat, ::Mat, B::LayoutArray) where Mat = error("Overload materialize!(::Ldiv{$(typeof(MemoryLayout(Mat))),$(typeof(MemoryLayout(B)))})")
__ldiv!(_, F, B) = LinearAlgebra.ldiv!(F, B)
@inline _ldiv!(A, B) = __ldiv!(A, factorize(A), B)
@inline _ldiv!(A::Factorization, B) = LinearAlgebra.ldiv!(A, B)
@inline _ldiv!(A::Factorization, B::LayoutArray) = error("Overload materialize!(::Ldiv{$(typeof(MemoryLayout(A))),$(typeof(MemoryLayout(B)))})")

@inline _ldiv!(dest, A, B) = ldiv!(dest, factorize(A), B)
@inline _ldiv!(dest, A::Factorization, B) = LinearAlgebra.ldiv!(dest, A, B)
@inline _ldiv!(dest, A::Transpose{<:Any,<:Factorization}, B) = LinearAlgebra.ldiv!(dest, A, B)
@inline _ldiv!(dest, A::Adjoint{<:Any,<:Factorization}, B) = LinearAlgebra.ldiv!(dest, A, B)

@inline ldiv(A, B) = materialize(Ldiv(A,B))
@inline rdiv(A, B) = materialize(Rdiv(A,B))

@inline ldiv!(A, B) = materialize!(Ldiv(A,B))
@inline rdiv!(A, B) = materialize!(Rdiv(A,B))

@inline ldiv!(C, A, B) = copyto!(C, Ldiv(A,B))
@inline rdiv!(C, A, B) = copyto!(C, Rdiv(A,B))

@inline materialize!(M::Ldiv) = _ldiv!(M.A, M.B)
@inline materialize!(M::Rdiv) = ldiv!(M.B', M.A')'
@inline copyto!(dest::AbstractArray, M::Rdiv) = copyto!(dest', Ldiv(M.B', M.A'))'
@inline copyto!(dest::AbstractArray, M::Ldiv) = _ldiv!(dest, M.A, copy(M.B))

const MatLdivVec{styleA, styleB, T, V} = Ldiv{styleA, styleB, <:AbstractMatrix{T}, <:AbstractVector{V}}
const MatLdivMat{styleA, styleB, T, V} = Ldiv{styleA, styleB, <:AbstractMatrix{T}, <:AbstractMatrix{V}}
const BlasMatLdivVec{styleA, styleB, T<:BlasFloat} = MatLdivVec{styleA, styleB, T, T}
const BlasMatLdivMat{styleA, styleB, T<:BlasFloat} = MatLdivMat{styleA, styleB, T, T}

const MatRdivMat{styleA, styleB, T, V} = Rdiv{styleA, styleB, <:AbstractMatrix{T}, <:AbstractMatrix{V}}
const BlasMatRdivMat{styleA, styleB, T<:BlasFloat} = MatRdivMat{styleA, styleB, T, T}

materialize!(M::Ldiv{ScalarLayout}) = Base.invoke(LinearAlgebra.ldiv!, Tuple{Number,AbstractArray}, M.A, M.B)
materialize!(M::Rdiv{<:Any,ScalarLayout}) = Base.invoke(LinearAlgebra.rdiv!, Tuple{AbstractArray,Number}, M.A, M.B)

function materialize!(M::Ldiv{ScalarLayout,<:SymmetricLayout})
    ldiv!(M.A, symmetricdata(M.B))
    M.B
end
function materialize!(M::Ldiv{ScalarLayout,<:HermitianLayout})
    ldiv!(M.A, hermitiandata(M.B))
    M.B
end
function materialize!(M::Rdiv{<:SymmetricLayout,ScalarLayout})
    rdiv!(symmetricdata(M.A), M.B)
    M.A
end
function materialize!(M::Rdiv{<:HermitianLayout,ScalarLayout})
    rdiv!(hermitiandata(M.A), M.B)
    M.A
end

macro _layoutldiv(Typ)
    ret = quote
        LinearAlgebra.ldiv!(A::$Typ, x::AbstractVector) = ArrayLayouts.ldiv!(A,x)
        LinearAlgebra.ldiv!(A::$Typ, x::AbstractMatrix) = ArrayLayouts.ldiv!(A,x)
        LinearAlgebra.ldiv!(A::$Typ, x::StridedVector) = ArrayLayouts.ldiv!(A,x)
        LinearAlgebra.ldiv!(A::$Typ, x::StridedMatrix) = ArrayLayouts.ldiv!(A,x)

        LinearAlgebra.ldiv!(A::Factorization, x::$Typ) = ArrayLayouts.ldiv!(A,x)

        LinearAlgebra.ldiv!(A::Bidiagonal, B::$Typ) = ArrayLayouts.ldiv!(A,B)

        Base.:\(A::$Typ, x::AbstractVector) = ArrayLayouts.ldiv(A,x)
        Base.:\(A::$Typ, x::AbstractMatrix) = ArrayLayouts.ldiv(A,x)

        Base.:\(x::AbstractMatrix, A::$Typ) = ArrayLayouts.ldiv(x,A)
        Base.:\(x::UpperTriangular, A::$Typ) = ArrayLayouts.ldiv(x,A)
        Base.:\(x::LowerTriangular, A::$Typ) = ArrayLayouts.ldiv(x,A)
        Base.:\(x::Diagonal, A::$Typ) = ArrayLayouts.ldiv(x,A)

        Base.:\(A::Bidiagonal{<:Number}, B::$Typ{<:Number}) = ArrayLayouts.ldiv(A,B)
        Base.:\(A::Bidiagonal, B::$Typ) = ArrayLayouts.ldiv(A,B)
        Base.:\(transA::Transpose{<:Number,<:Bidiagonal{<:Number}}, B::$Typ{<:Number}) = ArrayLayouts.ldiv(transA,B)
        Base.:\(transA::Transpose{<:Any,<:Bidiagonal}, B::$Typ) = ArrayLayouts.ldiv(transA,B)
        Base.:\(adjA::Adjoint{<:Number,<:Bidiagonal{<:Number}}, B::$Typ{<:Number}) = ArrayLayouts.ldiv(adjA,B)
        Base.:\(adjA::Adjoint{<:Any,<:Bidiagonal}, B::$Typ) = ArrayLayouts.ldiv(adjA,B)

        Base.:\(x::$Typ, A::$Typ) = ArrayLayouts.ldiv(x,A)

        Base.:/(A::$Typ, x::AbstractVector) = ArrayLayouts.rdiv(A,x)
        Base.:/(A::$Typ, x::AbstractMatrix) = ArrayLayouts.rdiv(A,x)

        Base.:/(x::AbstractMatrix, A::$Typ) = ArrayLayouts.rdiv(x,A)
        Base.:/(x::Diagonal, A::$Typ) = ArrayLayouts.rdiv(x,A)

        Base.:/(x::$Typ, A::$Typ) = ArrayLayouts.rdiv(x,A)
    end
    if Typ ≠ :LayoutVector
        ret = quote
            $ret
            Base.:\(A::$Typ, x::LayoutVector) = ArrayLayouts.ldiv(A,x)
        end
    end
    if Typ ≠ :LayoutMatrix
        ret = quote
            $ret
            Base.:\(A::$Typ, x::LayoutMatrix) = ArrayLayouts.ldiv(A,x)
        end
    end
    esc(ret)
end

macro layoutldiv(Typ)
    esc(quote
        ArrayLayouts.@_layoutldiv $Typ
        ArrayLayouts.@_layoutldiv UpperTriangular{T, <:$Typ{T}} where T
        ArrayLayouts.@_layoutldiv UnitUpperTriangular{T, <:$Typ{T}} where T
        ArrayLayouts.@_layoutldiv LowerTriangular{T, <:$Typ{T}} where T
        ArrayLayouts.@_layoutldiv UnitLowerTriangular{T, <:$Typ{T}} where T

        ArrayLayouts.@_layoutldiv UpperTriangular{T, <:SubArray{T,2,<:$Typ{T}}} where T
        ArrayLayouts.@_layoutldiv UnitUpperTriangular{T, <:SubArray{T,2,<:$Typ{T}}} where T
        ArrayLayouts.@_layoutldiv LowerTriangular{T, <:SubArray{T,2,<:$Typ{T}}} where T
        ArrayLayouts.@_layoutldiv UnitLowerTriangular{T, <:SubArray{T,2,<:$Typ{T}}} where T
    end)
end

@_layoutldiv LayoutVector


text:
@inline function throw_if_error(error_code::Cuint, str=nothing)
    if error_code == 0
        return
    else
        throw(get_error(error_code, str))
    end
end

function get_error(error_code::Cuint, str=nothing)
    if error_code == 1
        return ErrorException("Unknown SymEngine Exception")
    elseif error_code == 2
        return DivideError()
    elseif error_code == 3
        return ErrorException("Not implemented SymEngine feature")
    elseif error_code == 4
        return DomainError(str)
    elseif error_code == 5
        return Meta.ParseError(str)
    else
        return ErrorException("Unexpected SymEngine error code")
    end
end