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

## BitArray

# notes: bits are stored in contiguous chunks
#        unused bits must always be set to 0
"""
    BitArray{N} <: AbstractArray{Bool, N}

Space-efficient `N`-dimensional boolean array, using just one bit for each boolean value.

`BitArray`s pack up to 64 values into every 8 bytes, resulting in an 8x space efficiency
over `Array{Bool, N}` and allowing some operations to work on 64 values at once.

By default, Julia returns `BitArrays` from [broadcasting](@ref Broadcasting) operations
that generate boolean elements (including dotted-comparisons like `.==`) as well as from
the functions [`trues`](@ref) and [`falses`](@ref).

!!! note
    Due to its packed storage format, concurrent access to the elements of a `BitArray`
    where at least one of them is a write is not thread safe.

"""
mutable struct BitArray{N} <: AbstractArray{Bool, N}
    chunks::Vector{UInt64}
    len::Int
    dims::NTuple{N,Int}
    function BitArray{N}(::UndefInitializer, dims::Vararg{Int,N}) where N
        n = 1
        i = 1
        for d in dims
            d >= 0 || throw(ArgumentError("dimension size must be ≥ 0, got $d for dimension $i"))
            n *= d
            i += 1
        end
        nc = num_bit_chunks(n)
        chunks = Vector{UInt64}(undef, nc)
        nc > 0 && (chunks[end] = UInt64(0))
        b = new(chunks, n)
        N != 1 && (b.dims = dims)
        return b
    end
end

# note: the docs for the two signatures are unified, but only
# the first one is recognized by the help system; it would be nice
# to fix this.
"""
    BitArray(undef, dims::Integer...)
    BitArray{N}(undef, dims::NTuple{N,Int})

Construct an undef [`BitArray`](@ref) with the given dimensions.
Behaves identically to the [`Array`](@ref) constructor. See [`undef`](@ref).

# Examples
```julia-repl
julia> BitArray(undef, 2, 2)
2×2 BitMatrix:
 0  0
 0  0

julia> BitArray(undef, (3, 1))
3×1 BitMatrix:
 0
 0
 0
```
"""
BitArray(::UndefInitializer, dims::Integer...) = BitArray(undef, map(Int,dims))
BitArray{N}(::UndefInitializer, dims::Integer...) where {N} = BitArray{N}(undef, map(Int,dims))::BitArray{N}
BitArray(::UndefInitializer, dims::NTuple{N,Integer}) where {N} = BitArray{N}(undef, map(Int, dims)...)
BitArray{N}(::UndefInitializer, dims::NTuple{N,Integer}) where {N} = BitArray{N}(undef, map(Int, dims)...)

const BitVector = BitArray{1}
const BitMatrix = BitArray{2}

BitVector() = BitVector(undef, 0)

"""
    BitVector(nt::Tuple{Vararg{Bool}})

Construct a `BitVector` from a tuple of `Bool`.
# Examples
```julia-repl
julia> nt = (true, false, true, false)
(true, false, true, false)

julia> BitVector(nt)
4-element BitVector:
 1
 0
 1
 0
```
"""
function BitVector(nt::Tuple{Vararg{Bool}})
    bv = BitVector(undef, length(nt))
    bv .= nt
end

## utility functions ##

length(B::BitArray) = B.len
size(B::BitVector) = (B.len,)
size(B::BitArray) = B.dims

@inline function size(B::BitVector, d::Integer)
    d < 1 && throw_boundserror(size(B), d)
    ifelse(d == 1, B.len, 1)
end

isassigned(B::BitArray, i::Int) = 1 <= i <= length(B)

IndexStyle(::Type{<:BitArray}) = IndexLinear()

## aux functions ##

const _msk64 = ~UInt64(0)
@inline _div64(l) = l >> 6
@inline _mod64(l) = l & 63
@inline _blsr(x)= x & (x-1) #zeros the last set bit. Has native instruction on many archs. needed in multidimensional.jl
@inline _msk_end(l::Int) = _msk64 >>> _mod64(-l)
@inline _msk_end(B::BitArray) = _msk_end(length(B))
num_bit_chunks(n::Int) = _div64(n+63)

@inline get_chunks_id(i::Int) = _div64(i-1)+1, _mod64(i-1)

function glue_src_bitchunks(src::Vector{UInt64}, k::Int, ks1::Int, msk_s0::UInt64, ls0::Int)
    @inbounds begin
        chunk = ((src[k] & msk_s0) >>> ls0)
        if ks1 > k && ls0 > 0
            chunk_n = (src[k + 1] & ~msk_s0)
            chunk |= (chunk_n << (64 - ls0))
        end
    end
    return chunk
end

function copy_chunks!(dest::Vector{UInt64}, pos_d::Int, src::Vector{UInt64}, pos_s::Int, numbits::Int)
    numbits == 0 && return
    if dest === src && pos_d > pos_s
        return copy_chunks_rtol!(dest, pos_d, pos_s, numbits)
    end

    kd0, ld0 = get_chunks_id(pos_d)
    kd1, ld1 = get_chunks_id(pos_d + numbits - 1)
    ks0, ls0 = get_chunks_id(pos_s)
    ks1, ls1 = get_chunks_id(pos_s + numbits - 1)

    delta_kd = kd1 - kd0
    delta_ks = ks1 - ks0

    u = _msk64
    if delta_kd == 0
        msk_d0 = ~(u << ld0) | (u << (ld1+1))
    else
        msk_d0 = ~(u << ld0)
        msk_d1 = (u << (ld1+1))
    end
    if delta_ks == 0
        msk_s0 = (u << ls0) & ~(u << (ls1+1))
    else
        msk_s0 = (u << ls0)
    end

    chunk_s0 = glue_src_bitchunks(src, ks0, ks1, msk_s0, ls0)

    dest[kd0] = (dest[kd0] & msk_d0) | ((chunk_s0 << ld0) & ~msk_d0)

    delta_kd == 0 && return

    for i = 1 : kd1 - kd0 - 1
        chunk_s1 = glue_src_bitchunks(src, ks0 + i, ks1, msk_s0, ls0)

        chunk_s = (chunk_s0 >>> (64 - ld0)) | (chunk_s1 << ld0)

        dest[kd0 + i] = chunk_s

        chunk_s0 = chunk_s1
    end

    if ks1 >= ks0 + delta_kd
        chunk_s1 = glue_src_bitchunks(src, ks0 + delta_kd, ks1, msk_s0, ls0)
    else
        chunk_s1 = UInt64(0)
    end

    chunk_s = (chunk_s0 >>> (64 - ld0)) | (chunk_s1 << ld0)

    dest[kd1] = (dest[kd1] & msk_d1) | (chunk_s & ~msk_d1)

    return
end

function copy_chunks_rtol!(chunks::Vector{UInt64}, pos_d::Int, pos_s::Int, numbits::Int)
    pos_d == pos_s && return
    pos_d < pos_s && return copy_chunks!(chunks, pos_d, chunks, pos_s, numbits)

    left = numbits
    s = min(left, 64)
    b = left - s
    ps = pos_s + b
    pd = pos_d + b
    u = _msk64
    while left > 0
        kd0, ld0 = get_chunks_id(pd)
        kd1, ld1 = get_chunks_id(pd + s - 1)
        ks0, ls0 = get_chunks_id(ps)
        ks1, ls1 = get_chunks_id(ps + s - 1)

        delta_kd = kd1 - kd0
        delta_ks = ks1 - ks0

        if delta_kd == 0
            msk_d0 = ~(u << ld0) | (u << (ld1+1))
        else
            msk_d0 = ~(u << ld0)
            msk_d1 = (u << (ld1+1))
        end
        if delta_ks == 0
            msk_s0 = (u << ls0) & ~(u << (ls1+1))
        else
            msk_s0 = (u << ls0)
        end

        chunk_s0 = glue_src_bitchunks(chunks, ks0, ks1, msk_s0, ls0) & ~(u << s)
        chunks[kd0] = (chunks[kd0] & msk_d0) | ((chunk_s0 << ld0) & ~msk_d0)

        if delta_kd != 0
            chunk_s = (chunk_s0 >>> (64 - ld0))

            chunks[kd1] = (chunks[kd1] & msk_d1) | (chunk_s & ~msk_d1)
        end

        left -= s
        s = min(left, 64)
        b = left - s
        ps = pos_s + b
        pd = pos_d + b
    end
end

function fill_chunks!(Bc::Array{UInt64}, x::Bool, pos::Int, numbits::Int)
    numbits <= 0 && return
    k0, l0 = get_chunks_id(pos)
    k1, l1 = get_chunks_id(pos+numbits-1)

    u = _msk64
    if k1 == k0
        msk0 = (u << l0) & ~(u << (l1+1))
    else
        msk0 = (u << l0)
        msk1 = ~(u << (l1+1))
    end
    @inbounds if x
        Bc[k0] |= msk0
        for k = k0+1:k1-1
            Bc[k] = u
        end
        k1 > k0 && (Bc[k1] |= msk1)
    else
        Bc[k0] &= ~msk0
        for k = k0+1:k1-1
            Bc[k] = 0
        end
        k1 > k0 && (Bc[k1] &= ~msk1)
    end
end

copy_to_bitarray_chunks!(dest::Vector{UInt64}, pos_d::Int, src::BitArray, pos_s::Int, numbits::Int) =
    copy_chunks!(dest, pos_d, src.chunks, pos_s, numbits)

# pack 8 Bools encoded as one contiguous UIn64 into a single byte, e.g.:
# 0000001:0000001:00000000:00000000:00000001:00000000:00000000:00000001 → 11001001 → 0xc9
function pack8bools(z::UInt64)
    z |= z >>> 7
    z |= z >>> 14
    z |= z >>> 28
    z &= 0xFF
    return z
end

function copy_to_bitarray_chunks!(Bc::Vector{UInt64}, pos_d::Int, C::Array{Bool}, pos_s::Int, numbits::Int)
    kd0, ld0 = get_chunks_id(pos_d)
    kd1, ld1 = get_chunks_id(pos_d + numbits - 1)

    delta_kd = kd1 - kd0

    u = _msk64
    if delta_kd == 0
        msk_d0 = msk_d1 = ~(u << ld0) | (u << (ld1+1))
        lt0 = ld1
    else
        msk_d0 = ~(u << ld0)
        msk_d1 = (u << (ld1+1))
        lt0 = 63
    end

    bind = kd0
    ind = pos_s
    @inbounds if ld0 > 0
        c = UInt64(0)
        for j = ld0:lt0
            c |= (UInt64(C[ind]) << j)
            ind += 1
        end
        Bc[kd0] = (Bc[kd0] & msk_d0) | (c & ~msk_d0)
        bind += 1
    end

    nc = _div64(numbits - ind + pos_s)
    nc8 = (nc >>> 3) << 3
    if nc8 > 0
        ind8 = 1
        P8 = Ptr{UInt64}(pointer(C, ind)) # unaligned i64 pointer
        @inbounds for i = 1:nc8
            c = UInt64(0)
            for j = 0:7
                # unaligned load
                c |= (pack8bools(unsafe_load(P8, ind8)) << (j<<3))
                ind8 += 1
            end
            Bc[bind] = c
            bind += 1
        end
        ind += (ind8-1) << 3
    end
    @inbounds for i = (nc8+1):nc
        c = UInt64(0)
        for j = 0:63
            c |= (UInt64(C[ind]) << j)
            ind += 1
        end
        Bc[bind] = c
        bind += 1
    end
    @inbounds if bind ≤ kd1
        @assert bind == kd1
        c = UInt64(0)
        for j = 0:ld1
            c |= (UInt64(C[ind]) << j)
            ind += 1
        end
        Bc[kd1] = (Bc[kd1] & msk_d1) | (c & ~msk_d1)
    end
end

## More definitions in multidimensional.jl

# auxiliary definitions used when filling a BitArray via a Vector{Bool} cache
# (e.g. when constructing from an iterable, or in broadcast!)

const bitcache_chunks = 64 # this can be changed
const bitcache_size = 64 * bitcache_chunks # do not change this

dumpbitcache(Bc::Vector{UInt64}, bind::Int, C::Vector{Bool}) =
    copy_to_bitarray_chunks!(Bc, ((bind - 1) << 6) + 1, C, 1, min(bitcache_size, (length(Bc)-bind+1) << 6))

## custom iterator ##
function iterate(B::BitArray, i::Int=0)
    i >= length(B) && return nothing
    (B.chunks[_div64(i)+1] & (UInt64(1)<<_mod64(i)) != 0, i+1)
end

## similar, fill!, copy! etc ##

similar(B::BitArray) = BitArray(undef, size(B))
similar(B::BitArray, dims::Int...) = BitArray(undef, dims)
similar(B::BitArray, dims::Dims) = BitArray(undef, dims...)

similar(B::BitArray, T::Type{Bool}, dims::Dims) = BitArray(undef, dims)
# changing type to a non-Bool returns an Array
# (this triggers conversions like float(bitvector) etc.)
similar(B::BitArray, T::Type, dims::Dims) = Array{T}(undef, dims)

function fill!(B::BitArray, x)
    y = convert(Bool, x)
    isempty(B) && return B
    Bc = B.chunks
    if !y
        fill!(Bc, 0)
    else
        fill!(Bc, _msk64)
        Bc[end] &= _msk_end(B)
    end
    return B
end

"""
    falses(dims)

Create a `BitArray` with all values set to `false`.

# Examples
```jldoctest
julia> falses(2,3)
2×3 BitMatrix:
 0  0  0
 0  0  0
```
"""
falses(dims::DimOrInd...) = falses(dims)
falses(dims::NTuple{N, Union{Integer, OneTo}}) where {N} = falses(map(to_dim, dims))
falses(dims::NTuple{N, Integer}) where {N} = fill!(BitArray(undef, dims), false)
falses(dims::Tuple{}) = fill!(BitArray(undef, dims), false)

"""
    trues(dims)

Create a `BitArray` with all values set to `true`.

# Examples
```jldoctest
julia> trues(2,3)
2×3 BitMatrix:
 1  1  1
 1  1  1
```
"""
trues(dims::DimOrInd...) = trues(dims)
trues(dims::NTuple{N, Union{Integer, OneTo}}) where {N} = trues(map(to_dim, dims))
trues(dims::NTuple{N, Integer}) where {N} = fill!(BitArray(undef, dims), true)
trues(dims::Tuple{}) = fill!(BitArray(undef, dims), true)

function one(x::BitMatrix)
    m, n = size(x)
    m == n || throw(DimensionMismatch("multiplicative identity defined only for square matrices"))
    a = falses(n, n)
    for i = 1:n
        a[i,i] = true
    end
    return a
end

function copyto!(dest::BitArray, src::BitArray)
    length(src) > length(dest) && throw(BoundsError(dest, length(dest)+1))
    destc = dest.chunks; srcc = src.chunks
    nc = min(length(destc), length(srcc))
    nc == 0 && return dest
    @inbounds begin
        for i = 1 : nc - 1
            destc[i] = srcc[i]
        end
        if length(src) == length(dest)
            destc[nc] = srcc[nc]
        else
            msk_s = _msk_end(src)
            msk_d = ~msk_s
            destc[nc] = (msk_d & destc[nc]) | (msk_s & srcc[nc])
        end
    end
    return dest
end

function unsafe_copyto!(dest::BitArray, doffs::Integer, src::Union{BitArray,Array}, soffs::Integer, n::Integer)
    copy_to_bitarray_chunks!(dest.chunks, Int(doffs), src, Int(soffs), Int(n))
    return dest
end

copyto!(dest::BitArray, doffs::Integer, src::Array, soffs::Integer, n::Integer) =
    _copyto_int!(dest, Int(doffs), src, Int(soffs), Int(n))
function _copyto_int!(dest::BitArray, doffs::Int, src::Array, soffs::Int, n::Int)
    n == 0 && return dest
    soffs < 1 && throw(BoundsError(src, soffs))
    doffs < 1 && throw(BoundsError(dest, doffs))
    soffs+n-1 > length(src) && throw(BoundsError(src, length(src)+1))
    doffs+n-1 > length(dest) && throw(BoundsError(dest, length(dest)+1))
    return unsafe_copyto!(dest, doffs, src, soffs, n)
end

function copyto!(dest::BitArray, src::Array)
    length(src) > length(dest) && throw(BoundsError(dest, length(dest)+1))
    length(src) == 0 && return dest
    return unsafe_copyto!(dest, 1, src, 1, length(src))
end

function reshape(B::BitArray{N}, dims::NTuple{N,Int}) where N
    return dims == size(B) ? B : _bitreshape(B, dims)
end
reshape(B::BitArray, dims::Tuple{Vararg{Int}}) = _bitreshape(B, dims)
function _bitreshape(B::BitArray, dims::NTuple{N,Int}) where N
    prod(dims) == length(B) ||
        throw(DimensionMismatch("new dimensions $(dims) must be consistent with array size $(length(B))"))
    Br = BitArray{N}(undef, ntuple(i->0,Val(N))...)
    Br.chunks = B.chunks
    Br.len = prod(dims)
    N != 1 && (Br.dims = dims)
    return Br
end

## Constructors ##

function Array{T,N}(B::BitArray{N}) where {T,N}
    A = Array{T,N}(undef, size(B))
    Bc = B.chunks
    @inbounds for i = 1:length(A)
        A[i] = unsafe_bitgetindex(Bc, i)
    end
    return A
end

BitArray(A::AbstractArray{<:Any,N}) where {N} = BitArray{N}(A)
function BitArray{N}(A::AbstractArray{T,N}) where N where T
    B = BitArray(undef, convert(Dims{N}, size(A)::Dims{N}))
    Bc = B.chunks
    l = length(B)
    l == 0 && return B
    ind = 1
    @inbounds begin
        for i = 1:length(Bc)-1
            c = UInt64(0)
            for j = 0:63
                c |= (UInt64(convert(Bool, A[ind])::Bool) << j)
                ind += 1
            end
            Bc[i] = c
        end
        c = UInt64(0)
        for j = 0:_mod64(l-1)
            c |= (UInt64(convert(Bool, A[ind])::Bool) << j)
            ind += 1
        end
        Bc[end] = c
    end
    return B
end

function BitArray{N}(A::Array{Bool,N}) where N
    B = BitArray(undef, size(A))
    Bc = B.chunks
    l = length(B)
    l == 0 && return B
    copy_to_bitarray_chunks!(Bc, 1, A, 1, l)
    return B::BitArray{N}
end

reinterpret(::Type{Bool}, B::BitArray, dims::NTuple{N,Int}) where {N} = reinterpret(B, dims)
reinterpret(B::BitArray, dims::NTuple{N,Int}) where {N} = reshape(B, dims)

if nameof(@__MODULE__) === :Base  # avoid method overwrite
(::Type{T})(x::T) where {T<:BitArray} = copy(x)::T
BitArray(x::BitArray) = copy(x)
end

"""
    BitArray(itr)

Construct a [`BitArray`](@ref) generated by the given iterable object.
The shape is inferred from the `itr` object.

# Examples
```jldoctest
julia> BitArray([1 0; 0 1])
2×2 BitMatrix:
 1  0
 0  1

julia> BitArray(x+y == 3 for x = 1:2, y = 1:3)
2×3 BitMatrix:
 0  1  0
 1  0  0

julia> BitArray(x+y == 3 for x = 1:2 for y = 1:3)
6-element BitVector:
 0
 1
 0
 1
 0
 0
```
"""
BitArray(itr) = gen_bitarray(IteratorSize(itr), itr)
BitArray{N}(itr) where N = gen_bitarrayN(BitArray{N}, IteratorSize(itr), itr)

convert(T::Type{<:BitArray}, a::AbstractArray) = a isa T ? a : T(a)

# generic constructor from an iterable without compile-time info
# (we pass start(itr) explicitly to avoid a type-instability with filters)
gen_bitarray(isz::IteratorSize, itr) = gen_bitarray_from_itr(itr)

# generic iterable with known shape
function gen_bitarray(::HasShape, itr)
    B = BitArray(undef, size(itr))
    for (I,x) in zip(CartesianIndices(axes(itr)), itr)
        B[I] = x
    end
    return B
end

# generator with known shape or length
function gen_bitarray(::HasShape, itr::Generator)
    B = BitArray(undef, size(itr))
    return fill_bitarray_from_itr!(B, itr)
end
function gen_bitarray(::HasLength, itr)
    b = BitVector(undef, length(itr))
    return fill_bitarray_from_itr!(b, itr)
end

gen_bitarray(::IsInfinite, itr) =  throw(ArgumentError("infinite-size iterable used in BitArray constructor"))

gen_bitarrayN(::Type{BitVector}, itsz, itr)                        = gen_bitarray(itsz, itr)
gen_bitarrayN(::Type{BitVector}, itsz::HasShape{1}, itr)           = gen_bitarray(itsz, itr)
gen_bitarrayN(::Type{BitArray{N}}, itsz::HasShape{N}, itr) where N = gen_bitarray(itsz, itr)
# The first of these is just for ambiguity resolution
gen_bitarrayN(::Type{BitVector}, itsz::HasShape{N}, itr) where N      = throw(DimensionMismatch("cannot create a $T from a $N-dimensional iterator"))
gen_bitarrayN(@nospecialize(T::Type), itsz::HasShape{N}, itr) where N = throw(DimensionMismatch("cannot create a $T from a $N-dimensional iterator"))
gen_bitarrayN(@nospecialize(T::Type), itsz, itr) = throw(DimensionMismatch("cannot create a $T from a generic iterator"))

# The aux functions gen_bitarray_from_itr and fill_bitarray_from_itr! both
# use a Vector{Bool} cache for performance reasons

function gen_bitarray_from_itr(itr)
    B = empty!(BitVector(undef, bitcache_size))
    C = Vector{Bool}(undef, bitcache_size)
    Bc = B.chunks
    ind = 1
    cind = 1
    y = iterate(itr)
    while y !== nothing
        x, st = y
        @inbounds C[ind] = x
        ind += 1
        if ind > bitcache_size
            resize!(B, length(B) + bitcache_size)
            dumpbitcache(Bc, cind, C)
            cind += bitcache_chunks
            ind = 1
        end
        y = iterate(itr, st)
    end
    if ind > 1
        @inbounds C[ind:bitcache_size] .= false
        resize!(B, length(B) + ind - 1)
        dumpbitcache(Bc, cind, C)
    end
    return B
end

function fill_bitarray_from_itr!(B::BitArray, itr)
    n = length(B)
    C = Vector{Bool}(undef, bitcache_size)
    Bc = B.chunks
    ind = 1
    cind = 1
    y = iterate(itr)
    while y !== nothing
        x, st = y
        @inbounds C[ind] = x
        ind += 1
        if ind > bitcache_size
            dumpbitcache(Bc, cind, C)
            cind += bitcache_chunks
            ind = 1
        end
        y = iterate(itr, st)
    end
    if ind > 1
        @inbounds C[ind:bitcache_size] .= false
        dumpbitcache(Bc, cind, C)
    end
    return B
end

## Indexing: getindex ##

@inline function unsafe_bitgetindex(Bc::Vector{UInt64}, i::Int)
    i1, i2 = get_chunks_id(i)
    u = UInt64(1) << i2
    @inbounds r = (Bc[i1] & u) != 0
    return r
end

@inline function getindex(B::BitArray, i::Int)
    @boundscheck checkbounds(B, i)
    unsafe_bitgetindex(B.chunks, i)
end

## Indexing: setindex! ##

@inline function unsafe_bitsetindex!(Bc::Array{UInt64}, x::Bool, i::Int)
    i1, i2 = get_chunks_id(i)
    _unsafe_bitsetindex!(Bc, x, i1, i2)
end

@inline function _unsafe_bitsetindex!(Bc::Array{UInt64}, x::Bool, i1::Int, i2::Int)
    u = UInt64(1) << i2
    @inbounds begin
        c = Bc[i1]
        Bc[i1] = ifelse(x, c | u, c & ~u)
    end
end

@inline function setindex!(B::BitArray, x, i::Int)
    @boundscheck checkbounds(B, i)
    unsafe_bitsetindex!(B.chunks, convert(Bool, x), i)
    return B
end

indexoffset(i) = first(i)-1
indexoffset(::Colon) = 0

@propagate_inbounds function setindex!(B::BitArray, X::AbstractArray, J0::Union{Colon,UnitRange{Int}})
    _setindex!(IndexStyle(B), B, X, to_indices(B, (J0,))[1])
end

# Assigning an array of bools is more complicated, but we can still do some
# work on chunks by combining X and I 64 bits at a time to improve perf by ~40%
@inline function setindex!(B::BitArray, X::AbstractArray, I::BitArray)
    @boundscheck checkbounds(B, I)
    _unsafe_setindex!(B, X, I)
end
function _unsafe_setindex!(B::BitArray, X::AbstractArray, I::BitArray)
    Bc = B.chunks
    Ic = I.chunks
    length(Bc) == length(Ic) || throw_boundserror(B, I)
    lc = length(Bc)
    lx = length(X)
    last_chunk_len = _mod64(length(B)-1)+1

    c = 1
    for i = 1:lc
        @inbounds Imsk = Ic[i]
        @inbounds C = Bc[i]
        u = UInt64(1)
        for j = 1:(i < lc ? 64 : last_chunk_len)
            if Imsk & u != 0
                lx < c && throw_setindex_mismatch(X, c)
                @inbounds x = convert(Bool, X[c])
                C = ifelse(x, C | u, C & ~u)
                c += 1
            end
            u <<= 1
        end
        @inbounds Bc[i] = C
    end
    if length(X) != c-1
        throw_setindex_mismatch(X, c-1)
    end
    return B
end

## Dequeue functionality ##

function push!(B::BitVector, item)
    # convert first so we don't grow the bitarray if the assignment won't work
    item = convert(Bool, item)

    Bc = B.chunks

    l = _mod64(length(B))
    if l == 0
        _growend!(Bc, 1)
        Bc[end] = UInt64(0)
    end
    B.len += 1
    if item
        B[end] = true
    end
    return B
end

function append!(B::BitVector, items::BitVector)
    n0 = length(B)
    n1 = length(items)
    n1 == 0 && return B
    Bc = B.chunks
    k0 = length(Bc)
    k1 = num_bit_chunks(n0 + n1)
    if k1 > k0
        _growend!(Bc, k1 - k0)
        Bc[end] = UInt64(0)
    end
    B.len += n1
    copy_chunks!(Bc, n0+1, items.chunks, 1, n1)
    return B
end

append!(B::BitVector, items) = append!(B, BitVector(items))
append!(A::Vector{Bool}, items::BitVector) = append!(A, Array(items))

function prepend!(B::BitVector, items::BitVector)
    n0 = length(B)
    n1 = length(items)
    n1 == 0 && return B
    Bc = B.chunks
    k0 = length(Bc)
    k1 = num_bit_chunks(n0 + n1)
    if k1 > k0
        _growend!(Bc, k1 - k0)
        Bc[end] = UInt64(0)
    end
    B.len += n1
    copy_chunks!(Bc, 1 + n1, Bc, 1, n0)
    copy_chunks!(Bc, 1, items.chunks, 1, n1)
    return B
end

prepend!(B::BitVector, items) = prepend!(B, BitArray(items))
prepend!(A::Vector{Bool}, items::BitVector) = prepend!(A, Array(items))

function sizehint!(B::BitVector, sz::Integer)
    ccall(:jl_array_sizehint, Cvoid, (Any, UInt), B.chunks, num_bit_chunks(sz))
    return B
end

resize!(B::BitVector, n::Integer) = _resize_int!(B, Int(n))
function _resize_int!(B::BitVector, n::Int)
    n0 = length(B)
    n == n0 && return B
    n >= 0 || throw(BoundsError(B, n))
    if n < n0
        deleteat!(B, n+1:n0)
        return B
    end
    Bc = B.chunks
    k0 = length(Bc)
    k1 = num_bit_chunks(n)
    if k1 > k0
        _growend!(Bc, k1 - k0)
        Bc[end] = UInt64(0)
    end
    B.len = n
    return B
end

function pop!(B::BitVector)
    isempty(B) && throw(ArgumentError("argument must not be empty"))
    item = B[end]
    B[end] = false

    l = _mod64(length(B))
    l == 1 && _deleteend!(B.chunks, 1)
    B.len -= 1

    return item
end

function pushfirst!(B::BitVector, item)
    item = convert(Bool, item)

    Bc = B.chunks

    l = _mod64(length(B))
    if l == 0
        _growend!(Bc, 1)
        Bc[end] = UInt64(0)
    end
    B.len += 1
    if B.len == 1
        Bc[1] = item
        return B
    end
    for i = length(Bc) : -1 : 2
        Bc[i] = (Bc[i] << 1) | (Bc[i-1] >>> 63)
    end
    Bc[1] = UInt64(item) | (Bc[1] << 1)
    return B
end

function popfirst!(B::BitVector)
    isempty(B) && throw(ArgumentError("argument must not be empty"))
    @inbounds begin
        item = B[1]

        Bc = B.chunks

        for i = 1 : length(Bc) - 1
            Bc[i] = (Bc[i] >>> 1) | (Bc[i+1] << 63)
        end

        l = _mod64(length(B))
        if l == 1
            _deleteend!(Bc, 1)
        else
            Bc[end] >>>= 1
        end
        B.len -= 1
    end

    return item
end

insert!(B::BitVector, i::Integer, item) = _insert_int!(B, Int(i), item)
function _insert_int!(B::BitVector, i::Int, item)
    i = Int(i)
    n = length(B)
    1 <= i <= n+1 || throw(BoundsError(B, i))
    item = convert(Bool, item)

    Bc = B.chunks

    k, j = get_chunks_id(i)

    l = _mod64(length(B))
    if l == 0
        _growend!(Bc, 1)
        Bc[end] = UInt64(0)
    end
    B.len += 1

    for t = length(Bc) : -1 : k + 1
        Bc[t] = (Bc[t] << 1) | (Bc[t - 1] >>> 63)
    end

    msk_aft = (_msk64 << j)
    msk_bef = ~msk_aft
    Bc[k] = (msk_bef & Bc[k]) | ((msk_aft & Bc[k]) << 1)
    B[i] = item
    B
end

function _deleteat!(B::BitVector, i::Int)
    k, j = get_chunks_id(i)

    msk_bef = _msk64 >>> (63 - j)
    msk_aft = ~msk_bef
    msk_bef >>>= 1

    Bc = B.chunks

    @inbounds begin
        Bc[k] = (msk_bef & Bc[k]) | ((msk_aft & Bc[k]) >> 1)
        if length(Bc) > k
            Bc[k] |= (Bc[k + 1] << 63)
        end

        for t = k + 1 : length(Bc) - 1
            Bc[t] = (Bc[t] >>> 1) | (Bc[t + 1] << 63)
        end

        l = _mod64(length(B))

        if l == 1
            _deleteend!(Bc, 1)
        elseif length(Bc) > k
            Bc[end] >>>= 1
        end
    end

    B.len -= 1

    return B
end

function deleteat!(B::BitVector, i::Integer)
    i = Int(i)
    n = length(B)
    1 <= i <= n || throw(BoundsError(B, i))

    return _deleteat!(B, i)
end

function deleteat!(B::BitVector, r::UnitRange{Int})
    n = length(B)
    i_f = first(r)
    i_l = last(r)
    1 <= i_f || throw(BoundsError(B, i_f))
    i_l <= n || throw(BoundsError(B, n+1))

    Bc = B.chunks
    new_l = length(B) - length(r)
    delta_k = num_bit_chunks(new_l) - length(Bc)

    copy_chunks!(Bc, i_f, Bc, i_l+1, n-i_l)

    delta_k < 0 && _deleteend!(Bc, -delta_k)

    B.len = new_l

    if new_l > 0
        Bc[end] &= _msk_end(new_l)
    end

    return B
end

function deleteat!(B::BitVector, inds)
    n = new_l = length(B)
    y = iterate(inds)
    y === nothing && return B

    Bc = B.chunks

    (p, s) = y
    checkbounds(B, p)
    q = p+1
    new_l -= 1
    y = iterate(inds, s)
    while y !== nothing
        (i, s) = y
        if !(q <= i <= n)
            i < q && throw(ArgumentError("indices must be unique and sorted"))
            throw(BoundsError(B, i))
        end
        new_l -= 1
        if i > q
            copy_chunks!(Bc, p, Bc, Int(q), Int(i-q))
            p += i-q
        end
        q = i+1
        y = iterate(inds, s)
    end

    q <= n && copy_chunks!(Bc, p, Bc, Int(q), Int(n-q+1))

    delta_k = num_bit_chunks(new_l) - length(Bc)
    delta_k < 0 && _deleteend!(Bc, -delta_k)

    B.len = new_l

    if new_l > 0
        Bc[end] &= _msk_end(new_l)
    end

    return B
end

function splice!(B::BitVector, i::Integer)
    i = Int(i)
    n = length(B)
    1 <= i <= n || throw(BoundsError(B, i))

    v = B[i]   # TODO: change to a copy if/when subscripting becomes an ArrayView
    _deleteat!(B, i)
    return v
end

const _default_bit_splice = BitVector()

function splice!(B::BitVector, r::Union{UnitRange{Int}, Integer}, ins::AbstractArray = _default_bit_splice)
    _splice_int!(B, isa(r, UnitRange{Int}) ? r : Int(r), ins)
end
function _splice_int!(B::BitVector, r, ins)
    n = length(B)
    i_f, i_l = first(r), last(r)
    1 <= i_f <= n+1 || throw(BoundsError(B, i_f))
    i_l <= n || throw(BoundsError(B, n+1))

    Bins = convert(BitArray, ins)

    if (i_f > n)
        append!(B, Bins)
        return BitVector()
    end

    v = B[r]  # TODO: change to a copy if/when subscripting becomes an ArrayView

    Bc = B.chunks

    lins = length(Bins)
    ldel = length(r)

    new_l = length(B) + lins - ldel
    delta_k = num_bit_chunks(new_l) - length(Bc)

    delta_k > 0 && _growend!(Bc, delta_k)

    copy_chunks!(Bc, i_f+lins, Bc, i_l+1, n-i_l)
    copy_chunks!(Bc, i_f, Bins.chunks, 1, lins)

    delta_k < 0 && _deleteend!(Bc, -delta_k)

    B.len = new_l

    if new_l > 0
        Bc[end] &= _msk_end(new_l)
    end

    return v
end

function splice!(B::BitVector, r::Union{UnitRange{Int}, Integer}, ins)
    Bins = BitVector(undef, length(ins))
    i = 1
    for x in ins
        Bins[i] = Bool(x)
        i += 1
    end
    return splice!(B, r, Bins)
end

function empty!(B::BitVector)
    _deleteend!(B.chunks, length(B.chunks))
    B.len = 0
    return B
end

## Unary operators ##

function (-)(B::BitArray)
    A = zeros(Int, size(B))
    l = length(B)
    l == 0 && return A
    Bc = B.chunks
    ind = 1
    for i = 1:length(Bc)-1
        u = UInt64(1)
        c = Bc[i]
        for j = 1:64
            if c & u != 0
                A[ind] = -1
            end
            ind += 1
            u <<= 1
        end
    end
    u = UInt64(1)
    c = Bc[end]
    for j = 0:_mod64(l-1)
        if c & u != 0
            A[ind] = -1
        end
        ind += 1
        u <<= 1
    end
    return A
end

## Binary arithmetic operators ##

for f in (:+, :-)
    @eval function ($f)(A::BitArray, B::BitArray)
        r = Array{Int}(undef, promote_shape(size(A), size(B)))
        ay, by = iterate(A), iterate(B)
        ri = 1
        # promote_shape guarantees that A and B have the
        # same iteration space
        while ay !== nothing
            @inbounds r[ri] = ($f)(ay[1], by[1])
            ri += 1
            ay, by = iterate(A, ay[2]), iterate(B, by[2])
        end
        return r
    end
end

for f in (:/, :\)
    @eval begin
        ($f)(A::Union{BitMatrix,BitVector}, B::Union{BitMatrix,BitVector}) = ($f)(Array(A), Array(B))
    end
end
(/)(B::BitArray, x::Number) = (/)(Array(B), x)
(/)(x::Number, B::BitArray) = (/)(x, Array(B))

## promotion to complex ##

# TODO?

## comparison operators ##

function (==)(A::BitArray, B::BitArray)
    size(A) != size(B) && return false
    return A.chunks == B.chunks
end

## Data movement ##

# TODO some of this could be optimized

reverse(A::BitArray; dims::Integer) = _reverse_int(A, Int(dims))
function _reverse_int(A::BitArray, d::Int)
    nd = ndims(A)
    1 ≤ d ≤ nd || throw(ArgumentError("dimension $d is not 1 ≤ $d ≤ $nd"))
    sd = size(A, d)
    sd == 1 && return copy(A)

    B = similar(A)

    nnd = 0
    for i = 1:nd
        nnd += size(A,i)==1 || i==d
    end
    if nnd == nd
        # reverse along the only non-singleton dimension
        for i = 1:sd
            B[i] = A[sd+1-i]
        end
        return B
    end

    d_in = size(A)
    leading = d_in[1:(d-1)]
    M = prod(leading)
    N = length(A)
    stride = M * sd

    if M == 1
        for j = 0:stride:(N-stride)
            for i = 1:sd
                ri = sd+1-i
                B[j + ri] = A[j + i]
            end
        end
    else
        for i = 1:sd
            ri = sd+1-i
            for j=0:stride:(N-stride)
                offs = j + 1 + (i-1)*M
                boffs = j + 1 + (ri-1)*M
                copy_chunks!(B.chunks, boffs, A.chunks, offs, M)
            end
        end
    end
    return B
end

function reverse!(B::BitVector)
    # Basic idea: each chunk is divided into two blocks of size k = n % 64, and
    # h = 64 - k. Walk from either end (with indices i and j) reversing chunks
    # and separately ORing their two blocks into place.
    #
    #           chunk 3                  chunk 2                  chunk 1
    # ┌───────────────┬───────┐┌───────────────┬───────┐┌───────────────┬───────┐
    # │000000000000000│   E   ││       D       │   C   ││       B       │   A   │
    # └───────────────┴───────┘└───────────────┴───────┘└───────────────┴───────┘
    #                     k            h           k            h            k
    # yielding;
    # ┌───────────────┬───────┐┌───────────────┬───────┐┌───────────────┬───────┐
    # │000000000000000│  A'   ││      B'       │  C'   ││      D'       │  E'   │
    # └───────────────┴───────┘└───────────────┴───────┘└───────────────┴───────┘

    n = length(B)
    n == 0 && return B

    k = _mod64(n+63) + 1
    h = 64 - k

    i, j = 0, length(B.chunks)
    u = UInt64(0)
    v = bitreverse(B.chunks[j])
    B.chunks[j] = 0
    @inbounds while true
        i += 1
        if i == j
            break
        end
        u = bitreverse(B.chunks[i])
        B.chunks[i] = 0
        B.chunks[j] |= u >>> h
        B.chunks[i] |= v >>> h

        j -= 1
        if i == j
            break
        end
        v = bitreverse(B.chunks[j])
        B.chunks[j] = 0
        B.chunks[i] |= v << k
        B.chunks[j] |= u << k
    end

    if isodd(length(B.chunks))
        B.chunks[i] |= v >>> h
    else
        B.chunks[i] |= u << k
    end

    return B
end

reverse(v::BitVector) = reverse!(copy(v))

function (<<)(B::BitVector, i::UInt)
    n = length(B)
    i == 0 && return copy(B)
    A = falses(n)
    i < n && copy_chunks!(A.chunks, 1, B.chunks, Int(i+1), Int(n-i))
    return A
end

function (>>>)(B::BitVector, i::UInt)
    n = length(B)
    i == 0 && return copy(B)
    A = falses(n)
    i < n && copy_chunks!(A.chunks, Int(i+1), B.chunks, 1, Int(n-i))
    return A
end

"""
    >>(B::BitVector, n) -> BitVector

Right bit shift operator, `B >> n`. For `n >= 0`, the result is `B`
with elements shifted `n` positions forward, filling with `false`
values. If `n < 0`, elements are shifted backwards. Equivalent to
`B << -n`.

# Examples
```jldoctest
julia> B = BitVector([true, false, true, false, false])
5-element BitVector:
 1
 0
 1
 0
 0

julia> B >> 1
5-element BitVector:
 0
 1
 0
 1
 0

julia> B >> -1
5-element BitVector:
 0
 1
 0
 0
 0
```
"""
(>>)(B::BitVector, i::Union{Int, UInt}) = B >>> i

# signed integer version of shift operators with handling of negative values
"""
    <<(B::BitVector, n) -> BitVector

Left bit shift operator, `B << n`. For `n >= 0`, the result is `B`
with elements shifted `n` positions backwards, filling with `false`
values. If `n < 0`, elements are shifted forwards. Equivalent to
`B >> -n`.

# Examples
```jldoctest
julia> B = BitVector([true, false, true, false, false])
5-element BitVector:
 1
 0
 1
 0
 0

julia> B << 1
5-element BitVector:
 0
 1
 0
 0
 0

julia> B << -1
5-element BitVector:
 0
 1
 0
 1
 0
```
"""
(<<)(B::BitVector, i::Int) = (i >=0 ? B << unsigned(i) : B >> unsigned(-i))

"""
    >>>(B::BitVector, n) -> BitVector

Unsigned right bitshift operator, `B >>> n`. Equivalent to `B >> n`. See [`>>`](@ref) for
details and examples.
"""
(>>>)(B::BitVector, i::Int) = (i >=0 ? B >> unsigned(i) : B << unsigned(-i))

circshift!(dest::BitVector, src::BitVector, i::Integer) = _circshift_int!(dest, src, Int(i))
function _circshift_int!(dest::BitVector, src::BitVector, i::Int)
    i = Int(i)
    length(dest) == length(src) || throw(ArgumentError("destination and source should be of same size"))
    n = length(dest)
    i %= n
    i == 0 && return (src === dest ? src : copyto!(dest, src))
    Bc = (src === dest ? copy(src.chunks) : src.chunks)
    if i > 0 # right
        copy_chunks!(dest.chunks, i+1, Bc, 1, n-i)
        copy_chunks!(dest.chunks, 1, Bc, n-i+1, i)
    else # left
        i = -i
        copy_chunks!(dest.chunks, 1, Bc, i+1, n-i)
        copy_chunks!(dest.chunks, n-i+1, Bc, 1, i)
    end
    return dest
end

circshift!(B::BitVector, i::Integer) = circshift!(B, B, i)

## count & find ##

function bitcount(Bc::Vector{UInt64})
    n = 0
    @inbounds for i = 1:length(Bc)
        n += count_ones(Bc[i])
    end
    return n
end

count(B::BitArray) = bitcount(B.chunks)

function unsafe_bitfindnext(Bc::Vector{UInt64}, start::Int)
    chunk_start = _div64(start-1)+1
    within_chunk_start = _mod64(start-1)
    mask = _msk64 << within_chunk_start

    @inbounds begin
        if Bc[chunk_start] & mask != 0
            return (chunk_start-1) << 6 + trailing_zeros(Bc[chunk_start] & mask) + 1
        end

        for i = chunk_start+1:length(Bc)
            if Bc[i] != 0
                return (i-1) << 6 + trailing_zeros(Bc[i]) + 1
            end
        end
    end
    return nothing
end

# returns the index of the next true element, or nothing if all false
function findnext(B::BitArray, start::Integer)
    start = Int(start)
    start > 0 || throw(BoundsError(B, start))
    start > length(B) && return nothing
    unsafe_bitfindnext(B.chunks, start)
end

#findfirst(B::BitArray) = findnext(B, 1)  ## defined in array.jl

# aux function: same as findnext(~B, start), but performed without temporaries
function findnextnot(B::BitArray, start::Int)
    start > 0 || throw(BoundsError(B, start))
    start > length(B) && return nothing

    Bc = B.chunks
    l = length(Bc)
    l == 0 && return nothing

    chunk_start = _div64(start-1)+1
    within_chunk_start = _mod64(start-1)
    mask = ~(_msk64 << within_chunk_start)

    @inbounds if chunk_start < l
        if Bc[chunk_start] | mask != _msk64
            return (chunk_start-1) << 6 + trailing_ones(Bc[chunk_start] | mask) + 1
        end
        for i = chunk_start+1:l-1
            if Bc[i] != _msk64
                return (i-1) << 6 + trailing_ones(Bc[i]) + 1
            end
        end
        if Bc[l] != _msk_end(B)
            return (l-1) << 6 + trailing_ones(Bc[l]) + 1
        end
    elseif Bc[l] | mask != _msk_end(B)
        return (l-1) << 6 + trailing_ones(Bc[l] | mask) + 1
    end
    return nothing
end
findfirstnot(B::BitArray) = findnextnot(B,1)

# returns the index of the first matching element
function findnext(pred::Fix2{<:Union{typeof(isequal),typeof(==)},Bool},
                  B::BitArray, start::Integer)
    v = pred.x
    v == false && return findnextnot(B, Int(start))
    v == true && return findnext(B, start)
    return nothing
end
#findfirst(B::BitArray, v) = findnext(B, 1, v)  ## defined in array.jl

# returns the index of the first element for which the function returns true
findnext(testf::Function, B::BitArray, start::Integer) = _findnext_int(testf, B, Int(start))
function _findnext_int(testf::Function, B::BitArray, start::Int)
    f0::Bool = testf(false)
    f1::Bool = testf(true)
    !f0 && f1 && return findnext(B, start)
    f0 && !f1 && return findnextnot(B, start)

    start > 0 || throw(BoundsError(B, start))
    start > length(B) && return nothing
    f0 && f1 && return start
    return nothing # last case: !f0 && !f1
end
#findfirst(testf::Function, B::BitArray) = findnext(testf, B, 1)  ## defined in array.jl

function unsafe_bitfindprev(Bc::Vector{UInt64}, start::Int)
    chunk_start = _div64(start-1)+1
    mask = _msk_end(start)

    @inbounds begin
        if Bc[chunk_start] & mask != 0
            return (chunk_start-1) << 6 + (64 - leading_zeros(Bc[chunk_start] & mask))
        end

        for i = (chunk_start-1):-1:1
            if Bc[i] != 0
                return (i-1) << 6 + (64 - leading_zeros(Bc[i]))
            end
        end
    end
    return nothing
end

# returns the index of the previous true element, or nothing if all false
function findprev(B::BitArray, start::Integer)
    start = Int(start)
    start > 0 || return nothing
    start > length(B) && throw(BoundsError(B, start))
    unsafe_bitfindprev(B.chunks, start)
end

function findprevnot(B::BitArray, start::Int)
    start = Int(start)
    start > 0 || return nothing
    start > length(B) && throw(BoundsError(B, start))

    Bc = B.chunks

    chunk_start = _div64(start-1)+1
    mask = ~_msk_end(start)

    @inbounds begin
        if Bc[chunk_start] | mask != _msk64
            return (chunk_start-1) << 6 + (64 - leading_ones(Bc[chunk_start] | mask))
        end

        for i = chunk_start-1:-1:1
            if Bc[i] != _msk64
                return (i-1) << 6 + (64 - leading_ones(Bc[i]))
            end
        end
    end
    return nothing
end
findlastnot(B::BitArray) = findprevnot(B, length(B))

# returns the index of the previous matching element
function findprev(pred::Fix2{<:Union{typeof(isequal),typeof(==)},Bool},
                  B::BitArray, start::Integer)
    v = pred.x
    v == false && return findprevnot(B, Int(start))
    v == true && return findprev(B, start)
    return nothing
end
#findlast(B::BitArray, v) = findprev(B, 1, v)  ## defined in array.jl

# returns the index of the previous element for which the function returns true
findprev(testf::Function, B::BitArray, start::Integer) = _findprev_int(testf, B, Int(start))
function _findprev_int(testf::Function, B::BitArray, start::Int)
    f0::Bool = testf(false)
    f1::Bool = testf(true)
    !f0 && f1 && return findprev(B, start)
    f0 && !f1 && return findprevnot(B, start)

    start > 0 || return nothing
    start > length(B) && throw(BoundsError(B, start))
    f0 && f1 && return start
    return nothing # last case: !f0 && !f1
end
#findlast(testf::Function, B::BitArray) = findprev(testf, B, 1)  ## defined in array.jl

function findmax(a::BitArray)
    isempty(a) && throw(ArgumentError("BitArray must be non-empty"))
    m, mi = false, 1
    ti = 1
    ac = a.chunks
    for i = 1:length(ac)
        @inbounds k = trailing_zeros(ac[i])
        ti += k
        k == 64 || return (true, @inbounds keys(a)[ti])
    end
    return m, @inbounds keys(a)[mi]
end

function findmin(a::BitArray)
    isempty(a) && throw(ArgumentError("BitArray must be non-empty"))
    m, mi = true, 1
    ti = 1
    ac = a.chunks
    for i = 1:length(ac)-1
        @inbounds k = trailing_ones(ac[i])
        ti += k
        k == 64 || return (false, @inbounds keys(a)[ti])
    end
    l = Base._mod64(length(a)-1) + 1
    @inbounds k = trailing_ones(ac[end] & Base._msk_end(l))
    ti += k
    k == l || return (false, @inbounds keys(a)[ti])
    return (m, @inbounds keys(a)[mi])
end

# findall helper functions
# Generic case (>2 dimensions)
function allindices!(I, B::BitArray)
    ind = first(keys(B))
    for k = 1:length(B)
        I[k] = ind
        ind = nextind(B, ind)
    end
end

# Optimized case for vector
function allindices!(I, B::BitVector)
    I[:] .= 1:length(B)
end

# Optimized case for matrix
function allindices!(I, B::BitMatrix)
    k = 1
    for c = 1:size(B,2), r = 1:size(B,1)
        I[k] = CartesianIndex(r, c)
        k += 1
    end
end

@inline _overflowind(i1, irest::Tuple{}, size) = (i1, irest)
@inline function _overflowind(i1, irest, size)
    i2 = irest[1]
    while i1 > size[1]
        i1 -= size[1]
        i2 += 1
    end
    i2, irest = _overflowind(i2, tail(irest), tail(size))
    return (i1, (i2, irest...))
end

@inline _toind(i1, irest::Tuple{}) = i1
@inline _toind(i1, irest) = CartesianIndex(i1, irest...)

function findall(B::BitArray)
    nnzB = count(B)
    I = Vector{eltype(keys(B))}(undef, nnzB)
    nnzB == 0 && return I
    nnzB == length(B) && (allindices!(I, B); return I)
    Bc = B.chunks
    Bs = size(B)
    Bi = i1 = i = 1
    irest = ntuple(one, ndims(B) - 1)
    c = Bc[1]
    @inbounds while true
        while c == 0
            Bi == length(Bc) && return I
            i1 += 64
            Bi += 1
            c = Bc[Bi]
        end

        tz = trailing_zeros(c)
        c = _blsr(c)

        i1, irest = _overflowind(i1 + tz, irest, Bs)
        I[i] = _toind(i1, irest)
        i += 1
        i1 -= tz
    end
end

# For performance
findall(::typeof(!iszero), B::BitArray) = findall(B)

## Reductions ##

_sum(A::BitArray, dims)    = reduce(+, A, dims=dims)
_sum(B::BitArray, ::Colon) = count(B)

function all(B::BitArray)
    isempty(B) && return true
    Bc = B.chunks
    @inbounds begin
        for i = 1:length(Bc)-1
            Bc[i] == _msk64 || return false
        end
        Bc[end] == _msk_end(B) || return false
    end
    return true
end

function any(B::BitArray)
    isempty(B) && return false
    Bc = B.chunks
    @inbounds begin
        for i = 1:length(Bc)
            Bc[i] == 0 || return true
        end
    end
    return false
end

minimum(B::BitArray) = isempty(B) ? throw(ArgumentError("argument must be non-empty")) : all(B)
maximum(B::BitArray) = isempty(B) ? throw(ArgumentError("argument must be non-empty")) : any(B)

## map over bitarrays ##

# Specializing map is even more important for bitarrays than it is for generic
# arrays since there can be a 64x speedup by working at the level of Int64
# instead of looping bit-by-bit.

map(::Union{typeof(~), typeof(!)}, A::BitArray) = bit_map!(~, similar(A), A)
map(::typeof(zero), A::BitArray) = fill!(similar(A), false)
map(::typeof(one), A::BitArray) = fill!(similar(A), true)
map(::typeof(identity), A::BitArray) = copy(A)

map!(::Union{typeof(~), typeof(!)}, dest::BitArray, A::BitArray) = bit_map!(~, dest, A)
map!(::typeof(zero), dest::BitArray, A::BitArray) = fill!(dest, false)
map!(::typeof(one), dest::BitArray, A::BitArray) = fill!(dest, true)
map!(::typeof(identity), dest::BitArray, A::BitArray) = copyto!(dest, A)

for (T, f) in ((:(Union{typeof(&), typeof(*), typeof(min)}), :(&)),
               (:(Union{typeof(|), typeof(max)}),            :(|)),
               (:(Union{typeof(xor), typeof(!=)}),           :xor),
               (:(Union{typeof(>=), typeof(^)}),             :((p, q) -> p | ~q)),
               (:(typeof(<=)),                               :((p, q) -> ~p | q)),
               (:(typeof(==)),                               :((p, q) -> ~xor(p, q))),
               (:(typeof(<)),                                :((p, q) -> ~p & q)),
               (:(typeof(>)),                                :((p, q) -> p & ~q)))
    @eval map(::$T, A::BitArray, B::BitArray) = bit_map!($f, similar(A), A, B)
    @eval map!(::$T, dest::BitArray, A::BitArray, B::BitArray) = bit_map!($f, dest, A, B)
end

# If we were able to specialize the function to a known bitwise operation,
# map across the chunks. Otherwise, fall-back to the AbstractArray method that
# iterates bit-by-bit.
function bit_map!(f::F, dest::BitArray, A::BitArray) where F
    size(A) == size(dest) || throw(DimensionMismatch("sizes of dest and A must match"))
    isempty(A) && return dest
    destc = dest.chunks
    Ac = A.chunks
    for i = 1:(length(Ac)-1)
        destc[i] = f(Ac[i])
    end
    destc[end] = f(Ac[end]) & _msk_end(A)
    dest
end
function bit_map!(f::F, dest::BitArray, A::BitArray, B::BitArray) where F
    size(A) == size(B) == size(dest) || throw(DimensionMismatch("sizes of dest, A, and B must all match"))
    isempty(A) && return dest
    destc = dest.chunks
    Ac = A.chunks
    Bc = B.chunks
    for i = 1:(length(Ac)-1)
        destc[i] = f(Ac[i], Bc[i])
    end
    destc[end] = f(Ac[end], Bc[end]) & _msk_end(A)
    dest
end

## Filter ##

function filter(f, Bs::BitArray)
    boolmap::Array{Bool} = map(f, Bs)
    Bs[boolmap]
end

## Concatenation ##

function hcat(B::BitVector...)
    height = length(B[1])
    for j = 2:length(B)
        length(B[j]) == height ||
            throw(DimensionMismatch("dimensions must match"))
    end
    M = BitMatrix(undef, height, length(B))
    for j = 1:length(B)
        copy_chunks!(M.chunks, (height*(j-1))+1, B[j].chunks, 1, height)
    end
    return M
end

function vcat(V::BitVector...)
    n = 0
    for Vk in V
        n += length(Vk)
    end
    B = BitVector(undef, n)
    j = 1
    for Vk in V
        copy_chunks!(B.chunks, j, Vk.chunks, 1, length(Vk))
        j += length(Vk)
    end
    return B
end

function hcat(A::Union{BitMatrix,BitVector}...)
    nargs = length(A)
    nrows = size(A[1], 1)
    ncols = 0
    dense = true
    for j = 1:nargs
        Aj = A[j]
        nd = ndims(Aj)
        ncols += (nd==2 ? size(Aj,2) : 1)
        size(Aj, 1) == nrows ||
            throw(DimensionMismatch("row lengths must match"))
    end

    B = BitMatrix(undef, nrows, ncols)

    pos = 1
    for k = 1:nargs
        Ak = A[k]
        n = length(Ak)
        copy_chunks!(B.chunks, pos, Ak.chunks, 1, n)
        pos += n
    end
    return B
end

function vcat(A::BitMatrix...)
    nargs = length(A)
    nrows, nrowsA = 0, sizehint!(Int[], nargs)
    for a in A
        sz1 = size(a, 1)
        nrows += sz1
        push!(nrowsA, sz1)
    end
    ncols = size(A[1], 2)
    for j = 2:nargs
        size(A[j], 2) == ncols ||
            throw(DimensionMismatch("column lengths must match"))
    end
    B = BitMatrix(undef, nrows, ncols)
    Bc = B.chunks
    pos_d = 1
    pos_s = fill(1, nargs)
    for j = 1:ncols, k = 1:nargs
        copy_chunks!(Bc, pos_d, A[k].chunks, pos_s[k], nrowsA[k])
        pos_s[k] += nrowsA[k]
        pos_d += nrowsA[k]
    end
    return B
end

# general case, specialized for BitArrays and Integers
_cat(dims::Integer, X::Union{BitArray, Bool}...) = _cat(Int(dims)::Int, X...)
function _cat(dims::Int, X::Union{BitArray, Bool}...)
    dims = Int(dims)
    catdims = dims2cat(dims)
    shape = cat_shape(catdims, map(cat_size, X))
    A = falses(shape)
    return __cat(A, shape, catdims, X...)
end

# hvcat -> use fallbacks in abstractarray.jl

# BitArray I/O

write(s::IO, B::BitArray) = write(s, B.chunks)
function read!(s::IO, B::BitArray)
    n = length(B)
    Bc = B.chunks
    nc = length(read!(s, Bc))
    if length(Bc) > 0 && Bc[end] & _msk_end(n) ≠ Bc[end]
        Bc[end] &= _msk_end(n) # ensure that the BitArray is not broken
        throw(DimensionMismatch("read mismatch, found non-zero bits after BitArray length"))
    end
    return B
end

sizeof(B::BitArray) = sizeof(B.chunks)

max_stars_repo_name: ExpandingMan/julia

max_issues_repo_name: ExpandingMan/julia

max_forks_repo_name: ExpandingMan/julia


text:
function AbsoluteDate()
    return AbsoluteDate(())
end

function AbsoluteDate(arg0::AbsoluteDate, arg1::jdouble)
    return AbsoluteDate((AbsoluteDate, jdouble), arg0, arg1)
end

function AbsoluteDate(arg0::AbsoluteDate, arg1::jdouble, arg2::TimeScale)
    return AbsoluteDate((AbsoluteDate, jdouble, TimeScale), arg0, arg1, arg2)
end

function AbsoluteDate(arg0::Date, arg1::TimeScale)
    return AbsoluteDate((Date, TimeScale), arg0, arg1)
end

function AbsoluteDate(arg0::DateComponents, arg1::TimeComponents, arg2::TimeScale)
    return AbsoluteDate((DateComponents, TimeComponents, TimeScale), arg0, arg1, arg2)
end

function AbsoluteDate(arg0::DateComponents, arg1::TimeScale)
    return AbsoluteDate((DateComponents, TimeScale), arg0, arg1)
end

function AbsoluteDate(arg0::DateTimeComponents, arg1::TimeScale)
    return AbsoluteDate((DateTimeComponents, TimeScale), arg0, arg1)
end

function AbsoluteDate(arg0::JString, arg1::TimeScale)
    return AbsoluteDate((JString, TimeScale), arg0, arg1)
end

function AbsoluteDate(arg0::jint, arg1::Month, arg2::jint, arg3::TimeScale)
    return AbsoluteDate((jint, Month, jint, TimeScale), arg0, arg1, arg2, arg3)
end

function AbsoluteDate(arg0::jint, arg1::Month, arg2::jint, arg3::jint, arg4::jint, arg5::jdouble, arg6::TimeScale)
    return AbsoluteDate((jint, Month, jint, jint, jint, jdouble, TimeScale), arg0, arg1, arg2, arg3, arg4, arg5, arg6)
end

function AbsoluteDate(arg0::jint, arg1::jint, arg2::jint, arg3::TimeScale)
    return AbsoluteDate((jint, jint, jint, TimeScale), arg0, arg1, arg2, arg3)
end

function AbsoluteDate(arg0::jint, arg1::jint, arg2::jint, arg3::jint, arg4::jint, arg5::jdouble, arg6::TimeScale)
    return AbsoluteDate((jint, jint, jint, jint, jint, jdouble, TimeScale), arg0, arg1, arg2, arg3, arg4, arg5, arg6)
end

function compare_to(obj::AbsoluteDate, arg0::AbsoluteDate)
    return jcall(obj, "compareTo", jint, (AbsoluteDate,), arg0)
end

function create_besselian_epoch(::Type{AbsoluteDate}, arg0::jdouble)
    return jcall(AbsoluteDate, "createBesselianEpoch", AbsoluteDate, (jdouble,), arg0)
end

function create_jd_date(::Type{AbsoluteDate}, arg0::jint, arg1::jdouble, arg2::TimeScale)
    return jcall(AbsoluteDate, "createJDDate", AbsoluteDate, (jint, jdouble, TimeScale), arg0, arg1, arg2)
end

function create_julian_epoch(::Type{AbsoluteDate}, arg0::jdouble)
    return jcall(AbsoluteDate, "createJulianEpoch", AbsoluteDate, (jdouble,), arg0)
end

function create_mjd_date(::Type{AbsoluteDate}, arg0::jint, arg1::jdouble, arg2::TimeScale)
    return jcall(AbsoluteDate, "createMJDDate", AbsoluteDate, (jint, jdouble, TimeScale), arg0, arg1, arg2)
end

function duration_from(obj::AbsoluteDate, arg0::AbsoluteDate)
    return jcall(obj, "durationFrom", jdouble, (AbsoluteDate,), arg0)
end

function equals(obj::AbsoluteDate, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_class(obj::Object)
    return jcall(obj, "getClass", Class, ())
end

function get_components(obj::AbsoluteDate, arg0::TimeScale)
    return jcall(obj, "getComponents", DateTimeComponents, (TimeScale,), arg0)
end

function get_components(obj::AbsoluteDate, arg0::TimeZone)
    return jcall(obj, "getComponents", DateTimeComponents, (TimeZone,), arg0)
end

function get_components(obj::AbsoluteDate, arg0::TimeZone, arg1::TimeScale)
    return jcall(obj, "getComponents", DateTimeComponents, (TimeZone, TimeScale), arg0, arg1)
end

function get_components(obj::AbsoluteDate, arg0::jint)
    return jcall(obj, "getComponents", DateTimeComponents, (jint,), arg0)
end

function get_components(obj::AbsoluteDate, arg0::jint, arg1::TimeScale)
    return jcall(obj, "getComponents", DateTimeComponents, (jint, TimeScale), arg0, arg1)
end

function get_date(obj::AbsoluteDate)
    return jcall(obj, "getDate", AbsoluteDate, ())
end

function hash_code(obj::AbsoluteDate)
    return jcall(obj, "hashCode", jint, ())
end

function is_after(obj::AbsoluteDate, arg0::TimeStamped)
    return jcall(obj, "isAfter", jboolean, (TimeStamped,), arg0)
end

function is_after_or_equal_to(obj::AbsoluteDate, arg0::TimeStamped)
    return jcall(obj, "isAfterOrEqualTo", jboolean, (TimeStamped,), arg0)
end

function is_before(obj::AbsoluteDate, arg0::TimeStamped)
    return jcall(obj, "isBefore", jboolean, (TimeStamped,), arg0)
end

function is_before_or_equal_to(obj::AbsoluteDate, arg0::TimeStamped)
    return jcall(obj, "isBeforeOrEqualTo", jboolean, (TimeStamped,), arg0)
end

function is_between(obj::AbsoluteDate, arg0::TimeStamped, arg1::TimeStamped)
    return jcall(obj, "isBetween", jboolean, (TimeStamped, TimeStamped), arg0, arg1)
end

function is_between_or_equal_to(obj::AbsoluteDate, arg0::TimeStamped, arg1::TimeStamped)
    return jcall(obj, "isBetweenOrEqualTo", jboolean, (TimeStamped, TimeStamped), arg0, arg1)
end

function is_close_to(obj::AbsoluteDate, arg0::TimeStamped, arg1::jdouble)
    return jcall(obj, "isCloseTo", jboolean, (TimeStamped, jdouble), arg0, arg1)
end

function is_equal_to(obj::AbsoluteDate, arg0::TimeStamped)
    return jcall(obj, "isEqualTo", jboolean, (TimeStamped,), arg0)
end

function notify(obj::Object)
    return jcall(obj, "notify", void, ())
end

function notify_all(obj::Object)
    return jcall(obj, "notifyAll", void, ())
end

function offset_from(obj::AbsoluteDate, arg0::AbsoluteDate, arg1::TimeScale)
    return jcall(obj, "offsetFrom", jdouble, (AbsoluteDate, TimeScale), arg0, arg1)
end

function parse_ccsds_calendar_segmented_time_code(::Type{AbsoluteDate}, arg0::jbyte, arg1::Vector{jbyte})
    return jcall(AbsoluteDate, "parseCCSDSCalendarSegmentedTimeCode", AbsoluteDate, (jbyte, Vector{jbyte}), arg0, arg1)
end

function parse_ccsds_calendar_segmented_time_code(::Type{AbsoluteDate}, arg0::jbyte, arg1::Vector{jbyte}, arg2::TimeScale)
    return jcall(AbsoluteDate, "parseCCSDSCalendarSegmentedTimeCode", AbsoluteDate, (jbyte, Vector{jbyte}, TimeScale), arg0, arg1, arg2)
end

function parse_ccsds_day_segmented_time_code(::Type{AbsoluteDate}, arg0::jbyte, arg1::Vector{jbyte}, arg2::DateComponents)
    return jcall(AbsoluteDate, "parseCCSDSDaySegmentedTimeCode", AbsoluteDate, (jbyte, Vector{jbyte}, DateComponents), arg0, arg1, arg2)
end

function parse_ccsds_day_segmented_time_code(::Type{AbsoluteDate}, arg0::jbyte, arg1::Vector{jbyte}, arg2::DateComponents, arg3::TimeScale)
    return jcall(AbsoluteDate, "parseCCSDSDaySegmentedTimeCode", AbsoluteDate, (jbyte, Vector{jbyte}, DateComponents, TimeScale), arg0, arg1, arg2, arg3)
end

function parse_ccsds_unsegmented_time_code(::Type{AbsoluteDate}, arg0::jbyte, arg1::jbyte, arg2::Vector{jbyte}, arg3::AbsoluteDate)
    return jcall(AbsoluteDate, "parseCCSDSUnsegmentedTimeCode", AbsoluteDate, (jbyte, jbyte, Vector{jbyte}, AbsoluteDate), arg0, arg1, arg2, arg3)
end

function parse_ccsds_unsegmented_time_code(::Type{AbsoluteDate}, arg0::jbyte, arg1::jbyte, arg2::Vector{jbyte}, arg3::AbsoluteDate, arg4::AbsoluteDate)
    return jcall(AbsoluteDate, "parseCCSDSUnsegmentedTimeCode", AbsoluteDate, (jbyte, jbyte, Vector{jbyte}, AbsoluteDate, AbsoluteDate), arg0, arg1, arg2, arg3, arg4)
end

function shifted_by(obj::AbsoluteDate, arg0::jdouble)
    return jcall(obj, "shiftedBy", AbsoluteDate, (jdouble,), arg0)
end

function time_scales_offset(obj::AbsoluteDate, arg0::TimeScale, arg1::TimeScale)
    return jcall(obj, "timeScalesOffset", jdouble, (TimeScale, TimeScale), arg0, arg1)
end

function to_date(obj::AbsoluteDate, arg0::TimeScale)
    return jcall(obj, "toDate", Date, (TimeScale,), arg0)
end

function to_string(obj::AbsoluteDate)
    return jcall(obj, "toString", JString, ())
end

function to_string(obj::AbsoluteDate, arg0::TimeScale)
    return jcall(obj, "toString", JString, (TimeScale,), arg0)
end

function to_string(obj::AbsoluteDate, arg0::TimeZone)
    return jcall(obj, "toString", JString, (TimeZone,), arg0)
end

function to_string(obj::AbsoluteDate, arg0::TimeZone, arg1::TimeScale)
    return jcall(obj, "toString", JString, (TimeZone, TimeScale), arg0, arg1)
end

function to_string(obj::AbsoluteDate, arg0::jint)
    return jcall(obj, "toString", JString, (jint,), arg0)
end

function to_string(obj::AbsoluteDate, arg0::jint, arg1::TimeScale)
    return jcall(obj, "toString", JString, (jint, TimeScale), arg0, arg1)
end

function to_string_rfc3339(obj::AbsoluteDate, arg0::TimeScale)
    return jcall(obj, "toStringRfc3339", JString, (TimeScale,), arg0)
end

function wait(obj::Object)
    return jcall(obj, "wait", void, ())
end

function wait(obj::Object, arg0::jlong)
    return jcall(obj, "wait", void, (jlong,), arg0)
end

function wait(obj::Object, arg0::jlong, arg1::jint)
    return jcall(obj, "wait", void, (jlong, jint), arg0, arg1)
end

ext: jl


text:
using PDIPS

using LinearAlgebra
using Printf
using SparseArrays
using Test

@testset "Types" begin
    include("types.jl")
    include("algorithm.jl")
end


text:
# PS: Julia data structures are columnar
# type to operate on
type AD
  name::String
  value::Matrix{Float64}
  opNode::Bool
  grad::Matrix{Float64}
  gradOp::Function
  parents::Array{AD}
  level::Int64
  output::Matrix{Float64}
  outputOrder::Int64
  loss::Float64

  AD(val::Float64) = new(@getName(val), ones(1,1)*val, false, zeros(1,1), ad_constD, Array(AD,0), 1, zeros(1,1), -1, -1)
  AD(val::Int64) = new(@getName(val), ones(1,1)*val, false, zeros(1,1), ad_constD, Array(AD,0), 1, zeros(1,1), -1, -1)
  AD(val::Matrix{Float64}) = new(@getName(val), val, false, zeros(size(val)), ad_constD, Array(AD,0), 1, zeros(1,1), -1, -1)
  AD(name::String, val::Matrix{Float64}, opNode::Bool, grad::Matrix{Float64}, gradOp::Function, level::Int64) = new(name, val, opNode, grad, gradOp, Array(AD,0), level, zeros(1,1), -1, -1)
  AD(name::String, val::Matrix{Float64}, opNode::Bool, grad::Matrix{Float64}, gradOp::Function, level::Int64, output::Matrix{Float64}, outputOrder::Int64, loss::Float64) = new(name, val, opNode, grad, gradOp, Array(AD,0), level, output, outputOrder, loss)
end

# since we do not want to differentiate with respect to a constant
function ad_constD(this::AD)
  return zeros(1,1)
end


text:
local_dir = "/home/awise/Desktop/RvSpectML"
#local_dir = "/storage/work/afw5465/RvSpectML"

using Pkg
 Pkg.activate(joinpath(local_dir,"alex_sandbox/"))
 using Query, DataFrames, DataStructures, CSV, JLD, Statistics, StatsBase, Plots, JLD2, FileIO, ArgParse
 using RvSpectMLBase, EchelleInstruments, EchelleCCFs, RvSpectML, Scalpels

 include(joinpath(local_dir,"data_paths.jl"))
 include(joinpath(local_dir,"RvLineList/inputs/param.jl"))

function hasOneToOneMatch(X::A1,Y::A2) where { T1<:Real, A1<:AbstractArray{T1,1}, T2<:Real, A2<:AbstractArray{T2,1}}
   xTrue = zeros(Bool,size(X,1))
   for y in 1:length(Y)
      n_matches = 0
      match_i = -1
      for x in 1:length(X)
         if (abs(X[x]-Y[y])/X[x])<1e-5
            n_matches += 1
            match_i = x
         end
      end
      if n_matches == 1
         xTrue[match_i]=true
      end
   end
   return xTrue
end

"""
ccfs = Dict()
rvs = Dict()
times = Dict()
"""

targs = ["101501", "26965", "10700", "34411"]
norm_types = [:raw, :continuum, :blaze]
EXPRES_excalibur_onlys = [true, false]
mask_types = [:tophat, :gaussian]
mask_scale_factors = [1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0]
overlap_cutoffs = [2]
nbin=1
n=0

rms_rvs_espresso = DataFrame(targ=String[],norm_type=Symbol[],excalibur_only=Bool[],mask_type=Symbol[],mask_scale_factor=Float64[],overlap_cutoff=Int[], rv_rms=Float64[])
rms_rvs_espresso_nightly = DataFrame(targ=String[],norm_type=Symbol[],excalibur_only=Bool[],mask_type=Symbol[],mask_scale_factor=Float64[],overlap_cutoff=Int[], rv_rms=Float64[])
rms_rvs_espresso_within_night = DataFrame(targ=String[],norm_type=Symbol[],excalibur_only=Bool[],mask_type=Symbol[],mask_scale_factor=Float64[],overlap_cutoff=Int[], rv_rms=Float64[])
rms_rvs_empirical = DataFrame(targ=String[],norm_type=Symbol[],excalibur_only=Bool[],mask_type=Symbol[],mask_scale_factor=Float64[],overlap_cutoff=Int[], rv_rms=Float64[])
rms_rvs_empirical_nightly = DataFrame(targ=String[],norm_type=Symbol[],excalibur_only=Bool[],mask_type=Symbol[],mask_scale_factor=Float64[],overlap_cutoff=Int[], rv_rms=Float64[])
rms_rvs_empirical_within_night = DataFrame(targ=String[],norm_type=Symbol[],excalibur_only=Bool[],mask_type=Symbol[],mask_scale_factor=Float64[],overlap_cutoff=Int[], rv_rms=Float64[])

for norm_type in norm_types
   for EXPRES_excalibur_only in EXPRES_excalibur_onlys
      for mask_type in mask_types
         for mask_scale_factor in mask_scale_factors
            for overlap_cutoff in overlap_cutoffs
               results = load(joinpath(local_dir,"RvLineList","outputs",string("param_study_results","_norm_type=",norm_type,
               "_EXPRES_excalibur_only=",EXPRES_excalibur_only,"_mask_type=",mask_type,"_mask_scale_factor=",mask_scale_factor,
               "_overlap_cutoff=",overlap_cutoff,"_quant=",quant,"_nbin=",nbin,"_n=",n,".jld")))
               for targ in targs
                  #times = results["times"][targ]
                  #ccfs = results["ccfs"][targ]
                  rvs_espresso = results["rvs"][string("rvs_espresso_",targ)]
                  #rv_rms_espresso = std(rvs_espresso)
                  rv_rms_espresso = std(rvs_espresso .- mean(rvs_espresso))
                  rv_rms_espresso_nightly = std(bin_rvs_nightly(times=results["times"][targ],rvs=rvs_espresso .- mean(rvs_espresso)))
                  rv_rms_espresso_within_night = rms_rvs_within_night(times=results["times"][targ],rvs=rvs_espresso .- mean(rvs_espresso))
                  push!(rms_rvs_espresso,[targ,norm_type,EXPRES_excalibur_only,mask_type,mask_scale_factor,overlap_cutoff,rv_rms_espresso])
                  push!(rms_rvs_espresso_nightly,[targ,norm_type,EXPRES_excalibur_only,mask_type,mask_scale_factor,overlap_cutoff,rv_rms_espresso_nightly])
                  push!(rms_rvs_espresso_within_night,[targ,norm_type,EXPRES_excalibur_only,mask_type,mask_scale_factor,overlap_cutoff,rv_rms_espresso_within_night])
                  rvs_empirical = results["rvs"][string("rvs_empirical_",targ)]
                  #rv_rms_empirical = std(rvs_empirical)
                  rv_rms_empirical = std(rvs_empirical .- mean(rvs_empirical))
                  rv_rms_empirical_nightly = std(bin_rvs_nightly(times=results["times"][targ],rvs=rvs_empirical .- mean(rvs_empirical)))
                  rv_rms_empirical_within_night = rms_rvs_within_night(times=results["times"][targ],rvs=rvs_empirical .- mean(rvs_empirical))
                  push!(rms_rvs_empirical,[targ,norm_type,EXPRES_excalibur_only,mask_type,mask_scale_factor,overlap_cutoff,rv_rms_empirical])
                  push!(rms_rvs_empirical_nightly,[targ,norm_type,EXPRES_excalibur_only,mask_type,mask_scale_factor,overlap_cutoff,rv_rms_empirical_nightly])
                  push!(rms_rvs_empirical_within_night,[targ,norm_type,EXPRES_excalibur_only,mask_type,mask_scale_factor,overlap_cutoff,rv_rms_empirical_within_night])
               end
            end
         end
      end
   end
end
#commands used to visualize best RV params:
"""
i=1

sort(filter(row -> (row.targ == targs[i]) && (row.overlap_cutoff == 2), rms_rvs_espresso),:rv_rms)
sort(filter(row -> (row.targ == targs[i]) && (row.overlap_cutoff == 2), rms_rvs_espresso_nightly),:rv_rms)
sort(filter(row -> (row.targ == targs[i]) && (row.overlap_cutoff == 2), rms_rvs_espresso_within_night),:rv_rms)

sort(filter(row -> (row.targ == targs[i]) && (row.overlap_cutoff == 2), rms_rvs_empirical),:rv_rms)
sort(filter(row -> (row.targ == targs[i]) && (row.overlap_cutoff == 2), rms_rvs_empirical_nightly),:rv_rms)
sort(filter(row -> (row.targ == targs[i]) && (row.overlap_cutoff == 2), rms_rvs_empirical_within_night),:rv_rms)
"""
#index                     1        2        3        4
#best RMSs target:         101501,  26965,   10700,   34411
#best espresso RMSs:       3.44,    3.09,    1.94,    1.64
#espresso nightly:         3.36,    2.94,    1.19,    1.40
#espresso within-night:    0.51,    0.68,    1.32,    1.04

#best empirical RMSs:      4.14,    2.94,    1.69,    1.81
#empirical nightly:        4.08,    2.82,    0.91,    1.41
#empirical within-night:   0.60,    0.66,    1.33,    1.14

#HD101501 best ESPRESSO params: norm_type=:continuum, mask_type=:gaussian, mask_scale_factor=8.0
#     within-night, half the mask_scale_factor and norm_type=:raw

#all 3 other stars best ESPRESSO params: norm_type=:raw, mask_type=either, mask_scale_factor=1.0

#HD101501 empirical: continuum + excalibur_only is best, continum + not only excalibur is worst. gaussian > tophat, mask_scale_factor = max
#     within night, quarter the mask_scale_factor and norm_type=:raw
#HD26962 empirical: continuum + excalibur_only is best, mask_shape and mask_scale_factor have little effect
#     within night, raw is better than continuum
#HD10700 empirical: raw + excalibur_only is best, mask_shape has little effect and mask_scale_Factor = min
#     nightly, continuum now is better than raw.
#HD34411 empirical: continuum + excalibur_only is best, mask_shape has little effect and mask_scale_Factor has little effect but should be min
#     within night, raw is better than continuum

#usually continuum is better, but for HD10700 empirical and for ESPRESSO masks raw is better, need to keep testing both
#excalibur_only is best, mask type: gaussian >= tophat, mask_scale_factor =1.0
#within night, raw is better (also overall for 10700 since it is dominated by within-night), while nightly, continuum is better

#neid worst model in their grid for which we would still use the observations: air mass at 2 hrs from noon on dec 23rd, highest barometric pressure
#correlate water column (measured from spectra) with RVs to test effect

targs = ["101501", "26965", "10700", "34411"]
norm_types = [:raw, :continuum]
EXPRES_excalibur_onlys = [true]
mask_types = [:gaussian]
mask_scale_factors = [0.5,1.0,1.5]
overlap_cutoffs = [1,2,3,4,6]
quants = ["50", "70", "90"]
nbins = [1]
ns = [0]
fwtf = 2.0

rms_rvs_espresso = DataFrame(targ=String[],norm_type=Symbol[],excalibur_only=Bool[],mask_type=Symbol[],mask_scale_factor=Float64[],overlap_cutoff=Int[],quant=String[],nbin=Int[],n=Int[], rv_rms=Float64[])
rms_rvs_espresso_nightly = DataFrame(targ=String[],norm_type=Symbol[],excalibur_only=Bool[],mask_type=Symbol[],mask_scale_factor=Float64[],overlap_cutoff=Int[],quant=String[],nbin=Int[],n=Int[], rv_rms=Float64[])
rms_rvs_espresso_within_night = DataFrame(targ=String[],norm_type=Symbol[],excalibur_only=Bool[],mask_type=Symbol[],mask_scale_factor=Float64[],overlap_cutoff=Int[],quant=String[],nbin=Int[],n=Int[], rv_rms=Float64[])
rms_rvs_empirical = DataFrame(targ=String[],norm_type=Symbol[],excalibur_only=Bool[],mask_type=Symbol[],mask_scale_factor=Float64[],overlap_cutoff=Int[],quant=String[],nbin=Int[],n=Int[], rv_rms=Float64[])
rms_rvs_empirical_nightly = DataFrame(targ=String[],norm_type=Symbol[],excalibur_only=Bool[],mask_type=Symbol[],mask_scale_factor=Float64[],overlap_cutoff=Int[],quant=String[],nbin=Int[],n=Int[], rv_rms=Float64[])
rms_rvs_empirical_within_night = DataFrame(targ=String[],norm_type=Symbol[],excalibur_only=Bool[],mask_type=Symbol[],mask_scale_factor=Float64[],overlap_cutoff=Int[],quant=String[],nbin=Int[],n=Int[], rv_rms=Float64[])

for norm_type in norm_types
   for EXPRES_excalibur_only in EXPRES_excalibur_onlys
      for mask_type in mask_types
         for mask_scale_factor in mask_scale_factors
            for overlap_cutoff in overlap_cutoffs
               for quant in quants
                  for nbin in nbins
                     for n in ns
                        try
                          global results = load(joinpath(local_dir,"RvLineList","outputs",string("param_study_results","_norm_type=",norm_type,
                           "_EXPRES_excalibur_only=",EXPRES_excalibur_only,"_mask_type=",mask_type,"_mask_scale_factor=",mask_scale_factor,
                           "_overlap_cutoff=",overlap_cutoff,"_quant=",quant,"_nbin=",nbin,"_n=",n,".jld")))
                        catch ArgumentError
                           global results = NaN
                        end
                        for targ in targs
                           #if typeof(results)==Float64
                           #   rv_rms_espresso = NaN
                           #   rv_rms_empirical = NaN
                           #else
                           times = results["times"][targ]
                           ccfs_espresso = results["ccfs"][string("ccfs_espresso_",targ)]
                           v_grid_espresso = results["v_grids"][string("v_grid_espresso_",targ)]
                           pipeline_plan_espresso = results["pipeline_plans"][string("pipeline_plan_espresso_",targ)]
                           alg_fit_rv =  EchelleCCFs.MeasureRvFromCCFGaussian(frac_of_width_to_fit=fwtf)
                           #try
                           rvs_espresso = RvSpectML.calc_rvs_from_ccf_total(ccfs_espresso, pipeline_plan_espresso, v_grid=v_grid_espresso, times = times, recalc=true, bin_nightly=true, alg_fit_rv=alg_fit_rv)
                           #rv_rms_espresso = std(rvs_espresso)
                           rv_rms_espresso = std(rvs_espresso .- mean(rvs_espresso))
                           rv_rms_espresso_nightly = std(bin_rvs_nightly(times=results["times"][targ],rvs=rvs_espresso .- mean(rvs_espresso)))
                           rv_rms_espresso_within_night = rms_rvs_within_night(times=results["times"][targ],rvs=rvs_espresso .- mean(rvs_espresso))
                           #catch
                           #   rv_rms_espresso = NaN
                           #end
                           ccfs_empirical = results["ccfs"][string("ccfs_empirical_",targ)]
                           v_grid_empirical = results["v_grids"][string("v_grid_empirical_",targ)]
                           pipeline_plan_empirical = results["pipeline_plans"][string("pipeline_plan_empirical_",targ)]
                           alg_fit_rv =  EchelleCCFs.MeasureRvFromCCFGaussian(frac_of_width_to_fit=fwtf)
                           #try
                           rvs_empirical = RvSpectML.calc_rvs_from_ccf_total(ccfs_empirical, pipeline_plan_empirical, v_grid=v_grid_empirical, times = times, recalc=true, bin_nightly=true, alg_fit_rv=alg_fit_rv)
                           rv_rms_empirical = std(rvs_empirical .- mean(rvs_empirical))
                           rv_rms_empirical_nightly = std(bin_rvs_nightly(times=results["times"][targ],rvs=rvs_empirical .- mean(rvs_empirical)))
                           rv_rms_empirical_within_night = rms_rvs_within_night(times=results["times"][targ],rvs=rvs_empirical .- mean(rvs_empirical))
                           #catch
                           #   rv_rms_empirical = NaN
                           #end
                           #end
                           push!(rms_rvs_espresso,[targ,norm_type,EXPRES_excalibur_only,mask_type,mask_scale_factor,overlap_cutoff,quant,nbin,n,rv_rms_espresso])
                           push!(rms_rvs_espresso_nightly,[targ,norm_type,EXPRES_excalibur_only,mask_type,mask_scale_factor,overlap_cutoff,quant,nbin,n,rv_rms_espresso_nightly])
                           push!(rms_rvs_espresso_within_night,[targ,norm_type,EXPRES_excalibur_only,mask_type,mask_scale_factor,overlap_cutoff,quant,nbin,n,rv_rms_espresso_within_night])
                           push!(rms_rvs_empirical,[targ,norm_type,EXPRES_excalibur_only,mask_type,mask_scale_factor,overlap_cutoff,quant,nbin,n,rv_rms_empirical])
                           push!(rms_rvs_empirical_nightly,[targ,norm_type,EXPRES_excalibur_only,mask_type,mask_scale_factor,overlap_cutoff,quant,nbin,n,rv_rms_empirical_nightly])
                           push!(rms_rvs_empirical_within_night,[targ,norm_type,EXPRES_excalibur_only,mask_type,mask_scale_factor,overlap_cutoff,quant,nbin,n,rv_rms_empirical_within_night])
                        end
                     end
                  end
               end
            end
         end
      end
   end
end

#My hypothesis is that the best RVs to be in overlap_cutoff=1 and quant=90.
#Result: This is not generally true. Often overlap_cutoff=2 and/or quant=70 work better.
"""
i=1

sort(filter(row -> (row.targ == targs[i]), rms_rvs_empirical),:rv_rms)
sort(filter(row -> (row.targ == targs[i]), rms_rvs_empirical_nightly),:rv_rms)
sort(filter(row -> (row.targ == targs[i]), rms_rvs_empirical_within_night),:rv_rms)

"""

#index                     1        2        3        4
#best RMSs target:         101501,  26965,   10700,   34411
#best espresso RMSs:       3.44,    3.09,    1.94,    1.64
#espresso nightly:         3.36,    2.94,    1.19,    1.40
#espresso within-night:    0.51,    0.68,    1.32,    1.04

#round 1:
#best empirical RMSs:      4.14,    2.94,    1.69,    1.81
#empirical nightly:        4.08,    2.82,    0.91,    1.41
#empirical within-night:   0.60,    0.66,    1.33,    1.14

#round 2:
#best empirical RMSs:      4.2,     2.95,    1.60,    1.74
#empirical nightly:        3.95,    2.84,    0.86,    1.43
#empirical within-night:   0.50,    0.66,    1.33,    1.21

#HD101501 empirical: continuum, quant=50 is best; overlap_cutoff=6 helps a bit. Within-night, quant=70,overlap_cutoff=1 is best; raw helps a bit
#HD26962 empirical: continuum, quant=90, overlap_cutoff=1 is best. Within_night, same but raw is best
#HD10700 empirical: raw is best normally and within_night. overlap_cutoff = 2 to 5 normally and nightly, but 1 is okay for within_night. Quant always 90 or 70
#HD34411 empirical: overlap_cutoff=1 seems most important always. continuum is a bit better nightly, raw bit better within night

# Next, using overlap_cutoff=1 and quant=90, I compare summed empirical mask bins RV RMS to the same thing with
# 0, or 2 bins rejected
# problem: need to move the mask binning to AFTER mapping the mask onto the spectrum, so the bins are equally represented on the spectrum.

targs = ["101501", "26965", "10700", "34411"]
norm_types = [:raw, :continuum]
EXPRES_excalibur_onlys = [true]
mask_types = [:gaussian]
mask_scale_factors = [1.0, 4.0, 8.0]
overlap_cutoffs = [1, 6]
quants = ["90"]
nbins = [8]
ns = [1,2,3,4,5,6,7,8]
fwtf = 1.5

times = Dict()
ccfs = Dict()
v_grids = Dict()
pipeline_plans = Dict()

rms_rvs_empirical = DataFrame(targ=String[],norm_type=Symbol[],excalibur_only=Bool[],mask_type=Symbol[],mask_scale_factor=Float64[],overlap_cutoff=Int[],quant=String[],nbin=Int[],n=Int[], rv_rms=Float64[])

for norm_type in norm_types
   for EXPRES_excalibur_only in EXPRES_excalibur_onlys
      for mask_type in mask_types
         for mask_scale_factor in mask_scale_factors
            for overlap_cutoff in overlap_cutoffs
               for quant in quants
                  for nbin in nbins
                     for targ in targs
                        times[targ] = Dict()
                        ccfs[targ] = Dict()
                        v_grids[targ] = Dict()
                        pipeline_plans[targ] = Dict()
                        for n in ns
                           try
                              global results = load(joinpath(local_dir,"RvLineList","outputs",string("param_study_results","_norm_type=",norm_type,
                              "_EXPRES_excalibur_only=",EXPRES_excalibur_only,"_mask_type=",mask_type,"_mask_scale_factor=",mask_scale_factor,
                              "_overlap_cutoff=",overlap_cutoff,"_quant=",quant,"_nbin=",nbin,"_n=",n,".jld")))
                           catch ArgumentError
                              global results = NaN
                           end
                           if typeof(results)==Float64
                           else
                              times[targ][n] = results["times"][targ]
                              ccfs[targ][n] = results["ccfs"][string("ccfs_empirical_",targ)]
                              v_grids[targ][n] = results["v_grids"][string("v_grid_empirical_",targ)]
                              pipeline_plans[targ][n] = results["pipeline_plans"][string("pipeline_plan_empirical_",targ)]

                              alg_fit_rv =  EchelleCCFs.MeasureRvFromCCFGaussian(frac_of_width_to_fit=fwtf)
                              rvs_empirical = RvSpectML.calc_rvs_from_ccf_total(ccfs[targ][n], pipeline_plans[targ][n], v_grid=v_grids[targ][n], times = times[targ][n], recalc=true, bin_nightly=true, alg_fit_rv=alg_fit_rv)
                              rv_rms_empirical = std(rvs_empirical .- mean(rvs_empirical))
                              push!(rms_rvs_empirical,[targ,norm_type,EXPRES_excalibur_only,mask_type,mask_scale_factor,overlap_cutoff,quant,nbin,n,rv_rms_empirical])
                           end
                        end
                     end
                  end
               end
            end
         end
      end
   end
end

rms_rvs = DataFrame(targ=String[], number_of_bins_rejected=Int[],rv_rms=Float64[])
rms_rvs_nightly = DataFrame(targ=String[], number_of_bins_rejected=Int[],rv_rms=Float64[])
rms_rvs_within_night1 = DataFrame(targ=String[], number_of_bins_rejected=Int[],rv_rms=Float64[])

for targ in targs
   n_sorted = sort(filter(row -> (row.targ == targ) && (row.overlap_cutoff == 6) && (row.norm_type == :continuum) && (row.mask_scale_factor == 1.0) && (row.quant == "90"), rms_rvs_empirical),:rv_rms).n
   for n_reject in 0:(length(ns) -1 )
      sum_ccfs = ccfs[targ][n_sorted[1]]
      for i in 2:(length(ns)-n_reject)
         try
            sum_ccfs += ccfs[targ][i]
         catch end
      end
      v_grid = v_grids[targ][n_sorted[1]]
      pipeline_plan = pipeline_plans[targ][n_sorted[1]]
      alg_fit_rv =  EchelleCCFs.MeasureRvFromCCFGaussian(frac_of_width_to_fit=fwtf)
      try
         rvs = RvSpectML.calc_rvs_from_ccf_total(sum_ccfs, pipeline_plan, v_grid=v_grid, times = times[targ][n_sorted[1]], recalc=true, bin_nightly=true, alg_fit_rv=alg_fit_rv)
         rv_rms = std(rvs .- mean(rvs))
         rv_rms_nightly = std(bin_rvs_nightly(times=times[targ][1],rvs=rvs .- mean(rvs)))
         rv_rms_within_night = rms_rvs_within_night(times=times[targ][1],rvs=rvs .- mean(rvs))

         push!(rms_rvs,[targ,n_reject,rv_rms])
         push!(rms_rvs_nightly,[targ,n_reject,rv_rms_nightly])
         push!(rms_rvs_within_night1,[targ,n_reject,rv_rms_within_night])
      catch
         rv_rms = NaN
      end
   end
end

"""
i=1

sort(filter(row -> (row.targ == targs[i]), rms_rvs),:rv_rms)

"""

#best RMSs target:    101501,  26965,   10700,   34411
#best espresso RMSs:  3.43,    3.09,    1.94,    1.66
#best empirical RMSs: 4.48,    3.59,    4.03,    2.69

#HD101501 best ESPRESSO params: norm_type=:continuum, mask_type=:gaussian, mask_scale_factor=8.0

#all 3 other stars best ESPRESSO params: norm_type=:raw, mask_type=either, mask_scale_factor=1.0

outputdir = joinpath(local_dir,"RvLineList","outputs")
filelist = readdir(outputdir)
filelist = filter(fn -> occursin("orders_to_use",fn), filelist)

rv_rms = DataFrame(targ=String[], norm_type=Symbol[], orders_to_use = Vector{Int64}[],
   rv_rms=Float64[], rv_rms_nightly=Float64[], rv_rms_within_night=Float64[])

for i in 1:length(filelist)
   for targ in targs
      fn = filelist[i]
      results = load(joinpath(outputdir,fn))
      norm_type = occursin("norm_type=raw",fn) ? :raw : :continuum
      orders_to_use = results["orders_to_use"]

      rvs_empirical = results["rvs"][string("rvs_empirical_",targ)]
      #rv_rms_empirical = std(rvs_empirical)
      rv_rms_empirical = std(rvs_empirical .- mean(rvs_empirical))
      rv_rms_empirical_nightly = std(bin_rvs_nightly(times=results["times"][targ],rvs=rvs_empirical .- mean(rvs_empirical)))
      rv_rms_empirical_within_night = rms_rvs_within_night(times=results["times"][targ],rvs=rvs_empirical .- mean(rvs_empirical))
      push!(rv_rms,[targ,norm_type,orders_to_use,rv_rms_empirical,rv_rms_empirical_nightly,rv_rms_empirical_within_night])
   end
end

rv_rms_continuum = filter(row -> (row.norm_type == :continuum), rv_rms)
rv_rms_raw = filter(row -> (row.norm_type == :raw), rv_rms)

"""
i=1

sort(filter(row -> (row.targ == targs[i]) & (row.norm_type == :continuum), rv_rms), :rv_rms_nightly)

"""

# OLD RESULTS:
#index                     1        2        3        4
#best RMSs target:         101501,  26965,   10700,   34411
#best empirical RMSs:      4.2,     2.95,    1.60,    1.74
#empirical nightly:        3.95,    2.84,    0.86,    1.43
#empirical within-night:   0.50,    0.66,    1.33,    1.21

#for each order in 43:72, find its average RMS, averaged over all of runs it appeared in

order_rms = DataFrame(targ=String[], order=Int64[], count=Int64[], rv_rms=Float64[], rv_rms_nightly=Float64[], rv_rms_within_night=Float64[])
for targ in targs
   global rv_rms_targ = filter(row -> (row.targ == targ), rv_rms_continuum)
   for o in 43:72
      for i in 1:48
         if o in rv_rms_targ[i,:orders_to_use]
            #if length(filter(row -> (row.order == o) && (row.targ == targ), order_rms)) >= 1
            #   order_rms[]
            push!(order_rms,[targ,o,1,rv_rms_targ[i,:rv_rms],rv_rms_targ[i,:rv_rms_nightly],rv_rms_targ[i,:rv_rms_within_night]])
         end
      end
   end
end
order_rms_avg_targ = DataFrame(targ=String[], order=Int64[], count=Int64[], rv_rms=Float64[], rv_rms_nightly=Float64[], rv_rms_within_night=Float64[])
for targ in targs
   for o in 43:72
      df_local = filter(row -> (row.targ == targ) && (row.order==o), order_rms)
      push!(order_rms_avg_targ,[targ,o,sum(df_local[:,:count]),mean(df_local[:,:rv_rms]),mean(df_local[:,:rv_rms_nightly]),mean(df_local[:,:rv_rms_within_night])])
   end
end

order_rms_avg = DataFrame(targ=String[], order=Int64[], count=Int64[], rv_rms=Float64[], rv_rms_nightly=Float64[], rv_rms_within_night=Float64[])
for o in 43:72
   df_local = filter(row ->  (row.order==o), order_rms)
   push!(order_rms_avg,["all",o,sum(df_local[:,:count]),mean(df_local[:,:rv_rms]),mean(df_local[:,:rv_rms_nightly]),mean(df_local[:,:rv_rms_within_night])])
end

"""
i=1

println(sort(filter(row -> (row.targ == targs[i]), order_rms_avg_targ), :rv_rms_within_night))

i=2

println(sort(filter(row -> (row.targ == targs[i]), order_rms_avg_targ), :rv_rms_within_night))

i=3

println(sort(filter(row -> (row.targ == targs[i]), order_rms_avg_targ), :rv_rms_within_night))

i=4

println(sort(filter(row -> (row.targ == targs[i]), order_rms_avg_targ), :rv_rms_within_night))

"""

"""
i=1

println(sort( order_rms_avg, :rv_rms_within_night))

"""

max_stars_repo_name: alexander-wise/RvLineList

max_issues_repo_name: alexander-wise/RvLineList

max_forks_repo_name: alexander-wise/RvLineList


text:
#=
# Linear elasticity
=#

### If the Finch package has already been added, use this line #########
using Finch # Note: to add the package, first do: ]add "https://github.com/paralab/Finch.git"

### If not, use these four lines (working from the examples directory) ###
# if !@isdefined(Finch)
#     include("../Finch.jl");
#     using .Finch
# end
##########################################################################

init_finch("elasticity");

# Optionally generate a log
useLog("elasticitylog", level=3)

n = [20,4,4]; # number of elements in x,y,z
interval = [0,1,0,0.2,0,0.2]; # domain bounds

# Set up the configuration (order doesn't matter)
domain(3, grid=UNSTRUCTURED)
functionSpace(order=2)      # basis polynomial order

# Specify the problem
mesh(HEXMESH, elsperdim=n, bids=4, interval=interval)

u = variable("u", type=VECTOR)          # make a vector variable with symbol u
testSymbol("v", type=VECTOR)            # sets the symbol for a test function

boundary(u, 1, DIRICHLET, [0,0,0]) # x=0
boundary(u, 2, NEUMANN, [0,0,0])   # elsewhere
boundary(u, 3, NEUMANN, [0,0,0])
boundary(u, 4, NEUMANN, [0,0,0])

# Write the weak form
# coefficient("mu", "x>0.5 ? 0.2 : 10") # discontinuous mu
coefficient("mu", 1) # constant mu
coefficient("lambda", 1.25)
coefficient("f", [0,0,-10], type=VECTOR)
weakForm(u, "inner( (lambda * div(u) .* [1 0 0; 0 1 0; 0 0 1] + mu .* (grad(u) + transpose(grad(u)))), grad(v)) - dot(f,v)")

exportCode("elasticitycode")

println("solving")
t = @elapsed(solve(u));
println("solved")

println("time: "*string(t))
# Dump things to the log if desired
# log_dump_config();
# log_dump_prob();

# Write result to vtk file
#output_values(u, "elasticity", format="vtk");

finalize_finch() # Finish writing and close any files

println("u_z at end "*string(maximum(u.values[3,Finch.grid_data.bdry[2][3]])));

max_stars_repo_name: paralab/Finch

max_issues_repo_name: paralab/Finch

max_forks_repo_name: paralab/Finch


text:
# These functions are based on C functions in the CSparse library by Tim Davis.
# These are pure Julia implementations, and do not link to the CSparse library.
# CSparse can be downloaded from http://www.cise.ufl.edu/research/sparse/CSparse/CSparse.tar.gz
# CSparse is Copyright (c) 2006-2007, Timothy A. Davis and released under
# Lesser GNU Public License, version 2.1 or later.  A copy of the license can be
# downloaded from http://www.gnu.org/licenses/lgpl-2.1.html

# Because these functions are based on code covered by LGPL-2.1+ the same license
# must apply to the code in this file which is
# Copyright (c) 2013-2014 Viral Shah, Douglas Bates and other contributors

# Based on Direct Methods for Sparse Linear Systems, T. A. Davis, SIAM, Philadelphia, Sept. 2006.
# Section 2.4: Triplet form
# http://www.cise.ufl.edu/research/sparse/CSparse/
function sparse{Tv,Ti<:Integer}(I::AbstractVector{Ti}, J::AbstractVector{Ti},
                                V::AbstractVector{Tv},
                                nrow::Integer, ncol::Integer, combine::Function)

    if length(I) == 0; return spzeros(eltype(V),nrow,ncol); end

    # Work array
    Wj = Array(Ti, max(nrow,ncol)+1)

    # Allocate sparse matrix data structure
    # Count entries in each row
    Rnz = zeros(Ti, nrow+1)
    Rnz[1] = 1
    nz = 0
    for k=1:length(I)
        if V[k] != 0
            Rnz[I[k]+1] += 1
            nz += 1
        end
    end
    Rp = cumsum(Rnz)
    Ri = Array(Ti, nz)
    Rx = Array(Tv, nz)

    # Construct row form
    # place triplet (i,j,x) in column i of R
    # Use work array for temporary row pointers
    for i=1:nrow; Wj[i] = Rp[i]; end

    for k=1:length(I)
        ind = I[k]
        p = Wj[ind]
        Vk = V[k]
        if Vk != 0
            Wj[ind] += 1
            Rx[p] = Vk
            Ri[p] = J[k]
        end
    end

    # Reset work array for use in counting duplicates
    for j=1:ncol; Wj[j] = 0; end

    # Sum up duplicates and squeeze
    anz = 0
    for i=1:nrow
        p1 = Rp[i]
        p2 = Rp[i+1] - 1
        pdest = p1

        for p = p1:p2
            j = Ri[p]
            pj = Wj[j]
            if pj >= p1
                Rx[pj] = combine (Rx[pj], Rx[p])
            else
                Wj[j] = pdest
                if pdest != p
                    Ri[pdest] = j
                    Rx[pdest] = Rx[p]
                end
                pdest += 1
            end
        end

        Rnz[i] = pdest - p1
        anz += (pdest - p1)
    end

    # Transpose from row format to get the CSC format
    RiT = Array(Ti, anz)
    RxT = Array(Tv, anz)

    # Reset work array to build the final colptr
    Wj[1] = 1
    for i=2:(ncol+1); Wj[i] = 0; end
    for j = 1:nrow
        p1 = Rp[j]
        p2 = p1 + Rnz[j] - 1
        for p = p1:p2
            Wj[Ri[p]+1] += 1
        end
    end
    RpT = cumsum(Wj[1:(ncol+1)])

    # Transpose
    for i=1:length(RpT); Wj[i] = RpT[i]; end
    for j = 1:nrow
        p1 = Rp[j]
        p2 = p1 + Rnz[j] - 1
        for p = p1:p2
            ind = Ri[p]
            q = Wj[ind]
            Wj[ind] += 1
            RiT[q] = j
            RxT[q] = Rx[p]
        end
    end

    return SparseMatrixCSC(nrow, ncol, RpT, RiT, RxT)
end

## Transpose

# Based on Direct Methods for Sparse Linear Systems, T. A. Davis, SIAM, Philadelphia, Sept. 2006.
# Section 2.5: Transpose
# http://www.cise.ufl.edu/research/sparse/CSparse/

# NOTE: When calling transpose!(S,T), the colptr in the result matrix T must be set up
# by counting the nonzeros in every row of S.
function transpose!{Tv,Ti}(S::SparseMatrixCSC{Tv,Ti}, T::SparseMatrixCSC{Tv,Ti})
    (mT, nT) = size(T)
    colptr_T = T.colptr
    rowval_T = T.rowval
    nzval_T = T.nzval

    nnzS = nfilled(S)
    colptr_S = S.colptr
    rowval_S = S.rowval
    nzval_S = S.nzval

    w = copy(colptr_T)
    @inbounds for j = 1:mT, p = colptr_S[j]:(colptr_S[j+1]-1)
        ind = rowval_S[p]
        q = w[ind]
        w[ind] += 1
        rowval_T[q] = j
        nzval_T[q] = nzval_S[p]
    end

    return SparseMatrixCSC(mT, nT, colptr_T, rowval_T, nzval_T)
end

function transpose{Tv,Ti}(S::SparseMatrixCSC{Tv,Ti})
    (nT, mT) = size(S)
    nnzS = nfilled(S)
    rowval_S = S.rowval

    rowval_T = Array(Ti, nnzS)
    nzval_T = Array(Tv, nnzS)

    colptr_T = zeros(Ti, nT+1)
    colptr_T[1] = 1
    @inbounds for i=1:nfilled(S)
        colptr_T[rowval_S[i]+1] += 1
    end
    colptr_T = cumsum(colptr_T)

    T = SparseMatrixCSC(mT, nT, colptr_T, rowval_T, nzval_T)
    return transpose!(S, T)
end

# NOTE: When calling ctranspose!(S,T), the colptr in the result matrix T must be set up
# by counting the nonzeros in every row of S.
function ctranspose!{Tv,Ti}(S::SparseMatrixCSC{Tv,Ti}, T::SparseMatrixCSC{Tv,Ti})
    (mT, nT) = size(T)
    colptr_T = T.colptr
    rowval_T = T.rowval
    nzval_T = T.nzval

    nnzS = nfilled(S)
    colptr_S = S.colptr
    rowval_S = S.rowval
    nzval_S = S.nzval

    w = copy(colptr_T)
    @inbounds for j = 1:mT, p = colptr_S[j]:(colptr_S[j+1]-1)
        ind = rowval_S[p]
        q = w[ind]
        w[ind] += 1
        rowval_T[q] = j
        nzval_T[q] = conj(nzval_S[p])
    end

    return SparseMatrixCSC(mT, nT, colptr_T, rowval_T, nzval_T)
end

function ctranspose{Tv,Ti}(S::SparseMatrixCSC{Tv,Ti})
    (nT, mT) = size(S)
    nnzS = nfilled(S)
    rowval_S = S.rowval

    rowval_T = Array(Ti, nnzS)
    nzval_T = Array(Tv, nnzS)

    colptr_T = zeros(Ti, nT+1)
    colptr_T[1] = 1
    @inbounds for i=1:nfilled(S)
        colptr_T[rowval_S[i]+1] += 1
    end
    colptr_T = cumsum(colptr_T)

    T = SparseMatrixCSC(mT, nT, colptr_T, rowval_T, nzval_T)
    return ctranspose!(S, T)
end

# Compute the elimination tree of A using triu(A) returning the parent vector.
# A root node is indicated by 0. This tree may actually be a forest in that
# there may be more than one root, indicating complete separability.
# A trivial example is speye(n, n) in which every node is a root.
function etree{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}, postorder::Bool)
    m,n = size(A)
    Ap = A.colptr
    Ai = A.rowval
    parent = zeros(Ti, n)
    ancestor = zeros(Ti, n)
    for k in 1:n, p in Ap[k]:(Ap[k+1] - 1)
        i = Ai[p]
        while i != 0 && i < k
            inext = ancestor[i] # inext = ancestor of i
            ancestor[i] = k     # path compression
            if (inext == 0) parent[i] = k end # no anc., parent is k
            i = inext
        end
    end
    if !postorder return parent end
    head = zeros(Ti,n)                   # empty linked lists
    next = zeros(Ti,n)
    for j in n:-1:1                      # traverse in reverse order
        if (parent[j] == 0); continue; end # j is a root
        next[j] = head[parent[j]]        # add j to list of its parent
        head[parent[j]] = j
    end
    stack = Ti[]
    sizehint(stack, n)
    post = zeros(Ti,n)
    k = 1
    for j in 1:n
        if (parent[j] != 0) continue end # skip j if it is not a root
        push!(stack, j)                  # place j on the stack
        while (length(stack) > 0)        # while (stack is not empty)
            p = stack[end]               # p = top of stack
            i = head[p]                  # i = youngest child of p
            if (i == 0)
                pop!(stack)
                post[k] = p       # node p is the kth postordered node
                k += 1
            else
            head[p] = next[i]           # remove i from children of p
            push!(stack, i)
            end
        end
    end
    parent, post
end

etree(A::SparseMatrixCSC) = etree(A, false)

# find nonzero pattern of Cholesky L[k,1:k-1] using etree and triu(A[:,k])
# based on cs_ereach p. 43, "Direct Methods for Sparse Linear Systems"
function ereach{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}, k::Integer, parent::Vector{Ti})
    m,n = size(A); Ap = A.colptr; Ai = A.rowval
    s = Ti[]; sizehint(s, n)            # to be used as a stack
    visited = falses(n)
    visited[k] = true
    for p in Ap[k]:(Ap[k+1] - 1)
        i = Ai[p]                # A[i,k] is nonzero
        if i > k continue end    # only use upper triangular part of A
        while !visited[i]        # traverse up etree
            push!(s,i)           # L[k,i] is nonzero
            visited[i] = true
            i = parent[i]
        end
    end
    s
end

# based on cs_permute p. 21, "Direct Methods for Sparse Linear Systems"
function csc_permute{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}, pinv::Vector{Ti}, q::Vector{Ti})
    m,n = size(A); Ap = A.colptr; Ai = A.rowval; Ax = A.nzval
    if length(pinv) != m || length(q) !=  n
        error("dimension mismatch, size(A) = $(size(A)), length(pinv) = $(length(pinv)) and length(q) = $(length(q))")
    end
    if !isperm(pinv) || !isperm(q) error("both pinv and q must be permutations") end
    C = copy(A); Cp = C.colptr; Ci = C.rowval; Cx = C.nzval
    nz = zero(Ti)
    for k in 1:n
        Cp[k] = nz
        j = q[k]
        for t = Ap[j]:(Ap[j+1]-1)
            Cx[nz] = Ax[t]
            Ci[nz] = pinv[Ai[t]]
            nz += one(Ti)
        end
    end
    Cp[n + 1] = nz
    (C.').'                    # double transpose to order the columns
end

# based on cs_symperm p. 21, "Direct Methods for Sparse Linear Systems"
# form A[p,p] for a symmetric A stored in the upper triangle
function symperm{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}, pinv::Vector{Ti})
    m,n = size(A); Ap = A.colptr; Ai = A.rowval; Ax = A.nzval
    isperm(pinv) || error("perm must be a permutation")
    m == n == length(pinv) || error("dimension mismatch")
    C = copy(A); Cp = C.colptr; Ci = C.rowval; Cx = C.nzval
    w = zeros(Ti,n)
    for j in 1:n                   # count entries in each column of C
        j2 = pinv[j]
        for p in Ap[j]:(Ap[j+1]-1)
            (i = Ai[p]) > j || (w[max(pinv[i],j2)] += one(Ti))
        end
    end
    Cp[:] = cumsum(vcat(one(Ti),w))
    copy!(w,Cp[1:n])          # needed to be consistent with cs_cumsum
    for j in 1:n
        j2 = pinv[j]
        for p = Ap[j]:(Ap[j+1]-1)
            (i = Ai[p]) > j && continue
            i2 = pinv[i]
            ind = max(i2,j2)
            Ci[q = w[ind]] = min(i2,j2)
            w[ind] += 1
            Cx[q] = Ax[p]
        end
    end
    (C.').'                    # double transpose to order the columns
end

# Based on Direct Methods for Sparse Linear Systems, T. A. Davis, SIAM, Philadelphia, Sept. 2006.
# Section 2.7: Removing entries from a matrix
# http://www.cise.ufl.edu/research/sparse/CSparse/
function fkeep!{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}, f, other)
    nzorig = nfilled(A)
    nz = 1
    for j = 1:A.n
        p = A.colptr[j]                 # record current position
        A.colptr[j] = nz                # set new position
        while p < A.colptr[j+1]
            if f(A.rowval[p], j, A.nzval[p], other)
                A.nzval[nz] = A.nzval[p]
                A.rowval[nz] = A.rowval[p]
                nz += 1
            end
            p += 1
        end
    end
    A.colptr[A.n + 1] = nz
    nz -= 1
    if nz < nzorig
        resize!(A.nzval, nz)
        resize!(A.rowval, nz)
    end
    A
end

droptol!(A::SparseMatrixCSC, tol) = fkeep!(A, (i,j,x,other)->abs(x)>other, tol)
dropzeros!(A::SparseMatrixCSC) = fkeep!(A, (i,j,x,other)->x!=0, None)
triu!(A::SparseMatrixCSC) = fkeep!(A, (i,j,x,other)->(j>=i), None)
triu(A::SparseMatrixCSC) = triu!(copy(A))
tril!(A::SparseMatrixCSC) = fkeep!(A, (i,j,x,other)->(i>=j), None)
tril(A::SparseMatrixCSC) = tril!(copy(A))

max_stars_repo_name: tobig/julia

max_issues_repo_name: tobig/julia

max_forks_repo_name: tobig/julia