text:
type BlackKarasinski{TermStructureConsistentModelType, T <: TermStructure} <: OneFactorModel{TermStructureConsistentModelType}
  modT::TermStructureConsistentModelType
  a::ConstantParameter
  sigma::ConstantParameter
  ts::T
  privateConstraint::PrivateConstraint
  common::ModelCommon
end

function BlackKarasinski(ts::TermStructure, a::Float64 = 0.1, sigma = 0.1)
  a_const = ConstantParameter([a], PositiveConstraint())
  sigma_const = ConstantParameter([sigma], PositiveConstraint())

  privateConstraint = PrivateConstraint(ConstantParameter[a_const, sigma_const])

  return BlackKarasinski(TermStructureConsistentModelType(), a_const, sigma_const, ts, privateConstraint, ModelCommon()) # ShortRateModelCommon())
end

generate_arguments!(m::BlackKarasinski) = m # do nothing

type BlackKarasinskiDynamics{P <: Parameter} <: ShortRateDynamics
  process::OrnsteinUhlenbeckProcess
  fitting::P
  a::Float64
  sigma::Float64
end

# Dynamics #
BlackKarasinskiDynamics{P <: Parameter}(fitting::P, a::Float64, sigma::Float64) = BlackKarasinskiDynamics(OrnsteinUhlenbeckProcess(a, sigma), fitting, a, sigma)

short_rate(dynamic::BlackKarasinskiDynamics, t::Float64, x::Float64) = exp(x + dynamic.fitting(t))

immutable BlackKarasinskiHelper
  sz::Int
  xMin::Float64
  dx::Float64
  dt::Float64
  discountBondPrice::Float64
  statePrices::Vector{Float64}
end

function BlackKarasinskiHelper(i::Int, xMin::Float64, dx::Float64, discountBond::Float64, _tree::OneFactorShortRateTree)
  sz = get_size(_tree, i)
  statePrices = get_state_prices!(_tree, i)
  dt = _tree.tg.dt[i]

  return BlackKarasinskiHelper(sz, xMin, dx, dt, discountBond, statePrices)
end

function operator(bkh::BlackKarasinskiHelper)
  function _inner(theta::Float64)
    val = bkh.discountBondPrice
    x = bkh.xMin

    for j = 1:bkh.sz
      disc = exp(-exp(theta + x) * bkh.dt)
      val -= bkh.statePrices[j] * disc
      x += bkh.dx
    end

    return val
  end

  return _inner
end

# tree methods #
function tree(model::BlackKarasinski, grid::TimeGrid)
  phi = TermStructureFittingParameter(model.ts)

  numericDynamics = BlackKarasinskiDynamics(phi, get_a(model), get_sigma(model))
  trinomial = TrinomialTree(numericDynamics.process, grid)
  numericTree = OneFactorShortRateTree{BlackKarasinskiDynamics}(trinomial, numericDynamics, grid)

  reset_param_impl!(phi)

  val = 1.0
  vMin = -50.0
  vMax = 50.0

  @simd for i = 1:length(grid.times) - 1
    @inbounds discountBond = discount(model.ts, grid.times[i + 1])
    xMin = get_underlying(trinomial, i, 1)
    @inbounds dx = trinomial.dx[i]

    solverHelper = BlackKarasinskiHelper(i, xMin, dx, discountBond, numericTree)
    slvr = BrentSolver(1000) # max evals = 1000
    val = solve(slvr, operator(solverHelper), 1e-7, val, vMin, vMax)

    @inbounds set_params!(phi, grid.times[i], val)
  end

  return numericTree
end

ext: jl


text:
"""
    flatten(x::AbstractArray)

Reshape arbitrarly-shaped input into a matrix-shaped output
preserving the last dimension size.
Equivalent to `reshape(x, :, size(x)[end])`.
"""
function flatten(x::AbstractArray)
  return reshape(x, :, size(x)[end])
end

"""
    normalise(x; dims=ndims(x), ϵ=1e-5)

Normalise `x` to mean 0 and standard deviation 1 across the dimension(s) given by `dims`.
Per default, `dims` is the last dimension.
`ϵ` is a small additive factor added to the denominator for numerical stability.
"""
function normalise(x::AbstractArray; dims=ndims(x), ϵ=ofeltype(x, 1e-5))
  μ = mean(x, dims=dims)
    #   σ = std(x, dims=dims, mean=μ, corrected=false) # use this when Zygote#478 gets merged
  σ = std(x, dims=dims, corrected=false)
  return (x .- μ) ./ (σ .+ ϵ)
end


text:
using Javis

function ground(args...)
    background("white") # canvas background
    sethue("black") # pen color
end

function object(;p=O, color="black")
    sethue(color)
    s = ngon(p, 25, 4)
    poly(s, :fill, close=true)
    return p
end

myvideo = Video(500, 500)

actions = [Action((args...)-> object(;color="green"),
           Translation(Point(-250, 0), Point(250, 0))),
]

javis(
    myvideo,
    [
        BackgroundAction(1:500, ground), actions...
    ],
    pathname="particle.gif"
)


text:
using Oceananigans.Fields: cpudata

"""
    correct_field_size(arch, grid, FieldType, Tx, Ty, Tz)

Test that the field initialized by the FieldType constructor on `arch` and `grid`
has size `(Tx, Ty, Tz)`.
"""
correct_field_size(a, g, FieldType, Tx, Ty, Tz) = size(parent(FieldType(a, g))) == (Tx, Ty, Tz)

"""
    correct_reduced_field_size(loc, arch, grid, dims, Tx, Ty, Tz)

Test that the ReducedField at `loc`ation on `arch`itecture and `grid`
and reduced along `dims` has size `(Tx, Ty, Tz)`.
"""
correct_reduced_field_size(loc, arch, grid, dims, Tx, Ty, Tz) =
    size(parent(ReducedField(loc, arch, grid; dims=dims))) == (Tx, Ty, Tz)

"""
     correct_field_value_was_set(N, L, ftf, val)

Test that the field initialized by the field type function `ftf` on the grid g
can be correctly filled with the value `val` using the `set!(f::AbstractField, v)`
function.
"""
function correct_field_value_was_set(arch, grid, FieldType, val::Number)
    f = FieldType(arch, grid)
    set!(f, val)
    CUDA.@allowscalar return interior(f) ≈ val * ones(size(f))
end

function correct_reduced_field_value_was_set(arch, grid, loc, dims, val::Number)
    f = ReducedField(loc, arch, grid; dims=dims)
    set!(f, val)
    CUDA.@allowscalar return interior(f) ≈ val * ones(size(f))
end

function run_field_interpolation_tests(arch, FT)

    grid = RegularCartesianGrid(size=(4, 5, 7), x=(0, 1), y=(-π, π), z=(-5.3, 2.7))

    velocities = VelocityFields(arch, grid)
    tracers = TracerFields((:c,), arch, grid)

    (u, v, w), c = velocities, tracers.c

    # Choose a trilinear function so trilinear interpolation can return values that
    # are exactly correct.
    f(x, y, z) = exp(-1) + 3x - y/7 + z + 2x*y - 3x*z + 4y*z - 5x*y*z

    # Maximum expected rounding error is the unit in last place of the maximum value
    # of f over the domain of the grid.
    ε_max = f.(nodes((Face, Face, Face), grid, reshape=true)...) |> maximum |> eps

    set!(u, f)
    set!(v, f)
    set!(w, f)
    set!(c, f)

    # Check that interpolating to the field's own grid points returns
    # the same value as the field itself.

    ℑu = interpolate.(Ref(u), nodes(u, reshape=true)...)
    ℑv = interpolate.(Ref(v), nodes(v, reshape=true)...)
    ℑw = interpolate.(Ref(w), nodes(w, reshape=true)...)
    ℑc = interpolate.(Ref(c), nodes(c, reshape=true)...)

    @test all(isapprox.(ℑu, interior(u), atol=ε_max))
    @test all(isapprox.(ℑv, interior(v), atol=ε_max))
    @test all(isapprox.(ℑw, interior(w), atol=ε_max))
    @test all(isapprox.(ℑc, interior(c), atol=ε_max))

    # Check that interpolating between grid points works as expected.

    xs = reshape([0.3, 0.55, 0.73], (3, 1, 1))
    ys = reshape([-π/6, 0, 1+1e-7], (1, 3, 1))
    zs = reshape([-1.3, 1.23, 2.1], (1, 1, 3))

    F = f.(xs, ys, zs)

    ℑu = interpolate.(Ref(u), xs, ys, zs)
    ℑv = interpolate.(Ref(v), xs, ys, zs)
    ℑw = interpolate.(Ref(w), xs, ys, zs)
    ℑc = interpolate.(Ref(c), xs, ys, zs)

    @test all(isapprox.(ℑu, F, atol=ε_max))
    @test all(isapprox.(ℑv, F, atol=ε_max))
    @test all(isapprox.(ℑw, F, atol=ε_max))
    @test all(isapprox.(ℑc, F, atol=ε_max))

    return nothing
end

@testset "Fields" begin
    @info "Testing fields..."

    N = (4, 6, 8)
    L = (2π, 3π, 5π)
    H = (1, 1, 1)

    @testset "Field initialization" begin
        @info "  Testing field initialization..."
        for arch in archs, FT in float_types
            grid = RegularCartesianGrid(FT, size=N, extent=L, halo=H, topology=(Periodic, Periodic, Periodic))
            @test correct_field_size(arch, grid, CellField,  N[1] + 2 * H[1], N[2] + 2 * H[2], N[3] + 2 * H[3])
            @test correct_field_size(arch, grid, XFaceField, N[1] + 2 * H[1], N[2] + 2 * H[2], N[3] + 2 * H[3])
            @test correct_field_size(arch, grid, YFaceField, N[1] + 2 * H[1], N[2] + 2 * H[2], N[3] + 2 * H[3])
            @test correct_field_size(arch, grid, ZFaceField, N[1] + 2 * H[1], N[2] + 2 * H[2], N[3] + 2 * H[3])

            grid = RegularCartesianGrid(FT, size=N, extent=L, halo=H, topology=(Periodic, Periodic, Bounded))
            @test correct_field_size(arch, grid, CellField,  N[1] + 2 * H[1], N[2] + 2 * H[2], N[3] + 2 * H[3])
            @test correct_field_size(arch, grid, XFaceField, N[1] + 2 * H[1], N[2] + 2 * H[2], N[3] + 2 * H[3])
            @test correct_field_size(arch, grid, YFaceField, N[1] + 2 * H[1], N[2] + 2 * H[2], N[3] + 2 * H[3])
            @test correct_field_size(arch, grid, ZFaceField, N[1] + 2 * H[1], N[2] + 2 * H[2], N[3] + 2 * H[3] + 1)

            grid = RegularCartesianGrid(FT, size=N, extent=L, halo=H, topology=(Periodic, Bounded, Bounded))
            @test correct_field_size(arch, grid, CellField,  N[1] + 2 * H[1], N[2] + 2 * H[2], N[3] + 2 * H[3])
            @test correct_field_size(arch, grid, XFaceField, N[1] + 2 * H[1], N[2] + 2 * H[2], N[3] + 2 * H[3])
            @test correct_field_size(arch, grid, YFaceField, N[1] + 2 * H[1], N[2] + 1 + 2 * H[2], N[3] + 2 * H[3])
            @test correct_field_size(arch, grid, ZFaceField, N[1] + 2 * H[1], N[2] + 2 * H[2], N[3] + 1 + 2 * H[3])

            grid = RegularCartesianGrid(FT, size=N, extent=L, halo=H, topology=(Bounded, Bounded, Bounded))
            @test correct_field_size(arch, grid, CellField,  N[1] + 2 * H[1], N[2] + 2 * H[2], N[3] + 2 * H[3])
            @test correct_field_size(arch, grid, XFaceField, N[1] + 1 + 2 * H[1], N[2] + 2 * H[2], N[3] + 2 * H[3])
            @test correct_field_size(arch, grid, YFaceField, N[1] + 2 * H[1], N[2] + 1 + 2 * H[2], N[3] + 2 * H[3])
            @test correct_field_size(arch, grid, ZFaceField, N[1] + 2 * H[1], N[2] + 2 * H[2], N[3] + 1 + 2 * H[3])

            @test correct_reduced_field_size((Cell, Cell, Cell), arch, grid, 1,         1,               N[2] + 2 * H[2],     N[3] + 2 * H[3])
            @test correct_reduced_field_size((Face, Cell, Cell), arch, grid, 1,         1,               N[2] + 2 * H[2],     N[3] + 2 * H[3])
            @test correct_reduced_field_size((Cell, Face, Cell), arch, grid, 1,         1,               N[2] + 2 * H[2] + 1, N[3] + 2 * H[3])
            @test correct_reduced_field_size((Cell, Face, Face), arch, grid, 1,         1,               N[2] + 2 * H[2] + 1, N[3] + 2 * H[3] + 1)
            @test correct_reduced_field_size((Cell, Cell, Cell), arch, grid, 2,         N[1] + 2 * H[1], 1,                   N[3] + 2 * H[3])
            @test correct_reduced_field_size((Cell, Cell, Cell), arch, grid, 2,         N[1] + 2 * H[1], 1,                   N[3] + 2 * H[3])
            @test correct_reduced_field_size((Cell, Cell, Cell), arch, grid, 3,         N[1] + 2 * H[1], N[2] + 2 * H[2],     1)
            @test correct_reduced_field_size((Cell, Cell, Cell), arch, grid, (1, 2),    1,               1,                   N[3] + 2 * H[3])
            @test correct_reduced_field_size((Cell, Cell, Cell), arch, grid, (2, 3),    N[1] + 2 * H[1], 1,                   1)
            @test correct_reduced_field_size((Cell, Cell, Cell), arch, grid, (1, 2, 3), 1,               1,                   1)
        end
    end

    FieldTypes = (CellField, XFaceField, YFaceField, ZFaceField)
    reduced_dims = (1, 2, 3, (1, 2), (2, 3), (1, 3), (1, 2, 3))

    int_vals = Any[0, Int8(-1), Int16(2), Int32(-3), Int64(4), Int128(-5)]
    uint_vals = Any[6, UInt8(7), UInt16(8), UInt32(9), UInt64(10), UInt128(11)]
    float_vals = Any[0.0, -0.0, 6e-34, 1.0f10]
    rational_vals = Any[1//11, -23//7]
    other_vals = Any[π]
    vals = vcat(int_vals, uint_vals, float_vals, rational_vals, other_vals)

    @testset "Setting fields" begin
        @info "  Testing field setting..."

        for arch in archs, FT in float_types
	        ArrayType = array_type(arch)
            grid = RegularCartesianGrid(FT, size=N, extent=L, topology=(Periodic, Periodic, Bounded))

            for FieldType in FieldTypes, val in vals
                @test correct_field_value_was_set(arch, grid, FieldType, val)
            end

            for dims in reduced_dims, val in vals
                @test correct_reduced_field_value_was_set(arch, grid, (Cell, Cell, Cell), dims, val)
            end

            for FieldType in FieldTypes
                field = FieldType(arch, grid)
                sz = size(field)
                A = rand(FT, sz...) |> ArrayType
                set!(field, A)
                @test field.data[2, 4, 6] == A[2, 4, 6]
            end

            for dims in reduced_dims
                reduced_field = ReducedField((Cell, Cell, Cell), arch, grid, dims=dims)
                sz = size(reduced_field)
                A = rand(FT, sz...) |> ArrayType
                set!(reduced_field, A)
                @test reduced_field.data[1, 1, 1] == A[1, 1, 1]
            end
        end
    end

    @testset "Field utils" begin
        @info "  Testing field utils..."

        @test Fields.has_velocities(()) == false
        @test Fields.has_velocities((:u,)) == false
        @test Fields.has_velocities((:u, :v)) == false
        @test Fields.has_velocities((:u, :v, :w)) == true

		grid = RegularCartesianGrid(size=(4, 6, 8), extent=(1, 1, 1))
		ϕ = CellField(CPU(), grid)
		@test cpudata(ϕ).parent isa Array

		@hascuda begin
			ϕ = CellField(GPU(), grid)
			@test cpudata(ϕ).parent isa Array
		end
    end

    @testset "Field interpolation" begin
        @info "  Testing field interpolation..."

        for arch in archs, FT in float_types
            run_field_interpolation_tests(arch, FT)
        end
    end
end


text:
using Forge.DSL
using Colors

const MINIMUM = 1.0
const MAXIMUM = 20.
const STEP    = 2.0
const NELEMS  = (MAXIMUM-MINIMUM+1)/STEP
const DPOINTS = [5, 5, 5, 15, 15, 5, 15, 15]

points = []
dirs = []
for j in 0:NELEMS-1, i in 0:NELEMS-1
    y = MINIMUM + j*STEP
    x = MINIMUM + i*STEP
    push!(points, x, y)
    push!(dirs, sin(2pi*x/10), sin(2pi*y/10))
end

window = create_window(600, 480, "Vector Field Demo")
show_window(window)

chart = create_chart_2d()
set_axes_limits(chart, MINIMUM-1.0, MAXIMUM, MINIMUM-1.0, MAXIMUM)
set_axes_titles(chart, "x-axis", "y-axis")

plot = create_plot(chart, 4, Forge.FG_PLOT_SCATTER, Forge.FG_MARKER_CIRCLE)
set_color(plot, RGB(0.9, 0.9, 0))
set_legend(plot, "Convergence Points")
set_marker_size(plot, 24)

field = create_vector_field(chart, NELEMS*NELEMS)
set_color(field, RGB(0, 0.6, 0.3))

update_vertex_buffer(plot, DPOINTS)
update_vertex_buffer(field, points)
update_direction_buffer(field, dirs)

draw_chart(!isclosed, window, chart)
save_framebuffer(window, "field.png")


text:
# preliminary definitions: constants, macros
# and functions used throughout the code
const _msk64 = ~uint64(0)
macro _mskr(l) :(_msk64 >>> (63 & (64-$(esc(l))))) end
macro _div64(l) :($(esc(l)) >>> 6) end
macro _mod64(l) :($(esc(l)) & 63) end
macro _msk_end(l) :(@_mskr @_mod64 $(esc(l))) end
num_bit_chunks(n::Int) = @_div64 (n+63)

## BitArray

# notes: bits are stored in contiguous chunks
#        unused bits must always be set to 0
type BitArray{N} <: DenseArray{Bool, N}
    chunks::Vector{Uint64}
    len::Int
    dims::NTuple{N,Int}
    function BitArray(dims::Int...)
        if length(dims) != N
            error("number of dimensions must be $N (got $(length(dims)))")
        end
        n = 1
        for d in dims
            if d < 0
                error("dimension size must be nonnegative (got $d)")
            end
            n *= d
        end
        nc = num_bit_chunks(n)
        chunks = Array(Uint64, nc)
        if nc > 0
            chunks[end] = uint64(0)
        end
        b = new(chunks, n)
        if N != 1
            b.dims = dims
        end
        return b
    end
end

BitArray{N}(dims::NTuple{N,Int}) = BitArray{N}(dims...)
BitArray(dims::Int...) = BitArray(dims)

typealias BitVector BitArray{1}
typealias BitMatrix BitArray{2}

## utility functions ##

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

size(B::BitVector, d) = (d==1 ? B.len : d>1 ? 1 : error("dimensions should be positive (got $d)"))
size{N}(B::BitArray{N}, d) = (d>N ? 1 : B.dims[d])

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

## Aux functions ##

get_chunks_id(i::Integer) = @_div64(int(i)-1)+1, @_mod64(int(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::Integer, src::Vector{Uint64}, pos_s::Integer, numbits::Integer)
    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)

    if delta_kd == 0
        return
    end

    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::Integer, pos_s::Integer, numbits::Integer)
    if pos_d == pos_s
        return
    elseif pos_d < pos_s
        return copy_chunks!(chunks, pos_d, chunks, pos_s, numbits)
    end

    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::Integer, numbits::Integer)
    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

## custom iterator ##
start(B::BitArray) = 0
next(B::BitArray, i::Int) = (B.chunks[@_div64(i)+1] & (uint64(1)<<@_mod64(i)) != 0, i+1)
done(B::BitArray, i::Int) = i >= length(B)

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

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

similar(B::BitArray, T::Type{Bool}, dims::Dims) = BitArray(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, dims)

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

falses(args...) = fill!(BitArray(args...), false)
trues(args...) = fill!(BitArray(args...), 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 copy!(dest::BitArray, src::BitArray)
    destc = dest.chunks; srcc = src.chunks
    nc_d = length(destc)
    nc_s = length(srcc)
    nc = min(nc_s, nc_d)
    if nc == 0
        return dest
    end
    @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 length(src)
            msk_d = ~msk_s
            destc[nc] = (msk_d & destc[nc]) | (msk_s & srcc[nc])
        end
    end
    return dest
end

function copy!(dest::BitArray, pos_d::Integer, src::BitArray, pos_s::Integer, numbits::Integer)
    if pos_s+numbits-1 > length(src) || pos_d+numbits-1 > length(dest) || pos_d < 1 || pos_s < 1
        throw(BoundsError())
    end
    copy_chunks!(dest.chunks, pos_d, src.chunks, pos_s, numbits)
    return dest
end

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

## Conversions ##

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

convert{T,N}(::Type{BitArray}, A::AbstractArray{T,N}) = convert(BitArray{N},A)
function convert{T,N}(::Type{BitArray{N}}, A::AbstractArray{T,N})
    B = BitArray(size(A))
    Bc = B.chunks
    l = length(B)
    if l == 0
        return B
    end
    ind = 1
    @inbounds begin
        for i = 1:length(Bc)-1
            u = uint64(1)
            c = uint64(0)
            for j = 0:63
                if bool(A[ind])
                    c |= u
                end
                ind += 1
                u <<= 1
            end
            Bc[i] = c
        end
        u = uint64(1)
        c = uint64(0)
        for j = 0:@_mod64(l-1)
            if bool(A[ind])
                c |= u
            end
            ind += 1
            u <<= 1
        end
        Bc[end] = c
    end
    return B
end

convert{N}(::Type{BitArray{N}}, B::BitArray{N}) = B
convert{T,N}(::Type{AbstractArray{T,N}}, B::BitArray{N}) = convert(Array{T,N}, B)

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

# shorthand forms BitArray <-> Array
bitunpack{N}(B::BitArray{N}) = convert(Array{Bool,N}, B)
bitpack{T,N}(A::AbstractArray{T,N}) = convert(BitArray{N}, A)

## Random ##

function bitarray_rand_fill!(B::BitArray)
    if length(B) == 0
        return B
    end
    Bc = B.chunks
    for i = 1 : length(Bc) - 1
        Bc[i] = rand(Uint64)
    end
    msk = @_msk_end length(B)
    Bc[end] = msk & rand(Uint64)
    return B
end

## Indexing: getindex ##

function unsafe_bitgetindex(Bc::Vector{Uint64}, i::Int)
    return (Bc[@_div64(i-1)+1] & (uint64(1)<<@_mod64(i-1))) != 0
end

function getindex(B::BitArray, i::Int)
    1 <= i <= length(B) || throw(BoundsError())
    return unsafe_bitgetindex(B.chunks, i)
end

getindex(B::BitArray, i::Real) = getindex(B, to_index(i))

getindex(B::BitArray) = getindex(B, 1)

# 0d bitarray
getindex(B::BitArray{0}) = unsafe_bitgetindex(B.chunks, 1)

function getindex{T<:Real}(B::BitArray, I::AbstractVector{T})
    X = BitArray(length(I))
    lB = length(B)
    Xc = X.chunks
    Bc = B.chunks
    ind = 1
    for i in I
        # faster X[ind] = B[i]
        j = to_index(i)
        1 <= j <= lB || throw(BoundsError())
        unsafe_bitsetindex!(Xc, unsafe_bitgetindex(Bc, j), ind)
        ind += 1
    end
    return X
end

# logical indexing

# (multiple signatures for disambiguation)
for IT in [AbstractVector{Bool}, AbstractArray{Bool}]
    @eval function getindex(B::BitArray, I::$IT)
        checkbounds(B, I)
        n = sum(I)
        X = BitArray(n)
        Xc = X.chunks
        Bc = B.chunks
        ind = 1
        for i = 1:length(I)
            if I[i]
                # faster X[ind] = B[i]
                unsafe_bitsetindex!(Xc, unsafe_bitgetindex(Bc, i), ind)
                ind += 1
            end
        end
        return X
    end
end

## Indexing: setindex! ##

function unsafe_bitsetindex!(Bc::Array{Uint64}, x::Bool, i::Int)
    i1, i2 = get_chunks_id(i)
    u = uint64(1) << i2
    @inbounds begin
        if x
            Bc[i1] |= u
        else
            Bc[i1] &= ~u
        end
    end
end

setindex!(B::BitArray, x) = setindex!(B, convert(Bool,x), 1)

function setindex!(B::BitArray, x::Bool, i::Int)
    1 <= i <= length(B) || throw(BoundsError())
    unsafe_bitsetindex!(B.chunks, x, i)
    return B
end

# logical indexing

function setindex!(B::BitArray, x, I::AbstractArray{Bool})
    checkbounds(B, I)
    y = convert(Bool, x)
    Bc = B.chunks
    @inbounds for i = 1:length(I)
        # faster B[i] = y
        I[i] && unsafe_bitsetindex!(Bc, y, i)
    end
    return B
end

function setindex!(B::BitArray, X::AbstractArray, I::AbstractArray{Bool})
    checkbounds(B, I)
    Bc = B.chunks
    c = 1
    @inbounds for i = 1:length(I)
        if I[i]
            # faster B[i] = X[c]
            unsafe_bitsetindex!(Bc, convert(Bool, X[c]), i)
            c += 1
        end
    end
    if length(X) != c-1
        throw(DimensionMismatch("assigned $(length(X)) elements to length $(c-1) destination"))
    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
        ccall(:jl_array_grow_end, Void, (Any, Uint), 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)
    if n1 == 0
        return B
    end
    Bc = B.chunks
    k0 = length(Bc)
    k1 = num_bit_chunks(n0 + n1)
    if k1 > k0
        ccall(:jl_array_grow_end, Void, (Any, Uint), 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::AbstractVector{Bool}) = append!(B, bitpack(items))
append!(A::Vector{Bool}, items::BitVector) = append!(A, bitunpack(items))

function prepend!(B::BitVector, items::BitVector)
    n0 = length(B)
    n1 = length(items)
    if n1 == 0
        return B
    end
    Bc = B.chunks
    k0 = length(Bc)
    k1 = num_bit_chunks(n0 + n1)
    if k1 > k0
        ccall(:jl_array_grow_end, Void, (Any, Uint), 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::AbstractVector{Bool}) = prepend!(B, bitpack(items))
prepend!(A::Vector{Bool}, items::BitVector) = prepend!(A, bitunpack(items))

function sizehint(B::BitVector, sz::Integer)
    ccall(:jl_array_sizehint, Void, (Any, Uint), B.chunks, num_bit_chunks(sz))
    return B
end

function resize!(B::BitVector, n::Integer)
    if n < 0
        throw(BoundsError())
    end
    n0 = length(B)
    if n <= n0
        splice!(B, n+1:n0)
        return B
    end
    Bc = B.chunks
    k0 = length(Bc)
    k1 = num_bit_chunks(int(n))
    if k1 > k0
        ccall(:jl_array_grow_end, Void, (Any, Uint), Bc, k1 - k0)
        Bc[end] = uint64(0)
    end
    B.len = n
    return B
end

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

    l = @_mod64 length(B)
    if l == 1
        ccall(:jl_array_del_end, Void, (Any, Uint), B.chunks, 1)
    end
    B.len -= 1

    return item
end

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

    Bc = B.chunks

    l = @_mod64 length(B)
    if l == 0
        ccall(:jl_array_grow_end, Void, (Any, Uint), 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 shift!(B::BitVector)
    if isempty(B)
        error("argument must not be empty")
    end
    @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
            ccall(:jl_array_del_end, Void, (Any, Uint), Bc, 1)
        else
            Bc[end] >>>= 1
        end
        B.len -= 1
    end

    return item
end

function insert!(B::BitVector, i::Integer, item)
    if i < 1
        throw(BoundsError())
    end
    item = convert(Bool, item)

    n = length(B)
    if i > n
        x = falses(i - n)
        append!(B, x)
    else
        Bc = B.chunks

        k, j = get_chunks_id(i)

        l = @_mod64 length(B)
        if l == 0
            ccall(:jl_array_grow_end, Void, (Any, Uint), 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)
    end
    B[i] = item
end

function _deleteat!(B::BitVector, i::Integer)

    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
            ccall(:jl_array_del_end, Void, (Any, Uint), Bc, 1)
        elseif length(Bc) > k
            Bc[end] >>>= 1
        end
    end

    B.len -= 1

    return B
end

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

    return _deleteat!(B, i)
end

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

    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)

    if delta_k < 0
        ccall(:jl_array_del_end, Void, (Any, Uint), Bc, -delta_k)
    end

    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)
    s = start(inds)
    done(inds, s) && return B

    Bc = B.chunks

    (p, s) = next(inds, s)
    q = p+1
    new_l -= 1
    while !done(inds, s)
        (i,s) = next(inds, s)
        if !(q <= i <= n)
            i < q && error("indices must be unique and sorted")
            throw(BoundsError())
        end
        new_l -= 1
        if i > q
            copy_chunks!(Bc, p, Bc, q, i-q)
            p += i-q
        end
        q = i+1
    end

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

    delta_k = num_bit_chunks(new_l) - length(Bc)
    delta_k < 0 && ccall(:jl_array_del_end, Void, (Any, Uint), 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)
    n = length(B)
    if !(1 <= i <= n)
        throw(BoundsError())
    end

    v = B[i]   # TODO: change to a copy if/when subscripting becomes an ArrayView
    _deleteat!(B, i)
    return v
end
splice!(B::BitVector, i::Integer, ins::BitVector) = splice!(B, int(i):int(i), ins)
splice!(B::BitVector, i::Integer, ins::AbstractVector{Bool}) = splice!(B, i, bitpack(ins))

const _default_bit_splice = BitVector(0)

function splice!(B::BitVector, r::UnitRange{Int}, ins::BitVector = _default_bit_splice)
    n = length(B)
    i_f = first(r)
    i_l = last(r)
    if !(1 <= i_f <= n+1)
        throw(BoundsError())
    end
    if !(i_l <= n)
        throw(BoundsError())
    end
    if (i_f > n)
        append!(B, ins)
        return BitVector(0)
    end

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

    Bc = B.chunks

    lins = length(ins)
    ldel = length(r)

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

    if delta_k > 0
        ccall(:jl_array_grow_end, Void, (Any, Uint), Bc, delta_k)
    end
    copy_chunks!(Bc, i_f+lins, Bc, i_l+1, n-i_l)
    copy_chunks!(Bc, i_f, ins.chunks, 1, lins)
    if delta_k < 0
        ccall(:jl_array_del_end, Void, (Any, Uint), Bc, -delta_k)
    end

    B.len = new_l

    if new_l > 0
        Bc[end] &= @_msk_end new_l
    end

    return v
end
splice!(B::BitVector, r::UnitRange{Int}, ins::AbstractVector{Bool}) = splice!(B, r, bitpack(ins))

function empty!(B::BitVector)
    ccall(:jl_array_del_end, Void, (Any, Uint), B.chunks, length(B.chunks))
    B.len = 0
    return B
end

## Misc functions

for f in (:iround, :itrunc, :ifloor, :iceil, :abs)
    @eval ($f)(B::BitArray) = copy(B)
end

## Unary operators ##

function (-)(B::BitArray)
    A = zeros(Int, size(B))
    l = length(B)
    if l == 0
        return A
    end
    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
sign(B::BitArray) = copy(B)

function (~)(B::BitArray)
    C = similar(B)
    Bc = B.chunks
    if !isempty(Bc)
        Cc = C.chunks
        for i = 1:length(Bc)-1
            Cc[i] = ~Bc[i]
        end
        msk = @_msk_end length(B)
        Cc[end] = msk & (~Bc[end])
    end
    return C
end

function flipbits!(B::BitArray)
    Bc = B.chunks
    @inbounds if !isempty(Bc)
        for i = 1:length(B.chunks) - 1
            Bc[i] = ~Bc[i]
        end
        msk = @_msk_end length(B)
        Bc[end] = msk & (~Bc[end])
    end
    return B
end

!(B::BitArray) = ~B

## Binary arithmetic operators ##

for f in (:+, :-)
    @eval function ($f)(A::BitArray, B::BitArray)
        r = Array(Int, promote_shape(size(A), size(B)))
        ai = start(A)
        bi = start(B)
        ri = 1
        while !done(A, ai)
            a, ai = next(A, ai)
            b, bi = next(B, bi)
            @inbounds r[ri] = ($f)(a, b)
            ri += 1
        end
        return r
    end
end
for f in (:.+, :.-),
    (arg1, arg2, T, fargs) in ((:(B::BitArray), :(x::Bool)    , Int                                   , :(b, x)),
                               (:(B::BitArray), :(x::Number)  , :(promote_array_type(typeof(x), Bool)), :(b, x)),
                               (:(x::Bool)    , :(B::BitArray), Int                                   , :(x, b)),
                               (:(x::Number)  , :(B::BitArray), :(promote_array_type(typeof(x), Bool)), :(x, b)))
    @eval function ($f)($arg1, $arg2)
        r = Array($T, size(B))
        bi = start(B)
        ri = 1
        while !done(B, bi)
            b, bi = next(B, bi)
            @inbounds r[ri] = ($f)($fargs...)
            ri += 1
        end
        return r
    end
end

for f in (:/, :\)
    @eval begin
        ($f)(A::BitArray, B::BitArray) = ($f)(bitunpack(A), bitunpack(B))
    end
end
(/)(B::BitArray, x::Number) = (/)(bitunpack(B), x)
(/)(x::Number, B::BitArray) = (/)(x, bitunpack(B))

function div(A::BitArray, B::BitArray)
    shp = promote_shape(size(A), size(B))
    all(B) || throw(DivideError())
    return reshape(copy(A), shp)
end
div(A::BitArray, B::Array{Bool}) = div(A, bitpack(B))
div(A::Array{Bool}, B::BitArray) = div(bitpack(A), B)
function div(B::BitArray, x::Bool)
    return x ? copy(B) : throw(DivideError())
end
function div(x::Bool, B::BitArray)
    all(B) || throw(DivideError())
    return x ? trues(size(B)) : falses(size(B))
end
function div(x::Number, B::BitArray)
    all(B) || throw(DivideError())
    pt = promote_array_type(typeof(x), Bool)
    y = div(x, true)
    reshape(pt[ y for i = 1:length(B) ], size(B))
end

function mod(A::BitArray, B::BitArray)
    shp = promote_shape(size(A), size(B))
    all(B) || throw(DivideError())
    return falses(shp)
end
mod(A::BitArray, B::Array{Bool}) = mod(A, bitpack(B))
mod(A::Array{Bool}, B::BitArray) = mod(bitpack(A), B)
function mod(B::BitArray, x::Bool)
    return x ? falses(size(B)) : throw(DivideError())
end
function mod(x::Bool, B::BitArray)
    all(B) || throw(DivideError())
    return falses(size(B))
end
function mod(x::Number, B::BitArray)
    all(B) || throw(DivideError())
    pt = promote_array_type(typeof(x), Bool)
    y = mod(x, true)
    reshape(pt[ y for i = 1:length(B) ], size(B))
end

for f in (:div, :mod)
    @eval begin
        function ($f)(B::BitArray, x::Number)
            F = Array(promote_array_type(typeof(x), Bool), size(B))
            for i = 1:length(F)
                F[i] = ($f)(B[i], x)
            end
            return F
        end
    end
end

function (&)(B::BitArray, x::Bool)
    x ? copy(B) : falses(size(B))
end
(&)(x::Bool, B::BitArray) = B & x

function (|)(B::BitArray, x::Bool)
    x ? trues(size(B)) : copy(B)
end
(|)(x::Bool, B::BitArray) = B | x

function ($)(B::BitArray, x::Bool)
    x ? ~B : copy(B)
end
($)(x::Bool, B::BitArray) = B $ x

for f in (:&, :|, :$)
    @eval begin
        function ($f)(A::BitArray, B::BitArray)
            F = BitArray(promote_shape(size(A),size(B))...)
            Fc = F.chunks
            Ac = A.chunks
            Bc = B.chunks
            if !isempty(Ac) && !isempty(Bc)
                for i = 1:length(Fc) - 1
                    Fc[i] = ($f)(Ac[i], Bc[i])
                end
                msk = @_msk_end length(F)
                Fc[end] = msk & ($f)(Ac[end], Bc[end])
            end
            return F
        end
        ($f)(A::Array{Bool}, B::BitArray) = ($f)(bitpack(A), B)
        ($f)(B::BitArray, A::Array{Bool}) = ($f)(B, bitpack(A))
        ($f)(x::Number, B::BitArray) = ($f)(x, bitunpack(B))
        ($f)(B::BitArray, x::Number) = ($f)(bitunpack(B), x)
    end
end

function (.^)(B::BitArray, x::Bool)
    x ? copy(B) : trues(size(B))
end
function (.^)(x::Bool, B::BitArray)
    x ? trues(size(B)) : ~B
end
function (.^)(x::Number, B::BitArray)
    z = x ^ false
    u = x ^ true
    reshape([ B[i] ? u : z for i = 1:length(B) ], size(B))
end
function (.^)(B::BitArray, x::Integer)
    if x == 0
        return trues(size(B))
    elseif x < 0
        throw(DomainError())
    else
        return copy(B)
    end
end
function (.^){T<:Number}(B::BitArray, x::T)
    if x == 0
        return ones(typeof(true ^ x), size(B))
    elseif T <: Real && x > 0
        return convert(Array{T}, B)
    else
        z = nothing
        u = nothing
        zerr = nothing
        uerr = nothing
        try
            z = false^x
        catch err
            zerr = err
        end
        try
            u = true^x
        catch err
            uerr = err
        end
        if zerr == nothing && uerr == nothing
            t = promote_type(typeof(z), typeof(u))
        elseif zerr == nothing
            t = typeof(z)
        else
            t = typeof(u)
        end
        F = Array(t, size(B))
        for i = 1:length(B)
            if B[i]
                if uerr == nothing
                    F[i] = u
                else
                    throw(uerr)
                end
            else
                if zerr == nothing
                    F[i] = z
                else
                    throw(zerr)
                end
            end
        end
        return F
    end
end

(.*)(x::Bool, B::BitArray) = x & B
(.*)(B::BitArray, x::Bool) = B & x
(.*)(x::Number, B::BitArray) = x .* bitunpack(B)
(.*)(B::BitArray, x::Number) = bitunpack(B) .* x

## promotion to complex ##

# TODO?

## comparison operators ##

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

function (!=)(A::BitArray, B::BitArray)
    if size(A) != size(B)
        return true
    end
    return A.chunks != B.chunks
end

## Data movement ##

# TODO some of this could be optimized

function slicedim(A::BitArray, d::Integer, i::Integer)
    d_in = size(A)
    leading = d_in[1:(d-1)]
    d_out = tuple(leading..., 1, d_in[(d+1):end]...)

    M = prod(leading)
    N = length(A)
    stride = M * d_in[d]

    B = BitArray(d_out)
    index_offset = 1 + (i-1)*M

    l = 1

    if M==1
        for j=0:stride:(N-stride)
            B[l] = A[j + index_offset]
            l += 1
        end
    else
        for j=0:stride:(N-stride)
            offs = j + index_offset
            for k=0:(M-1)
                B[l] = A[offs + k]
                l += 1
            end
        end
    end
    return B
end

function flipdim(A::BitArray, d::Integer)
    nd = ndims(A)
    sd = d > nd ? 1 : size(A, d)
    if sd == 1
        return copy(A)
    end

    B = similar(A)

    nnd = 0
    for i = 1:nd
        nnd += int(size(A,i)==1 || i==d)
    end
    if nnd==nd
        # flip 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_bits(src::Uint64)
    z    = src
    z    = ((z >>>  1) & 0x5555555555555555) | ((z <<  1) & 0xaaaaaaaaaaaaaaaa)
    z    = ((z >>>  2) & 0x3333333333333333) | ((z <<  2) & 0xcccccccccccccccc)
    z    = ((z >>>  4) & 0x0f0f0f0f0f0f0f0f) | ((z <<  4) & 0xf0f0f0f0f0f0f0f0)
    z    = ((z >>>  8) & 0x00ff00ff00ff00ff) | ((z <<  8) & 0xff00ff00ff00ff00)
    z    = ((z >>> 16) & 0x0000ffff0000ffff) | ((z << 16) & 0xffff0000ffff0000)
    return ((z >>> 32) & 0x00000000ffffffff) | ((z << 32) & 0xffffffff00000000)
end

function reverse!(B::BitVector)
    n = length(B)
    if n == 0
        return B
    end
    pnc = length(B.chunks) & 1
    hnc = (length(B.chunks) >>> 1)

    aux_chunks = Array(Uint64, 1)

    for i = 1 : hnc
        j = ((i - 1) << 6)
        aux_chunks[1] = reverse_bits(B.chunks[i])
        copy_chunks!(B.chunks, j+1, B.chunks, n-63-j, 64)
        B.chunks[i] = reverse_bits(B.chunks[i])
        copy_chunks!(B.chunks, n-63-j, aux_chunks, 1, 64)
    end

    if pnc == 0
        return B
    end

    i = hnc + 1
    j = hnc << 6
    l = (@_mod64 (n+63)) + 1
    msk = @_mskr l

    aux_chunks[1] = reverse_bits(B.chunks[i] & msk)
    aux_chunks[1] >>>= (64 - l)
    copy_chunks!(B.chunks, j+1, aux_chunks, 1, l)

    return B
end

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

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

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

(>>)(B::BitVector, i::Int32) = B >>> i
(>>)(B::BitVector, i::Integer) = B >>> i

function rol(B::BitVector, i::Integer)
    n = length(B)
    i %= n
    i == 0 && return copy(B)
    i < 0 && return ror(B, -i)
    A = BitArray(n)
    copy_chunks!(A.chunks, 1, B.chunks, i+1, n-i)
    copy_chunks!(A.chunks, n-i+1, B.chunks, 1, i)
    return A
end

function ror(B::BitVector, i::Integer)
    n = length(B)
    i %= n
    i == 0 && return copy(B)
    i < 0 && return rol(B, -i)
    A = BitArray(n)
    copy_chunks!(A.chunks, i+1, B.chunks, 1, n-i)
    copy_chunks!(A.chunks, 1, B.chunks, n-i+1, i)
    return A
end

#TODO: rol!, ror!

## countnz & find ##

function countnz(B::BitArray)
    n = 0
    Bc = B.chunks
    @inbounds for i = 1:length(Bc)
        n += count_ones(Bc[i])
    end
    return n
end

# returns the index of the next non-zero element, or 0 if all zeros
function findnext(B::BitArray, start::Integer)
    if start <= 0
        throw(BoundsError())
    elseif start > length(B)
        return 0
    end
    Bc = B.chunks

    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 0
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::Integer)
    if start <= 0
        throw(BoundsError())
    elseif start > length(B)
        return 0
    end

    Bc = B.chunks
    l = length(Bc)
    if l == 0
        return 0
    end

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

    @inbounds begin
        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
        ce = Bc[end]
    end
    if ce != @_msk_end length(B)
        return (l-1) << 6 + trailing_ones(ce) + 1
    end
    return 0
end
findfirstnot(B::BitArray) = findnextnot(B,1)

# returns the index of the first matching element
function findnext(B::BitArray, v, start::Integer)
    if v == false
        return findnextnot(B, start)
    elseif v == true
        return findnext(B, start)
    else
        return 0
    end
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
function findnext(testf::Function, B::BitArray, start::Integer)
    f0::Bool = testf(false)
    f1::Bool = testf(true)
    if length(B) == 0 || !(f0 || f1)
        return 0
    elseif f0 && f1
        return 1
    elseif !f0 && f1
        return findnext(B, start)
    else
        return findnextnot(B, start)
    end
end
#findfirst(testf::Function, B::BitArray) = findnext(testf, B, 1)  ## defined in array.jl

function find(B::BitArray)
    l = length(B)
    nnzB = countnz(B)
    I = Array(Int, nnzB)
    if nnzB == 0
        return I
    end
    Bc = B.chunks
    Bcount = 1
    Icount = 1
    for i = 1:length(Bc)-1
        u = uint64(1)
        c = Bc[i]
        for j = 1:64
            if c & u != 0
                I[Icount] = Bcount
                Icount += 1
            end
            Bcount += 1
            u <<= 1
        end
    end
    u = uint64(1)
    c = Bc[end]
    for j = 0:@_mod64(l-1)
        if c & u != 0
            I[Icount] = Bcount
            Icount += 1
        end
        Bcount += 1
        u <<= 1
    end
    return I
end

findn(B::BitVector) = find(B)

function findn(B::BitMatrix)
    nnzB = countnz(B)
    I = Array(Int, nnzB)
    J = Array(Int, nnzB)
    count = 1
    for j=1:size(B,2), i=1:size(B,1)
        if B[i,j]
            I[count] = i
            J[count] = j
            count += 1
        end
    end
    return (I, J)
end

function findnz(B::BitMatrix)
    I, J = findn(B)
    return (I, J, trues(length(I)))
end

nonzeros(B::BitArray) = trues(countnz(B))

## Reductions ##

sum(A::BitArray, region) = reducedim(+,A,region,0,Array(Int,reduced_dims(A,region)))

sum(B::BitArray) = countnz(B)

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

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

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

## map over bitarrays ##

function map!(f::Callable, A::Union(StridedArray,BitArray))
    for i=1:length(A)
        A[i] = f(A[i])
    end
    return A
end

function map!(f::Callable, dest::Union(StridedArray,BitArray), A::Union(StridedArray,BitArray))
    for i=1:length(A)
        dest[i] = f(A[i])
    end
    return dest
end

function map!(f::Callable, dest::Union(StridedArray,BitArray), A::Union(StridedArray,BitArray), B::Union(StridedArray,BitArray))
    for i=1:length(A)
        dest[i] = f(A[i], B[i])
    end
    return dest
end

function map!(f::Callable, dest::Union(StridedArray,BitArray), A::Union(StridedArray,BitArray), B::Number)
    for i=1:length(A)
        dest[i] = f(A[i], B)
    end
    return dest
end

function map!(f::Callable, dest::Union(StridedArray,BitArray), A::Number, B::Union(StridedArray,BitArray))
    for i=1:length(B)
        dest[i] = f(A, B[i])
    end
    return dest
end

function map!(f::Callable, dest::Union(StridedArray,BitArray), As::Union(StridedArray,BitArray)...)
    n = length(As[1])
    i = 1
    ith = a->a[i]
    for i=1:n
        dest[i] = f(map(ith, As)...)
    end
    return dest
end

## Filter ##

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

## Transpose ##

transpose(B::BitVector) = reshape(copy(B), 1, length(B))

# fast 8x8 bit transpose from Henry S. Warrens's "Hacker's Delight"
# http://www.hackersdelight.org/HDcode/transpose8.c.txt
function transpose8x8(x::Uint64)
    y = x
    t = (y $ (y >>> 7)) & 0x00aa00aa00aa00aa
    y = y $ t $ (t << 7)
    t = (y $ (y >>> 14)) & 0x0000cccc0000cccc
    y = y $ t $ (t << 14)
    t = (y $ (y >>> 28)) & 0x00000000f0f0f0f0
    return y $ t $ (t << 28)
end

function form_8x8_chunk(Bc::Vector{Uint64}, i1::Int, i2::Int, m::Int, cgap::Int, cinc::Int, nc::Int, msk8::Uint64)
    x = uint64(0)

    k, l = get_chunks_id(i1 + (i2 - 1) * m)
    r = 0
    for j = 1 : 8
        if k > nc
            break
        end
        x |= ((Bc[k] >>> l) & msk8) << r
        if l + 8 >= 64 && nc > k
            r0 = 8 - (@_mod64 (l + 8))
            x |= (Bc[k + 1] & (msk8 >>> r0)) << (r + r0)
        end
        k += cgap + (l + cinc >= 64 ? 1 : 0)
        l = @_mod64 (l + cinc)
        r += 8
    end
    return x
end

# note: assumes B is filled with 0's
function put_8x8_chunk(Bc::Vector{Uint64}, i1::Int, i2::Int, x::Uint64, m::Int, cgap::Int, cinc::Int, nc::Int, msk8::Uint64)
    k, l = get_chunks_id(i1 + (i2 - 1) * m)
    r = 0
    for j = 1 : 8
        if k > nc
            break
        end
        Bc[k] |= ((x >>> r) & msk8) << l
        if l + 8 >= 64 && nc > k
            r0 = 8 - (@_mod64 (l + 8))
            Bc[k + 1] |= ((x >>> (r + r0)) & (msk8 >>> r0))
        end
        k += cgap + (l + cinc >= 64 ? 1 : 0)
        l = @_mod64 (l + cinc)
        r += 8
    end
    return
end

function transpose(B::BitMatrix)
    l1 = size(B, 1)
    l2 = size(B, 2)
    Bt = falses(l2, l1)

    cgap1 = @_div64 l1
    cinc1 = @_mod64 l1

    cgap2 = @_div64 l2
    cinc2 = @_mod64 l2

    Bc = B.chunks
    Btc = Bt.chunks

    nc = length(Bc)

    for i = 1 : 8 : l1

        msk8_1 = uint64(0xff)
        if (l1 < i + 7)
            msk8_1 >>>= i + 7 - l1
        end

        for j = 1 : 8 : l2
            x = form_8x8_chunk(Bc, i, j, l1, cgap1, cinc1, nc, msk8_1)
            x = transpose8x8(x)

            msk8_2 = uint64(0xff)
            if (l2 < j + 7)
                msk8_2 >>>= j + 7 - l2
            end

            put_8x8_chunk(Btc, j, i, x, l2, cgap2, cinc2, nc, msk8_2)
        end
    end
    return Bt
end

ctranspose(B::BitArray) = transpose(B)

## Permute array dims ##

function permutedims(B::Union(BitArray,StridedArray), perm)
    dimsB = size(B)
    ndimsB = length(dimsB)
    (ndimsB == length(perm) && isperm(perm)) || error("no valid permutation of dimensions")
    dimsP = ntuple(ndimsB, i->dimsB[perm[i]])::typeof(dimsB)
    P = similar(B, dimsP)
    permutedims!(P, B, perm)
end

## Concatenation ##

function hcat(B::BitVector...)
    height = length(B[1])
    for j = 2:length(B)
        length(B[j]) == height || error("dimensions must match")
    end
    M = BitArray(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 = BitArray(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)
        if size(Aj, 1) != nrows; error("rows must match"); end
    end

    B = BitArray(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 = sum(a->size(a, 1), A)::Int
    ncols = size(A[1], 2)
    for j = 2:nargs
        if size(A[j], 2) != ncols; error("columns must match"); end
    end
    B = BitArray(nrows, ncols)
    Bc = B.chunks
    nrowsA = [size(a, 1) for a in A]
    Ac = [a.chunks for a in A]
    pos_d = 1
    pos_s = ones(Int, nargs)
    for j = 1:ncols
        for k = 1:nargs
            copy_chunks!(Bc, pos_d, Ac[k], pos_s[k], nrowsA[k])
            pos_s[k] += nrowsA[k]
            pos_d += nrowsA[k]
        end
    end
    return B
end

# general case, specialized for BitArrays and Integers
function cat(catdim::Integer, X::Union(BitArray, Integer)...)
    nargs = length(X)
    # using integers results in conversion to Array{Int}
    # (except in the all-Bool case)
    has_bitarray = false
    has_integer = false
    for a in X
        if isa(a, BitArray)
            has_bitarray = true
        else
            has_integer = true
        end
    end
    # just integers and no BitArrays -> general case
    has_bitarray || return invoke(cat, (Integer, Any...), catdim, X...)
    dimsX = map((a->isa(a,BitArray) ? size(a) : (1,)), X)
    ndimsX = map((a->isa(a,BitArray) ? ndims(a) : 1), X)
    d_max = maximum(ndimsX)

    if catdim > d_max + 1
        for i=1:nargs
            if dimsX[1] != dimsX[i]
                error("all inputs must have same dimensions when concatenating along a higher dimension");
            end
        end
    elseif nargs >= 2
        for d=1:d_max
            if d == catdim; continue; end
            len = d <= ndimsX[1] ? dimsX[1][d] : 1
            for i = 2:nargs
                if len != (d <= ndimsX[i] ? dimsX[i][d] : 1)
                    error("mismatch in dimension ", d)
                end
            end
        end
    end

    cat_ranges = ntuple(nargs, i->(catdim <= ndimsX[i] ? dimsX[i][catdim] : 1))

    function compute_dims(d)
        if d == catdim
            if catdim <= d_max
                return sum(cat_ranges)
            else
                return nargs
            end
        else
            if d <= ndimsX[1]
                return dimsX[1][d]
            else
                return 1
            end
        end
    end

    ndimsC = max(catdim, d_max)
    dimsC = ntuple(ndimsC, compute_dims)::(Int...)
    typeC = promote_type(map(x->isa(x,BitArray) ? eltype(x) : typeof(x), X)...)
    if !has_integer || typeC == Bool
        C = BitArray(dimsC)
    else
        C = Array(typeC, dimsC)
    end

    range = 1
    for k = 1:nargs
        nextrange = range+cat_ranges[k]
        cat_one = ntuple(ndimsC, i->(i != catdim ?
                                     (1:dimsC[i]) : (range:nextrange-1) ))
        # note: when C and X are BitArrays, this calls
        #       the special assign with ranges
        C[cat_one...] = X[k]
        range = nextrange
    end
    return C
end

# hvcat -> use fallbacks in abstractarray.jl

isequal(A::BitArray, B::BitArray) = (A == B)

# Hashing

hash(B::BitArray) = hash((size(B), B.chunks))

max_stars_repo_name: alyst/julia

max_issues_repo_name: alyst/julia

max_forks_repo_name: alyst/julia


text:
"""
```
default_settings!(m::AbstractDSGEModel)
```

Default Settings are constructed, initialized and added to `m.settings`.
"""
function default_settings!(m::AbstractDSGEModel)

    settings = m.settings

    # I/O File locations
    saveroot = normpath(joinpath(dirname(@__FILE__), "..","save"))
    datapath = normpath(joinpath(dirname(@__FILE__), "..","save","input_data"))

    settings[:saveroot] = Setting(:saveroot, saveroot, "Root of data directory structure")
    settings[:dataroot] = Setting(:dataroot, datapath, "Input data directory path")

    # Data settings for released and conditional data. Default behavior is to set vintage
    # of data to today's date.
    vint = Dates.format(now(), DSGE_DATE_FORMAT)
    settings[:data_vintage] = Setting(:data_vintage, vint, true, "vint",
        "Data vintage")
    settings[:data_id] = Setting(:data_id, 3,
        "Dataset identifier")
    settings[:cond_vintage] = Setting(:cond_vintage, vint,
        "Conditional data vintage")
    settings[:cond_id] = Setting(:cond_id, 2,
        "Conditional dataset identifier")
    settings[:cond_full_names] = Setting(:cond_full_names, [:obs_gdp, :obs_corepce, :obs_spread,
                                                            :obs_nominalrate, :obs_longrate],
        "Observables used in conditional forecasts")
    settings[:cond_semi_names] = Setting(:cond_semi_names, [:obs_spread, :obs_nominalrate,
                                                            :obs_longrate],
        "Observables used in semiconditional forecasts")
    settings[:use_population_forecast] = Setting(:use_population_forecast, false,
        "Whether to use population forecasts as data")
    settings[:population_mnemonic] = Setting(:population_mnemonic, Nullable(:CNP16OV__FRED),
        "Mnemonic of FRED data series for computing per-capita values (a Nullable{Symbol})")
    settings[:hpfilter_population] = Setting(:hpfilter_population, true,
        "Whether to HP filter combined population and forecast")
    settings[:rate_expectations_source] = Setting(:rate_expectations_source, :ois,
        "Data source of expectations for anticipated shocks")

    # Dates
    settings[:date_presample_start] = Setting(:date_presample_start, quartertodate("1959-Q3"),
        "Start date of pre-sample")
    settings[:date_mainsample_start] = Setting(:date_mainsample_start, quartertodate("1960-Q1"),
        "Start date of main sample")
    settings[:date_zlb_start] = Setting(:date_zlb_start, quartertodate("2008-Q4"),
        "Start date of first zero lower bound regime")
    settings[:date_forecast_start] = Setting(:date_forecast_start,
                                             Dates.lastdayofquarter(Dates.today()),
        "Start date of forecast period")
    settings[:date_conditional_end] = Setting(:date_conditional_end,
                                              Dates.lastdayofquarter(Dates.today()),
        "End date of conditional data period")

    # Anticipated shocks
    settings[:n_mon_anticipated_shocks] = Setting(:n_mon_anticipated_shocks, 0,
        "Number of anticipated policy shocks")
    settings[:n_mon_anticipated_shocks_padding] = Setting(:n_mon_anticipated_shocks_padding, 20,
        "Padding for anticipated policy shocks")

    # Monetary policy shock name and observable
    settings[:monetary_policy_shock] = Setting(:monetary_policy_shock, :rm_sh,
                                               "Name of the monetary policy shock in field `exogenous_shocks` " *
                                               "of a concrete subtype of AbstractDSGEModel")
    settings[:nominal_rate_observable] = Setting(:nominal_rate_observable, :obs_nominalrate,
                                                 "Name of the observable for monetary policy's nominal interest rate.")

    # General computation
    settings[:use_parallel_workers] = Setting(:use_parallel_workers, true,
        "Use available parallel workers in computations")

    # Solution
    settings[:solution_method] = Setting(:solution_method, :gensys, "What method the model uses to solve for the
                                         state space model form. The options are :gensys and :klein")

    # Optimization
    settings[:reoptimize] = Setting(:reoptimize, true,
        "Optimize the posterior mode. If false, reads in mode from a file.")
    settings[:optimization_method] = Setting(:optimization_method, :csminwel,
        "Method for finding the posterior mode")
    settings[:optimization_iterations] = Setting(:optimization_iterations, 100,
        "Number of iterations the optimizer should run for")
    settings[:optimization_step_size] = Setting(:optimization_step_size, 0.01,
        "Step size scaling factor for optimization")
	settings[:simulated_annealing_temperature] = Setting(:simulated_annealing_temperature,
                                                         Optim.log_temperature,
        "Temperature function for simulated annealing")
   settings[:simulated_annealing_block_proportion] = Setting(:simulated_annealing_block_proportion,
                                                             0.3,
        "Fraction of parameters to vary for each proposed move in simulated annealing")
   settings[:optimization_ftol] = Setting(:optimization_ftol, 1e-10,
        "Relative function difference threshold for optimization")
    settings[:optimization_xtol] = Setting(:optimization_xtol, 1e-10,
        "Relative input vector difference threshold for optimization")
    settings[:optimization_gtol] = Setting(:optimization_gtol, 1e-10,
        "Relative gradient difference threshold for optimization")
    settings[:combined_optimizer_max_cycles] = Setting(:combined_optimizer_max_cycles, 4,
        "Total number of cycles to use in the combined optimization routine")
    settings[:optimization_attempts] = Setting(:optimization_attempts, 4,
        "Number of times to attempt optimization in estimate()")

    # Estimation
    settings[:sampling_method] = Setting(:sampling_method, :MH,
        "The sampling method used, either :MH or :SMC")
    settings[:calculate_hessian] = Setting(:calculate_hessian, true,
        "Calculate the hessian at the mode")
    settings[:n_hessian_test_params] = Setting(:n_hessian_test_params, typemax(Int),
        "Max number of free params for which to calculate Hessian")
    settings[:use_chand_recursion] = Setting(:use_chand_recursion, false,
        "Use Chandrasekhar Recursions instead of standard Kalman filter")

    # Metropolis-Hastings
    settings[:n_mh_simulations] = Setting(:n_mh_simulations, 5000,
        "Number of draws saved (after thinning) per block in Metropolis-Hastings")
    settings[:n_mh_blocks] = Setting(:n_mh_blocks, 22,
        "Number of blocks for Metropolis-Hastings")
    settings[:n_mh_param_blocks] = Setting(:n_mh_param_blocks, 1,
        "Number of parameter blocks for Metropolis-Hastings")
    settings[:mh_adaptive_accept] = Setting(:mh_adaptive_accept, false,
        "Whether to adaptively solve for acceptance rate in Metropolis-Hastings")
    settings[:mh_target_accept] = Setting(:mh_target_accept, .25,
        "Target acceptance rate when adaptively adjusting the acceptance probability in Metropolis-Hastings")
    settings[:mh_c] = Setting(:mh_c, 0.5,
        "Initial step size used for adaptive acceptance rate in Metropolis-Hastings")
    settings[:mh_α] = Setting(:mh_α, 1.0,
        "Mixture proportion for adaptive acceptance rate in Metropolis-Hastings")
    settings[:n_mh_burn] = Setting(:n_mh_burn, 2,
        "Number of blocks to use as burn-in in Metropolis-Hastings")
    settings[:mh_thin] = Setting(:mh_thin, 5,
        "Metropolis-Hastings thinning step")
    settings[:mh_cc] = Setting(:mh_cc, 0.09,
        "Jump size for Metropolis-Hastings (after initialization)")
    settings[:mh_cc0] = Setting(:mh_cc0, 0.01,
        "Jump size for initialization of Metropolis-Hastings")

    # Forecast
    settings[:forecast_block_size] = Setting(:forecast_block_size, 5000,
        "Number of draws in each forecast block (before thinning by forecast_jstep)")
    settings[:forecast_start_block] = Setting(:forecast_start_block, 1,
        "Block at which to resume forecasting")
    settings[:forecast_input_file_overrides] = Setting(:forecast_input_file_overrides,
        Dict{Symbol, String}())
    settings[:forecast_jstep] = Setting(:forecast_jstep, 5,
        "Forecast thinning step (in addition to MH thinning step")
    settings[:forecast_uncertainty_override] = Setting(:forecast_uncertainty_override,
                                                       Nullable{Bool}(),
        "If non-null, overrides default drawing states/shocks behavior in smoother and forecast")
    settings[:forecast_smoother] = Setting(:forecast_smoother, :durbin_koopman,
        "Choice of smoother to use during forecasting. Can be :hamilton, :koopman, " *
        ":carter_kohn, or :durbin_koopman")
    settings[:forecast_horizons] = Setting(:forecast_horizons, 60,
        "Number of periods to forecast ahead")
    settings[:forecast_tdist_shocks] = Setting(:forecast_tdist_shocks, false,
        "Draw Students-t distributed shocks in forecast")
    settings[:forecast_tdist_df_val] = Setting(:forecast_tdist_df_val, 15,
        "Students-t degrees of freedom fixed value")
    settings[:forecast_zlb_value] = Setting(:forecast_zlb_value, 0.13/4,
        "Value of the zero lower bound in forecast periods, if we choose to enforce it")
    settings[:shockdec_startdate] = Setting(:shockdec_startdate, Nullable{Date}(),
        "Date of start of shock decomposition output period. If null, then shockdec starts " *
        "at date_mainsample_start")
    settings[:shockdec_enddate] = Setting(:shockdec_enddate, Nullable{Date}(),
        "Date of end of shock decomposition output period. If null, then shockdec ends at " *
        "date_forecast_end")
    settings[:impulse_response_horizons] = Setting(:impulse_response_horizons, 40,
        "Number of periods for which to calculate an impulse response")
    settings[:compute_shockdec_bands] = Setting(:compute_shockdec_bands, false,
        "Whether or not to compute bands for shock decomposition. Setting to false saves " *
        "signficant storage space.")
    settings[:time_varying_trends] = Setting(:time_varying_trends, true,
        "Indicator for whether the trends (CCC and DD) in the state space system are potentially time-varying. If false, " *
        "some functions will use shortcuts for non-time-varying trends.")

	# Sequential Monte Carlo
    settings[:n_particles] = Setting(:n_particles, 2000,
        "Number of particles for use in SMC")
    settings[:n_Φ] = Setting(:n_Φ, 300,
        "Number of stages in the tempering schedule")
  	settings[:λ] = Setting(:λ, 2.1,
        "The 'bending coefficient' λ in Φ(n) = (n/N(Φ))^λ")
    settings[:n_smc_blocks] = Setting(:n_smc_blocks, 1,
        "The number of parameter blocks in SMC")
    settings[:step_size_smc] = Setting(:step_size_smc, .5,
        "Scaling factor for covariance of the particles. Controls size of steps in mutation step")
    settings[:n_mh_steps_smc] = Setting(:n_mh_steps_smc, 1,
        "Number of Metropolis Hastings steps to attempt during the mutation step.")
    settings[:target_accept] = Setting(:target_accept, .25,
        "The initial target acceptance rate for new particles during mutation")
    settings[:resampler_smc] = Setting(:resampler_smc, :multinomial,
        "Which resampling method to use in SMC")
    settings[:mixture_proportion] = Setting(:mixture_proportion, 1.,
        "The mixture proportion for the mutation step's proposal distribution")
    settings[:debug_assertion] = Setting(:debug_assertion, false,
        "Print output (if applicable) when encountering an assertion error during SMC to help with debugging.")

    # Endogenous ϕ Schedule
    settings[:use_fixed_schedule] = Setting(:use_fixed_schedule, true,
        "Boolean indicating whether or not to use a fixed tempering (ϕ) schedule")
    settings[:resampling_threshold] = Setting(:resampling_threshold, 0.5,
        "Threshold s.t. particles will be resampled when the population drops below threshold * N")
    settings[:adaptive_tempering_target_smc] = Setting(:adaptive_tempering_target_smc, 0.97, false,
         "adpt", "Coefficient of the sample size to be targeted when solving for an endogenous ϕ " *
                 "or 0.0 if using a fixed schedule")

    # Generalized Tempering Settings
    settings[:tempered_update_prior_weight] = Setting(:tempered_update_prior_weight, 0.0, false,
         "prior_weight", "When bridging from old estimation, how much weight to put on prior.")
    settings[:smc_iteration] = Setting(:smc_iteration, 1, false, "iter",
        "The iteration index for the number of times smc has been run on the same data vintage. " *
        "Primarily for numerical accuracy/testing purposes.")
    settings[:previous_data_vintage] = Setting(:previous_data_vintage, vint,
        "The old data vintage from which to start SMC when using a tempered update.")

    # Alternative policy
    ## TODO: delete this once alternative policies have been completely ported to using regime_eqcond_info
    # baseline_policy = AltPolicy(:historical, eqcond, solve, forecast_init = identity)
    # settings[:alternative_policy] = Setting(:alternative_policy, baseline_policy)

    # Temporary alternative policy defaults to zero rate
    settings[:temporary_altpolicy_names] = Setting(:temporary_altpolicy_names, Symbol[:zlb_rule])

    return settings
end

"""
```
default_test_settings!(m::AbstractDSGEModel)
```

The following Settings are constructed, initialized and added to
`m.test_settings`. Their purposes are identical to those in
`m.settings`, but these values are used to test DSGE.jl.

### I/O Locations and identifiers
- `saveroot::Setting{String}`: A temporary directory in /tmp/
- `dataroot::Setting{String}`: dsgeroot/test/reference/
- `data_vintage::Setting{String}`: \"_REF\"

### Metropolis-Hastings
- `n_mh_simulations::Setting{Int}`: 100
- `n_mh_blocks::Setting{Int}`: 1
- `n_mh_burn::Setting{Int}`: 0
- `mh_thin::Setting{Int}`: 1
"""
function default_test_settings!(m::AbstractDSGEModel)

    test = m.test_settings

    # I/O
    dataroot = normpath(joinpath(dirname(@__FILE__), "..", "test", "reference", "input_data"))
    saveroot = mktempdir()

    # General
    test[:saveroot] = Setting(:saveroot, saveroot,
        "Where to write files when in test mode")
    test[:dataroot] = Setting(:dataroot, dataroot,
        "Location of input files when in test mode" )
    test[:data_vintage] = Setting(:data_vintage, "REF", true, "vint",
        "Reference data identifier")
    test[:cond_vintage] = Setting(:cond_vintage, "REF",
        "Vintage of conditional data")
    test[:use_parallel_workers] = Setting(:use_parallel_workers, false, false, "parw",
        "Use available parallel workers in computations")
    test[:n_hessian_test_params] = Setting(:n_hessian_test_params, 3, false, "mhfp",
        "Max number of free params for which to calculate Hessian")

    # Metropolis-Hastings
    test[:n_mh_simulations] = Setting(:n_mh_simulations, 100, false, "nsim",
        "Number of parameter draws per block for testing Metropolis-Hastings")
    test[:n_mh_blocks] = Setting(:n_mh_blocks, 1, false, "nblc",
        "Number of blocks to draw parameters for testing Metropolis-Hastings")
    test[:n_mh_burn] = Setting(:n_mh_burn, 0, false, "nbrn",
        "Number of burn-in blocks for testing Metropolis-Hastings")
    test[:mh_thin] = Setting(:mh_thin, 1, false, "thin",
        "Thinning step for testing Metropolis-Hastings")

    # Forecast
    test[:date_forecast_start] = Setting(:date_forecast_start, quartertodate("2015-Q4"),
        "Start date of forecast period")
    test[:forecast_horizons] = Setting(:forecast_horizons, 2,
        "Number of periods to forecast ahead")
    test[:forecast_jstep] = Setting(:forecast_jstep, 1,
        "Forecast thinning step (in addition to MH thinning step")
    test[:impulse_response_horizons] = Setting(:impulse_response_horizons, 2,
        "Number of periods for which to calculate an impulse response")

    return test
end


text:
"""
Module specifying the structure of a two-player cribbage game.

The code here represents the reference implementation for other versions of the
    game (e.g.: three-player, four-player).
"""
module TwoPlayer

import Cribbage: CribbageGame, GameState, Player, nextplayer, UnexpectedPlayerException

export incrementscore, isgameover, players, reportscore, rotatedealer, TwoPlayerGame, TwoPlayerGameState

struct TwoPlayerGame <: CribbageGame
    p₁ :: Player
    p₂ :: Player
    winningscore :: Int8
    function TwoPlayerGame(
        p₁::Player, p₂::Player, winningscore::Integer=121)
        @assert(winningscore == 121 || winningscore == 61, "winningscore is either 61 or 121")
        @assert(p₁ ≢ p₂)
        new(p₁, p₂, winningscore)
    end
end

struct TwoPlayerGameState <: GameState
    game :: TwoPlayerGame
    dealer :: Player
    s₁ :: Int8
    s₂ :: Int8
    function TwoPlayerGameState(
        game::TwoPlayerGame, dealer::Player, s₁::Integer=0, s₂::Integer=0)
        @assert(dealer ∈ players(game), "dealer is not one of the players of the game")
        @assert(0 ≤ s₁ ≤ game.winningscore, "s₁ is an invalid score")
        @assert(0 ≤ s₂ ≤ game.winningscore, "s₂ is an invalid score")
        new(game, dealer, s₁, s₂)
    end
end

"""
    numplayers(game)
Returns the number of players in `game`.
"""
function numplayers(game::TwoPlayerGame)::Int8
    return 2
end

"""
    cardsdealt(game)
Returns the number of cards to be dealt to each player in `game`.
"""
function cardsdealt(game::TwoPlayerGame)::Int8
    return 6
end

"""
    cardstocrib(game)
Returns the number of cards to seed the crib with during each hand of `game`.
"""
function cardstocrib(game::TwoPlayerGame)::Int8
    return 0
end

"""
    players(game)
Returns an array of players in `game`.
"""
function players(game::TwoPlayerGame)::Array{Player,1}
    return [ game.p₁, game.p₂ ]
end

"""
    incrementscore(gamestate, player, amount)
Returns a new game state, starting from `gamestate`, in which the score for
    `player` has been incremented by `amount`.
"""
function incrementscore(
        gamestate::TwoPlayerGameState,
        player::Player,
        amount::Integer
    )
    @assert(amount ≥ 0)
    @debug "Incrementing the score for $(player.name) by $(amount)"
    if player === gamestate.game.p₁
        return TwoPlayerGameState(
            gamestate.game,
            gamestate.dealer,
            min(gamestate.s₁ + amount, gamestate.game.winningscore),
            gamestate.s₂
        )
    elseif player === gamestate.game.p₂
        return TwoPlayerGameState(
            gamestate.game,
            gamestate.dealer,
            gamestate.s₁,
            min(gamestate.s₂ + amount, gamestate.game.winningscore)
        )
    else
        throw(UnexpectedPlayerException("$(player.name) is not one of the players in this game"))
    end
end

"""
    isgameover(gamestate)
Returns true if any player has a winning score in `gamestate`.
"""
function isgameover(gamestate::TwoPlayerGameState)::Bool
    return gamestate.s₁ == gamestate.game.winningscore || gamestate.s₂ == gamestate.game.winningscore
end

"""
    reportscore(gamestate)
Returns a string representing the score in `gamestate`.
"""
function reportscore(gamestate::TwoPlayerGameState)::String
    return "$(gamestate.s₁) - $(gamestate.s₂)"
end

"""
    rotatedealer(gamestate)
Returns a new game state with the next dealer selected.
"""
function rotatedealer(gamestate::TwoPlayerGameState)::TwoPlayerGameState
    return TwoPlayerGameState(
        gamestate.game,
        nextplayer(players(gamestate.game), gamestate.dealer),
        gamestate.s₁,
        gamestate.s₂
    )
end

end  # module TwoPlayer


text:
using CUDA
using Adapt

mutable struct Test_Cfg
    A::CuArray
    B::CuArray
end
function Adapt.adapt_structure(to, test::Test_Cfg)
    A = Adapt.adapt_structure(to, test.A)
    B = Adapt.adapt_structure(to, test.B)
    Test_Cfg(A, B)
end
struct Kernel_test_Cfg
    A::Float32
    B::Int32
end
function Adapt.adapt_structure(to, test::Kernel_test_Cfg)
    A = Adapt.adapt_structure(to, test.A)
    B = Adapt.adapt_structure(to, test.B)
    Kernel_test_Cfg(A, B)
end

function kernel_test(test)

    @cuprintln(test.A)
    @cuprintln(test.B)
    return
end
#test = Kernel_test_Cfg(2.0f0,3)
# @cuda kernel_test(test)

function maze_kernel(maze,offset,maze_columns,maze_rows)
    total_amount_of_cells = maze_columns * maze_rows

    x_coordinate_stack = @cuDynamicSharedMem(Int32,total_amount_of_cells,offset)
    y_coordinate_stack = @cuDynamicSharedMem(Int32,total_amount_of_cells,offset+sizeof(x_coordinate_stack))
    neighbours = @cuDynamicSharedMem(Int32,4,offset+sizeof(x_coordinate_stack)+sizeof(y_coordinate_stack))
    @inbounds fill!(maze,1)
    cell_x_coordinate = 1
    cell_y_coordinate = 1
    amount_of_cells_visited = 1
    cell_stack_index = 1

    while amount_of_cells_visited < total_amount_of_cells

        @inbounds fill!(neighbours,0)

        if cell_y_coordinate-1 >= 1
            for i in 1:4
                if maze[cell_y_coordinate-1,cell_x_coordinate,i,blockIdx().x] == 0
                    break
                end
                if i == 4
                    neighbours[1] = 1
                end
            end
        end
        if cell_x_coordinate+1 <= maze_columns
            for i in 1:4
                if maze[cell_y_coordinate,cell_x_coordinate+1,i,blockIdx().x] == 0
                    break
                end
                if i == 4
                    neighbours[2] = 1
                end
            end
        end
        if cell_y_coordinate+1 <= maze_rows
            for i in 1:4
                if maze[cell_y_coordinate+1,cell_x_coordinate,i,blockIdx().x] == 0
                    break
                end
                if i == 4
                    neighbours[3] = 1
                end
            end
        end
        if cell_x_coordinate-1 >= 1
            for i in 1:4
                if maze[cell_y_coordinate,cell_x_coordinate-1,i,blockIdx().x] == 0
                    break
                end
                if i == 4
                    neighbours[4] = 1
                end
            end
        end
        if neighbours[1] == 0 && neighbours[2] == 0 && neighbours[3] == 0 && neighbours[4] == 0
            cell_stack_index = cell_stack_index - 1
            cell_x_coordinate = x_coordinate_stack[cell_stack_index]
            cell_y_coordinate = y_coordinate_stack[cell_stack_index]
            continue
        end

        move_x_coordinate = 0
        move_y_coordinate = 0
        rand_index = (abs(rand(Int32)) % 4)
        for i in 1:4
            index = ((rand_index+i) % 4) + 1
            if neighbours[index] == 1
                if index == 1
                    move_y_coordinate = -1
                    break
                end
                if index == 2
                    move_x_coordinate = 1
                    break
                end
                if index == 3
                    move_y_coordinate = 1
                    break
                end
                if index == 4
                    move_x_coordinate = -1
                    break
                end
            end
        end
        if move_y_coordinate == -1
            maze[cell_y_coordinate,cell_x_coordinate,1,blockIdx().x] = 0
            maze[cell_y_coordinate+move_y_coordinate,cell_x_coordinate,3,blockIdx().x] = 0
        end
        if move_x_coordinate == 1
            maze[cell_y_coordinate,cell_x_coordinate,2,blockIdx().x] = 0
            maze[cell_y_coordinate,cell_x_coordinate+move_x_coordinate,4,blockIdx().x] = 0
        end
        if move_y_coordinate == 1
            maze[cell_y_coordinate,cell_x_coordinate,3,blockIdx().x] = 0
            maze[cell_y_coordinate+move_y_coordinate,cell_x_coordinate,1,blockIdx().x] = 0
        end
        if move_x_coordinate == -1
            maze[cell_y_coordinate,cell_x_coordinate,4,blockIdx().x] = 0
            maze[cell_y_coordinate,cell_x_coordinate+move_x_coordinate,2,blockIdx().x] = 0
        end

        #step5: add origin cell to stack
        x_coordinate_stack[cell_stack_index] = cell_x_coordinate
        y_coordinate_stack[cell_stack_index] = cell_y_coordinate
        cell_stack_index = cell_stack_index +1
        #step6: set coordinates to new cell
        cell_x_coordinate = cell_x_coordinate + move_x_coordinate
        cell_y_coordinate = cell_y_coordinate + move_y_coordinate
        amount_of_cells_visited = amount_of_cells_visited +1
    end
    return
end

function kernel_env_step(maze,action,input,environment_config_array)
    maze_columns = 5
    maze_rows = 5
    maze_cell_size = 80
    agent_radius = 12
    agent_movement_range = 10.0f0
    point_radius = 8
    reward_per_collected_positive_point = 500.00f0
    reward_per_collected_negative_point = -700.00f0

    screen_height = maze_rows * maze_cell_size
    screen_width = maze_columns * maze_cell_size

    agent_x_coordinate = environment_config_array[1] + clamp(floor(action[1] * env_cfg.agent_movement_range),-env_cfg.agent_movement_range,env_cfg.agent_movement_range)
    agent_x_coordinate = min(agent_x_coordinate,screen_width - env_cfg.agent_radius) # Eastern border
    agent_x_coordinate = max(agent_x_coordinate,env_cfg.agent_radius) # Western Border

    agent_y_coordinate = environment_config_array[2] + clamp(floor(action[2] * env_cfg.agent_movement_range),-env_cfg.agent_movement_range,env_cfg.agent_movement_range)
    agent_y_coordinate = min(agent_y_coordinate,screen_height - env_cfg.agent_radius) # Eastern border
    agent_y_coordinate = max(agent_y_coordinate,env_cfg.agent_radius) # Western Border

    # Get cell indizes of agents current position
    cell_x = convert(Int32,ceil(agent_x_coordinate / env_cfg.maze_cell_size))
    overlay_y = convert(Int32,ceil(agent_y_coordinate / env_cfg.maze_cell_size))
    cell_y = maze_rows - (overlay_y-1)

    # Get edge coordinates of current cell
    x_left = env_cfg.maze_cell_size * (cell_x - 1)
    x_right = env_cfg.maze_cell_size * cell_x
    y_bottom = env_cfg.maze_cell_size * (overlay_y - 1)
    y_top = env_cfg.maze_cell_size * overlay_y

    if maze[cell_y,cell_x,1] == 1 #check for Northern Wall
        agent_y_coordinate = min(agent_y_coordinate,y_top - env_cfg.agent_radius)
    end
    if maze[cell_y,cell_x,3] == 1 #check for Southern Wall
        agent_y_coordinate = max(agent_y_coordinate,y_bottom + env_cfg.agent_radius)
    end
    if maze[cell_y,cell_x,2] == 1 #check for Eastern Wall
        agent_x_coordinate = min(agent_x_coordinate,x_right - env_cfg.agent_radius)
    end
    if maze[cell_y,cell_x,4] == 1 #check for Western Wall
        agent_x_coordinate = max(agent_x_coordinate,x_left + env_cfg.agent_radius)
    end
    # Check agent collision with top-left edge (prevents sneaking through the edge)
    if (agent_x_coordinate - x_left < env_cfg.agent_radius) && (y_top - agent_y_coordinate < env_cfg.agent_radius)
        agent_x_coordinate = x_left + env_cfg.agent_radius
        agent_y_coordinate = y_top - env_cfg.agent_radius
    end
    # Check agent collision with top-right edge (prevents sneaking through the edge)
    if (x_right - agent_x_coordinate < env_cfg.agent_radius) && (y_top - agent_y_coordinate < env_cfg.agent_radius)
        agent_x_coordinate = x_right - env_cfg.agent_radius
        agent_y_coordinate = y_top - env_cfg.agent_radius
    end
    # Check agent collision with bottom-right edge (prevents sneaking through the edge)
    if (x_right - agent_x_coordinate < env_cfg.agent_radius) && (agent_y_coordinate - y_bottom < env_cfg.agent_radius)
        agent_x_coordinate = x_right - env_cfg.agent_radius
        agent_y_coordinate = y_bottom + env_cfg.agent_radius
    end
    # Check agent collision with bottom-left edge (prevents sneaking through the edge)
    if (agent_x_coordinate - x_left < env_cfg.agent_radius) && (agent_y_coordinate - y_bottom < env_cfg.agent_radius)
        agent_x_coordinate = x_left + env_cfg.agent_radius
        agent_y_coordinate = y_bottom + env_cfg.agent_radius
    end
end

maze_columns = 5
maze_rows = 5
number_individuals = 10
maze = fill(0,(maze_rows,maze_columns,4,number_individuals))

maze_gpu = CuArray(maze)
@cuda blocks=number_individuals shmem=2*maze_columns*maze_rows+4 maze_kernel(maze_gpu,0,maze_columns,maze_rows)
CUDA.synchronize()
display(Array(maze_gpu))


text:
using JuMP
using Gurobi

""" generate an array or matrix full of 1 or -1 independently with prob p
size is either an integer (if we want to generate an array) or a tuple (if we want
a matrix)
"""
function pm_one(size, p::Float64)
    bernoulli_array = rand(size...) .< p  # pointwise compare random numbers to p
    2 .* bernoulli_array .- 1
end

""" generate a dataset of n people, with approximately pn of them having +1 and the rest having -1 """
function make_data(n::Int, p::Float64)
    pm_one(n, p)
end

""" generate num_queries random queries whose coefficients are +1 or -1
Compute noisy answers for them by adding uniform noise, where the noise can be at most bound.
"""
function noisy_queries(num_queries::Int, data::Array, bound::Float64)
    n = length(data)
    size = (num_queries, n)
    query_matrix = pm_one(size, 0.5)
    answers = query_matrix * data
    noisy_answers = answers + bound .* (2 .* rand(num_queries) .- 1)
    (query_matrix, noisy_answers)
end

"""
Attempt to reconstruct the database using the noisy queries
"""
function attack_db(query_matrix, noisy_answers, bound)
    model = Model(Gurobi.Optimizer)
    num_queries, n = size(query_matrix)
    @variable(model, -1 <= x[1:n] <= 1)
    @constraint(model, query_matrix * x .<= noisy_answers .+ bound)
    @constraint(model, noisy_answers .- bound .<= query_matrix * x)
    @objective(model, Min, 1)
    optimize!(model)
    2 .* (value.(x) .>= 0) .- 1
end

"""  When you call this function, you must name the parameters you want to provide
e.g., example(n=100)
"""
function example(;n=100, bound=5.0, num_queries=200)
    data = make_data(n, 0.5)
    queries, answers = noisy_queries(num_queries, data, bound)
    x = attack_db(queries, answers, bound)
    num_correct = sum(x .== data)
    println("\n*****\n Got $(num_correct) correct out of $(n)\n****")
    queries, answers, data, x
end

max_stars_repo_name: dkifer/dpcourse

max_issues_repo_name: dkifer/dpcourse

max_forks_repo_name: dkifer/dpcourse


text:
@doc """
    init_jsi(modeldirs, noises, commands)

!!! note "Advanced"
        init_jsi(modeldirs, noises, commands; <keyword arguments>)

Initialize the package JustSayIt.

# Arguments
- `modeldirs::Dict{String, String}`: the directories where the speech recognition models to be used are located.
- `noises::Dict{String, <:AbstractArray{String}}`: for each model, an array of strings with noises (tokens that are to be ignored in the speech as, e.g., "huh").
- `commands::Dict{String, Function}`: the commands to be recognized with their mapping to a function.
!!! note "Advanced keyword arguments"
    - `vosk_log_level::Integer=-1`: the vosk log level (see in the Vosk documentation for details).

See also: [`finalize_jsi`](@ref)
"""
init_jsi

let
    global init_jsi, command, command_names, model, noises, noises_names, recognizer, controller, set_controller
    _commands::Dict{String, Function}                                           = Dict{String, Function}()
    _models::Dict{String, PyObject}                                             = Dict{String, PyObject}()
    _noises::Dict{String, <:AbstractArray{String}}                              = Dict{String, Array{String}}()
    _recognizers::Dict{String, PyObject}                                        = Dict{String, PyObject}()
    _controllers::Dict{String, PyObject}                                        = Dict{String, PyObject}()
    command(name::AbstractString)::Function                                     = _commands[name]
    command_names()::Base.KeySet{String, Dict{String, Function}}                = keys(_commands)
    model(name::AbstractString=DEFAULT_MODEL_NAME)::PyObject                    = _models[name]
    noises(modelname::AbstractString)                                           = _noises[modelname]  # NOTE: no return value declaration as would be <:AbstractArray{String} which is not possible.
    noises_names()::Base.KeySet{String, Dict{String, <:AbstractArray{String}}}  = keys(_noises)
    recognizer(id::AbstractString)::PyObject                                    = _recognizers[id]
    controller(name::AbstractString)::PyObject                                  = if (name in keys(_controllers)) return _controllers[name] else error("The controller for $name is not available as it has not been set up in init_jsi.") end
    set_controller(name::AbstractString, c::PyObject)                           = (_controllers[name] = c; return)

    function init_jsi(modeldirs::Dict{String, String}, noises::Dict{String, <:AbstractArray{String}}, commands::Dict{String, Function}; vosk_log_level::Integer=-1)
        Vosk.SetLogLevel(vosk_log_level)

        # Validate and store the commands, adding the help command to it.
        if haskey(commands, COMMAND_NAME_EXIT) error("the command name $COMMAND_NAME_EXIT is reserved for exiting command interpretation. Please choose another command name for your command.") end
        if haskey(commands, COMMAND_NAME_SLEEP) error("the command name $COMMAND_NAME_SLEEP is reserved for putting JustSayIt to sleep. Please choose another command name for your command.") end
        if haskey(commands, COMMAND_NAME_AWAKE) error("the command name $COMMAND_NAME_AWAKE is reserved for awaking JustSayIt. Please choose another command name for your command.") end
        _commands = commands

        # Verify that DEFAULT_MODEL_NAME is listed in modeldirs and noises. Set the values for the  other surely required models (i.e. which are used in the Commands submodule) to the same as the default if not available.
        if !haskey(modeldirs, DEFAULT_MODEL_NAME) error("a directory for 'default' model is mandatory.") end
        if !haskey(noises, DEFAULT_MODEL_NAME) error("a noises list for 'default' model is mandatory (e.g. `NOISES_ENGLISH`).") end
        if !haskey(modeldirs, TYPE_MODEL_NAME)
            modeldirs[TYPE_MODEL_NAME] = modeldirs[DEFAULT_MODEL_NAME]
            @warn("no modeldir given for model \"type\" - falling back to model 'default' for typing.")
        end
        if !haskey(noises, TYPE_MODEL_NAME)
            noises[TYPE_MODEL_NAME] = noises[DEFAULT_MODEL_NAME]
            @warn("no noises given for model \"type\" - falling back to noises configuration of model 'default' when typing.")
        end
        if haskey(modeldirs, "") error("an empty string is not valid as model identifier.") end
        _noises = noises

        # Set up a default recognizer for each model as well as the command name recognizer.
        for modelname in keys(modeldirs)
            if !isdir(modeldirs[modelname])
                tilde_errmsg = ""
                if startswith(modeldirs[modelname], "~") tilde_errmsg = " Tildes ('~') in strings are not expanded for portability reasons (consider using 'homedir()')." end
                error("directory $(modeldirs[modelname]) does not exist.$tilde_errmsg")
            end
            _models[modelname] = Vosk.Model(modeldirs[modelname])
            _recognizers[modelname] = Vosk.KaldiRecognizer(model(modelname), SAMPLERATE)
            if modelname == DEFAULT_MODEL_NAME
                grammar = json([keys(commands)..., COMMAND_NAME_EXIT, COMMAND_NAME_SLEEP, COMMAND_NAME_AWAKE, noises[modelname]..., UNKNOWN_TOKEN])
                _recognizers[COMMAND_RECOGNIZER_ID] = Vosk.KaldiRecognizer(model(modelname), SAMPLERATE, grammar)
            end
        end

        # Set up the special recognizers defined by voiceargs.
        for f_name in voicearg_f_names()
            for voicearg in keys(voiceargs(f_name))
                kwargs = voiceargs(f_name)[voicearg]
                if haskey(kwargs, :model) && !haskey(modeldirs, kwargs[:model]) error("no directory was given for the model $(kwargs[:model]) required for voicearg $voicearg in function $f_name.") end # NOTE: this error should only ever occur for user defined @voicearg functions; all funtions in submodule Commands must use
                if haskey(kwargs, :valid_input)
                    modelname = haskey(kwargs, :model) ? kwargs[:model] : DEFAULT_MODEL_NAME
                    grammar = json([kwargs[:valid_input]..., noises[modelname]..., UNKNOWN_TOKEN])
                    set_recognizer(f_name, voicearg, Vosk.KaldiRecognizer(model(modelname), SAMPLERATE, grammar))
                end
            end
            @debug "Voiceargs of function `$f_name`:" voiceargs(f_name)
        end

        # Call the garbage collector to avoid having it running right after starting to use JustSayIt.
        GC.gc()
    end

end

_noises(args...) = noises(args...)