text:
using ArgParse
using BioSequences
using CodecZlib
using FASTX
## Command Line Interface ##
cli = ArgParseSettings()
@add_arg_table cli begin
	"--window"
		help = "window size"
		arg_type = Int
		default = 100
	"fasta"
		help = "fasta file, required"
		required = true
end
arg = parse_args(ARGS, cli)
## FASTA reader setup ##

reader = nothing
if splitext(arg["fasta"])[2] == ".gz"
	reader = FASTA.Reader(GzipDecompressorStream(open(arg["fasta"])))
else
	reader = FASTA.Reader(open(arg["fasta"]))
end

## Main loop ##
w = arg["window"]
for record in reader
	#println(">", FASTA.identifier(record))
	chrom = FASTA.sequence(record)
	println(length(dna_seq)) #tell length of every one of the sequenes: run skew.jl test.fa
	for i = 1:length(dna_seq)-w
		g = 0
		c = 0
		for j = i:i+w
			if     dna_seq[j] == DNA_C c+=1
			elseif dna_seq[j] == DNA_G g+=1
			end
		end
		println("$chrome\t$i\t",(g-c)/(g+c))
	end

end

"""
chromsome, coordinate, skew file
what does /dev/null do? get rid of output
produces skew coordinates over whole thing
22 seconds with 100 size window
NOW
GRAB THIS SHIT FROM MEGHANA'S Code
lets try to improve it!!!!! and compare the time (see screengrab)

"""

max_stars_repo_name: KorfLab/julie

max_issues_repo_name: KorfLab/julie

max_forks_repo_name: KorfLab/julie


text:
# Test of SpaceAnalysis

using SpaceAnalysis, Test, Random, DelimitedFiles

"Calculates the magnetic field from a dipole with magnetic moment `M` at `r`."
function dipole(rIn, M)
   x, y, z = rIn
   r = sqrt(x^2 + y^2 + z^2)
   Coef = 1e-7 / r^5 # μ₀/(4*π*r^5)

   B = [3*x^2-r^2 3*x*y     3*x*z;
        3*y*x     3*y^2-r^2 3*y*z;
        3*z*x     3*z*y     3*z^2-r^2
       ] * M * Coef
end

@testset "SpaceAnalysis.jl" begin
   @testset "Spectrum" begin
      Fs = 1000.     # Sampling frequency
      T = 1/Fs       # Sampling period
      L = 1500       # Length of signal
      t = (0:L-1)*T  # Time vector
      # Form a signal containing a 50 Hz sinusoid of amplitude 0.7 and
      # a 120 Hz sinusoid of amplitude 1.
      S = @. 0.7*sinpi(2*50*t) + sinpi(2*120*t)
      # Corrupt the signal with zero-mean white noise with a variance of 4.
      X = S + 2*randn(MersenneTwister(1234), length(t))

      f, P1 = spectrum(X, Fs)
      @test f[end] == 500. && maximum(P1) == 0.9513531477953538

      @test mag2db(f)[end] == 26.989700043360187
   end

   @testset "MVA" begin # minimum variance analysis
      include("MVA_galileo_G8.jl")
      F = mva_analysis("data/Galileo_G8_flyby_MAG.dat.gz", 2656:2775)
      eigenRef = [2079.360,  78.142,   34.309]
      @test F.values ≈ eigenRef atol=1e-3
   end

   @testset "MFA" begin # mean field-aligned coordinate
      # Generate artificial signal
      f = 7e-3 # frequency, [Hz]
      h = 1    # sample period, [s]
      t = range(0.0, 286.0, step=h)
      y = @. sinpi(2 * f * t)
      ynoise = y + 0.1*randn(MersenneTwister(1234), length(t))

      fmax = 1e-2
      fsample = 1 / h
      ybk = emd_average(ynoise, fmax, fsample)

      @test ybk[100] ≈ -0.9506051118519016

      # Transfroming an ideal Earth dipole field along a straight line in x
      Rₑ = 6.378e6
      x = range(3, 20, step=0.1) * Rₑ
      y = ones(size(x))
      z = 2 .* y
      Bx, By, Bz = similar(z), similar(z), similar(z)
      BMoment_Earth = [0.0, 0.0, 7.94e22] # [V*s/(A*m)]
      for i in axes(x, 1)
         Bx[i], By[i], Bz[i] = dipole([x[i], y[i], z[i]], BMoment_Earth)
      end
      # Bμ ∥ B should be dominant.
      Bμ, Bϕ, Bν = mfa(x, y, z, Bx, By, Bz; method="MAVG", nw=15)

      @test Bμ[end] == 3.8254010321827695e-9

      # Introduce 2 oscillations with different freqs s.t. there will be multiple IMFs.
      Bx += 1e-9*sinpi.(2 * f * t[1:length(x)])
      Bx += 1e-10*cospi.(2 * 10f * t[1:length(x)])

      # Bμ ∥ B should be dominant.
      Bμ, Bϕ, Bν = mfa(x, y, z, Bx, By, Bz; method="EMD", fsample, fmax, verbose=true)

      @test Bμ[end] == 3.956605392869394e-9
   end

   @testset "signal" begin # signal generation
      # Alfvénic perturbation
      n0 = 1e6 # density, [amu/m^3]
      T0 = 1e6 # temperature, [K]
      V0 = [-700., 0.0, 0.0]*1e3 # velocity, [m/s]
      B0 = [-3., 3., 0.]*1e-9 # magnetic field, [T]

      varBG = BackgroundVariable(n0, T0, V0..., B0...)

      s = generate_signal(varBG; fsignal=0.1, signal=:alfven, dv=1e5, fsample=2, tend=10,
         dir="y")

      save_signal(s)

      d = readdlm("sw.dat", '\t', Float64, header=true)

      @test d[1][10, end-1] == 1.4167e-9
      rm("sw.dat")

      # Density pulse
      n0 = 2e6 # density, [amu/m^3]
      T0 = 5e5 # temperature, [K]
      V0 = [-600., 0.0, 0.0]*1e3 # [m/s]
      B0 = [5., 0., 0.]*1e-9 # [T]

      varBG = BackgroundVariable(n0, T0, V0..., B0...)

      fsignal = 1/150 # perturbation frequency, [Hz]
      signal  = :fast # perturbation type
      dv      = 5e4   # perturbation magnitude in velocity, [m/s]
      dn      = 1     # perturbation magnitude in density, [amu/cc]
      dir     = "y"   # perturbed direction(s)
      fsample = 2     # sampling frequency, [Hz]
      tstart  = 300   # [s]
      tend    = 450   # [s]

      s = generate_signal(varBG; fsignal, signal, dv, fsample, tstart, tend, dir)
      @test s.n[end] ≈ n0

   end
end

ext: jl


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

using Base: findprevnot, findnextnot
using Random, LinearAlgebra, Test

tc(r1::NTuple{N,Any}, r2::NTuple{N,Any}) where {N} = all(x->tc(x...), [zip(r1,r2)...])
tc(r1::BitArray{N}, r2::Union{BitArray{N},Array{Bool,N}}) where {N} = true
tc(r1::SubArray{Bool,N1,BitArray{N2}}, r2::SubArray{Bool,N1,<:Union{BitArray{N2},Array{Bool,N2}}}) where {N1,N2} = true
tc(r1::Transpose{Bool,BitVector}, r2::Union{Transpose{Bool,BitVector},Transpose{Bool,Vector{Bool}}}) = true
tc(r1::T, r2::T) where {T} = true
tc(r1,r2) = false

bitcheck(b::BitArray) = Test._check_bitarray_consistency(b)
bitcheck(x) = true
bcast_setindex!(b, x, I...) = (b[I...] .= x; b)

function check_bitop_call(ret_type, func, args...; kwargs...)
    r1 = func(args...; kwargs...)
    r2 = func(map(x->(isa(x, BitArray) ? Array(x) : x), args)...; kwargs...)
    ret_type ≢ nothing && !isa(r1, ret_type) && @show ret_type, typeof(r1)
    ret_type ≢ nothing && @test isa(r1, ret_type)
    @test tc(r1, r2)
    @test isequal(r1, ret_type ≡ nothing ? r2 : r2)
    @test bitcheck(r1)
end
macro check_bit_operation(ex, ret_type)
    @assert Meta.isexpr(ex, :call)
    Expr(:call, :check_bitop_call, esc(ret_type), map(esc, ex.args)...)
end
macro check_bit_operation(ex)
    @assert Meta.isexpr(ex, :call)
    Expr(:call, :check_bitop_call, nothing, map(esc, ex.args)...)
end

let t0 = time_ns()
    global timesofar
    function timesofar(str)
        return # no-op, comment to see timings
        t1 = time_ns()
        println(str, ": ", (t1-t0)/1e9, " seconds")
        t0 = t1
    end
end

@testset "empty bitvector" begin
    @test BitVector() == BitVector(undef, 0)
end

# vectors size
v1 = 260
# matrices size
n1, n2 = 17, 20
# arrays size
s1, s2, s3, s4 = 5, 8, 3, 7

allsizes = [((), BitArray{0}), ((v1,), BitVector),
            ((n1,n2), BitMatrix), ((s1,s2,s3,s4), BitArray{4})]

@testset "trues and falses for size $sz" for (sz,T) in allsizes
    a = falses(sz...)
    @test a == falses(sz)
    @test !any(a)
    @test sz == size(a)
    b = trues(sz...)
    @test b == trues(sz)
    @test all(b)
    @test sz == size(b)
    c = trues(size(a))
    @test all(c)
    @test !any(a)
    @test sz == size(c)
    d = falses(size(b))
    @test !any(d)
    @test all(b)
    @test sz == size(d)
    @test !isassigned(a, 0)
    @test !isassigned(b, 0)
    for ii in 1:prod(sz)
        @test isassigned(a, ii)
        @test isassigned(b, ii)
    end
    @test !isassigned(a, length(a) + 1)
    @test !isassigned(b, length(b) + 1)
end

@testset "Conversions for size $sz" for (sz, T) in allsizes
    b1 = rand!(falses(sz...))
    @test isequal(BitArray(Array(b1)), b1)
    @test isequal(convert(Array{Float64,ndims(b1)}, b1),
                  convert(Array{Float64,ndims(b1)}, Array(b1)))
    @test isequal(convert(AbstractArray{Float64,ndims(b1)}, b1),
                  convert(AbstractArray{Float64,ndims(b1)}, Array(b1)))

    i1 = rand!(zeros(Bool, sz...), false:true)
    @test isequal(Array(BitArray(i1)), i1)
end

timesofar("conversions")

@testset "utility functions" begin
    b1 = bitrand(v1)
    @test isequal(fill!(b1, true), trues(size(b1)))
    @test isequal(fill!(b1, false), falses(size(b1)))

    for (sz,T) in allsizes
        @test isequal(Array(trues(sz...)), fill(true, sz...))
        @test isequal(Array(falses(sz...)), fill(false, sz...))

        b1 = rand!(falses(sz...))
        @test isa(b1, T)

        @check_bit_operation length(b1) Int
        @check_bit_operation ndims(b1)  Int
        @check_bit_operation size(b1)   Tuple{Vararg{Int}}

        b2 = similar(b1)
        u1 = Array(b1)
        @check_bit_operation copyto!(b2, b1) T
        @check_bit_operation copyto!(b2, u1) T
    end

    @testset "copyto!" begin
        let b1 = trues(1)
            @test all(copyto!(b1, []))
        end
        for n in [1; 1023:1025]
            b1 = falses(n)
            for m in [1; 10; 1023:1025]
                u1 = fill(true, m)
                for fu! in [u->fill!(u, true), u->rand!(u)]
                    fu!(u1)
                    c1 = convert(Vector{Int}, u1)
                    for i1 in [1; 10; 53:65; 1013:1015; 1020:1025], i2 in [1; 3; 10; 511:513], l in [1; 5; 10; 511:513; 1023:1025]
                        for fb! in [b->fill!(b, false), b->rand!(b)]
                            fb!(b1)
                            if i1 < 1 || i1 > n || (i2 + l - 1 > m) || (i1 + l - 1 > n)
                                @test_throws BoundsError copyto!(b1, i1, u1, i2, l)
                            else
                                @check_bit_operation copyto!(b1, i1, u1, i2, l) BitArray
                                @check_bit_operation copyto!(b1, i1, c1, i2, l) BitArray
                            end
                        end
                    end
                end
            end
        end
    end

    @test_throws BoundsError size(trues(5), 0)

    @testset "reshape and resize!" begin
        b1 = bitrand(n1, n2)
        @check_bit_operation reshape(b1, (n2,n1)) BitMatrix
        @test_throws DimensionMismatch reshape(b1, (1,n1))

        @test @inferred(reshape(b1, n1*n2)) == @inferred(reshape(b1, (n1*n2,))) == @inferred(reshape(b1, Val(1))) == @inferred(reshape(b1, :))
        @test @inferred(reshape(b1, n1, n2)) === @inferred(reshape(b1, Val(2))) === b1
        @test @inferred(reshape(b1, n2, :)) == @inferred(reshape(b1, (n2, n1))) != @inferred(reshape(b1, Val(2)))

        b1 = bitrand(s1, s2, s3, s4)
        @check_bit_operation reshape(b1, (s3,s1,s2,s4)) BitArray{4}
        @test_throws DimensionMismatch reshape(b1, (1,n1))

        b1 = bitrand(v1)
        @test_throws BoundsError resize!(b1, -1)
        @check_bit_operation resize!(b1, v1 ÷ 2) BitVector
        gr(b) = (resize!(b, v1)[(v1÷2):end] .= 1; b)
        @check_bit_operation gr(b1) BitVector
    end

    @testset "sizeof (issue #7515)" begin
        @test sizeof(BitVector(undef, 64)) == 8
        @test sizeof(BitVector(undef, 65)) == 16
    end
end

timesofar("utils")

@testset "Constructors" begin
    @testset "non-Int dims constructors" begin
        b1 = BitVector(undef, Int32(v1))
        b2 = BitVector(undef, Int64(v1))
        @test size(b1) == size(b2)

        for c in [trues, falses]
            b1 = c(Int32(v1))
            b2 = c(Int64(v1))
            @test b1 == b2
        end
    end

    @testset "constructors from iterables" begin
        for g in ((x%7==3 for x = 1:v1),
                  (x%7==3 for x = 1:v1 if x>5),
                  ((x+y)%5==2 for x = 1:n1, y = 1:n2),
                  ((x+y+z+t)%5==2 for x = 1:s2, y = 1:s2, z = 1:s3, t = 1:s4),
                  ((x+y)%5==2 for x = 1:n1 for y = 1:n2))
            @test BitArray(g) == BitArray(collect(g))
        end
        @test_throws DimensionMismatch BitVector(false)
        @test_throws DimensionMismatch BitVector((iszero(i%4) for i in 1:n1, j in 1:n2))
        @test_throws DimensionMismatch BitMatrix((isodd(i) for i in 1:3))
    end

    @testset "constructor from NTuple" begin
        for nt in ((true, false, false), NTuple{0,Bool}(), (false,), (true,))
            @test BitVector(nt) == BitVector(collect(nt))
        end
    end

    @testset "one" begin
        @test Array(one(BitMatrix(undef, 2,2))) == Matrix(I, 2, 2)
        @test_throws DimensionMismatch one(BitMatrix(undef, 2,3))
    end

    # constructors should copy
    a = trues(3)
    @test BitArray(a) !== a
    @test BitVector(a) !== a

    # issue #24062
    @test_throws InexactError BitArray([0, 1, 2, 3])
    @test_throws MethodError BitArray([0, ""])

    # construction with poor inference
    f(c) = BitVector(c[1])
    @test @inferred(f(AbstractVector[[0,1]])) == [false, true]
end

timesofar("constructors")

@testset "Indexing" begin
    @testset "0d for size $sz" for (sz,T) in allsizes
        b1 = rand!(falses(sz...))
        if length(b1) == 1
            @check_bit_operation getindex(b1)         Bool
            @check_bit_operation setindex!(b1, true)  T
            @check_bit_operation setindex!(b1, false) T
        else
            @test_throws BoundsError getindex(b1)
            @test_throws BoundsError setindex!(b1, true)
            @test_throws BoundsError setindex!(b1, false)
        end
    end

    @testset "linear for size $sz" for (sz,T) in allsizes[2:end]
        l = *(sz...)
        b1 = rand!(falses(sz...))
        for j = 1:l
            @check_bit_operation getindex(b1, j) Bool
        end
        for j in [0, 1, 63, 64, 65, 127, 128, 129, 191, 192, 193, l-1, l]
            @check_bit_operation getindex(b1, 1:j)   BitVector
            @check_bit_operation getindex(b1, j+1:l) BitVector
        end
        for j in [1, 63, 64, 65, 127, 128, 129, div(l,2)]
            m1 = j:(l-j)
            @check_bit_operation getindex(b1, m1) BitVector
        end

        t1 = findall(bitrand(l))
        @check_bit_operation getindex(b1, t1)        BitVector

        for j = 1:l
            x = rand(Bool)
            @check_bit_operation setindex!(b1, x, j) T
        end

        y = rand(0.0:1.0)
        @check_bit_operation setindex!(b1, y, 100) T

        for j in [1, 63, 64, 65, 127, 128, 129, 191, 192, 193, l-1]
            x = rand(Bool)
            @check_bit_operation fill!(b1, x) T
            rand!(b1)
            @check_bit_operation bcast_setindex!(b1, x, 1:j)
            b2 = bitrand(j)
            for bb in (b2, view(b2, 1:j), view(Array{Any}(b2), :))
                @check_bit_operation setindex!(b1, bb, 1:j) T
            end
            x = rand(Bool)
            @check_bit_operation bcast_setindex!(b1, x, j+1:l) T
            b2 = bitrand(l-j)
            @check_bit_operation setindex!(b1, b2, j+1:l) T
        end
        for j in [1, 63, 64, 65, 127, 128, 129, div(l,2)]
            m1 = j:(l-j)
            x = rand(Bool)
            @check_bit_operation bcast_setindex!(b1, x, m1) T
            b2 = bitrand(length(m1))
            @check_bit_operation setindex!(b1, b2, m1) T
        end
        x = rand(Bool)
        @check_bit_operation bcast_setindex!(b1, x, 1:100) T
        b2 = bitrand(100)
        @check_bit_operation setindex!(b1, b2, 1:100) T

        y = rand(0.0:1.0)
        @check_bit_operation bcast_setindex!(b1, y, 1:100) T

        t1 = findall(bitrand(l))
        x = rand(Bool)
        @check_bit_operation bcast_setindex!(b1, x, t1) T
        b2 = bitrand(length(t1))
        @check_bit_operation setindex!(b1, b2, t1) T

        y = rand(0.0:1.0)
        @check_bit_operation bcast_setindex!(b1, y, t1) T
    end

    @testset "multidimensional" begin
        rand_m1m2() = rand(1:n1), rand(1:n2)
        b1 = bitrand(n1, n2)
        m1, m2 = rand_m1m2()
        b2 = bitrand(m1, m2)
        @check_bit_operation copyto!(b1, b2) BitMatrix

        function gen_getindex_data(c)
            m1, m2 = rand_m1m2()
            put!(c, (m1, m2, Bool))
            m1, m2 = rand_m1m2()
            put!(c, (m1, 1:m2, BitVector))
            put!(c, (m1, :, BitVector))
            m1, m2 = rand_m1m2()
            put!(c, (m1, randperm(m2), BitVector))
            m1, m2 = rand_m1m2()
            put!(c, (1:m1, m2, BitVector))
            put!(c, (:, m2, BitVector))
            m1, m2 = rand_m1m2()
            put!(c, (1:m1, 1:m2, BitMatrix))
            put!(c, (:, :, BitMatrix))
            m1, m2 = rand_m1m2()
            put!(c, (1:m1, randperm(m2), BitMatrix))
            put!(c, (:, randperm(m2), BitMatrix))
            m1, m2 = rand_m1m2()
            put!(c, (randperm(m1), m2, BitVector))
            m1, m2 = rand_m1m2()
            put!(c, (randperm(m1), 1:m2, BitMatrix))
            put!(c, (randperm(m1), :, BitMatrix))
            m1, m2 = rand_m1m2()
            put!(c, (randperm(m1), randperm(m2), BitMatrix))
        end

        for (k1, k2, T) in Channel(gen_getindex_data)
            # println(typeof(k1), " ", typeof(k2), " ", T) # uncomment to debug
            @check_bit_operation getindex(b1, k1, k2) T
        end

        m1, m2 = rand_m1m2()
        @check_bit_operation getindex(b1, 1:m1, m2, 1) BitVector
        @check_bit_operation getindex(b1, :, randperm(m2), 1) BitMatrix

        b1 = bitrand(s1, s2, s3, s4)
        function gen_getindex_data4(c)
            m1, m2, m3, m4 = (:, :, :, :)
            put!(c, (m1, m2, m3, m4, BitArray{4}))

            m1, m2, m3, m4 = (2, :, :, :)
            put!(c, (m1, m2, m3, m4, BitArray{3}))
            m1, m2, m3, m4 = (:, :, 2, :)
            put!(c, (m1, m2, m3, m4, BitArray{3}))
            m1, m2, m3, m4 = (:, :, :, 2)
            put!(c, (m1, m2, m3, m4, BitArray{3}))

            m1, m2, m3, m4 = (2, :, :, 2)
            put!(c, (m1, m2, m3, m4, BitArray{2}))
            m1, m2, m3, m4 = (:, 2, :, 2)
            put!(c, (m1, m2, m3, m4, BitArray{2}))
            m1, m2, m3, m4 = (2, :, 2, :)
            put!(c, (m1, m2, m3, m4, BitArray{2}))
            m1, m2, m3, m4 = (2, 2, :, :)
            put!(c, (m1, m2, m3, m4, BitArray{2}))

            m1, m2, m3, m4 = (:, 2, 2, 2)
            put!(c, (m1, m2, m3, m4, BitArray{1}))
            m1, m2, m3, m4 = (2, 2, :, 2)
            put!(c, (m1, m2, m3, m4, BitArray{1}))

            m1, m2, m3, m4 = (4, 3:7, 2:2, 2)
            put!(c, (m1, m2, m3, m4, BitArray{2}))
            m1, m2, m3, m4 = (1:2, 5, 1, 2:7)
            put!(c, (m1, m2, m3, m4, BitArray{2}))

            m1, m2, m3, m4 = (2:3, 2:7, 1:2, 4:6)
            put!(c, (m1, m2, m3, m4, BitArray{4}))
        end

        for (k1, k2, k3, k4, T) in Channel(gen_getindex_data4)
            #println(typeof(k1), " ", typeof(k2), " ", typeof(k3), " ", typeof(k4), " ", T) # uncomment to debug
            @check_bit_operation getindex(b1, k1, k2, k3, k4) T
        end

        b1 = bitrand(n1, n2)
        function gen_setindex_data(c)
            m1, m2 = rand_m1m2()
            put!(c, (rand(Bool), m1, m2))
            m1, m2 = rand_m1m2()
            put!(c, (rand(Bool), m1, 1:m2))
            put!(c, (rand(Bool), m1, :))
            put!(c, (bitrand(m2), m1, 1:m2))
            m1, m2 = rand_m1m2()
            put!(c, (rand(Bool), m1, randperm(m2)))
            put!(c, (bitrand(m2), m1, randperm(m2)))
            m1, m2 = rand_m1m2()
            put!(c, (rand(Bool), 1:m1, m2))
            put!(c, (rand(Bool), :, m2))
            put!(c, (bitrand(m1), 1:m1, m2))
            m1, m2 = rand_m1m2()
            put!(c, (rand(Bool), 1:m1, 1:m2))
            put!(c, (rand(Bool), :, :))
            put!(c, (bitrand(m1, m2), 1:m1, 1:m2))
            m1, m2 = rand_m1m2()
            put!(c, (rand(Bool), 1:m1, randperm(m2)))
            put!(c, (rand(Bool), :, randperm(m2)))
            put!(c, (bitrand(m1, m2), 1:m1, randperm(m2)))
            m1, m2 = rand_m1m2()
            put!(c, (rand(Bool), randperm(m1), m2))
            put!(c, (bitrand(m1), randperm(m1), m2))
            m1, m2 = rand_m1m2()
            put!(c, (rand(Bool), randperm(m1), 1:m2))
            put!(c, (rand(Bool), randperm(m1), :))
            put!(c, (bitrand(m1,m2), randperm(m1), 1:m2))
            m1, m2 = rand_m1m2()
            put!(c, (rand(Bool), randperm(m1), randperm(m2)))
            put!(c, (bitrand(m1,m2), randperm(m1), randperm(m2)))
        end

        for (b2, k1, k2) in Channel(gen_setindex_data)
            # println(typeof(b2), " ", typeof(k1), " ", typeof(k2)) # uncomment to debug
            if b2 isa AbstractArray
                for bb in (b2, view(b2, :), view(Array{Any}(b2), :))
                    @check_bit_operation setindex!(b1, bb, k1, k2) BitMatrix
                end
            else
                if k1 isa Integer && k2 isa Integer
                    @check_bit_operation setindex!(b1, b2, k1, k2) BitMatrix
                else
                    @check_bit_operation bcast_setindex!(b1, b2, k1, k2) BitMatrix
                end
            end
        end

        m1, m2 = rand_m1m2()
        b2 = bitrand(1, 1, m2)
        @check_bit_operation setindex!(b1, b2, m1, 1:m2) BitMatrix
        x = rand(Bool)
        b2 = bitrand(1, m2, 1)
        @check_bit_operation bcast_setindex!(b1, x, m1, 1:m2, 1)  BitMatrix
        @check_bit_operation setindex!(b1, b2, m1, 1:m2, 1) BitMatrix

        b1 = bitrand(s1, s2, s3, s4)
        function gen_setindex_data4(c)
            m1, m2, m3, m4 = (:, :, :, :)
            put!(c, (rand(Bool), m1, m2, m3, m4))
            put!(c, (bitrand(s1, s2, s3, s4), m1, m2, m3, m4))

            m1, m2, m3, m4 = (2, :, :, :)
            put!(c, (rand(Bool), m1, m2, m3, m4))
            put!(c, (bitrand(s2, s3, s4), m1, m2, m3, m4))
            m1, m2, m3, m4 = (:, :, 2, :)
            put!(c, (bitrand(s1, s2, s4), m1, m2, m3, m4))
            m1, m2, m3, m4 = (:, :, :, 2)
            put!(c, (rand(Bool), m1, m2, m3, m4))
            put!(c, (bitrand(s1, s2, s3), m1, m2, m3, m4))

            m1, m2, m3, m4 = (2, :, :, 2)
            put!(c, (rand(Bool), m1, m2, m3, m4))
            put!(c, (bitrand(s2, s3), m1, m2, m3, m4))
            m1, m2, m3, m4 = (:, 2, :, 2)
            put!(c, (rand(Bool), m1, m2, m3, m4))
            put!(c, (bitrand(s1, s3), m1, m2, m3, m4))
            m1, m2, m3, m4 = (2, :, 2, :)
            put!(c, (bitrand(s2, s4), m1, m2, m3, m4))
            m1, m2, m3, m4 = (:, 2, 2, :)
            put!(c, (bitrand(s1, s4), m1, m2, m3, m4))

            m1, m2, m3, m4 = (:, 2, 2, 2)
            put!(c, (bitrand(s1), m1, m2, m3, m4))
            m1, m2, m3, m4 = (2, 2, :, 2)
            put!(c, (bitrand(s3), m1, m2, m3, m4))

            m1, m2, m3, m4 = (4, 3:7, 2:2, 2)
            put!(c, (bitrand(5, 1), m1, m2, m3, m4))
            m1, m2, m3, m4 = (1:2, 5, 1, 2:7)
            put!(c, (rand(Bool), m1, m2, m3, m4))
            put!(c, (bitrand(2, 6), m1, m2, m3, m4))

            m1, m2, m3, m4 = (2:3, 2:7, 1:2, 4:6)
            put!(c, (bitrand(2, 6, 2, 3), m1, m2, m3, m4))
        end

        for (b2, k1, k2, k3, k4) in Channel(gen_setindex_data4)
            # println(typeof(b2), " ", typeof(k1), " ", typeof(k2), " ", typeof(k3), " ", typeof(k4)) # uncomment to debug
            if b2 isa Bool
                @check_bit_operation bcast_setindex!(b1, b2, k1, k2, k3, k4) BitArray{4}
            else
                @check_bit_operation setindex!(b1, b2, k1, k2, k3, k4) BitArray{4}
            end
        end

        for p1 = [rand(1:v1) 1 63 64 65 191 192 193]
            for p2 = [rand(1:v1) 1 63 64 65 191 192 193]
                for n = 0 : min(v1 - p1 + 1, v1 - p2 + 1)
                    b1 = bitrand(v1)
                    b2 = bitrand(v1)
                    @check_bit_operation copyto!(b1, p1, b2, p2, n) BitVector
                end
            end
        end
    end

    @testset "logical indexing" begin
        b1 = bitrand(n1, n2)

        t1 = bitrand(n1, n2)
        @test isequal(Array(b1[t1]), Array(b1)[t1])
        @test isequal(Array(b1[t1]), Array(b1)[Array(t1)])

        t1 = bitrand(n1)
        t2 = bitrand(n2)
        @test isequal(Array(b1[t1, t2]), Array(b1)[t1, t2])
        @test isequal(Array(b1[t1, t2]), Array(b1)[Array(t1), Array(t2)])

        b1 = bitrand(n1, n2)
        t1 = bitrand(n1, n2)
        @check_bit_operation bcast_setindex!(b1, true, t1) BitMatrix

        t1 = bitrand(n1, n2)
        b2 = bitrand(count(t1))
        @check_bit_operation setindex!(b1, b2, t1) BitMatrix

        m1 = rand(1:n1)
        m2 = rand(1:n2)
        t1 = bitrand(n1)
        b2 = bitrand(count(t1), m2)
        k2 = randperm(m2)
        @check_bit_operation setindex!(b1, b2, t1, 1:m2)       BitMatrix
        @check_bit_operation setindex!(b1, b2, t1, n2-m2+1:n2) BitMatrix
        @check_bit_operation setindex!(b1, b2, t1, k2)         BitMatrix

        t2 = bitrand(n2)
        b2 = bitrand(m1, count(t2))
        k1 = randperm(m1)
        @check_bit_operation setindex!(b1, b2, 1:m1, t2)       BitMatrix
        @check_bit_operation setindex!(b1, b2, n1-m1+1:n1, t2) BitMatrix
        @check_bit_operation setindex!(b1, b2, k1, t2)         BitMatrix
    end
end

timesofar("indexing")

@testset "Deque Functionality" begin
    b1 = BitVector()
    i1 = Bool[]
    for m = 1:v1
        x = rand(Bool)
        push!(b1, x)
        push!(i1, x)
        @test isequal(Array(b1), i1)
        @test bitcheck(b1)
    end

    for m1 = 0:v1, m2 = [0, 1, 63, 64, 65, 127, 128, 129]
        b1 = bitrand(m1)
        b2 = bitrand(m2)
        i1 = Array(b1)
        i2 = Array(b2)
        # Append from array
        @test isequal(Array(append!(b1, b2)), append!(i1, i2))
        @test isequal(Array(append!(b1, i2)), append!(i1, b2))
        @test bitcheck(b1)
        # Append from HasLength iterator
        @test isequal(Array(append!(b1, (v for v in b2))), append!(i1, i2))
        @test isequal(Array(append!(b1, (v for v in i2))), append!(i1, b2))
        @test bitcheck(b1)
        # Append from SizeUnknown iterator
        @test isequal(Array(append!(b1, (v for v in b2 if true))), append!(i1, i2))
        @test isequal(Array(append!(b1, (v for v in i2 if true))), append!(i1, b2))
        @test bitcheck(b1)
    end

    for m1 = 0:v1, m2 = [0, 1, 63, 64, 65, 127, 128, 129]
        b1 = bitrand(m1)
        b2 = bitrand(m2)
        i1 = Array(b1)
        i2 = Array(b2)
        @test isequal(Array(prepend!(b1, b2)), prepend!(i1, i2))
        @test isequal(Array(prepend!(b1, i2)), prepend!(i1, b2))
        @test bitcheck(b1)
    end

    b1 = bitrand(v1)
    i1 = Array(b1)
    for m = 1:v1
        jb = pop!(b1)
        ji = pop!(i1)
        @test jb == ji
        @test isequal(Array(b1), i1)
        @test bitcheck(b1)
    end

    @test length(b1) == 0

    b1 = BitVector()
    i1 = Bool[]
    for m = 1:v1
        x = rand(Bool)
        pushfirst!(b1, x)
        pushfirst!(i1, x)
        @test isequal(Array(b1), i1)
        @test bitcheck(b1)
    end

    b1 = bitrand(v1)
    i1 = Array(b1)
    for m = 1:v1
        jb = popfirst!(b1)
        ji = popfirst!(i1)
        @test jb == ji
        @test isequal(Array(b1), i1)
        @test bitcheck(b1)
    end
    @test length(b1) == 0

    b1 = BitVector()
    @test_throws BoundsError insert!(b1, 2, false)
    @test_throws BoundsError insert!(b1, 0, false)
    i1 = Array(b1)
    for m = 1:v1
        j = rand(1:m)
        x = rand(Bool)
        @test insert!(b1, j, x) === b1
        insert!(i1, j, x)
        @test isequal(Array(b1), i1)
        @test bitcheck(b1)
    end

    b1 = bitrand(v1)
    i1 = Array(b1)
    for j in [63, 64, 65, 127, 128, 129, 191, 192, 193]
        x = rand(0:1)
        @test insert!(b1, j, x) === b1
        insert!(i1, j, x)
        @test isequal(Array(b1), i1)
        @test bitcheck(b1)
    end

    b1 = bitrand(v1)
    i1 = Array(b1)
    for m = v1:-1:1
        j = rand(1:m)
        b = splice!(b1, j)
        i = splice!(i1, j)
        @test isequal(Array(b1), i1)
        @test bitcheck(b1)
        @test b == i
        @test bitcheck(b)
    end
    @test length(b1) == 0

    b1 = bitrand(v1)
    i1 = Array(b1)
    for m = v1:-1:1
        j = rand(1:m)
        deleteat!(b1, j)
        deleteat!(i1, j)
        @test isequal(Array(b1), i1)
        @test bitcheck(b1)
    end
    @test length(b1) == 0

    b1 = bitrand(v1)
    @test_throws ArgumentError deleteat!(b1, [1, 1, 2])
    @test_throws BoundsError deleteat!(b1, [1, length(b1)+1])
    @test_throws BoundsError deleteat!(b1, [length(b1)+rand(1:100)])
    @test_throws BoundsError deleteat!(bitrand(1), [-1, 0, 1])

    @test_throws BoundsError deleteat!(BitVector(), 1)
    @test_throws BoundsError deleteat!(BitVector(), [1])
    @test_throws BoundsError deleteat!(BitVector(), [2])
    @test deleteat!(BitVector(), []) == BitVector()
    @test deleteat!(BitVector(), Bool[]) == BitVector()

    b1 = bitrand(v1)
    i1 = Array(b1)
    for j in [63, 64, 65, 127, 128, 129, 191, 192, 193]
        b = splice!(b1, j)
        i = splice!(i1, j)
        @test isequal(Array(b1), i1)
        @test bitcheck(b1)
        @test b == i
        @test bitcheck(b)
    end

    b1 = bitrand(v1)
    i1 = Array(b1)
    for j in [63, 64, 65, 127, 128, 129, 191, 192, 193]
        deleteat!(b1, j)
        deleteat!(i1, j)
        @test isequal(Array(b1), i1)
        @test bitcheck(b1)
    end

    b1 = bitrand(v1)
    i1 = Array(b1)
    for m1 = 1:v1, m2 = m1:v1
        b2 = copy(b1)
        i2 = copy(i1)
        b = splice!(b2, m1:m2)
        i = splice!(i2, m1:m2)
        @test isequal(Array(b2), i2)
        @test bitcheck(b2)
        @test b == i
        @test bitcheck(b)
    end

    b1 = bitrand(v1)
    i1 = Array(b1)
    for m1 = 1:v1, m2 = m1:v1
        b2 = copy(b1)
        i2 = copy(i1)
        deleteat!(b2, m1:m2)
        deleteat!(i2, m1:m2)
        @test isequal(Array(b2), i2)
        @test bitcheck(b2)
    end

    b1 = bitrand(v1)
    i1 = Array(b1)
    for m1 = [1:7:v1; v1+1], m2 = [(m1-1):5:(v1-1); v1], v2 = [0, 1, 63, 64, 65, 127, 128, 129, 191, 192, 193, rand(1:v1)]
        b2 = copy(b1)
        i2 = copy(i1)
        b3 = bitrand(v2)
        i3 = Array(b3)
        b = splice!(b2, m1:m2, b3)
        i = splice!(i2, m1:m2, i3)
        @test isequal(Array(b2), i2)
        @test bitcheck(b2)
        @test b == i
        @test bitcheck(b)

        b2 = copy(b1)
        i2 = copy(i1)
        i3 = map(Int, bitrand(v2))
        b = splice!(b2, m1:m2, i3)
        i = splice!(i2, m1:m2, i3)
        @test isequal(Array(b2), i2)
        @test bitcheck(b2)
        @test b == i
        @test bitcheck(b)
        b2 = copy(b1)
        i2 = copy(i1)
        i3 = Dict(j => rand(0:1) for j = 1:v2)
        b = splice!(b2, m1:m2, values(i3))
        i = splice!(i2, m1:m2, values(i3))
        @test isequal(Array(b2), i2)
        @test bitcheck(b2)
        @test b == i
        @test bitcheck(b)
    end

    b1 = bitrand(v1)
    i1 = Array(b1)
    for m1 = 1:v1, v2 = [0, 1, 63, 64, 65, 127, 128, 129, 191, 192, 193, rand(1:v1)]
        b2 = copy(b1)
        i2 = copy(i1)
        b3 = bitrand(v2)
        i3 = Array(b3)
        b = splice!(b2, m1, b3)
        i = splice!(i2, m1, i3)
        @test isequal(Array(b2), i2)
        @test bitcheck(b2)
        @test b == i
        @test bitcheck(b)
    end

    b1 = bitrand(v1)
    i1 = Array(b1)
    for m1 = 1:(v1-1), m2 = (m1+1):v1
        locs = bitrand(m2-m1+1)
        m = [m1:m2...][locs]
        b2 = copy(b1)
        i2 = copy(i1)
        deleteat!(b2, m)
        deleteat!(i2, m)
        @test isequal(Array(b2), i2)
        @test bitcheck(b2)
    end

    b1 = bitrand(v1)
    i1 = Array(b1)
    empty!(b1)
    empty!(i1)
    @test isequal(Array(b1), i1)
    @test bitcheck(b1)
end

timesofar("dequeue")

@testset "Unary operators" begin
    b1 = bitrand(n1, n2)
    @check_bit_operation broadcast(~, b1)  BitMatrix
    @check_bit_operation broadcast(!, b1)  BitMatrix
    @check_bit_operation (-)(b1)  Matrix{Int}
    @check_bit_operation broadcast(sign, b1) BitMatrix
    @check_bit_operation real(b1) BitMatrix
    @check_bit_operation imag(b1) BitMatrix
    @check_bit_operation conj(b1) BitMatrix

    b0 = falses(0)
    @check_bit_operation broadcast(~, b0)  BitVector
    @check_bit_operation broadcast(!, b0)  BitVector
    @check_bit_operation (-)(b0)  Vector{Int}
    @check_bit_operation broadcast(sign, b0) BitVector

    @testset "in-place .!" begin
        b1 = bitrand(n1, n2)
        i1 = Array(b1)
        b1 .= .!b1
        @test b1 == .~i1
        @test bitcheck(b1)
    end
end

timesofar("unary arithmetic")

@testset "Binary arithmetic operators" begin
    @testset "Matrix{Bool}/Matrix{Bool}" begin
        b1 = bitrand(n1, n2)
        b2 = bitrand(n1, n2)
        @check_bit_operation broadcast(&, b1, b2)  BitMatrix
        @check_bit_operation broadcast(|, b1, b2)  BitMatrix
        @check_bit_operation broadcast(xor, b1, b2)  BitMatrix
        @check_bit_operation broadcast(nand, b1, b2)  BitMatrix
        @check_bit_operation broadcast(nor, b1, b2)  BitMatrix
        @check_bit_operation (+)(b1, b2)  Matrix{Int}
        @check_bit_operation (-)(b1, b2)  Matrix{Int}
        @check_bit_operation broadcast(*, b1, b2) BitMatrix
        @check_bit_operation broadcast(/, b1, b2) Matrix{Float64}
        @check_bit_operation broadcast(^, b1, b2) BitMatrix
        @check_bit_operation (/)(b1,1) Matrix{Float64}

        b2 = trues(n1, n2)
        @check_bit_operation broadcast(div, b1, b2) BitMatrix
        @check_bit_operation broadcast(mod, b1, b2) BitMatrix
        @check_bit_operation broadcast(div, b1, Array(b2)) BitMatrix
        @check_bit_operation broadcast(mod, b1, Array(b2)) BitMatrix
        @check_bit_operation broadcast(div, Array(b1), b2) BitMatrix
        @check_bit_operation broadcast(mod, Array(b1), b2) BitMatrix

        b1 = bitrand(n1, n1)
        while abs(det(Array{Float64}(b1))) ≤ 1e-6
            b1 = bitrand(n1, n1)
        end
        b2 = bitrand(n1, n1)

        @check_bit_operation (*)(b1, b2) Matrix{Int}
        @check_bit_operation (/)(b1, b1) Matrix{Float64}
        @check_bit_operation (\)(b1, b1) Matrix{Float64}

        b0 = falses(0)
        @check_bit_operation broadcast(&, b0, b0)  BitVector
        @check_bit_operation broadcast(|, b0, b0)  BitVector
        @check_bit_operation broadcast(xor, b0, b0)  BitVector
        @check_bit_operation broadcast(nand, b0, b0)  BitVector
        @check_bit_operation broadcast(nor, b0, b0)  BitVector
        @check_bit_operation broadcast(*, b0, b0) BitVector
        @check_bit_operation (*)(b0, b0') BitMatrix
    end

    @testset "Matrix{Bool}/Matrix{Int}" begin
        b1 = bitrand(n1, n2)
        i2 = rand(1:10, n1, n2)
        @check_bit_operation broadcast(&, b1, i2)  Matrix{Int}
        @check_bit_operation broadcast(|, b1, i2)  Matrix{Int}
        @check_bit_operation broadcast(xor, b1, i2)  Matrix{Int}
        @check_bit_operation broadcast(nand, b1, i2)  Matrix{Int}
        @check_bit_operation broadcast(nor, b1, i2)  Matrix{Int}
        @check_bit_operation (+)(b1, i2)  Matrix{Int}
        @check_bit_operation (-)(b1, i2)  Matrix{Int}
        @check_bit_operation broadcast(*, b1, i2) Matrix{Int}
        @check_bit_operation broadcast(/, b1, i2) Matrix{Float64}
        @check_bit_operation broadcast(^, b1, i2) BitMatrix
        @check_bit_operation broadcast(div, b1, i2)  Matrix{Int}
        @check_bit_operation broadcast(mod, b1, i2)  Matrix{Int}
    end

    @testset "Matrix{Bool}/Matrix{Float64}" begin
        b1 = bitrand(n1, n2)
        f2 = 1.0 .+ rand(n1, n2)
        @check_bit_operation broadcast(*, b1, f2) Matrix{Float64}
        @check_bit_operation broadcast(/, b1, f2) Matrix{Float64}
        @check_bit_operation broadcast(^, b1, f2) Matrix{Float64}
        @check_bit_operation broadcast(div, b1, f2)  Matrix{Float64}
        @check_bit_operation broadcast(mod, b1, f2)  Matrix{Float64}
    end

    @testset "Number/Matrix" begin
        b2 = bitrand(n1, n2)
        i1 = rand(1:10)
        u1 = UInt8(i1)
        f1 = Float64(i1)
        ci1 = complex(i1)
        cu1 = complex(u1)
        cf1 = complex(f1)

        @check_bit_operation broadcast(&, i1, b2)  Matrix{Int}
        @check_bit_operation broadcast(|, i1, b2)  Matrix{Int}
        @check_bit_operation broadcast(xor, i1, b2)  Matrix{Int}
        @check_bit_operation broadcast(nand, i1, b2)  Matrix{Int}
        @check_bit_operation broadcast(nor, i1, b2)  Matrix{Int}
        @check_bit_operation broadcast(+, i1, b2)  Matrix{Int}
        @check_bit_operation broadcast(-, i1, b2)  Matrix{Int}
        @check_bit_operation broadcast(*, i1, b2) Matrix{Int}

        @check_bit_operation broadcast(&, u1, b2)  Matrix{UInt8}
        @check_bit_operation broadcast(|, u1, b2)  Matrix{UInt8}
        @check_bit_operation broadcast(xor, u1, b2)  Matrix{UInt8}
        @check_bit_operation broadcast(nand, u1, b2)  Matrix{UInt8}
        @check_bit_operation broadcast(nor, u1, b2)  Matrix{UInt8}
        @check_bit_operation broadcast(+, u1, b2)  Matrix{UInt8}
        @check_bit_operation broadcast(-, u1, b2)  Matrix{UInt8}
        @check_bit_operation broadcast(*, u1, b2) Matrix{UInt8}

        for (x1,t1) = [(f1, Float64),
                       (ci1, Complex{Int}),
                       (cu1, Complex{UInt8}),
                       (cf1, ComplexF64)]
            @check_bit_operation broadcast(+, x1, b2)  Matrix{t1}
            @check_bit_operation broadcast(-, x1, b2)  Matrix{t1}
            @check_bit_operation broadcast(*, x1, b2) Matrix{t1}
        end

        b2 = trues(n1, n2)
        @check_bit_operation broadcast(/, true, b2)  Matrix{Float64}
        @check_bit_operation broadcast(div, true, b2)   BitMatrix
        @check_bit_operation broadcast(mod, true, b2)   BitMatrix
        @check_bit_operation broadcast(/, false, b2) Matrix{Float64}
        @check_bit_operation broadcast(div, false, b2)  BitMatrix
        @check_bit_operation broadcast(mod, false, b2)  BitMatrix

        @check_bit_operation broadcast(/, i1, b2) Matrix{Float64}
        @check_bit_operation broadcast(div, i1, b2)  Matrix{Int}
        @check_bit_operation broadcast(mod, i1, b2)  Matrix{Int}

        @check_bit_operation broadcast(/, u1, b2) Matrix{Float64}
        @check_bit_operation broadcast(div, u1, b2)  Matrix{UInt8}
        @check_bit_operation broadcast(mod, u1, b2)  Matrix{UInt8}

        @check_bit_operation broadcast(/, f1, b2) Matrix{Float64}
        @check_bit_operation broadcast(div, f1, b2)  Matrix{Float64}
        @check_bit_operation broadcast(mod, f1, b2)  Matrix{Float64}

        @check_bit_operation broadcast(/, ci1, b2) Matrix{ComplexF64}
        @check_bit_operation broadcast(/, cu1, b2) Matrix{ComplexF64}
        @check_bit_operation broadcast(/, cf1, b2) Matrix{ComplexF64}

        b2 = bitrand(n1, n2)
        @check_bit_operation broadcast(^, false, b2) BitMatrix
        @check_bit_operation broadcast(^, true, b2)  BitMatrix
        @check_bit_operation broadcast(^, 0x0, b2)   Matrix{UInt8}
        @check_bit_operation broadcast(^, 0x1, b2)   Matrix{UInt8}
        @check_bit_operation broadcast(^, -1, b2)    Matrix{Int}
        @check_bit_operation broadcast(^, 0, b2)     Matrix{Int}
        @check_bit_operation broadcast(^, 1, b2)     Matrix{Int}
        @check_bit_operation broadcast(^, 0.0, b2)   Matrix{Float64}
        @check_bit_operation broadcast(^, 1.0, b2)   Matrix{Float64}
        @check_bit_operation broadcast(^, 0.0im, b2) Matrix{ComplexF64}
        @check_bit_operation broadcast(^, 1.0im, b2) Matrix{ComplexF64}
        @check_bit_operation broadcast(^, 0im, b2)   Matrix{Complex{Int}}
        @check_bit_operation broadcast(^, 1im, b2)   Matrix{Complex{Int}}
        @check_bit_operation broadcast(^, 0x0*im, b2) Matrix{Complex{UInt8}}
        @check_bit_operation broadcast(^, 0x1*im, b2) Matrix{Complex{UInt8}}
    end

    @testset "Matrix/Number" begin
        b1 = bitrand(n1, n2)
        i2 = rand(1:10)
        u2 = UInt8(i2)
        f2 = Float64(i2)
        ci2 = complex(i2)
        cu2 = complex(u2)
        cf2 = complex(f2)
        b2 = Array(bitrand(n1,n2))

        @check_bit_operation broadcast(&, b1, true)   BitMatrix
        @check_bit_operation broadcast(&, b1, false)  BitMatrix
        @check_bit_operation broadcast(&, true, b1)   BitMatrix
        @check_bit_operation broadcast(&, false, b1)  BitMatrix
        @check_bit_operation broadcast(|, b1, true)   BitMatrix
        @check_bit_operation broadcast(|, b1, false)  BitMatrix
        @check_bit_operation broadcast(|, true, b1)   BitMatrix
        @check_bit_operation broadcast(|, false, b1)  BitMatrix
        @check_bit_operation broadcast(xor, b1, true)   BitMatrix
        @check_bit_operation broadcast(xor, b1, false)  BitMatrix
        @check_bit_operation broadcast(xor, true, b1)   BitMatrix
        @check_bit_operation broadcast(xor, false, b1)  BitMatrix
        @check_bit_operation broadcast(nand, b1, true)   BitMatrix
        @check_bit_operation broadcast(nand, b1, false)  BitMatrix
        @check_bit_operation broadcast(nand, true, b1)   BitMatrix
        @check_bit_operation broadcast(nand, false, b1)  BitMatrix
        @check_bit_operation broadcast(nor, b1, true)   BitMatrix
        @check_bit_operation broadcast(nor, b1, false)  BitMatrix
        @check_bit_operation broadcast(nor, true, b1)   BitMatrix
        @check_bit_operation broadcast(nor, false, b1)  BitMatrix
        @check_bit_operation broadcast(+, b1, true)   Matrix{Int}
        @check_bit_operation broadcast(+, b1, false)  Matrix{Int}
        @check_bit_operation broadcast(-, b1, true)   Matrix{Int}
        @check_bit_operation broadcast(-, b1, false)  Matrix{Int}
        @check_bit_operation broadcast(*, b1, true)  BitMatrix
        @check_bit_operation broadcast(*, b1, false) BitMatrix
        @check_bit_operation broadcast(*, true, b1)  BitMatrix
        @check_bit_operation broadcast(*, false, b1) BitMatrix
        @check_bit_operation broadcast(/, b1, true)  Matrix{Float64}
        @check_bit_operation broadcast(/, b1, false) Matrix{Float64}
        @check_bit_operation broadcast(div, b1, true)   BitMatrix
        @check_bit_operation broadcast(mod,b1, true)    BitMatrix

        @check_bit_operation broadcast(&, b1, b2)  BitMatrix
        @check_bit_operation broadcast(|, b1, b2)  BitMatrix
        @check_bit_operation broadcast(xor, b1, b2)  BitMatrix
        @check_bit_operation broadcast(nand, b1, b2)  BitMatrix
        @check_bit_operation broadcast(nor, b1, b2)  BitMatrix
        @check_bit_operation broadcast(&, b2, b1)  BitMatrix
        @check_bit_operation broadcast(|, b2, b1)  BitMatrix
        @check_bit_operation broadcast(xor, b2, b1)  BitMatrix
        @check_bit_operation broadcast(nand, b2, b1)  BitMatrix
        @check_bit_operation broadcast(nor, b2, b1)  BitMatrix
        @check_bit_operation broadcast(&, b1, i2)  Matrix{Int}
        @check_bit_operation broadcast(|, b1, i2)  Matrix{Int}
        @check_bit_operation broadcast(xor, b1, i2)  Matrix{Int}
        @check_bit_operation broadcast(nand, b1, i2)  Matrix{Int}
        @check_bit_operation broadcast(nor, b1, i2)  Matrix{Int}
        @check_bit_operation broadcast(+, b1, i2)  Matrix{Int}
        @check_bit_operation broadcast(-, b1, i2)  Matrix{Int}
        @check_bit_operation broadcast(*, b1, i2) Matrix{Int}
        @check_bit_operation broadcast(/, b1, i2) Matrix{Float64}
        @check_bit_operation broadcast(div, b1, i2)  Matrix{Int}
        @check_bit_operation broadcast(mod, b1, i2)  Matrix{Int}

        @check_bit_operation broadcast(&, b1, u2)  Matrix{UInt8}
        @check_bit_operation broadcast(|, b1, u2)  Matrix{UInt8}
        @check_bit_operation broadcast(xor, b1, u2)  Matrix{UInt8}
        @check_bit_operation broadcast(nand, b1, u2)  Matrix{UInt8}
        @check_bit_operation broadcast(nor, b1, u2)  Matrix{UInt8}
        @check_bit_operation broadcast(+, b1, u2)  Matrix{UInt8}
        @check_bit_operation broadcast(-, b1, u2)  Matrix{UInt8}
        @check_bit_operation broadcast(*, b1, u2) Matrix{UInt8}
        @check_bit_operation broadcast(/, b1, u2) Matrix{Float64}
        @check_bit_operation broadcast(div, b1, u2)  Matrix{UInt8}
        @check_bit_operation broadcast(mod, b1, u2)  Matrix{UInt8}

        @check_bit_operation broadcast(+, b1, f2)  Matrix{Float64}
        @check_bit_operation broadcast(-, b1, f2)  Matrix{Float64}
        @check_bit_operation broadcast(*, b1, f2) Matrix{Float64}
        @check_bit_operation broadcast(/, b1, f2) Matrix{Float64}
        @check_bit_operation broadcast(div, b1, f2)  Matrix{Float64}
        @check_bit_operation broadcast(mod, b1, f2)  Matrix{Float64}

        @check_bit_operation broadcast(+, b1, ci2)  Matrix{Complex{Int}}
        @check_bit_operation broadcast(-, b1, ci2)  Matrix{Complex{Int}}
        @check_bit_operation broadcast(*, b1, ci2) Matrix{Complex{Int}}
        @check_bit_operation broadcast(/, b1, ci2) Matrix{ComplexF64}

        @check_bit_operation broadcast(+, b1, cu2)  Matrix{Complex{UInt8}}
        @check_bit_operation broadcast(-, b1, cu2)  Matrix{Complex{UInt8}}
        @check_bit_operation broadcast(*, b1, cu2) Matrix{Complex{UInt8}}
        @check_bit_operation broadcast(/, b1, cu2) Matrix{ComplexF64}

        @check_bit_operation broadcast(+, b1, cf2)  Matrix{ComplexF64}
        @check_bit_operation broadcast(-, b1, cf2)  Matrix{ComplexF64}
        @check_bit_operation broadcast(*, b1, cf2) Matrix{ComplexF64}
        @check_bit_operation broadcast(/, b1, cf2) Matrix{ComplexF64}

        @check_bit_operation broadcast(^, b1, false) BitMatrix
        @check_bit_operation broadcast(^, b1, true)  BitMatrix
        @check_bit_operation broadcast(^, b1, 0x0)   BitMatrix
        @check_bit_operation broadcast(^, b1, 0x1)   BitMatrix
        @check_bit_operation broadcast(^, b1, 0)     BitMatrix
        @check_bit_operation broadcast(^, b1, 1)     BitMatrix
        @check_bit_operation broadcast(^, b1, -1.0)  Matrix{Float64}
        @check_bit_operation broadcast(^, b1, 0.0)   Matrix{Float64}
        @check_bit_operation broadcast(^, b1, 1.0)   Matrix{Float64}
        @check_bit_operation broadcast(^, b1, 0.0im) Matrix{ComplexF64}
        @check_bit_operation broadcast(^, b1, 0x0*im) Matrix{ComplexF64}
        @check_bit_operation broadcast(^, b1, 0im)   Matrix{ComplexF64}
        @test_throws DomainError broadcast(^, b1, -1)

        b1 = trues(n1, n2)
        @check_bit_operation broadcast(^, b1, -1.0im) Matrix{ComplexF64}
        @check_bit_operation broadcast(^, b1, 1.0im)  Matrix{ComplexF64}
        @check_bit_operation broadcast(^, b1, -1im)   Matrix{ComplexF64}
        @check_bit_operation broadcast(^, b1, 1im)    Matrix{ComplexF64}
        @check_bit_operation broadcast(^, b1, 0x1*im)  Matrix{ComplexF64}
    end

    @testset "Matrix/Vector" begin
        b1 = bitrand(n1, n2)
        b2 = bitrand(n1)
        b3 = bitrand(n2)

        @check_bit_operation broadcast(&, b1, b2)             BitMatrix
        @check_bit_operation broadcast(&, b1, transpose(b3))  BitMatrix
        @check_bit_operation broadcast(&, b2, b1)             BitMatrix
        @check_bit_operation broadcast(&, transpose(b3), b1)  BitMatrix
        @check_bit_operation broadcast(|, b1, b2)             BitMatrix
        @check_bit_operation broadcast(|, b1, transpose(b3))  BitMatrix
        @check_bit_operation broadcast(|, b2, b1)             BitMatrix
        @check_bit_operation broadcast(|, transpose(b3), b1)  BitMatrix
        @check_bit_operation broadcast(xor, b1, b2)             BitMatrix
        @check_bit_operation broadcast(xor, b1, transpose(b3))  BitMatrix
        @check_bit_operation broadcast(xor, b2, b1)             BitMatrix
        @check_bit_operation broadcast(xor, transpose(b3), b1)  BitMatrix
        @check_bit_operation broadcast(nand, b1, b2)             BitMatrix
        @check_bit_operation broadcast(nand, b1, transpose(b3))  BitMatrix
        @check_bit_operation broadcast(nand, b2, b1)             BitMatrix
        @check_bit_operation broadcast(nand, transpose(b3), b1)  BitMatrix
        @check_bit_operation broadcast(nor, b1, b2)             BitMatrix
        @check_bit_operation broadcast(nor, b1, transpose(b3))  BitMatrix
        @check_bit_operation broadcast(nor, b2, b1)             BitMatrix
        @check_bit_operation broadcast(nor, transpose(b3), b1)  BitMatrix
        @check_bit_operation broadcast(+, b1, b2)             Matrix{Int}
        @check_bit_operation broadcast(+, b1, transpose(b3))  Matrix{Int}
        @check_bit_operation broadcast(+, b2, b1)             Matrix{Int}
        @check_bit_operation broadcast(+, transpose(b3), b1)  Matrix{Int}
        @check_bit_operation broadcast(-, b1, b2)             Matrix{Int}
        @check_bit_operation broadcast(-, b1, transpose(b3))  Matrix{Int}
        @check_bit_operation broadcast(-, b2, b1)             Matrix{Int}
        @check_bit_operation broadcast(-, transpose(b3), b1)  Matrix{Int}
        @check_bit_operation broadcast(*, b1, b2)             BitMatrix
        @check_bit_operation broadcast(*, b1, transpose(b3))  BitMatrix
        @check_bit_operation broadcast(*, b2, b1)             BitMatrix
        @check_bit_operation broadcast(*, transpose(b3), b1)  BitMatrix
        @check_bit_operation broadcast(/, b1, b2)             Matrix{Float64}
        @check_bit_operation broadcast(/, b1, transpose(b3))  Matrix{Float64}
        @check_bit_operation broadcast(/, b2, b1)             Matrix{Float64}
        @check_bit_operation broadcast(/, transpose(b3), b1)  Matrix{Float64}
    end
end

timesofar("binary arithmetic")

@testset "Binary comparison operators" begin
    b1 = bitrand(n1, n2)
    b2 = bitrand(n1, n2)
    @check_bit_operation broadcast(==, b1, b2) BitMatrix
    @check_bit_operation broadcast(!=, b1, b2) BitMatrix
    @check_bit_operation broadcast(<, b1, b2) BitMatrix
    @check_bit_operation broadcast(<=, b1, b2) BitMatrix
end

timesofar("binary comparison")

@testset "Data movement" begin
    b1 = bitrand(s1, s2, s3, s4)
    for d = 1:4
        j = rand(1:size(b1, d))
        #for j = 1 : size(b1, d)
            @check_bit_operation selectdim(b1, d, j) SubArray{Bool, 3, BitArray{4}}
        #end
        @check_bit_operation reverse(b1, dims=d) BitArray{4}
    end
    @test_throws ArgumentError reverse(b1, dims=5)

    b1 = bitrand(n1, n2)
    for k = 1:4
        @check_bit_operation rotl90(b1, k) BitMatrix
    end

    for m = 0:v1
        b1 = bitrand(m)
        @check_bit_operation reverse(b1) BitVector
    end

    b1 = bitrand(v1)
    for m = [rand(1:v1)-1, 0, 1, 63, 64, 65, 191, 192, 193, v1-1]
        @test isequal(b1 << m, [ b1[m+1:end]; falses(m) ])
        @test isequal(b1 >>> m, [ falses(m); b1[1:end-m] ])
        @test isequal(b1 << -m, b1 >> m)
        @test isequal(b1 >>> -m, b1 << m)
        @test isequal(circshift(b1, -m), [ b1[m+1:end]; b1[1:m] ])
        @test isequal(circshift(b1, m), [ b1[end-m+1:end]; b1[1:end-m] ])
        @test isequal(circshift(b1, m), circshift(b1, m - length(b1)))
    end

    b = bitrand(v1)
    i = bitrand(v1)
    for m = [rand(1:v1), 63, 64, 65, 191, 192, 193, v1-1]
        j = rand(1:m)
        b1 = circshift!(i, b, j)
        i1 = circshift!(b, j)
        @test b1 == i1
        b2 = circshift!(i1, b1, -j)
        i2 = circshift!(b1, -j)
        @test b2 == i2

        @check_bit_operation selectdim(b1, 1, m) SubArray{Bool, 0}
    end
    @check_bit_operation selectdim(b1, 1, :) SubArray{Bool, 1}
end

timesofar("datamove")

@testset "count & find" begin
    for m = 0:v1, b1 in Any[bitrand(m), trues(m), falses(m)]
        @check_bit_operation count(b1) Int

        @check_bit_operation findfirst(b1) Union{Int,Nothing}

        @check_bit_operation findfirst(!iszero, b1)    Union{Int,Nothing}
        @check_bit_operation findfirst(iszero, b1)     Union{Int,Nothing}
        @check_bit_operation findfirst(isequal(3), b1) Union{Int,Nothing}

        @check_bit_operation findfirst(x->x, b1)     Union{Int,Nothing}
        @check_bit_operation findfirst(x->!x, b1)    Union{Int,Nothing}
        @check_bit_operation findfirst(Returns(true ), b1)  Union{Int,Nothing}
        @check_bit_operation findfirst(Returns(false), b1) Union{Int,Nothing}

        @check_bit_operation findall(b1) Vector{Int}
    end

    b1 = trues(v1)
    for i = 0:(v1-1)
        @test findfirst(b1 >> i) == i+1
        @test Base.findfirstnot(.~(b1 >> i)) == i+1
    end

    for i = 3:(v1-1), j = 2:i
        submask = b1 << (v1-j+1)
        @test findnext((b1 >> i) .| submask, j) == i+1
        @test findnextnot((.~(b1 >> i)) .⊻ submask, j) == i+1
    end

    # Do a few more thorough tests for findall
    b1 = bitrand(n1, n2)
    @check_bit_operation findall(b1) Vector{CartesianIndex{2}}
    @check_bit_operation findall(!iszero, b1) Vector{CartesianIndex{2}}

    # tall-and-skinny (test index overflow logic in findall)
    @check_bit_operation findall(bitrand(1, 1, 1, 250)) Vector{CartesianIndex{4}}

    # empty dimensions
    @check_bit_operation findall(bitrand(0, 0, 10)) Vector{CartesianIndex{3}}

    # sparse (test empty 64-bit chunks in findall)
    b1 = falses(8, 8, 8)
    b1[3,3,3] = b1[6,6,6] = true
    @check_bit_operation findall(b1) Vector{CartesianIndex{3}}

    # BitArrays of various dimensions
    for dims = 0:8
        t = Tuple(fill(2, dims))
        ret_type = Vector{dims == 1 ? Int : CartesianIndex{dims}}
        @check_bit_operation findall(trues(t)) ret_type
        @check_bit_operation findall(falses(t)) ret_type
        @check_bit_operation findall(bitrand(t)) ret_type
    end

    @test count(trues(2, 2), init=0x03) === 0x07
    @test count(trues(2, 2, 2), dims=2) == fill(2, 2, 1, 2)
end

timesofar("find")

@testset "Findnext/findprev" begin
    b1 = trues(v1)
    b2 = falses(v1)
    for i = 1:v1
        @test findprev(b1, i) == findprev(isequal(true), b1, i) == findprev(identity, b1, i)
        @test findprevnot(b2, i) == findprev(!, b2, i) == i
    end

    odds = broadcast(isodd, 1:2000)
    evens = broadcast(iseven, 1:2000)

    for i = 1:2:2000
        @test findprev(odds,i)  == findprevnot(evens,i) == i
        @test findnext(odds,i)  == findnextnot(evens,i) == i
        @test findprev(evens,i) == findprevnot(odds,i)  == (i > 1    ? i-1 : nothing)
        @test findnext(evens,i) == findnextnot(odds,i)  == (i < 2000 ? i+1 : nothing)
    end
    for i = 2:2:2000
        @test findprev(odds,i)  == findprevnot(evens,i) == i-1
        @test findprev(evens,i) == findprevnot(odds,i)  == i
        @test findnext(evens,i) == findnextnot(odds,i)  == i
        @test findnext(odds,i)  == findnextnot(evens,i) == (i < 2000 ? i+1 : nothing)
    end

    elts = (1:64:(64*64+1)) .+ (0:64)
    n = maximum(elts)
    for c = [falses, trues]
        b1 = c(n)
        b1[elts] = .!b1[elts]
        b2 = .~b1
        i1 = Array(b1)
        for i = 1:n
            @test findprev(b1, i) == findprev(i1, i) == findprevnot(b2, i) == findprev(!, b2, i)
            @test findnext(b1, i) == findnext(i1, i) == findnextnot(b2, i) == findnext(!, b2, i)
        end
    end

    b1 = falses(1000)
    b1[77] = true
    b1[777] = true
    b2 = .~b1
    @test_throws BoundsError findprev(b1, 1001)
    @test_throws BoundsError findprevnot(b2, 1001)
    @test_throws BoundsError findprev(!, b2, 1001)
    @test_throws BoundsError findprev(identity, b1, 1001)
    @test_throws BoundsError findprev(Returns(false), b1, 1001)
    @test_throws BoundsError findprev(Returns(true ), b1, 1001)
    @test findprev(b1, 1000) == findprevnot(b2, 1000) == findprev(!, b2, 1000) == 777
    @test findprev(b1, 777)  == findprevnot(b2, 777)  == findprev(!, b2, 777)  == 777
    @test findprev(b1, 776)  == findprevnot(b2, 776)  == findprev(!, b2, 776)  == 77
    @test findprev(b1, 77)   == findprevnot(b2, 77)   == findprev(!, b2, 77)   == 77
    @test findprev(b1, 76)   == findprevnot(b2, 76)   == findprev(!, b2, 76)   == nothing
    @test findprev(b1, -1)   == findprevnot(b2, -1)   == findprev(!, b2, -1)   == nothing
    @test findprev(identity, b1, -1) == nothing
    @test findprev(Returns(false), b1, -1) == nothing
    @test findprev(Returns(true), b1, -1) == nothing
    @test_throws BoundsError findnext(b1, -1)
    @test_throws BoundsError findnextnot(b2, -1)
    @test_throws BoundsError findnext(!, b2, -1)
    @test_throws BoundsError findnext(identity, b1, -1)
    @test_throws BoundsError findnext(Returns(false), b1, -1)
    @test_throws BoundsError findnext(Returns(true), b1, -1)
    @test findnext(b1, 1)    == findnextnot(b2, 1)    == findnext(!, b2, 1)    == 77
    @test findnext(b1, 77)   == findnextnot(b2, 77)   == findnext(!, b2, 77)   == 77
    @test findnext(b1, 78)   == findnextnot(b2, 78)   == findnext(!, b2, 78)   == 777
    @test findnext(b1, 777)  == findnextnot(b2, 777)  == findnext(!, b2, 777)  == 777
    @test findnext(b1, 778)  == findnextnot(b2, 778)  == findnext(!, b2, 778)  == nothing
    @test findnext(b1, 1001) == findnextnot(b2, 1001) == findnext(!, b2, 1001) == nothing
    @test findnext(identity, b1, 1001) == findnext(Returns(false), b1, 1001) == findnext(Returns(true), b1, 1001) == nothing

    @test findlast(b1) == Base.findlastnot(b2) == 777
    @test findfirst(b1) == Base.findfirstnot(b2) == 77

    b0 = BitVector()
    @test findprev(Returns(true), b0, -1) == nothing
    @test_throws BoundsError findprev(Returns(true), b0, 1)
    @test_throws BoundsError findnext(Returns(true), b0, -1)
    @test findnext(Returns(true), b0, 1) == nothing

    b1 = falses(10)
    @test findprev(Returns(true), b1, 5) == 5
    @test findnext(Returns(true), b1, 5) == 5
    @test findprev(Returns(true), b1, -1) == nothing
    @test findnext(Returns(true), b1, 11) == nothing
    @test findprev(Returns(false), b1, 5) == nothing
    @test findnext(Returns(false), b1, 5) == nothing
    @test findprev(Returns(false), b1, -1) == nothing
    @test findnext(Returns(false), b1, 11) == nothing
    @test_throws BoundsError findprev(Returns(true), b1, 11)
    @test_throws BoundsError findnext(Returns(true), b1, -1)

    @testset "issue 32568" for T = (UInt, BigInt)
        for x = (1, 2)
            @test findnext(evens, T(x)) isa keytype(evens)
            @test findnext(iseven, evens, T(x)) isa keytype(evens)
            @test findnext(isequal(true), evens, T(x)) isa keytype(evens)
            @test findnext(isequal(false), evens, T(x)) isa keytype(evens)
        end
        for x = (3, 4)
            @test findprev(evens, T(x)) isa keytype(evens)
            @test findprev(iseven, evens, T(x)) isa keytype(evens)
            @test findprev(isequal(true), evens, T(x)) isa keytype(evens)
            @test findprev(isequal(false), evens, T(x)) isa keytype(evens)
        end
    end

    for l = [1, 63, 64, 65, 127, 128, 129]
        f = falses(l)
        t = trues(l)
        @test findprev(f, l) == findprevnot(t, l) == nothing
        @test findprev(t, l) == findprevnot(f, l) == l
        b1 = falses(l)
        b1[end] = true
        b2 = .~b1
        @test findprev(b1, l) == findprevnot(b2, l) == l
        @test findprevnot(b1, l) == findprev(b2, l) == (l == 1 ? nothing : l-1)
        if l > 1
            b1 = falses(l)
            b1[end-1] = true
            b2 = .~b1
            @test findprev(b1, l) == findprevnot(b2, l) == l-1
            @test findprevnot(b1, l) == findprev(b2, l) == l
        end
    end
end

@testset "Reductions" begin
    b1 = bitrand(s1, s2, s3, s4)
    m1 = 1
    m2 = 3
    @check_bit_operation maximum(b1, dims=(m1, m2)) BitArray{4}
    @check_bit_operation minimum(b1, dims=(m1, m2)) BitArray{4}
    @check_bit_operation sum(b1, dims=(m1, m2)) Array{Int,4}

    @check_bit_operation maximum(b1) Bool
    @check_bit_operation minimum(b1) Bool
    @check_bit_operation any(b1) Bool
    @check_bit_operation all(b1) Bool
    @check_bit_operation sum(b1) Int

    b0 = falses(0)
    @check_bit_operation any(b0) Bool
    @check_bit_operation all(b0) Bool
    @check_bit_operation sum(b0) Int
end

timesofar("reductions")

@testset "map over bitarrays" begin
    for l = [0, 1, 63, 64, 65, 127, 128, 129, 255, 256, 257, 6399, 6400, 6401]
        b1 = bitrand(l)
        b2 = bitrand(l)
        @test map(~, b1) == map(x->~x, b1) == broadcast(~, b1)
        @test map(identity, b1) == map(x->x, b1) == b1
        @test map(zero, b1) == map(Returns(false), b1) == falses(l)
        @test map(one, b1) == map(Returns(true), b1) == trues(l)

        @test map(&, b1, b2) == map((x,y)->x&y, b1, b2) == broadcast(&, b1, b2)
        @test map(|, b1, b2) == map((x,y)->x|y, b1, b2) == broadcast(|, b1, b2)
        @test map(⊻, b1, b2) == map((x,y)->x⊻y, b1, b2) == broadcast(⊻, b1, b2) == broadcast(xor, b1, b2)
        @test map(⊼, b1, b2) == map((x,y)->x⊼y, b1, b2) == broadcast(⊼, b1, b2) == broadcast(nand, b1, b2)
        @test map(⊽, b1, b2) == map((x,y)->x⊽y, b1, b2) == broadcast(⊽, b1, b2) == broadcast(nor, b1, b2)

        @test map(^, b1, b2) == map((x,y)->x^y, b1, b2) == b1 .^ b2
        @test map(*, b1, b2) == map((x,y)->x*y, b1, b2) == b1 .* b2

        @test map(min, b1, b2) == map((x,y)->min(x,y), b1, b2) == min.(b1, b2)
        @test map(max, b1, b2) == map((x,y)->max(x,y), b1, b2) == max.(b1, b2)

        @test map(<, b1, b2)  == map((x,y)->x<y, b1, b2)  == (b1 .< b2)
        @test map(<=, b1, b2) == map((x,y)->x<=y, b1, b2) == (b1 .<= b2)
        @test map(==, b1, b2) == map((x,y)->x==y, b1, b2) == (b1 .== b2)
        @test map(>=, b1, b2) == map((x,y)->x>=y, b1, b2) == (b1 .>= b2)
        @test map(>, b1, b2)  == map((x,y)->x>y, b1, b2)  == (b1 .> b2)
        @test map(!=, b1, b2) == map((x,y)->x!=y, b1, b2) == (b1 .!= b2)

        @testset "map! for length $l" begin
            b = BitVector(undef, l)
            @test map!(~, b, b1) == map!(x->~x, b, b1) == broadcast(~, b1) == b
            @test map!(!, b, b1) == map!(x->!x, b, b1) == broadcast(~, b1) == b
            @test map!(identity, b, b1) == map!(x->x, b, b1) == b1 == b
            @test map!(zero, b, b1) == map!(Returns(false), b, b1) == falses(l) == b
            @test map!(one, b, b1) == map!(Returns(true), b, b1) == trues(l) == b

            @test map!(&, b, b1, b2) == map!((x,y)->x&y, b, b1, b2) == broadcast(&, b1, b2) == b
            @test map!(|, b, b1, b2) == map!((x,y)->x|y, b, b1, b2) == broadcast(|, b1, b2) == b
            @test map!(⊻, b, b1, b2) == map!((x,y)->x⊻y, b, b1, b2) == broadcast(⊻, b1, b2) == broadcast(xor, b1, b2) == b

            @test map!(^, b, b1, b2) == map!((x,y)->x^y, b, b1, b2) == b1 .^ b2 == b
            @test map!(*, b, b1, b2) == map!((x,y)->x*y, b, b1, b2) == b1 .* b2 == b

            @test map!(min, b, b1, b2) == map!((x,y)->min(x,y), b, b1, b2) == min.(b1, b2) == b
            @test map!(max, b, b1, b2) == map!((x,y)->max(x,y), b, b1, b2) == max.(b1, b2) == b

            @test map!(<, b, b1, b2)  == map!((x,y)->x<y, b, b1, b2)  == (b1 .< b2)  == b
            @test map!(<=, b, b1, b2) == map!((x,y)->x<=y, b, b1, b2) == (b1 .<= b2) == b
            @test map!(==, b, b1, b2) == map!((x,y)->x==y, b, b1, b2) == (b1 .== b2) == b
            @test map!(>=, b, b1, b2) == map!((x,y)->x>=y, b, b1, b2) == (b1 .>= b2) == b
            @test map!(>, b, b1, b2)  == map!((x,y)->x>y, b, b1, b2)  == (b1 .> b2)  == b
            @test map!(!=, b, b1, b2) == map!((x,y)->x!=y, b, b1, b2) == (b1 .!= b2) == b
        end
    end

    @testset "Issue #17970" begin
        A17970 = [1,2,3] .== [3,2,1]
        B17970 = map(x -> x ? 1 : 2, A17970)
        @test B17970::Array{Int,1} == [2,1,2]
        C17970 = map(x -> x ? false : true, A17970)
        @test C17970::BitArray{1} == map(~, A17970)
    end
end

## Filter ##

# TODO

@testset "transpose" begin
    b1 = bitrand(v1)
    @check_bit_operation transpose(b1) Transpose{Bool,BitVector}

    for m1 = 0:n1, m2 = 0:n2
        b1 = bitrand(m1, m2)
        @check_bit_operation copy(b1') BitMatrix
    end
end

timesofar("transpose")

@testset "Permutedims" begin
    b1 = bitrand(s1, s2, s3, s4)
    p = randperm(4)
    @check_bit_operation permutedims(b1, p) BitArray{4}
    @check_bit_operation permutedims(b1, tuple(p...)) BitArray{4}
end

timesofar("permutedims")

@testset "Concatenation" begin
    b1 = bitrand(v1)
    b2 = bitrand(v1)
    @check_bit_operation hcat(b1, b2) BitMatrix
    for m = 1:(v1-1)
        @check_bit_operation vcat(b1[1:m], b1[m+1:end]) BitVector
    end
    @test_throws DimensionMismatch hcat(b1,trues(n1+1))
    @test_throws DimensionMismatch hcat(hcat(b1, b2),trues(n1+1))

    b1 = bitrand(n1, n2)
    b2 = bitrand(n1)
    b3 = bitrand(n1, n2)
    b4 = bitrand(1, n2)
    @check_bit_operation hcat(b1, b2, b3) BitMatrix
    @check_bit_operation vcat(b1, b4, b3) BitMatrix
    @test_throws DimensionMismatch vcat(b1, b4, trues(n1,n2+1))

    b1 = bitrand(s1, s2, s3, s4)
    b2 = bitrand(s1, s3, s3, s4)
    b3 = bitrand(s1, s2, s3, s1)
    @check_bit_operation cat(b1, b2, dims=2) BitArray{4}
    @check_bit_operation cat(b1, b3, dims=4) BitArray{4}
    @check_bit_operation cat(b1, b1, dims=6) BitArray{6}

    b1 = bitrand(1, v1, 1)
    @check_bit_operation cat(0, b1, 1, 1, b1, dims=2) Array{Int,3}
    @check_bit_operation cat(3, b1, 4, 5, b1, dims=2) Array{Int,3}
    @check_bit_operation cat(false, b1, true, true, b1, dims=2) BitArray{3}

    b1 = bitrand(n1, n2)
    for m1 = 1:(n1-1), m2 = 1:(n2-1)
        @test isequal([b1[1:m1,1:m2] b1[1:m1,m2+1:end]; b1[m1+1:end,1:m2] b1[m1+1:end,m2+1:end]], b1)
    end
end

timesofar("cat")

@testset "Linear algebra" begin
    b1 = bitrand(v1)
    b2 = bitrand(v1)
    @check_bit_operation dot(b1, b2) Int

    b1 = bitrand(n1, n2)
    @test_throws ArgumentError tril(b1, -n1 - 2)
    @test_throws ArgumentError tril(b1, n2)
    @test_throws ArgumentError triu(b1, -n1)
    @test_throws ArgumentError triu(b1, n2 + 2)
    for k in (-n1 - 1):(n2 - 1)
        @check_bit_operation tril(b1, k) BitMatrix
    end
    for k in (-n1 + 1):(n2 + 1)
        @check_bit_operation triu(b1, k) BitMatrix
    end

    for sz = [(n1,n1), (n1,n2), (n2,n1)], (f,isf) = [(tril,istril), (triu,istriu)]
        b1 = bitrand(sz...)
        @check_bit_operation isf(b1) Bool
        b1 = f(bitrand(sz...))
        @check_bit_operation isf(b1) Bool
    end

    b1 = bitrand(n1,n1)
    b1 .|= copy(b1')
    @check_bit_operation issymmetric(b1) Bool
    @check_bit_operation ishermitian(b1) Bool

    b1 = bitrand(n1)
    b2 = bitrand(n2)
    @check_bit_operation kron(b1, b2) BitVector

    b1 = bitrand(s1, s2)
    b2 = bitrand(s3, s4)
    @check_bit_operation kron(b1, b2) BitMatrix

    b1 = bitrand(v1)
    @check_bit_operation diff(b1) Vector{Int}

    b1 = bitrand(n1, n2)
    @check_bit_operation diff(b1, dims=1) Matrix{Int}
    @check_bit_operation diff(b1, dims=2) Matrix{Int}

    b1 = bitrand(n1, n1)
    @test ((svdb1, svdb1A) = (svd(b1), svd(Array(b1)));
            svdb1.U == svdb1A.U && svdb1.S == svdb1A.S && svdb1.V == svdb1A.V)
    @test ((qrb1, qrb1A) = (qr(b1), qr(Array(b1)));
            qrb1.Q == qrb1A.Q && qrb1.R == qrb1A.R)

    b1 = bitrand(v1)
    @check_bit_operation diagm(0 => b1) BitMatrix

    b1 = bitrand(v1)
    b2 = bitrand(v1)
    @check_bit_operation diagm(-1 => b1, 1 => b2) BitMatrix

    b1 = bitrand(n1, n1)
    @check_bit_operation diag(b1)
end

timesofar("linalg")

@testset "findmax, findmin" begin
    b1 = trues(0)
    @test_throws ArgumentError findmax(b1)
    @test_throws ArgumentError findmin(b1)

    for b1 in [falses(v1), trues(v1),
               BitArray([1,0,1,1,0]),
               BitArray([0,0,1,1,0]),
               BitArray([1 0; 1 1]),
               BitArray([0 0; 1 1]),
               bitrand(v1)]
        @check_bit_operation findmin(b1)
        @check_bit_operation findmax(b1)
    end
end

@testset "I/O" begin
    b1 = bitrand(v1)
    fname = ""
    b1[v1 ÷ 2 + 1] = true
    b1[end] = true
    try
        fname = tempname()
        open(fname, "w") do f
            write(f, b1)
        end
        b2 = falses(v1)
        read!(fname, b2)
        @test bitcheck(b2)
        @test b1 == b2
        b2 = falses(v1 ÷ 10, 10)
        read!(fname, b2)
        @test bitcheck(b2)
        @test reshape(b1, v1 ÷ 10, 10) == b2
        b2 = falses(v1 + 65)
        @test bitcheck(b2)
        @test_throws EOFError read!(fname, b2)
        @test bitcheck(b2)
        b2 = falses(v1 ÷ 2)
        @test_throws DimensionMismatch read!(fname, b2)
        @test bitcheck(b2)
        b2 = falses(v1 - 1)
        @test_throws DimensionMismatch read!(fname, b2)
        @test bitcheck(b2)

        b1 = BitVector()
        open(fname, "w") do f
            write(f, b1)
        end
        b2 = BitVector()
        read!(fname, b2)
        @test b1 == b2
        @test bitcheck(b2)
        b2 = trues(1)
        @test_throws EOFError read!(fname, b2)
        @test bitcheck(b2)
    finally
         isfile(fname) && rm(fname)
    end
end

timesofar("I/O")

@testset "not strided" begin
    @test_throws ErrorException pointer(trues(1))
    @test_throws ErrorException pointer(trues(1),1)
    b = falses(3)
    b[:] = view(trues(10), [1,3,7])
    @test b == trues(3)
end

@testset "chunked broadcast" begin
    for (f,g,h) in ((&,|,!),(*,xor,identity),(|,xor,sign),(&,&,~),(|,|,!))
        fg = (A, B, C)->f.(A, g.(B, C))
        fgh = (A, B, C)->f.(A, g.(B, h.(C)))
        for n in (1, 63, 64, 65, 127, 128, 129)
            for ((A,B,C),T) in ((bitrand.((n,n,n)), BitVector), (bitrand.((n,n,n), 2), BitMatrix))
                @check_bit_operation broadcast(f, A) T
                @check_bit_operation broadcast(g, A) T
                @check_bit_operation broadcast(h, A) T
                @check_bit_operation fg(A, B, C) T
                @check_bit_operation fg(true, B, C) T
                @check_bit_operation fg(A, false, C) T
                @check_bit_operation fg(A, B, true) T
                @check_bit_operation fgh(A, B, C) T
                @check_bit_operation fgh(true, B, C) T
                @check_bit_operation fgh(A, false, C) T
                @check_bit_operation fgh(A, B, true) T
            end
        end
    end
end

@testset "SIMD violations (issue #27482)" begin
    @test all(any!(falses(10), trues(10, 10)))
    @check_bit_operation any!(falses(10), trues(10, 10))
    @check_bit_operation any!(falses(100), trues(100, 100))
    @check_bit_operation any!(falses(1000), trues(1000, 100))
    @check_bit_operation all!(falses(10), trues(10, 10))
    @check_bit_operation all!(falses(100), trues(100, 100))
    @check_bit_operation all!(falses(1000), trues(1000, 100))
end

@testset "multidimensional concatenation returns BitArrays" begin
    a = BitVector(ones(5))
    @test typeof([a ;;; a]) <: BitArray
    @test typeof([a a ;;; a a]) <: BitArray
    @test typeof([a a ;;; [a a]]) <: BitArray
end

@testset "deleteat! additional tests" begin
    for v in ([1, 2, 3], [true, true, true], trues(3))
        @test_throws BoundsError deleteat!(v, true:true)
    end

    for v in ([1], [true], trues(1))
        @test length(deleteat!(v, false:false)) == 1
        @test isempty(deleteat!(v, true:true))
    end

    x = trues(3)
    x[3] = false
    @test deleteat!(x, [UInt8(2)]) == [true, false]
    @test_throws ArgumentError deleteat!(x, Any[true])
    @test_throws ArgumentError deleteat!(x, Any[1, true])
    @test_throws ArgumentError deleteat!(x, Any[2, 1])
    @test_throws BoundsError deleteat!(x, Any[4])
    @test_throws BoundsError deleteat!(x, Any[2, 4])

    function test_equivalence(n::Int)
        x1 = rand(Bool, n)
        x2 = BitVector(x1)
        inds1 = rand(Bool, n)
        inds2 = BitVector(inds1)
        return deleteat!(copy(x1), findall(inds1)) ==
               deleteat!(copy(x1), inds1) ==
               deleteat!(copy(x2), inds1) ==
               deleteat!(copy(x1), inds2) ==
               deleteat!(copy(x2), inds2)
    end

    Random.seed!(1234)
    for n in 1:20, _ in 1:100
        @test test_equivalence(n)
    end
end

max_stars_repo_name: digital-carver/julia

max_issues_repo_name: digital-carver/julia

max_forks_repo_name: digital-carver/julia


text:
# set voltage
# SV
export psp_setvoltage

"""
    psp_setvoltage(io, voltage)

Set output voltage.
"""
function psp_setvoltage(io_psp::IO,voltage::Real)
  @assert voltage<100 "voltage must be < 100V"
  vstring = @sprintf("%02.2f",voltage)
  if length(vstring) == 4
    set_voltage_string = "SV "*"0"*vstring*"\r"
  else
    set_voltage_string = "SV "*vstring*"\r"
  end
  write(io_psp, set_voltage_string)
  return nothing
end


text:
# ps3000setadvtriggerchannelconditions

export ps3000setadvtriggerchannelconditions, PS3000TriggerState, TriggerConditions

baremodule PS3000TriggerState
  const CONDITION_DONT_CARE = 0x0000
  const CONDITION_TRUE = 0x0001
  const CONDITION_FALSE = 0x0002
end

immutable TriggerConditions
    a :: Int16
    b :: Int16
    c :: Int16
    d :: Int16
    pulsewidthqualifier :: Int16
end

function ps3000setadvtriggerchannelconditions(handle     :: Int16,
                                              conditions :: TriggerConditions)
  nconditions = 1
  ps_status = ccall((:ps3000SetAdvTriggerChannelConditions, ps3000driver), Int16,
    (Int16, Ref{TriggerConditions}, Int16),
    handle, Ref{TriggerConditions}(conditions), nconditions)
  if ps_status == 0
    error("ps3000 Error: one of the parameters is out of range")
  end
  return nothing
end

function ps3000setadvtriggerchannelconditions(handle     :: Int16)
  # call with no conditions to turn trigger off
  nconditions = 0
  ps_status = ccall((:ps3000SetAdvTriggerChannelConditions, ps3000driver), Int16,
    (Int16, Ptr{Void}, Int16),
    handle, Ptr{Void}(C_NULL), nconditions)
  if ps_status == 0
    error("ps3000 Error: one of the parameters is out of range")
  end
  return nothing
end


text:
module Adaptation

import Base: string
import LinearAlgebra, Statistics
using ..AdvancedHMC: DEBUG

abstract type AbstractAdaptor end

include("stepsize.jl")
include("precond.jl")

abstract type AbstractCompositeAdaptor <: AbstractAdaptor end

# TODO: generalise this to a list of adaptors
struct NaiveCompAdaptor <: AbstractCompositeAdaptor
    pc  :: AbstractPreConditioner
    ssa :: StepSizeAdaptor
end

function adapt!(tp::NaiveCompAdaptor, θ::AbstractVector{<:Real}, α::AbstractFloat)
    adapt!(tp.ssa, θ, α)
    adapt!(tp.pc, θ, α)
end

function getM⁻¹(ca::AbstractCompositeAdaptor)
    return getM⁻¹(ca.pc)
end

function getϵ(ca::AbstractCompositeAdaptor)
    return getϵ(ca.ssa)
end

include("stan_adaption.jl")

export adapt!, getϵ, getM⁻¹,
       NesterovDualAveraging,
       UnitPreConditioner, DiagPreConditioner, DensePreConditioner,
       NaiveCompAdaptor, StanNUTSAdaptor

end # module


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

"""
    RCM(adjacency::Dict{Int, Vector{Int}}, degrees::Dict{Int, Int}, P::Int)

Calculate the Reverse Cuthill-McKee Algorithm for the adjacency graph.

# Examples

```jldoctest
julia> adjacency = Dict(
                   1 => [3, 8, 9],
                   2 => [6, 7, 3, 8],
                   3 => [1, 2],
                   4 => [12, 13],
                   5 => [13, 16],
                   6 => [2, 7],
                   7 => [2, 6, 10],
                   8 => [1, 2, 10, 15],
                   9 => [1, 15],
                  10 => [8, 7],
                  11 => [12, 13, 14],
                  12 => [4, 11],
                  13 => [5, 16, 4, 11],
                  14 => [11, 16],
                  15 => [9, 8],
                  16 => [13, 5, 14]);

julia> degrees = Dict(
                 1 => 3,
                 2 => 4,
                 3 => 2,
                 4 => 2,
                 5 => 2,
                 6 => 2,
                 7 => 3,
                 8 => 4,
                 9 => 2,
                10 => 2,
                11 => 3,
                12 => 2,
                13 => 4,
                14 => 2,
                15 => 2,
                16 => 3);

julia> P = 9;

julia> Q = 12;

julia> RCM(adjacency, degrees, P)
16-element Array{Int64,1}:
 16
 5
 14
 13
 11
 4
 12
 6
 7
 2
 10
 3
 8
 1
 15
 9

```

# References

* Wikipedia contributors. "Cuthill–McKee algorithm". Wikipedia, The Free Encyclopedia. Wikipedia, The Free Encyclopedia, 11 Jul. 2017. Web. 17 Jul. 2017. https://en.wikipedia.org/wiki/Cuthill%E2%80%93McKee_algorithm
"""
function RCM(adjacency::Dict{Int, Vector{Int}}, degrees::Dict{Int, Int}, P::Int, Q::Int)
    n = length(adjacency)
    R = Int[P]
    for i=1:n
        if i > length(R)
            push!(R, Q)
            A = sort(adjacency[Q], by=j->degrees[j])
            for b in A
                if !(b in R)
                    push!(R, b)
                end
            end
        else
            t = sort(adjacency[R[i]], by=j->degrees[j])
            for T in t
                if !(T in R)
                    push!(R, T)
                end
            end
        end
    end
    new_order = reverse(R)
    return new_order
end


text:
optimizer_glpk = () -> POI.Optimizer(GLPK.Optimizer())
optimizer_ipopt = () -> POI.Optimizer(Ipopt.Optimizer())
const OPTIMIZER_GLPK =
    MOI.Bridges.full_bridge_optimizer(optimizer_glpk(), Float64)
const OPTIMIZER_CACHED_IPOPT = MOI.Utilities.CachingOptimizer(
    MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
    MOI.Bridges.full_bridge_optimizer(optimizer_ipopt(), Float64),
)
const CONFIG = MOIT.Config()

@testset "GLPK tests" begin
    # TODO see why tests error or fail
    MOIT.runtests(
        OPTIMIZER_GLPK,
        CONFIG;
        exclude = [
            # GLPK returns INVALID_MODEL instead of INFEASIBLE
            "test_constraint_ZeroOne_bounds_3",
            # Upstream issue: https://github.com/jump-dev/MathOptInterface.jl/issues/1431
            "test_model_LowerBoundAlreadySet",
            "test_model_UpperBoundAlreadySet",
        ],
    )
end

@testset "Ipopt tests" begin
    MOI.set(OPTIMIZER_CACHED_IPOPT, MOI.Silent(), true)
    # Without fixed_variable_treatment set, duals are not computed for variables
    # that have lower_bound == upper_bound.
    MOI.set(
        OPTIMIZER_CACHED_IPOPT,
        MOI.RawOptimizerAttribute("fixed_variable_treatment"),
        "make_constraint",
    )
    MOIT.runtests(
        OPTIMIZER_CACHED_IPOPT,
        MOIT.Config(
            atol = 1e-4,
            rtol = 1e-4,
            optimal_status = MOI.LOCALLY_SOLVED,
            exclude = Any[
                MOI.ConstraintBasisStatus,
                MOI.DualObjectiveValue,
                MOI.ObjectiveBound,
            ],
        );
        exclude = String[
            # Tests purposefully excluded:
            #  - Upstream: ZeroBridge does not support ConstraintDual
            "test_conic_linear_VectorOfVariables_2",
            #  - Excluded because this test is optional
            "test_model_ScalarFunctionConstantNotZero",
            #  - Excluded because Ipopt returns NORM_LIMIT instead of
            #    DUAL_INFEASIBLE
            "test_solve_TerminationStatus_DUAL_INFEASIBLE",
            #  - Excluded because Ipopt returns INVALID_MODEL instead of
            #    LOCALLY_SOLVED
            "test_linear_VectorAffineFunction_empty_row",
            #  - Excluded because Ipopt returns LOCALLY_INFEASIBLE instead of
            #    INFEASIBLE
            "INFEASIBLE",
            "test_conic_linear_INFEASIBLE",
            "test_conic_linear_INFEASIBLE_2",
            "test_solve_DualStatus_INFEASIBILITY_CERTIFICATE_",
            #  - Excluded due to upstream issue
            "test_model_LowerBoundAlreadySet",
            "test_model_UpperBoundAlreadySet",
            #  - CachingOptimizer does not throw if optimizer not attached
            "test_model_copy_to_UnsupportedAttribute",
            "test_model_copy_to_UnsupportedConstraint",
        ],
    )
end


text:
module HeisenbergImproper
__precompile__(false)

export He, heis, hn, left, right

using Nemo
using Luxor
import Luxor.transform
using Colors
using StructEquality

QQ = Rational{Int}

using ...CoreDefs

@def_structequal struct He <: DGroupElem
    a::QQ
    b::QQ
    c::QQ
end
@enum HePTile begin
    Left=1
    Right=2
end

function Base.:*(x :: He, y :: He) where {T}
    return He( x.a + y.a, x.b + y.b + x.a*y.c, x.c + y.c)
end

function CoreDefs.dilate(λ, x :: He)
    return He(x.a, λ * x.b, λ * x.c )
end

function Base.inv(x :: He)
    return He(-x.a, -x.b, x.a*x.b - x.c)
end

function left(a,b,c)
    return (He(a,b,c), Left)
end
function left()
    return left(0,0,0)
end

function right(a,b,c)
    return (He(a,b,c), Right)
end
function right()
    return right(0,0,0)
end

function heis()
    heisenberg_subst = Dict([
        (Left, [
            left(0,-2,-1),
            right(0,2,-1),
            left(-1,1,-1),
            left(0,2,3),
        ]),
        (Right, [
            right(0,-2,3),
            left(0,2,3),
            right(1,-1,3),
            right(0,2,-1),
        ])
    ])

    return SubSystem(heisenberg_subst, 2//1)
end

function CoreDefs.embed_aff(x :: He)
    return float([1, 0, x.a, 1, x.b, x.c])
end

function CoreDefs.draw(p::HePTile, action)
    if p == Left
        Luxor.poly([
            Point(-2,-1),
            Point(2,3),
            Point(2,-1)
        ], action, close=true)
    else
        Luxor.poly([
            Point(-2,3),
            Point(2,3),
            Point(2,-1)
        ], action, close=true)
    end

end

function hn(n)
    substitute(heis(), [left()], n)
end

function color(tile)
    return Int(tile[2])
end

end


text:
################################# sample #######################################
# IDEA: pre-allocation?
# IDEA: lazy GP sampling?
using LinearAlgebraExtensions: vecofvec # converts a matrix into a vector of vectors
sample(G::Univariate{T}) where T = mean(G) + std(G) * randn(T)

function sample(G::Univariate{T}, k::Int) where T
    x = randn(T, k)
    x .= mean(G) .+ std(G) .* x
end

# sample(G::AbstractGaussian{<:Real}) = sample(G, 1)[]
function sample(G::Multivariate; tol::Real = tol)
    vec(sample(G, 1, tol = tol)[1])
end

# TODO: sample(G, k) should return a vector with k elements, each of which is a sample
# multi-variate Gaussian sample, falls back to cholesky
function sample(G::Multivariate{<:AbstractVector, <:AbstractMatrix}, k::Int; tol::Real = tol)
    C = cholesky(G, tol = tol)
    return sample(C, k, tol = tol)
end
# requires that Σ.L*Σ.L' = Σ of the original multivariate Gaussian
function sample(G::Gaussian{<:AbstractVector, <:Cholesky}, k::Int; tol::Real = tol)
    vecofvec(mean(G) .+ cov(G).L * randn(length(G), k))
end
function sample(G::Gaussian{<:AbstractVector, <:CholeskyPivoted}, k::Int; tol::Real = tol)
    C = cov(G)
    C.tol ≤ tol || throw("CholeskyPivoted does not have sufficient error tolerance: $(C.tol) ≰ $tol")
    L = LowRank(C)
    sample(Gaussian(mean(G), L), k, tol = tol)
end
using LinearAlgebraExtensions: LowRank
function sample(G::Gaussian{<:AbstractVector, <:LowRank}, k::Int; tol::Real = tol)
    L = cov(G)
    L.tol ≤ tol || throw("LowRank does not have sufficient error tolerance: $(L.tol) ≰ $tol")
    L.U ≡ L.V' || throw("sampling with LowRank only supported if U ≡ V'")
    vecofvec(mean(G) .+ L.U * randn(rank(L), k))
end

function sample(G::Gaussian{<:AbstractVector, <:Union{Diagonal, UniformScaling}},
                k::Int; tol::Real = tol)
    vecofvec(mean(G) .+ sqrt(cov(G)) * randn((length(G), k)))
end

function sample(G::Gaussian{<:AbstractVector, <:Woodbury}, k::Int; tol::Real = tol)
    W = cov(G)
    W.U ≡ W.V' || throw("sampling with Woodbury only supported if U ≡ V'")
    s(A) = sample(Gaussian(A), k, tol = tol)
    s(W.A) .+ .*((W.U,), s(W.C)) .+ (mean(G),)
end


text:
using Base64: Base64EncodePipe

"""
    isjsonmime(m::MIME)

Determine whether or not a MIME type is JSON data that should be included
literally in a MIMEBundle (otherwise they will be included as a JSON-encoded
string).
"""
isjsonmime(m::MIME) = istextmime(m) && endswith(string(m), "+json")
isjsonmime(m::AbstractString) = isjsonmime(MIME(m))

# """
#     maybejson(m::MIME, s::AbstractString)
#
# Wrap `s` in a `JSONText` if `m` is a JSON MIME type. This prevents the JSON from
# being doubly-encoded when `MIMEBundles` are JSON encoded.
# """
# function maybejson(m::MIME, s::AbstractString)
#     if isjsonmime(m)
#         return JSONText(s)
#     end
#     return String(s)
# end

"""
    israwtext(m::MIME, obj)

Determine whether or not `obj` should be displayed literally for the specified
MIME type. This is useful to display data that's already encoded properly with
a string MIME type and is used by `mimestring`.
"""
israwtext(m::MIME, x::AbstractString) = istextmime(m)
israwtext(::MIME"text/plain", x::AbstractString) = false
israwtext(::MIME, x) = false

"""
    mimestring(mime::MIME, obj[; limit=true, color=true])

Generate a string using the specified MIME type in a IO context with `:limit`
set to `true`. Binary data is base64-encoded.
"""
function mimestring(mime::MIME, obj; limit=true, color=true)::String
    # If the MIME type is textual and we're given some kind of string, we
    # assume that the data is already represented as that type and there's
    # no need to do anything else to it.
    # e.g. `mimestring("text/html", "<html>...</html>")`
    if israwtext(mime, obj)
        return obj
    end

    # Otherwise, we need to "coerce" obj to a correct MIME representation via
    # `Base.show`.
    buf = IOBuffer()
    if istextmime(mime)
        # Otherwise, we're trying to show some non-string object (or the MIME
        # type should never be considered raw, like text/plain) using the
        # specified MIME, so we defer to `Base.show`.
        # e.g. `mimestring("text/markdown", md"**hello!**")`
        show(IOContext(buf, :limit=>limit, :color=>color), mime, obj)
    else
        b64 = Base64EncodePipe(buf)

        # We were called with a binary MIME and a bunch of binary data (i.e. a
        # vector of bytes), so we should just pass it directly through.
        # e.g. `mimestring("image/png", read("/path/to/catparty.gif")`
        if isa(obj, Vector{UInt8})
            write(b64, obj)
        # We were called with some non-byte-vector object so we need to "coerce"
        # it to the correct MIME type representation.
        else
            show(IOContext(b64, :limit => true, :color => true), mime, obj)
        end

        close(b64)
    end

    return String(take!(buf))
end

function mimerepr(mime::MIME, obj)
    if !showable(mime, obj)
        return nothing, nothing
    end
    try
        return mime, mimestring(mime, obj)
    catch exc
        @warn "Exception while trying to generate MIME repr:" mime obj
        return nothing, nothing
    end
end

function mimerepr(mimes::AbstractVector{MIME}, obj)
    for mime in mimes
        result = mimerepr(mime, obj)
        if result[1] !== nothing
            return result
        end
    end

    # No MIMEs returned anything.
    return nothing, nothing
end


text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

#
# package typedb.client.api;
#
# import typedb.client.api.concept.ConceptManager;
# import typedb.client.api.logic.LogicManager;
# import typedb.client.api.query.QueryFuture;
# import typedb.client.api.query.QueryManager;
# import typedb.protocol.TransactionProto;
#
# import javax.annotation.CheckReturnValue;
# import java.util.stream.Stream;
#
# public interface TypeDBTransaction extends AutoCloseable {
#
#     @CheckReturnValue
#     boolean isOpen();
#
#     @CheckReturnValue
#     Type type();
#
#     @CheckReturnValue
#     TypeDBOptions options();
#
#     @CheckReturnValue
#     ConceptManager concepts();
#
#     @CheckReturnValue
#     LogicManager logic();
#
#     @CheckReturnValue
#     QueryManager query();
#
#     void commit();
#
#     void rollback();
#
#     void close();
#
#     enum Type {
#         READ(0),
#         WRITE(1);
#
#         private final int id;
#         private final boolean isWrite;
#
#         Type(int id) {
#             this.id = id;
#             this.isWrite = id == 1;
#         }
#
#         public static Type of(int value) {
#             for (Type t : values()) {
#                 if (t.id == value) return t;
#             }
#             return null;
#         }
#
#         public int id() {
#             return id;
#         }
#
#         public boolean isRead() { return !isWrite; }
#
#         public boolean isWrite() { return isWrite; }
#
#         public TransactionProto.Transaction.Type proto() {
#             return TransactionProto.Transaction.Type.forNumber(id);
#         }
#     }
#
#     interface Extended extends TypeDBTransaction {
#
#         TransactionProto.Transaction.Res execute(TransactionProto.Transaction.Req.Builder request);
#
#         QueryFuture<TransactionProto.Transaction.Res> query(TransactionProto.Transaction.Req.Builder request);
#
#         Stream<TransactionProto.Transaction.ResPart> stream(TransactionProto.Transaction.Req.Builder request);
#     }
# }


text:
module ColorFunctions
	using Colors
	import ..Lattices: intersect_lattice_with_plane, dimension
	using LightGraphs

	is_leaf(g::SimpleDiGraph, s) = !mapreduce(|, vertices(g)) do in; has_edge(g, in, s) end

	mutable struct ColorMapping
		inner::Dict{Any, Color}
		full::Dict{Any, Color}
	end
	ColorMapping(seed::Dict{<:Any, <:Color}) = ColorMapping(seed, copy(seed))
	min_colors(palette) = ColorMapping(Dict([0 =>colorant"white", 1 => palette[1]]))

	default_palette = distinguishable_colors(256, colorant"midnightblue")

	"""
		lighten(C, p)

		Lighten a color `C` by `p` percent. Internally converts to HSV colorspace,
		reduces saturation, increases value (brightness) and converts back to
		the orginal colorspace.
	"""
	function lighten(C,p)
	    T = typeof(C)
	    hsv = HSV(C)
	    convert(T, HSV(hsv.h, min(hsv.s * (1-0*p), 1.0), min(hsv.v * (1+p), 1.0) ))
	end

	function color_phylo!(C::ColorMapping, state; depth=1,
		 palette=default_palette, transition=(x::Color, d)->x)

		P = state.Phylogeny

		## Find all vertices which are at most `depth` removed from the root.
		## Each of which is going to be assigned its own color.
		## They form the 'inner' part.
		if depth == -1
			depth = nv(P)
		end
		ring = neighborhood(P, 1, depth, dir=:in)
		# How many colors are currently in use?
		j = length(unique(values(C.inner))) - 1

		### Prune inner colors ###
		# Assign genotypes to indices.
		inner_gs = map(x->state.meta.genotypes[x], ring)
		# If a genotype is no longer present, remove it from the color-map.
		for k in keys(C.inner)
			if k!=0 && !(k in inner_gs)
				delete!(C.inner, k)
				j -= 1
				# @info "$k pruned"
			end
		end
		# Assign colors to new genotypes within the inner ring.
		for ringv in ring
			g = state.meta.genotypes[ringv]
			if !haskey(C.inner ,g)
				push!(C.inner, g => palette[j%length(palette)+1])
			end
			j += 1
		end

		### Build a new _full_ color-map from the inner colors.
		C.full = copy(C.inner)
		# @info C.inner
		for ringv in ring
			nb = neighborhood(P, ringv, nv(P), dir=:in)
			g_parent = state.meta.genotypes[ringv]
			for (d,n) in enumerate(nb)
				g = state.meta.genotypes[n]
				newcolor = range(C.inner[g_parent], colorant"white", length=round(Int, 1.5*length(nb)))[d]
				#oc = RGBA(C.inner[g_parent])
				#newcolor = RGBA(oc.r, oc.g, oc.b, clamp( oc.alpha*(1-0.3), 0., 1.))
				push!(C.full, g => newcolor)
			end
		end

		map(reshape(state.lattice.data, length(state.lattice.data))) do s
			C.full[s]
		end
	end
	color_phylo(state; depth=1, palette=default_palette, kwargs...) = color_phylo!(min_colors(palette), state; depth=depth, kwargs...)

	function color_phylo_by_genomic_distance!(C::ColorMapping, state; depth=1,
		 palette=reverse(colormap("Blues", 20))
		)
		P = state.Phylogeny

		C.full = copy(C.inner)

		for j in 2:length(state.meta.genotypes)
			l = length(state.meta.snps[j])
			# push!(C.full, state.meta.genotypes[j] => palette[l%length(palette)+1])
			push!(C.full, state.meta.genotypes[j] => palette[min(l, 20)] )
		end

		map(reshape(state.lattice.data, state.lattice.Na^2)) do s
			C.full[s]
		end
	end
	function color_generic(state;palette=rand(RGB,maximum(state.lattice.data)))
		v0 = state.lattice.data
		map(v0) do x
			if x==0
				RGBA{Float32}(1.0,1.0,1.0,0.0)
			else
				RGBA{Float32}(palette[x])
			end
		end |> X->reshape(X, length(v0))
	end

	function color_genotype(state, g; color=colorant"blue")
		v0 = state.lattice.data
		color = RGBA(color)
		map(v0) do x
			if x!=g
				RGBA{Float32}(1.0,1.0,1.0,0.0)
			else
				color
			end
		end |> X->reshape(X, length(v0))
	end

	function color_neighbors(state, I::CartesianIndex)
		nn = neighbors(state.lattice, I)
		C = similar(state.lattice.data, RGBA{Float32})
		fill!(C, colorant"gray")

		C[nn] .= colorant"cornflowerblue"
		return reshape(C, length(C))
	end

	function color_column(state, q)
		C = similar(state.lattice.data, RGBA{Float32})
		fill!(C, colorant"gray")

		C[:,q] .= colorant"cornflowerblue"
		return reshape(C, length(C))
	end

	function color_row(state, r)
		C = similar(state.lattice.data, RGBA{Float32})
		fill!(C, colorant"gray")

		C[r,:] .= colorant"cornflowerblue"
		return reshape(C, length(C))
	end

	function color_rand(state)
		C = similar(state.lattice.data, RGBA{Float32})
		reshape(rand(RGB{Float32}, size(state.lattice.data)), length(state.lattice.data))
	end

	function hide!(v, A::BitArray)
		for i in eachindex(v)
			if !A[i]
				v[i] = RGBA{Float32}(0., 0., 0., 0.)
			end
		end
		v
	end

#= End module =#
end


text:
module Tokenization

include("bioes.jl")
include("nn.jl")
include("decoder.jl")

end


text:
# The Computer Language Benchmarks Game
# https://salsa.debian.org/benchmarksgame-team/benchmarksgame/
#
# contributed by Daniel Jones
# fixed by David Campbell
# modified by Jarrett Revels, Alex Arslan, Yichao Yu

using Printf

const variants = [
      "agggtaaa|tttaccct",
      "[cgt]gggtaaa|tttaccc[acg]",
      "a[act]ggtaaa|tttacc[agt]t",
      "ag[act]gtaaa|tttac[agt]ct",
      "agg[act]taaa|ttta[agt]cct",
      "aggg[acg]aaa|ttt[cgt]ccct",
      "agggt[cgt]aa|tt[acg]accct",
      "agggta[cgt]a|t[acg]taccct",
      "agggtaa[cgt]|[acg]ttaccct"
]

const subs = [
    (r"tHa[Nt]", "<4>"),
    (r"aND|caN|Ha[DS]|WaS", "<3>"),
    (r"a[NSt]|BY", "<2>"),
    (r"<[^>]*>", "|"),
    (r"\|[^|][^|]*\|", "-")
]

function perf_regex_dna()
    seq = read(stdin, String)
    l1 = length(seq)

    seq = replace(seq, r">.*\n|\n" => "")
    l2 = length(seq)

    for v in variants
        k = 0
        for m in eachmatch(Regex(v), seq)
            k += 1
        end
        @printf("%s %d\n", v, k)
    end

    for (u, v) in subs
        seq = replace(seq, u => v)
    end

    println()
    println(l1 + 1)  # why + 1??
    println(l2)
    println(length(seq))
end

perf_regex_dna()


text:
using DimensionalData
using DimensionalData: @dim, hasdim, Dimension, IndependentDim
using DimensionalData: basetypeof
using Dates

Time = DimensionalData.Ti

AbDimArray = DimensionalData.AbstractDimArray

export At, Between, Near # Selectors from DimensionalArrays.jl
export hasdim, DimensionalArray
export get_var_as_dimarray
export Time, Lon, Lat, dims, Coord, Hei, Ti
export EqArea, Grid, spacestructure

@dim Lon IndependentDim "Longitude" "lon"
@dim Lat IndependentDim "Latitude" "lat"
@dim Coord IndependentDim "Coordinates (spatial)"
@dim Hei IndependentDim "Height" "height"
@dim Pre IndependentDim "Pressure" "pressure"

STANDARD_DIMS = (Lon, Lat, Time, Hei, Pre, Coord)

const COMMONNAMES = Dict(
    "lat" => Lat,
    "latitude" => Lat,
    "lon" => Lon,
    "long" => Lon,
    "longitude" => Lon,
    "time" => Time,
    "height" => Hei,
    "altitude" => Hei,
    "pressure" => Pre,
)

# the trait EqArea is for equal area grids. Functions can use the `spacestructure` and
# dispatch on `EqArea` or other types while still being type-stable
struct EqArea end
struct Grid end
spacestructure(a::AbDimArray) = spacestructure(dims(a))
function spacestructure(dims)
    if hasdim(dims, Coord)
        EqArea()
    elseif hasdim(dims, Lon) || hasdim(dims, Lat)
        Grid()
    else
        nothing
    end
end

export ClimArray
struct ClimArray{T,N,D<:Tuple,R<:Tuple,A<:AbstractArray{T,N},Na<:AbstractString,Me} <: AbstractDimensionalArray{T,N,D,A}
    data::A
    dims::D
    refdims::R
    name::Na
    attrib::Me
end
ClimArray(A::DimensionalArray) = ClimArray(A.data, A.dims, A.refdims, A.name, A.metadata)

"""
    ClimArray(A::Array, dims::Tuple; name = "", attrib = nothing)

`ClimArray` is a structure that contains numerical array data bundled with dimensional
information, a name and an `attrib` field (typically a dictionary) that holds general
attributes.
You can think of `ClimArray` as a in-memory representation of a CFVariable.

At the moment, a `ClimArray` is using `DimensionalArray` from DimensionalData.jl, and
all basic handling of `ClimArray` is offered by `DimensionalData` (see below).

`ClimArray` is created by passing in standard array data `A` and a tuple of dimensions `dims`.

## Example
```julia
using ClimateBase, Dates
Time = ClimateBase.Ti # more intuitive name for time dimension
lats = -90:5:90
lons = 0:10:359
t = Date(2000, 3, 15):Month(1):Date(2020, 3, 15)
# dimensional information:
dimensions = (Lon(lons), Lat(lats), Time(t))
data = rand(36, 37, 241) # numeric data
A = ClimArray(data, dimensions)
```
"""
ClimArray(A::AbstractArray, dims::Tuple; refdims=(), name="", attrib=nothing) =
ClimArray(A, DimensionalData.formatdims(A, dims), refdims, name, attrib)
ClimArray(A::AbstractArray, dims::Tuple, name::String; refdims=(), attrib=nothing) =
ClimArray(A, DimensionalData.formatdims(A, dims), refdims, name, attrib)

Base.parent(A::ClimArray) = A.data
Base.@propagate_inbounds Base.setindex!(A::ClimArray, x, I::Vararg{DimensionalData.StandardIndices}) =
    setindex!(A.data, x, I...)

DimensionalData.metadata(A::ClimArray) = A.attrib
DimensionalData.rebuild(A::ClimArray, data::Any, dims::Tuple=dims(A), refdims=DimensionalData.refdims(A),
name="", attrib=nothing) = ClimArray(data, dims, refdims, name, attrib)
DimensionalData.basetypeof(::ClimArray) = ClimArray

# The following basic methods allow indexing with tuples, (Time(5), Lon(3))
Base.getindex(A::ClimArray, i::Tuple) = A[i...]
Base.setindex!(A::ClimArray, x, i::Tuple) = setindex!(A, x, i...)
Base.view(A::ClimArray, i::Tuple) = view(A, i...)

# Remove reference dims from printing, and show attributes if any
function Base.show(io::IO, A::ClimArray)
    summary(io, A)
    print(io, "and")
    printstyled(io, " data: "; color=:green)
    dataA = data(A)
    print(io, summary(dataA), "\n")
    DimensionalData.custom_show(io, data(A))
end

# Define summary
function Base.summary(io::IO, A::ClimArray)
    l = nameof(typeof(A))
    printstyled(io, nameof(typeof(A)); color=:blue)
    if A.name != ""
        print(io, " (named ")
        printstyled(io, A.name; color=:blue)
        print(io, ")")
    end

    print(io, " with dimensions:\n")
    for d in dims(A)
        print(io, " ", d, "\n")
    end
    if !isnothing(A.attrib)
        printstyled(io, "attributes: "; color=:magenta)
        show(io, MIME"text/plain"(), A.attrib)
        print(io, '\n')
    end
end

Base.ones(A::AbDimArray) = basetypeof(A)(ones(size(A)), dims(A))