text:
# ref: https://github.com/apache/parquet-format/blob/master/Encodings.md

macro bitwidth(i)
    quote
        ceil(Int, log(2, $(esc(i))+1))
    end
end
macro bit2bytewidth(i)
    quote
        ceil(Int, $(esc(i))/8)
    end
end
macro byt2itype(i)
    quote
        ($(esc(i)) <= 4) ? Int32 : ($(esc(i)) <= 8) ? Int64 : Int128
    end
end
macro byt2uitype(i)
    quote
        ($(esc(i)) <= 4) ? UInt32 : ($(esc(i)) <= 8) ? UInt64 : UInt128
    end
end
macro byt2uitype_small(i)
    quote
        ($(esc(i)) <= 1) ? UInt8 : ($(esc(i)) <= 2) ? UInt16 : ($(esc(i)) <= 4) ? UInt32 : ($(esc(i)) <= 8) ? UInt64 : UInt128
    end
end

const MSB = 0x80
const MASK7 = 0x7f
const MASK8 = 0xff
const MASK3 = 0x07
function MASKN(nbits::UInt8)
    byte_width = @bit2bytewidth(nbits)
    type_small = @byt2uitype_small(byte_width)
    MASKN(nbits, type_small)
end
function MASKN(nbits::UInt8, ::Type{T}) where {T}
    O = convert(T, 0x1)
    (O << nbits) - O
end

#read_fixed(io::IO, typ::Type{UInt32}) = _read_fixed(io, convert(UInt32,0), 4)
#read_fixed(io::IO, typ::Type{UInt64}) = _read_fixed(io, convert(UInt64,0), 8)
read_fixed(io::IO, typ::Type{Int32}) = reinterpret(Int32, _read_fixed(io, convert(UInt32,0), 4))
#read_fixed(io::IO, typ::Type{Int64}) = reinterpret(Int64, _read_fixed(io, convert(UInt64,0), 8))
#read_fixed(io::IO, typ::Type{Int128}) = reinterpret(Int128, _read_fixed(io, convert(UInt128, 0), 12))   # INT96: 12 bytes little endian
#read_fixed(io::IO, typ::Type{Float32}) = reinterpret(Float32, _read_fixed(io, convert(UInt32,0), 4))
#read_fixed(io::IO, typ::Type{Float64}) = reinterpret(Float64, _read_fixed(io, convert(UInt64,0), 8))
function _read_fixed(io::IO, ret::T, N::Int) where {T <: Unsigned}
    for n in 0:(N-1)
        byte = convert(T, read(io, UInt8))
        ret |= (byte << *(8,n))
    end
    ret
end

function _read_fixed_bigendian(io::IO, ret::T, N::Int) where {T <: Unsigned}
    for n in (N-1):0
        byte = convert(T, read(io, UInt8))
        ret |= (byte << *(8,n))
    end
end

mutable struct InputState
    data::Vector{UInt8}
    offset::Int
end

mutable struct OutputState{T}
    data::Vector{T}
    offset::Int
end

function OutputState(::Type{T}, size) where {T}
    arr = Array{T}(undef, size)
    OutputState{T}(arr, 0)
end

function ensure_additional_size(iostate, additional_size)
    needed_size = iostate.offset + additional_size
    if length(iostate.data) < needed_size
        resize!(iostate.data, needed_size)
    end
    nothing
end

function reset_to_size(iostate, size)
    iostate.offset = 0
    (length(iostate.data) == size) || resize!(iostate.data, size)
    nothing
end

read_fixed(inp::InputState, typ::Type{UInt32}) = _read_fixed(inp, convert(UInt32,0), 4)
read_fixed(inp::InputState, typ::Type{UInt64}) = _read_fixed(inp, convert(UInt64,0), 8)
read_fixed(inp::InputState, typ::Type{Int32}) = reinterpret(Int32, _read_fixed(inp, convert(UInt32,0), 4))
read_fixed(inp::InputState, typ::Type{Int64}) = reinterpret(Int64, _read_fixed(inp, convert(UInt64,0), 8))
read_fixed(inp::InputState, typ::Type{Int128}) = reinterpret(Int128, _read_fixed(inp, convert(UInt128, 0), 12))
read_fixed(inp::InputState, typ::Type{Float32}) = reinterpret(Float32, _read_fixed(inp, convert(UInt32,0), 4))
read_fixed(inp::InputState, typ::Type{Float64}) = reinterpret(Float64, _read_fixed(inp, convert(UInt64,0), 8))
function _read_fixed(inp::InputState, ret::T, N::Int) where {T <: Unsigned}
    data = inp.data
    offset = inp.offset
    for n in 0:(N-1)
        byte = convert(T, data[1+offset])
        offset += 1
        ret |= (byte << *(8,n))
    end
    inp.offset = offset
    ret
end

function _read_varint(inp::InputState, ::Type{T}) where {T <: Integer}
    data = inp.data
    offset = inp.offset
    res = zero(T)
    n = 0
    byte = UInt8(MSB)
    while (byte & MSB) != 0
        byte = data[1+offset]
        offset += 1
        res |= (convert(T, byte & MASK7) << (7*n))
        n += 1
    end
    inp.offset = offset

    # in case of overflow, consider it as missing field and return default value
    if (n-1) > sizeof(T)
        #@debug("overflow reading $T. returning 0")
        return zero(T)
    end
    res
end

# parquet types:   BOOLEAN, INT32, INT64,  INT96,   FLOAT,  DOUBLE, BYTE_ARRAY, FIXED_LEN_BYTE_ARRAY
# enum values:           0,     1,     2,      3,       4,       5,          6,                    7
const PLAIN_JTYPES = (Bool, Int32, Int64, Int128, Float32, Float64,      UInt8,               UInt8)

# read plain encoding (PLAIN = 0)
read_plain_byte_array(inp::InputState) = read_plain_byte_array(inp::InputState, read_fixed(inp, Int32))
function read_plain_byte_array(inp::InputState, count::Int32)
    arr = inp.data[(1+inp.offset):(inp.offset+count)]  # TODO: Return subarr?
    inp.offset += count
    arr
end

# read plain values _Type.BOOLEAN
function read_plain_values(inp::InputState, out::OutputState{Bool}, count::Int32)
    #@debug("reading plain values", type=Bool, count=count)
    read_bitpacked_booleans(inp, out, count)
end
# read plain values _Type.BYTE_ARRAY
function read_plain_values(inp::InputState, out::OutputState{Vector{UInt8}}, count::Int32)
    # _Type.FIXED_LEN_BYTE_ARRAY is most likely same as byte array
    #@debug("reading plain values", type=Vector{UInt8}, count=count)
    arr = out.data
    offset = out.offset
    @assert (offset + count) <= length(arr)
    @inbounds for i in 1:count
        arr[i+offset] = read_plain_byte_array(inp)
    end
    out.offset += count
    nothing
end
function read_plain_values(inp::InputState, out::OutputState{T}, count::Int32, converter_fn::Function, storage_type::Int32) where T <: Union{Decimal,Float64,Int16,Int32,Int64,Int128}
    arr = out.data
    offset = out.offset
    @assert (offset + count) <= length(arr)
    if storage_type === _Type.FIXED_LEN_BYTE_ARRAY
        elem_bytes_len = Int32((length(inp.data) - inp.offset) / length(arr))
        #@debug("reading decimal plain values", offset, count, length(arr), length(inp.data), inp.offset, elem_bytes_len)
        @inbounds for i in 1:count
            arr[i+offset] = converter_fn(read_plain_byte_array(inp, elem_bytes_len))
        end
    elseif storage_type === _Type.INT64
        @inbounds for i in 1:count
            arr[i+offset] = converter_fn(read_fixed(inp, Int64))
        end
    elseif storage_type === _Type.INT32
        @inbounds for i in 1:count
            arr[i+offset] = converter_fn(read_fixed(inp, Int32))
        end
    else
        error("unsupported storage type $(storage_type) for $T")
    end
    out.offset += count
    nothing
end
function read_plain_values(inp::InputState, out::OutputState{String}, count::Int32, converter_fn::Function, storage_type::Int32)
    #@debug("reading plain values", type=Vector{UInt8}, count=count)
    arr = out.data
    offset = out.offset
    @assert (offset + count) <= length(arr)
    if storage_type === _Type.BYTE_ARRAY
        @inbounds for i in 1:count
            arr[i+offset] = converter_fn(read_plain_byte_array(inp))
        end
    elseif storage_type === _Type.FIXED_LEN_BYTE_ARRAY
        elem_bytes_len = Int32((length(inp.data) - inp.offset) / length(arr))
        @inbounds for i in 1:count
            arr[i+offset] = converter_fn(read_plain_byte_array(inp, elem_bytes_len))
        end
    else
        error("unsupported storage type $(storage_type) for String")
    end
    out.offset += count
    nothing
end
# read_plain_values of type T
function read_plain_values(inp::InputState, out::OutputState{T}, count::Int32) where {T}
    #@debug("reading plain values", type=T, count=count)
    arr = out.data
    offset = out.offset
    @assert (offset + count) <= length(arr)
    @inbounds for i in 1:count
        arr[i+offset] = read_fixed(inp, T)
    end
    #@debug("read $(length(arr)) plain values")
    out.offset += count
    nothing
end
function read_plain_values(inp::InputState, out::OutputState{DateTime}, count::Int32, converter_fn::Function, storage_type::Int32)
    #@debug("reading plain values", type=T, count=count)
    arr = out.data
    offset = out.offset
    @assert (offset + count) <= length(arr)
    if storage_type === _Type.INT96
        @inbounds for i in 1:count
            arr[i+offset] = converter_fn(read_fixed(inp, Int128))
        end
    else
        error("unsupported storage type $(storage_type) for DateTime")
    end
    #@debug("read $(length(arr)) plain values")
    out.offset += count
    nothing
end

function read_bitpacked_booleans(inp::InputState, out::OutputState{Bool}, count::Int32)
    #@debug("reading bitpacked booleans", count)
    arrpos = 1
    bits = UInt8(0)
    bitpos = 9

    data = inp.data
    data_offset = inp.offset

    arr = out.data
    offset = out.offset
    @assert (offset+count) <= length(arr)

    @inbounds while arrpos <= count
        if bitpos > 8
            bits = data[1+data_offset]
            data_offset += 1
            #@debug("bits", bits, bitstring(bits))
            bitpos = 1
        end
        arr[arrpos+offset] = Bool(bits & 0x1)
        arrpos += 1
        bits >>= 1
        bitpos += 1
    end
    out.offset += count
    inp.offset = data_offset
    nothing
end

# read data dictionary (RLE_DICTIONARY = 8, or PLAIN_DICTIONARY = 2 in a data page)
function read_data_dict(inp::InputState, count::Int32)
    bits = inp.data[inp.offset+=1]
    #@info("bits", bits)
    byte_width = @bit2bytewidth(bits)
    typ = @byt2itype(byte_width)
    out = OutputState(typ, count)

    #@info("reading read_hybrid", count, read_len)
    read_hybrid(inp, out, count, bits, byte_width; read_len=false)
    out.data
end

# read RLE or bit backed format (RLE = 3)
function read_hybrid(inp::InputState, out::OutputState{T}, count::Int32, bits::UInt8, byt::Int; read_len::Bool=true) where {T}
    len = Int32(0)
    if read_len
        len = read_fixed(inp, Int32)
    end
    #@debug("reading hybrid data", len, count, bits)
    mask = MASKN(bits)
    arr = out.data
    arrpos = 1
    offset = out.offset - 1    # to counter arrpos starting at 1
    while arrpos <= count
        runhdr = _read_varint(inp, Int)
        isbitpack = ((runhdr & 0x1) == 0x1)
        runhdr >>= 1
        nrunhdrbits = runhdr * 8
        nitems = min(isbitpack ? nrunhdrbits : runhdr, count - arrpos + 1)

        if isbitpack
            runcount = min(nrunhdrbits, length(arr)-offset-arrpos)
            read_bitpacked_run(inp, out, runcount, bits, byt, mask)
            if nrunhdrbits > runcount
                nbytes_to_skip = div(nrunhdrbits - runcount, 8)
                @debug("skipping trailing bytes in bitpacked run", nbytes_to_skip)
                inp.offset += nbytes_to_skip
            end
            #out.offset += runcount
        else # rle
            read_type = @byt2uitype(byt)
            read_rle_run(inp, out, nitems, bits, byt, read_type)
            #out.offset += nitems
        end
        arrpos += nitems
    end
    nothing
end

function read_rle_run(inp::InputState, out::OutputState{T}, count::Int, bits::UInt8, byt::Int, read_type::Type{V}) where {T,V}
    #@debug("read_rle_run", count, T, bits, byt)
    rawval = _read_fixed(inp, zero(V), byt)
    val = reinterpret(T, rawval)
    arr = out.data
    offset = out.offset
    @assert length(arr) >= (count+offset)
    @inbounds for idx in 1:count
        arr[idx+offset] = val
    end
    out.offset += count
    nothing
end

function read_bitpacked_run(inp::InputState, out::OutputState{T}, count::Int, bits::UInt8, byt::Int, mask::V=MASKN(bits)) where {T,V}
    bitbuff = zero(V)
    nbitsbuff = UInt8(0)
    shift = UInt8(0)

    data = inp.data
    arr = out.data
    offset = out.offset
    arridx = 1
    dataidx = 1 + inp.offset
    while arridx <= count
        #@debug("arridx:$arridx nbitsbuff:$nbitsbuff shift:$shift bits:$bits")
        if nbitsbuff > 0
            # we have leftover bits, which must be appended
            if nbitsbuff < bits
                # but only append if we need to read more in this cycle
                @inbounds arr[arridx+offset] = bitbuff & MASKN(nbitsbuff, V)
                shift = nbitsbuff
                nbitsbuff = UInt8(0)
                bitbuff = zero(V)
            end
        end

        # fill buffer
        while (nbitsbuff + shift) < bits
             # shift 8 bits and read directly into bitbuff
            bitbuff |= (V(data[dataidx]) << nbitsbuff)
            dataidx += 1
            nbitsbuff += UInt8(8)
        end

        # set values
        while ((nbitsbuff + shift) >= bits) && (arridx <= count)
            if shift > 0
                remshift = bits - shift
                #@debug("setting part from bitbuff nbitsbuff:$nbitsbuff, shift:$shift, remshift:$remshift")
                arr[arridx+offset] |= convert(T, (bitbuff << shift) & mask)
                bitbuff >>= remshift
                nbitsbuff -= remshift
                shift = UInt8(0)
            else
                #@debug("setting all from bitbuff nbitsbuff:$nbitsbuff")
                arr[arridx+offset] = convert(T, bitbuff & mask)
                bitbuff >>= bits
                nbitsbuff -= bits
            end
            arridx += 1
        end
    end
    inp.offset = dataidx - 1
    out.offset += count
    nothing
end

# read bit packed in deprecated format (BIT_PACKED = 4)
function read_bitpacked_run_old(inp::InputState, out::OutputState{T}, count::Int32, bits::UInt8, mask::V=MASKN(bits)) where {T <: Integer, V <: Integer}
    # the mask is of the smallest bounding type for bits
    # T is one of the types that map on to the appropriate Julia type in Parquet (which may be larger than the mask type)
    bitbuff = zero(V)
    nbitsbuff = 0

    data = inp.data
    arr = out.data
    offset = out.offset
    arridx = Int32(1)
    dataidx = Int32(1 + inp.offset)
    while arridx <= count
        diffnbits = bits - nbitsbuff
        while diffnbits > 8
            # shift 8 bits and read directly into bitbuff
            bitbuff <<= 8
            bitbuff |= data[dataidx]
            dataidx += Int32(1)
            nbitsbuff += 8
            diffnbits -= 8
        end

        if diffnbits > 0
            # read next byte from input
            nxtdata = data[dataidx]
            dataidx += Int32(1)
            # shift bitbuff by diffnbits, add diffnbits and set result
            nbitsbuff = 8 - diffnbits
            arr[arridx+offset] = convert(T, ((bitbuff << diffnbits) | (nxtdata >> nbitsbuff)) & mask)
            arridx += Int32(1)
            # keep remaining bits in bitbuff
            bitbuff <<= 8
            bitbuff |= nxtdata
        else
            # set result
            arr[arridx+offset] = convert(T, (bitbuff >> abs(diffnbits)) & mask)
            arridx += Int32(1)
            nbitsbuff -= bits
        end
    end
    inp.offset = dataidx - 1
    out.offset += count
    nothing
end

function logical_timestamp(barr; offset::Dates.Period=Dates.Second(0))
    nanos = read(IOBuffer(barr[1:8]), Int64)
    julian_days = read(IOBuffer(barr[9:12]), Int32)
    Dates.julian2datetime(julian_days) + Dates.Nanosecond(nanos) + offset
end

function logical_timestamp(i128::Int128; offset::Dates.Period=Dates.Second(0))
    iob = IOBuffer()
    write(iob, i128)
    logical_timestamp(take!(iob); offset=offset)
end

logical_string(bytes::Vector{UInt8}) = String(bytes)

function logical_decimal(bytes::Union{Int64,Int32,Vector{UInt8}}, precision::Integer, scale::Integer; use_float::Bool=false)
    T = logical_decimal_unscaled_type(Int32(precision))
    if scale == 0
        logical_decimal_integer(bytes, T)
    elseif use_float
        logical_decimal_float64(bytes, T, Int32(scale))
    else
        logical_decimal_scaled(bytes, T, Int32(scale))
    end
end

function logical_decimal_integer(intval::Union{Int64,Int32}, ::Type{T}) where T <: Union{UInt16,UInt32,UInt64,UInt128}
    signed(T)(intval)
end

function logical_decimal_integer(bytes::Vector{UInt8}, ::Type{T}) where T <: Union{UInt16,UInt32,UInt64,UInt128}
    N = length(bytes)
    uintval = T(0)
    for idx in 1:N
        uintval |= (T(bytes[idx]) << *(8,N-idx))
    end
    reinterpret(signed(T), uintval)
end

function logical_decimal_float64(bytes::Union{Int64,Int32,Vector{UInt8}}, ::Type{T}, scale::Int32) where T <: Union{UInt16,UInt32,UInt64,UInt128}
    if scale > 0
        logical_decimal_integer(bytes, T) / 10^scale
    else
        logical_decimal_integer(bytes, T) * 10^abs(scale)
    end
end

function logical_decimal_scaled(bytes::Union{Int64,Int32,Vector{UInt8}}, ::Type{T}, scale::Int32) where T <: Union{UInt16,UInt32,UInt64,UInt128}
    intval = logical_decimal_integer(bytes, T)
    sign = (intval < 0) ? 1 : 0
    intval = abs(intval)
    Decimal(sign, intval, -scale)
end


text:
isinplace(W::AbstractNoiseProcess{T,N,S,inplace}) where {T,N,S,inplace} = inplace

mutable struct NoiseProcess{T,N,Tt,T2,T3,ZType,F,F2,inplace,S1,S2,RSWM,C,RNGType} <: AbstractNoiseProcess{T,N,Vector{T2},inplace}
  dist::F
  bridge::F2
  t::Vector{Tt}
  u::Vector{T2} # Aliased pointer to W for the AbstractVectorOfArray interface
  W::Vector{T2}
  Z::ZType
  curt::Tt
  curW::T2
  curZ::T3
  dt::Tt
  dW::T2
  dZ::T3
  dWtilde::T2
  dZtilde::T3
  dWtmp::T2
  dZtmp::T3
  S₁::S1
  S₂::S2
  rswm::RSWM
  maxstacksize::Int
  maxstacksize2::Int
  save_everystep::Bool
  iter::Int
  rng::RNGType
  reset::Bool
  reseed::Bool
  continuous::Bool
  cache::C

  function NoiseProcess{iip}(t0,W0,Z0,dist,bridge;
                         rswm = RSWM(),save_everystep=true,
                         rng = Xorshifts.Xoroshiro128Plus(rand(UInt64)),
                         reset = true, reseed = true, continuous = true,
                         cache = nothing) where iip
    S₁ = ResettableStacks.ResettableStack{iip}(
                          Tuple{typeof(t0),typeof(W0),typeof(Z0)})
    S₂ = ResettableStacks.ResettableStack{iip}(
                          Tuple{typeof(t0),typeof(W0),typeof(Z0)})
    if Z0==nothing
      Z=nothing
      curZ = nothing
      dZ = nothing
      dZtilde= nothing
      dZtmp = nothing
    else
      Z=[copy(Z0)]
      curZ = copy(Z0)
      dZ = copy(Z0)
      dZtilde= copy(Z0)
      dZtmp = copy(Z0)
    end
    W = [copy(W0)]
    N = length((size(W0)..., length(W)))
    new{eltype(eltype(W0)),N,typeof(t0),typeof(W0),typeof(dZ),typeof(Z),
                  typeof(dist),typeof(bridge),
                  iip,typeof(S₁),typeof(S₂),typeof(rswm),typeof(cache),typeof(rng)}(
                  dist,bridge,[t0],W,W,Z,t0,
                  copy(W0),curZ,t0,copy(W0),dZ,copy(W0),dZtilde,copy(W0),dZtmp,
                  S₁,S₂,rswm,0,0,save_everystep,0,rng,reset,reseed,continuous,cache)
  end

end
(W::NoiseProcess)(t) = interpolate!(W,nothing,nothing,t)
(W::NoiseProcess)(u,p,t) = interpolate!(W,u,p,t)
(W::NoiseProcess)(out1,out2,u,p,t) = interpolate!(out1,out2,W,u,p,t)
adaptive_alg(W::NoiseProcess) = adaptive_alg(W.rswm)

function NoiseProcess(t0,W0,Z0,dist,bridge;kwargs...)
  iip=DiffEqBase.isinplace(dist,7)
  NoiseProcess{iip}(t0,W0,Z0,dist,bridge;kwargs...)
end

mutable struct SimpleNoiseProcess{T,N,Tt,T2,T3,ZType,F,F2,inplace,RNGType} <: AbstractNoiseProcess{T,N,Vector{T2},inplace}
  dist::F
  bridge::F2
  t::Vector{Tt}
  u::Vector{T2} # Aliased pointer to W for the AbstractVectorOfArray interface
  W::Vector{T2}
  Z::ZType
  curt::Tt
  curW::T2
  curZ::T3
  dt::Tt
  dW::T2
  dZ::T3
  dWtilde::T2
  dZtilde::T3
  dWtmp::T2
  dZtmp::T3
  save_everystep::Bool
  iter::Int
  rng::RNGType
  reset::Bool
  reseed::Bool

  function SimpleNoiseProcess{iip}(t0,W0,Z0,dist,bridge;
                         save_everystep=true,
                         rng = Xorshifts.Xoroshiro128Plus(rand(UInt64)),
                         reset = true, reseed = true) where iip
    if Z0==nothing
      Z=nothing
      curZ = nothing
      dZ = nothing
      dZtilde= nothing
      dZtmp = nothing
    else
      Z=[copy(Z0)]
      curZ = copy(Z0)
      dZ = copy(Z0)
      dZtilde= copy(Z0)
      dZtmp = copy(Z0)
    end
    W = [copy(W0)]
    N = length((size(W0)..., length(W)))
    new{eltype(eltype(W0)),N,typeof(t0),typeof(W0),typeof(dZ),typeof(Z),
                  typeof(dist),typeof(bridge),iip,typeof(rng)}(
                  dist,bridge,[t0],W,W,Z,t0,
                  copy(W0),curZ,t0,copy(W0),dZ,copy(W0),dZtilde,copy(W0),dZtmp,
                  save_everystep,0,rng,reset,reseed)
  end

end
(W::SimpleNoiseProcess)(t) = interpolate!(W,nothing,nothing,t)
(W::SimpleNoiseProcess)(u,p,t) = interpolate!(W,u,p,t)
(W::SimpleNoiseProcess)(out1,out2,u,p,t) = interpolate!(out1,out2,W,u,p,t)

function SimpleNoiseProcess(t0,W0,Z0,dist,bridge;kwargs...)
  iip=DiffEqBase.isinplace(dist,7)
  SimpleNoiseProcess{iip}(t0,W0,Z0,dist,bridge;kwargs...)
end

mutable struct NoiseWrapper{T,N,Tt,T2,T3,T4,ZType,inplace} <: AbstractNoiseProcess{T,N,Vector{T2},inplace}
  t::Vector{Tt}
  u::Vector{T2}
  W::Vector{T2}
  Z::ZType
  curt::Tt
  curW::T2
  curZ::T3
  dt::Tt
  dW::T2
  dZ::T3
  source::T4
  reset::Bool
  reverse::Bool
end

function NoiseWrapper(source::AbstractNoiseProcess{T,N,Vector{T2},inplace};
                      reset=true,reverse=false,indx=nothing) where {T,N,T2,inplace}

  if indx === nothing
    if reverse
      indx = length(source.t)
    else
      indx = 1
    end
  end

  if source.Z==nothing
    Z=nothing
    curZ = nothing
    dZ = nothing
  else
    Z=[copy(source.Z[indx])]
    curZ = copy(source.Z[indx])
    dZ = copy(source.Z[indx])
  end
  W = [copy(source.W[indx])]

  NoiseWrapper{T,N,typeof(source.t[1]),typeof(source.W[1]),typeof(dZ),typeof(source),typeof(Z),inplace}(
                [source.t[indx]],W,W,Z,source.t[indx],copy(source.W[indx]),curZ,source.t[indx],copy(source.W[indx]),dZ,source,reset,reverse)
end

(W::NoiseWrapper)(t) = interpolate!(W,nothing,nothing,t)
(W::NoiseWrapper)(u,p,t) = interpolate!(W,u,p,t)
(W::NoiseWrapper)(out1,out2,u,p,t) = interpolate!(out1,out2,W,u,p,t,reverse=W.reverse)
adaptive_alg(W::NoiseWrapper) = adaptive_alg(W.source)

mutable struct NoiseFunction{T,N,wType,zType,Tt,T2,T3,inplace} <: AbstractNoiseProcess{T,N,nothing,inplace}
  W::wType
  Z::zType
  curt::Tt
  curW::T2
  curZ::T3
  dt::Tt
  dW::T2
  dZ::T3
  reset::Bool

  function NoiseFunction{iip}(t0,W,Z=nothing;
                         noise_prototype=W(nothing,nothing,t0),reset=true) where iip
    curt = t0
    dt = t0
    curW = copy(noise_prototype)
    dW = copy(noise_prototype)
    if Z==nothing
      curZ = nothing
      dZ = nothing
    else
      curZ = copy(noise_prototype)
      dZ = copy(noise_prototype)
    end
    new{typeof(noise_prototype),ndims(noise_prototype),typeof(W),typeof(Z),
                  typeof(curt),typeof(curW),typeof(curZ),iip}(W,Z,curt,curW,curZ,
                  dt,dW,dZ,reset)
  end

end

(W::NoiseFunction)(t) = W(nothing,nothing,t)
function (W::NoiseFunction)(u,p,t)
  if W.Z != nothing
    if isinplace(W)
      out2 = similar(W.dZ)
      W.Z(out2,u,p,t)
    else
      out2 = W.Z(u,p,t)
    end
  else
    out2 = nothing
  end
  if isinplace(W)
    out1 = similar(W.dW)
    W.W(out1,u,p,t)
  else
    out1 = W.W(u,p,t)
  end
  out1,out2
end
function (W::NoiseFunction)(out1,out2,u,p,t)
  W.W(out1,u,p,t)
  W.Z != nothing && W.Z(out2,u,p,t)
end

function NoiseFunction(t0,W,Z=nothing;kwargs...)
  iip=DiffEqBase.isinplace(W,4)
  NoiseFunction{iip}(t0,W,Z;kwargs...)
end

mutable struct NoiseGrid{T,N,Tt,T2,T3,ZType,inplace} <: AbstractNoiseProcess{T,N,Vector{T2},inplace}
  t::Vector{Tt}
  u::Vector{T2}
  W::Vector{T2}
  Z::ZType
  curt::Tt
  curW::T2
  curZ::T3
  dt::Tt
  dW::T2
  dZ::T3
  step_setup::Bool
  reset::Bool
end

function NoiseGrid(t,W,Z=nothing;reset=true)
  val = W[1]
  curt = t[1]
  dt = t[1]
  curW = copy(val)
  dW = copy(val)
  if Z==nothing
    curZ = nothing
    dZ = nothing
  else
    curZ = copy(Z[1])
    dZ = copy(Z[1])
  end
  typeof(val) <: AbstractArray ? iip = true : iip = false
  NoiseGrid{typeof(val),ndims(val),typeof(dt),typeof(dW),typeof(dZ),typeof(Z),iip}(
            t,W,W,Z,curt,curW,curZ,dt,dW,dZ,true,reset)
end

(W::NoiseGrid)(t) = interpolate!(W,t)
(W::NoiseGrid)(u,p,t) = interpolate!(W,t)
(W::NoiseGrid)(out1,out2,u,p,t) = interpolate!(out1,out2,W,t)

mutable struct NoiseApproximation{T,N,Tt,T2,T3,S1,S2,ZType,inplace} <: AbstractNoiseProcess{T,N,Vector{T2},inplace}
  t::Vector{Tt}
  u::Vector{T2}
  W::Vector{T2}
  Z::ZType
  curt::Tt
  curW::T2
  curZ::T3
  dt::Tt
  dW::T2
  dZ::T3
  source1::S1
  source2::S2
  reset::Bool
end

function NoiseApproximation(source1::DEIntegrator,source2::Union{DEIntegrator,Nothing}=nothing;
                   reset=true)
  _source1 = deepcopy(source1)
  _source2 = deepcopy(source2)
  if _source2==nothing
    Z=nothing
    curZ = nothing
    dZ = nothing
  else
    Z=_source2.sol.u
    curZ = copy(_source2.u)
    dZ = copy(_source2.u)
    _source2.opts.advance_to_tstop = true
  end
  val = copy(_source1.u)
  t = _source1.sol.t
  W = _source1.sol.u
  curW = copy(_source1.u)
  dW = copy(_source1.u)
  dt = _source1.dt
  curt = _source1.t
  _source1.opts.advance_to_tstop = true
  NoiseApproximation{typeof(val),ndims(val),typeof(curt),typeof(curW),typeof(curZ),
                     typeof(_source1),typeof(_source2),typeof(Z),
                     isinplace(_source1.sol.prob)}(
                     t,W,W,Z,curt,curW,curZ,dt,dW,dZ,_source1,_source2,reset)
end

(W::NoiseApproximation)(t) = interpolate!(W,t)
(W::NoiseApproximation)(u,p,t) = interpolate!(W,t)
(W::NoiseApproximation)(out1,out2,u,p,t) = interpolate!(out1,out2,W,t)

mutable struct VirtualBrownianTree{T,N,F,F2,Tt,T2,T3,T2tmp,T3tmp,seedType,tolType,RNGType,inplace} <: AbstractNoiseProcess{T,N,Vector{T2},inplace}
  dist::F
  bridge::F2
  t::Vector{Tt} # tstart::Tt to tend::Tt
  u::Vector{T2}
  W::Vector{T2}
  Z::Vector{T3}
  curt::Tt
  curW::T2
  curZ::T3
  dt::Tt
  dW::T2
  dZ::T3
  W0tmp::T2tmp
  W1tmp::T2tmp
  Z0tmp::T3tmp
  Z1tmp::T3tmp
  seeds::Vector{seedType}
  atol::tolType
  rng::RNGType
  tree_depth::Int
  search_depth::Int
  step_setup::Bool
  reset::Bool
end

function VirtualBrownianTree{iip}(t0,W0,Z0,dist,bridge;
                      tend=nothing,Wend=nothing,Zend=nothing,
                      atol=1e-10,tree_depth::Int=4,
                      search_depth=nothing,
                      rng = RandomNumbers.Random123.Threefry4x(),
                      reset=true) where iip

  if search_depth==nothing
    if atol<1e-10
      search_depth = 50 # maximum search depth
    else
      search_depth = floor(Int,(-log2(atol)+2))
    end
  end

  (tree_depth >= search_depth) && error("search depth of the VBT must be greater than the depth of the cached tree")

  curt = t0
  dt = t0
  curW = copy(W0)
  dW = copy(W0)
  if Z0==nothing
    curZ = nothing
    dZ = nothing
  else
    curZ = copy(Z0)
    dZ = copy(Z0)
  end

  # assign final time and state
  if tend==nothing
    tend = t0 + one(t0)
  end

  if Wend==nothing
     Wend = curW + dist(dW,nothing,tend-t0,nothing,nothing,nothing,rng)
  end

  if Zend==nothing && Z0!=nothing
     Zend = curW + bridge(dZ,nothing,tend-t0,nothing,nothing,nothing,rng)
  end

  t, W, Z, seeds = create_VBT_cache(bridge,t0,W0,Z0,tend,Wend,Zend,rng,tree_depth,search_depth)

  if iip
    W0tmp, W1tmp = copy(W0), copy(Wend)
    if Z0!=nothing
      Z0tmp,Z1tmp = copy(Z0), copy(Zend)
    else
      Z0tmp,Z1tmp = nothing, nothing
    end
  else
    W0tmp,W1tmp,Z0tmp,Z1tmp = nothing,nothing,nothing,nothing
  end

  VirtualBrownianTree{
   typeof(W0),ndims(W0),
   typeof(dist),typeof(bridge),typeof(curt),typeof(curW),typeof(curZ),
   typeof(W0tmp),typeof(Z0tmp),typeof(seeds[1]),typeof(atol),typeof(rng),
   iip}(dist,bridge,t,W,W,Z,curt,curW,curZ,dt,dW,dZ,W0tmp,W1tmp,Z0tmp,Z1tmp,seeds,atol,rng,tree_depth,
   search_depth,true,reset)
end

(W::VirtualBrownianTree)(t) = interpolate!(W,nothing, nothing, t)
(W::VirtualBrownianTree)(u,p,t) = interpolate!(W,u,p,t)
(W::VirtualBrownianTree)(out1,out2,u,p,t) = interpolate!(out1,out2,W,u,p,t)

function VirtualBrownianTree(t0,W0,Z0=nothing,dist=WHITE_NOISE_DIST,bridge=VBT_BRIDGE;kwargs...)
  VirtualBrownianTree{false}(t0,W0,Z0,dist,bridge;kwargs...)
end

function VirtualBrownianTree!(t0,W0,Z0=nothing,dist=INPLACE_WHITE_NOISE_DIST,bridge=INPLACE_VBT_BRIDGE;kwargs...)
  VirtualBrownianTree{true}(t0,W0,Z0,dist,bridge;kwargs...)
end

ext: jl


text:
# Point conversion routines
using Test
import CryptoGroups: bits2octet, octet2bits, spec, <|, AffinePoint, specialize, generator

let
    x = UInt8[0, 129]
    @test bits2octet(octet2bits(x)) == x

    a = octet2bits(x)
    # First leftmostbits are padded to zero to make a full octet
    @test bits2octet(a[4:end]) == bits2octet(a)

    @test octet2bits(x, 13) == octet2bits(x)[4:end]
end

curve_spec = spec(:P_256)

P = specialize(AffinePoint, curve_spec)

p = P <| generator(curve_spec)

@test P <| Vector{UInt8} <| p == p

@test P <| (Vector{UInt8}, :uncompressed) <| p == p


text:
"""
A local fixture that creates a temporary directory and returns its name;
the directory and all its contents is removed during the teardown.
"""
const temporary_dir = @local_fixture begin
    dir = mktempdir()
    @produce dir
    rm(dir, recursive=true)
end

"""
A global fixture that returns the dictionary with the current run options
(see [`Run options`](@ref run_options_manual) in the manual for the full list.
"""
const run_options = run_options_fixture()


text:
#---------------------------------------------------------
# # [NUFFT Overview](@id 1-nufft)
#---------------------------------------------------------

#=
This example illustrates how to use Nonuniform FFT (NUFFT)
for image reconstruction in MRI
using the Julia language.

This entire page was generated using a single Julia file:
[1-nufft.jl](@__REPO_ROOT_URL__/mri/1-nufft.jl).
=#
#md # In any such Julia documentation,
#md # you can access the source code
#md # using the "Edit on GitHub" link in the top right.

#md # The corresponding notebook can be viewed in
#md # [nbviewer](http://nbviewer.jupyter.org/) here:
#md # [`1-nufft.ipynb`](@__NBVIEWER_ROOT_URL__/mri/1-nufft.ipynb),
#md # and opened in [binder](https://mybinder.org/) here:
#md # [`1-nufft.ipynb`](@__BINDER_ROOT_URL__/mri/1-nufft.ipynb).

#=
Some MRI scans use non-Cartesian sampling patterns
(like radial and spiral k-space trajectories, among others),
often in the interest of acquisition speed.

Image reconstruction from fully sampled Cartesian k-space data
typically uses simple inverse FFT operations,
whereas non-Cartesian sampling requires more complicated methods.

The examples here illustrate both non-iterative (aka "gridding")
and iterative methods for non-Cartesian MRI reconstruction.
For simplicity the examples consider the case of single-coil data
and ignore the effects of B0 field inhomogeneity.

First we add the Julia packages that are need for these examples.
Change `false` to `true` in the following code block
if you are using any of the following packages for the first time.
=#

if false
    import Pkg
    Pkg.add([
        "ImagePhantoms"
        "Plots"
        "Unitful"
        "UnitfulRecipes"
        "LaTeXStrings"
        "MIRTjim"
        "MIRT"
        "InteractiveUtils"
    ])
end

# Now tell this Julia session to use the following packages for this example.
# Run `Pkg.add()` in the preceding code block first, if needed.

using ImagePhantoms: shepp_logan, SheppLoganEmis, spectrum, phantom #, Gauss2
using Plots; default(label="", markerstrokecolor=:auto)
using Unitful: mm # Allows use of physical units (mm here)
using UnitfulRecipes # Ensures that plots show appropriate units when applicable
using LaTeXStrings # for LaTeX in plot labels, e.g., L"\alpha_n"
using MIRTjim: jim, prompt # jiffy image display
using MIRT: Anufft, diffl_map, ncg
using InteractiveUtils: versioninfo

# The following line is helpful when running this jl-file as a script;
# this way it will prompt user to hit a key after each image is displayed.

isinteractive() && jim(:prompt, true);

# ## Radial k-space sampling

#=
We focus on radial sampling as a simple representative non-Cartesian case.
Consider imaging a 256mm × 256mm field of FOV
with the goal of reconstructing a 128 × 128 pixel image.
The following radial and angular k-space sampling is reasonable.
=#

N = 128
FOV = 256mm # physical units!
Δx = FOV / N # pixel size
kmax = 1 / 2Δx
kr = ((-N÷2):(N÷2)) / (N÷2) * kmax # radial sampling in k-space
Nϕ = 3N÷2 # theoretically should be about π/2 * N
kϕ = (0:Nϕ-1)/Nϕ * π # angular samples
νx = kr * cos.(kϕ)' # N × Nϕ k-space sampling in cycles/mm
νy = kr * sin.(kϕ)'
plot(νx, νy,
    xlabel=L"\nu_x", ylabel=L"\nu_y",
    aspect_ratio = 1,
    title = "Radial k-space sampling",
)

#
isinteractive() && prompt();

#=
For the NUFFT routines considered here,
the sampling arguments must be "Ω" values:
digital frequencies have pseudo-units of radians / pixel.
=#

Ωx = (2π * Δx) * νx # N × Nϕ grid of k-space sample locations
Ωy = (2π * Δx) * νy # in pseudo-units of radians / sample

scatter(Ωx, Ωy,
    xlabel=L"\Omega_x", ylabel=L"\Omega_y",
    xticks=((-1:1)*π, ["-π", "0", "π"]),
    yticks=((-1:1)*π, ["-π", "0", "π"]),
    xlims=(-π,π) .* 1.1,
    ylims=(-π,π) .* 1.1,
    aspect_ratio = 1, markersize = 0.5,
    title = "Radial k-space sampling",
)

#
isinteractive() && prompt();

# ## Radial k-space data for Shepp-Logan phantom

#=
Get the ellipse parameters for a MRI-suitable version of the Shepp-Logan phantom
and calculate (analytically) the radial k-space data.
Then display in polar coordinates.
=#

object = shepp_logan(SheppLoganEmis(); fovs=(FOV,FOV))
#object = [Gauss2(18mm, 0mm, 100mm, 70mm, 0, 1)] # useful for validating DCF
data = spectrum(object).(νx,νy)
data = data / oneunit(eltype(data)) # abandon units at this point
jim(kr, kϕ, abs.(data), title="k-space data magnitude",
    xlabel=L"k_r",
    ylabel=L"k_{\phi}",
    xticks = (-1:1) .* maximum(abs, kr),
    yticks = (0,π),
    ylims = (0,π),
    aspect_ratio = :none,
)

# ## Non-iterative gridding image reconstruction

#=
It would be impossible for a radiologist
to diagnose a patient
from the k-space data in polar coordinates,
so image reconstruction is needed.

The simplest approach
is to (nearest-neighbor) interpolate the k-space data
onto a Cartesian grid,
and then reconstruction with an inverse FFT.

One way to do that interpolation
is to use a `Histogram` method
in Julia's statistics package.
=#

using StatsBase: fit, Histogram, weights

#=
The following function is a work-around because `weights` in StatsBase is
[limited to Real data](https://github.com/JuliaStats/StatsBase.jl/issues/745),
so here we bin the real and imaginary k-space data separately,
and handle the units.
=#

function histogram(coord, vals::AbstractArray{<:Number}, edges)
    u = oneunit(eltype(vals))
    wr = weights(real(vec(vals / u)))
    wi = weights(imag(vec(vals / u)))
    tmp1 = fit(Histogram, coord, wr, edges)
    tmp2 = fit(Histogram, coord, wi, edges)
    return u * complex.(tmp1.weights, tmp2.weights)
end

kx = N * Δx * νx # N × Nϕ k-space sampling in cycles/mm
ky = N * Δx * νy
bin = (-(N÷2):(N÷2)) .- 0.5 # (N+1,) histogram bin edges
gridded1 = histogram((vec(kx), vec(ky)), data, (bin,bin))

using FFTW: ifft, fftshift
tmp = fftshift(ifft(fftshift(gridded1)))
x = (-(N÷2):(N÷2-1)) * Δx
y = (-(N÷2):(N÷2-1)) * Δx
jim(x, y, tmp, title="Elementary gridding reconstruction")

# That crummy gridding method does not work well.
# Here's what the phantom should look like:

ideal = phantom(x, y, object, 2)
clim = (0, 9)
p0 = jim(x, y, ideal, title="True Shepp-Logan phantom"; clim)

# ## NUFFT approach to gridding

#=
Basic nearest-neighbor gridding
does not provide acceptable image quality in MRI,
so now we turn to using the NUFFT.
For MRI purposes,
the NUFFT is a function that maps Cartesian spaced image data
into non-Cartesian k-space data.
The
[NFFT.jl](https://github.com/tknopp/NFFT.jl)
package has functions
for computing the NUFFT and its adjoint.
These are linear mappings (generalizations of matrices),
so instead of calling those functions directly,
here we use the NUFFT linear map object
defined in MIRT
that provides a non-Cartesian Fourier encoding "matrix".
=#

A = Anufft([vec(Ωx) vec(Ωy)], (N,N); n_shift = [N/2,N/2])

#=
This linear map is constructed to map a N × N image into `length(Ωx)` k-space samples.
So its
[adjoint](https://en.wikipedia.org/wiki/Adjoint)
goes the other direction.
However, an adjoint is *not* an inverse!
=#

gridded2 = A' * vec(data)
jim(x, y, gridded2, title="NUFFT gridding without DCF")

# ## Density compensation

#=
To get a decent image with NUFFT-based gridding of non-Cartesian data,
one must compensate for the k-space sampling density.
See
[this book chapter](https://web.eecs.umich.edu/~fessler/book/c-four.pdf)
for details.

Because this example uses radial sampling,
we can borrow ideas from tomography,
especially the
[ramp filter](https://en.wikipedia.org/wiki/Radon_transform#Radon_inversion_formula),
to define a reasonable density compensation function (DCF).

Here is a basic DCF version that uses the ramp filter in a simple way,
corresponding to the areas of annular segments
(Voronoi cells in polar coordinates).
The `dcf .* data` line uses Julia's
[broadcast](https://docs.julialang.org/en/v1/manual/arrays/#Broadcasting)
feature
to apply the 1D DCF to each radial spoke.
=#

dν = 1/FOV # k-space radial sample spacing
dcf = pi / Nϕ * dν * abs.(kr) # (N+1) weights along k-space polar coordinate
dcf = dcf / oneunit(eltype(dcf)) # kludge: units not working with LinearMap now
gridded3 = A' * vec(dcf .* data)
p3 = jim(x, y, gridded3, title="NUFFT gridding with simple ramp-filter DCF"; clim)

#=
The image is more reasonable than without any DCF,
but we can do better (quantitatively) using the correction of
[Lauzon&Rutt, 1996](http://doi.org/10.1002/mrm.1910360617)
and
[Joseph, 1998](http://doi.org/10.1002/mrm.1910400317).
=#

dcf = pi / Nϕ * dν * abs.(kr) # see lauzon:96:eop, joseph:98:sei
dcf[kr .== 0/mm] .= pi * (dν/2)^2 / Nϕ # area of center disk
dcf = dcf / oneunit(eltype(dcf)) # kludge: units not working with LinearMap now
gridded4 = A' * vec(dcf .* data)
p4 = jim(x, y, gridded4, title="NUFFT gridding with better ramp-filter DCF"; clim)

# A profile helps illustrate the improvement.

pp = plot(x, real(gridded4[:,N÷2]), label="Modified ramp DCF")
plot!(x, real(gridded3[:,N÷2]), label="Basic ramp DCF")
plot!(x, real(ideal[:,N÷2]), label="Ideal", xlabel=L"x", ylabel="middle profile")

#
isinteractive() && prompt();

#=
Finally we have made a NUFFT gridded image with DCF
that has the appropriate range of values,
but it still looks less than ideal.
So next we try an iterative approach.

## Iterative MR image reconstruction using NUFFT

As an initial iterative approach,
let's apply a few iterations of conjugate gradient (CG)
to seek the minimizer of the least-squares cost function:
```math
\arg \min_{x} \frac{1}{2} \| A x - y \|_2^2.
```
CG is well-suited to minimizing quadratic cost functions,
but we do not expect the image to be great quality
because radial sampling omits the "corners" of k-space
so the NUFFT operator ``A`` is badly conditioned.
=#

dx = FOV / N # pixel size
dx = dx / oneunit(dx) # abandon units for now
gradf = u -> u - vec(data/dx^2) # gradient of f(u) = 1/2 \| u - y \|^2
curvf = u -> 1 # curvature of f(u)
x0 = gridded4 # initial guess: best gridding reconstruction
xls, _ = ncg([A], [gradf], [curvf], x0; niter = 20)
p5 = jim(x, y, xls, "LS-CG reconstruction"; clim)

#=
To improve the results, we include regularization.
Here we would like to reconstruct an image
by finding the minimizer of a regularized LS cost function
such as the following:
```math
\arg \min_{x} \frac{1}{2} \| A x - y \|_2^2 + \beta R(x)
, \qquad
R(x) = 1' \psi.(T x).
```
We focus here on the case of edge-preserving regularization
where ``T`` is a 2D finite-differencing operator
and ``\psi`` is a potential function.
This operator maps a N×N image into a N×N×2 array
with the horizontal and vertical finite differences.
=#

T = diffl_map((N,N), [1,2] ; T = ComplexF32)

# Applying this operator to the ideal image illustrated its action:
jim(x, y, T * ideal; ncol=1, title="Horizontal and vertical finite differences")

#=
## Edge-preserving regularization

We use the Fair potential function:
a rounded corner version of absolute value,
an approximation to anisotropic total variation (TV).
=#

β = 2^13 # regularization parameter
δ = 0.05 # edge-preserving parameter
wpot = z -> 1 / (1 + abs(z)/δ) # weighting function

#=
## Nonlinear CG algorithm

We apply a (nonlinear) CG algorithm
to seek the minimizer of the cost function.
Nonlinear CG is well suited to convex problems
that are locally quadratic
like the regularized cost function considered here.
See
[this survey paper](http://doi.org/10.1109/MSP.2019.2943645)
for an overview of optimization methods for MRI.
=#

B = [A, T] # see MIRT.ncg
gradf = [u -> u - vec(data/dx^2), # data-term gradient, correct for pixel area
         u -> β * (u .* wpot.(u))] # regularizer gradient
curvf = [u -> 1, u -> β] # curvature of quadratic majorizers
x0 = gridded4 # initial guess is best gridding reconstruction
xhat, _ = ncg(B, gradf, curvf, x0; niter = 90)
p6 = jim(x, y, xhat, "Iterative reconstruction"; clim)

# Here is a comparison of the profiles.

plot!(pp, x, real(xls[:,N÷2]), label="LS-CG")
plot!(pp, x, real(xhat[:,N÷2]), label="Iterative edge-preserving", color=:black)

#
isinteractive() && prompt();

#=
In this case, iterative image reconstruction provides the best looking image.
One might argue this simulation was doomed to succeed,
because the phantom is piece-wise constant,
which is the best case for edge-preserving regularization.
On the other hand,
this was not an
[inverse crime](http://doi.org/10.1016/j.cam.2005.09.027)
([see also here](http://arxiv.org/abs/2109.08237))
because the k-space data came from the analytical spectrum of ellipses,
rather than from a discrete image.

### Caveats

* The phantom used here was real-valued, which is unrealistic
  (although the reconstruction methods did not "know" it was real).
* This simulation is for a single-coil scan
  whereas modern MRI scanners generally use multiple receive coils.
* There was no statistical noise in this simulation.
=#

# ## Reproducibility

# This page was generated with the following version of Julia:

io = IOBuffer(); versioninfo(io); split(String(take!(io)), '\n')

# And with the following package versions

import Pkg; Pkg.status()

max_stars_repo_name: JuliaImageRecon/Examples

max_issues_repo_name: JuliaImageRecon/Examples

max_forks_repo_name: JuliaImageRecon/Examples


text:
const PENDING = 0
const APPROVED = 1
const REJECTED = 2

const IN = 0
const OUT = 1
const ERR = 2
const ENV = 4

const WAITING = 0
const RUNNING = 1
const SUCCESS = 2
const FAILED = 3

struct Colony
  colonyid::String
  name::String
end

struct Runtime
  runtimeid::String
  runtimetype::String
  name::String
  colonyid::String
  cpu::String
  cores::Int64
  mem::Int64
  gpu::String
  gpus::Int64
  state::Int64

  function Runtime(runtimeid::String, runtimetype::String, name::String, colonyid::String, cpu::String, cores::Int64, mem::Int64, gpu::String,  gpus::Int64, state::Int64)
    new(runtimeid, runtimetype, name, colonyid, cpu, cores, mem, gpu, gpus, state)
  end

  function Runtime(runtimeid::String, runtimetype::String, name::String, colonyid::String, cpu::String, cores::Int64, mem::Int64, gpu::String,  gpus::Int64)
    new(runtimeid, runtimetype, name, colonyid, cpu, cores, mem, gpu, gpus, PENDING)
  end

  function Runtime(runtimeid::String, runtimetype::String, name::String, colonyid::String)
    new(runtimeid, runtimetype, name, colonyid, "", 1, 0, "", 0, PENDING)
  end

end

struct Conditions
  colonyid::String
  runtimeids::Array{String,1}
  runtimetype::String
  mem::UInt16
  cores::UInt16
  gpus::UInt16
end

struct ProcessSpec
  timeout::Int16
  maxretries::Int16
  conditions::Conditions
  env::Dict{String,String}
end

struct Attribute
  attributeid::String
  targetid::String
  attributetype::UInt16
  key::String
  value::String

  function Attribute(targetid::String, key::String, value::String)
    new("", targetid, OUT, key, value)
  end

  function Attribute(attributeid::String, targetid::String,attributetype::UInt16, key::String, value::String)
    new(attributeid, targetid, attributetype, key, value)
  end
end

struct Process
  processid::String
  assignedruntimeid::String
  isassigned::Bool
  state::UInt16
  submissiontime::String
  starttime::String
  endtime::String
  deadline::String
  retries::UInt16
  attributes::Array{Attribute,1}
  spec::ProcessSpec
end


text:
type BasicContUnvParameterNState{N<:Real} <: ParameterNState{Continuous, Univariate}
  value::Vector{N}
  loglikelihood::Vector{N}
  logprior::Vector{N}
  logtarget::Vector{N}
  gradloglikelihood::Vector{N}
  gradlogprior::Vector{N}
  gradlogtarget::Vector{N}
  tensorloglikelihood::Vector{N}
  tensorlogprior::Vector{N}
  tensorlogtarget::Vector{N}
  dtensorloglikelihood::Vector{N}
  dtensorlogprior::Vector{N}
  dtensorlogtarget::Vector{N}
  diagnosticvalues::Matrix
  n::Integer
  diagnostickeys::Vector{Symbol}
  copy::Function

  function BasicContUnvParameterNState(
    n::Integer,
    monitor::Vector{Bool}=[true; fill(false, 12)],
    diagnostickeys::Vector{Symbol}=Symbol[],
    ::Type{N}=Float64,
    diagnosticvalues::Matrix=Array(Any, length(diagnostickeys), isempty(diagnostickeys) ? 0 : n)
  )
    instance = new()

    l = Array(Integer, 13)
    for i in 1:13
      l[i] = (monitor[i] == false ? zero(Integer) : n)
    end

    fnames = fieldnames(BasicContUnvParameterNState)
    for i in 1:13
      setfield!(instance, fnames[i], Array(N, l[i]))
    end

    instance.diagnosticvalues = diagnosticvalues

    instance.n = n
    instance.diagnostickeys = diagnostickeys

    instance.copy = eval(codegen(:copy, instance, monitor))

    instance
  end
end

BasicContUnvParameterNState{N<:Real}(
  n::Integer,
  monitor::Vector{Bool}=[true; fill(false, 12)],
  diagnostickeys::Vector{Symbol}=Symbol[],
  ::Type{N}=Float64,
  diagnosticvalues::Matrix=Array(Any, length(diagnostickeys), isempty(diagnostickeys) ? 0 : n)
) =
  BasicContUnvParameterNState{N}(n, monitor, diagnostickeys, N, diagnosticvalues)

function BasicContUnvParameterNState{N<:Real}(
  n::Integer,
  monitor::Vector{Symbol},
  diagnostickeys::Vector{Symbol}=Symbol[],
  ::Type{N}=Float64,
  diagnosticvalues::Matrix=Array(Any, length(diagnostickeys), isempty(diagnostickeys) ? 0 : n)
)
  fnames = fieldnames(BasicContUnvParameterNState)
  BasicContUnvParameterNState(n, [fnames[i] in monitor ? true : false for i in 1:13], diagnostickeys, N, diagnosticvalues)
end

typealias ContUnvMarkovChain BasicContUnvParameterNState

codegen(f::Symbol, nstate::BasicContUnvParameterNState, monitor::Vector{Bool}) = codegen(Val{f}, nstate, monitor)

function codegen(::Type{Val{:copy}}, nstate::BasicContUnvParameterNState, monitor::Vector{Bool})
  body = []
  fnames = fieldnames(BasicContUnvParameterNState)
  local f::Symbol # f must be local to avoid compiler errors. Alternatively, this variable declaration can be omitted

  for j in 1:13
    if monitor[j]
      f = fnames[j]
      push!(body, :(getfield(_nstate, $(QuoteNode(f)))[_i] = getfield(_state, $(QuoteNode(f)))))
    end
  end

  if !isempty(nstate.diagnosticvalues)
    push!(body, :(_nstate.diagnosticvalues[:, _i] = _state.diagnosticvalues))
  end

  @gensym _copy

  quote
    function $_copy(_nstate::BasicContUnvParameterNState, _state::BasicContUnvParameterState, _i::Integer)
      $(body...)
    end
  end
end

Base.eltype{N<:Real}(::Type{BasicContUnvParameterNState{N}}) = N
Base.eltype{N<:Real}(::BasicContUnvParameterNState{N}) = N

Base.(:(==)){S<:BasicContUnvParameterNState}(z::S, w::S) =
  reduce(&, [getfield(z, n) == getfield(w, n) for n in fieldnames(S)[1:16]])

Base.isequal{S<:BasicContUnvParameterNState}(z::S, w::S) =
  reduce(&, [isequal(getfield(z, n), getfield(w, n)) for n in fieldnames(S)[1:16]])

function Base.show{N<:Real}(io::IO, nstate::BasicContUnvParameterNState{N})
  fnames = fieldnames(BasicContUnvParameterNState)
  fbool = map(n -> !isempty(getfield(nstate, n)), fnames[1:13])
  indentation = "  "

  println(io, "BasicContUnvParameterNState:")

  println(io, indentation*"eltype: $(eltype(nstate))")
  println(io, indentation*"number of states = $(nstate.n)")

  print(io, indentation*"monitored components:")
  if !any(fbool)
    println(io, " none")
  else
    print(io, "\n")
    for i in 1:13
      if fbool[i]
        println(io, string(indentation^2, fnames[i]))
      end
    end
  end

  print(io, indentation*"diagnostics:")
  if isempty(nstate.diagnostickeys)
    print(io, " none")
  else
    for k in nstate.diagnostickeys
      print(io, "\n")
      print(io, string(indentation^2, k))
    end
  end
end


text:
using BinDeps
using Compat

@BinDeps.setup

@static if is_apple()
  ostype = "apple"
  elseif is_windows()
    ostype = "windows"
  elseif is_linux()
    ostype = "linux"
  else
    error("No wkhtmltox library available for this OS")
  end

ostype *= Int==Int64 ? "64" : "32"

libwkhtml = library_dependency("libwkhtml",
                               aliases=["wkhtmltox.dll",
                                        "libwkhtmltox.so.0.12.4"])

archivemap = Dict(
 "apple32"   => "_osx-carbon-i386.pkg",
 "windows32" => "_mingw-w64-cross-win32.exe",
 "linux32"   => "_linux-generic-i386.tar.xz",
 "apple64"   => "_osx-cocoa-x86-64.pkg",
 "windows64" => "_mingw-w64-cross-win64.exe",
 "linux64"   => "_linux-generic-amd64.tar.xz",
 )

url = "https://downloads.wkhtmltopdf.org/0.12/0.12.4/" *
      "wkhtmltox-0.12.4" *
      archivemap[ostype]

# if Windows .exe file installer, trick Unpacker into seeing a zip file
downloadname = basename(url) * ( is_windows() ? ".zip" : "" )

libmap = Dict(
 "apple32"   => "",
 "windows32" => "bin/wkhtmltox.dll",
 "linux32"   => "wkhtmltox/lib/libwkhtmltox.so.0.12.4",
 "apple64"   => "",
 "windows64" => "bin/wkhtmltox.dll",
 "linux64"   => "wkhtmltox/lib/libwkhtmltox.so.0.12.4",
 )

libpath = libmap[ostype]
libfile = basename(libpath)

libdir = BinDeps.libdir(libwkhtml)
srcdir = BinDeps.srcdir(libwkhtml)
downloadsdir = BinDeps.downloadsdir(libwkhtml)

type FileCopyRule <: BinDeps.BuildStep
    src::AbstractString
    dest::AbstractString
end
Base.run(fc::FileCopyRule) = isfile(fc.dest) || cp(fc.src, fc.dest)

provides(BuildProcess,
	(@build_steps begin
    CreateDirectory(downloadsdir, true)
		FileDownloader(url, joinpath(downloadsdir, downloadname))
		CreateDirectory(srcdir, true)
    FileUnpacker(joinpath(downloadsdir, downloadname), srcdir, libpath)
		CreateDirectory(libdir, true)
    FileCopyRule(joinpath(srcdir,libpath), joinpath(libdir,libfile))
	end), libwkhtml)

# push!(BinDeps.defaults, BuildProcess)

@BinDeps.install Dict(:libwkhtml => :libwkhtml)


text:
function boundary_data(shape1::Shape{T}, inside_medium::Acoustic{T,Dim}, outside_medium::Acoustic{T,Dim},
        sim::FrequencySimulation{T,Dim,P}, ωs::Vector{T};
        dr = T(10)^(3*Dim) * eps(T), kws...) where {T<:AbstractFloat,Dim,P<:Acoustic{T,Dim}}

    in_m = inside_medium
    out_m = outside_medium

    # points just inside particles
    inside1_points = boundary_points(shape1; dr = - dr - 10*eps(T))[:]
    inside2_points = boundary_points(shape1; dr = - 10*eps(T))[:]
    inside_points = mean([inside1_points,inside2_points])

    # points just outside particles
    outside1_points = boundary_points(shape1; dr = 10*eps(T))[:]
    outside2_points = boundary_points(shape1; dr = dr + 10*eps(T))[:]
    outside_points  = mean([outside1_points,outside2_points])

    # results from just inside particles
    L = length(inside1_points)
    N = 6; #number of sets of points
    all_points = [inside1_points; inside2_points; inside_points; outside1_points; outside2_points; outside_points]

    results = run(sim, all_points, ωs; kws...)
    in1 = field(results)[1:L,:]
    in2 = field(results)[L+1:2L,:]
    in  = field(results)[2L+1:3L,:]

    out1 = field(results)[3L+1:4L,:]
    out2 = field(results)[4L+1:5L,:]
    out  = field(results)[5L+1:6L,:]

    in_pressure = FrequencySimulationResult(in, inside_points, Vector(ωs))
    out_pressure = FrequencySimulationResult(out, outside_points, Vector(ωs))

    fields = (in2 - in1)/(dr * in_m.ρ)
    in_displace = FrequencySimulationResult(fields, (inside1_points + inside2_points) ./ 2, Vector(ωs))

    fields = (out2 - out1)/(dr * out_m.ρ)
    out_displace = FrequencySimulationResult(fields, (outside1_points + outside2_points) ./ 2, Vector(ωs))

    return ([in_pressure, out_pressure], [in_displace, out_displace])
end


text:
# -*- coding: utf-8 -*-

type Board
    bb::Array{BitBoard,1} # bb[MJKI] のように使う
    WhitePieces::BitBoard
    BlackPieces::BitBoard
    OccupiedSquares::BitBoard

    WhitePiecesInHands::Array{Int,1}
    BlackPiecesInHands::Array{Int,1}

    nextMove::Int # SENTE or GOTE

    Material::Int                # incrementally updated, total material on board,
                                 # in centipawns, from white’s side of view
    square::Array{Int,1}         # incrementally updated, this array is usefull if we want to
                                 # probe what kind of piece is on a particular square.
    viewRotated::Bool            # only used for displaying the board. true or false.

    kingposW::Int
    kingposB::Int

    Board() = begin # 空っぽにする
        t = new([uint128(0) for x=1:MJGONUM]::Array{BitBoard,1},
                uint128(0),uint128(0),uint128(0),
                [0,0,0,0,0,0,0,0]::Array{Int,1},[0,0,0,0,0,0,0,0]::Array{Int,1},
                SENTE,0.0,
                [0 for x=1:NumSQ], false, 0, 0)
        t
    end
end

const BitSet = [(uint128(1) << i) for i=0:80]::Array{BitBoard,1}
#const BitSetRL90 = [BitSet[RL90[sq]] for sq=A9:I1]::Array{BitBoard,1}
#const BitSetRL45 = [BitSet[RL45[sq]] for sq=A9:I1]::Array{BitBoard,1}
#const BitSetRR45 = [BitSet[RR45[sq]] for sq=A9:I1]::Array{BitBoard,1}

function DisplayBitBoard(bb::BitBoard,isRotate::Bool)
    boardc = Array(String,NumSQ)

    for sq = 1:NumSQ
        boardc[sq] = ((bb & BitSet[sq]) != uint128(0)) ? "1": "."
    end

    println("as binary integer:")

    for i = 81:-1:1
        print("$(boardc[i])");
    end

    println("")

    println( "firstOne=", leading_zeros(bb),
            ", LastOne=", trailing_zeros(bb),
            ", BitCnt=", count_ones(bb))

    println("")

    if isRotate

        println("   ihgfedcba")
        for rank = 1:NumRank
            print("   ")
            for file = NumFile:-1:1
	        print(boardc[FR2IDX(file,rank)])
            end
            println(" $rank")
        end
        println("")
    else
        println("   abcdefghi")
        for rank = NumRank:-1:1
            print(" $rank ")
            for file = 1:NumFile
	        print(boardc[FR2IDX(file,rank)])
            end
            println("")
        end
        println("")
    end
end

function inBoard( irank::Int, ifile::Int)
    ( irank >= rank1 ) && (irank <= rank9) && (ifile >= file1) && (ifile <= file9 )
end

function setAttack( irank::Int, ifile::Int, bb::BitBoard )
    if inBoard( irank, ifile)
        bb $= BitSet[FR2IDX(ifile,irank)]
    end
    bb
end

function SquareInit(sq::Array{Int,1}, mochi_sente::Array{Int,1}, mochi_gote::Array{Int,1}, sengo::Int, bo::Board)
    bo.bb = [uint128(0) for x=1:MJGONUM]::Array{BitBoard,1}
    bo.WhitePieces = uint128(0)::BitBoard
    bo.BlackPieces = uint128(0)::BitBoard
    bo.OccupiedSquares = uint128(0)::BitBoard

    bo.WhitePiecesInHands = mochi_sente::Array{Int,1}
    bo.BlackPiecesInHands = mochi_gote::Array{Int,1}

    bo.Material = 0::Int

    bo.nextMove = sengo::Int

    bo.square = [0 for x=1:NumSQ]::Array{Int,1}
    bo.viewRotated = false::Bool

    for idx = 1:NumSQ
        piece = sq[idx]
        bo.square[idx] = piece
        if piece == MJNONE
            continue
        end
        if piece == MJOU
            bo.kingposW = idx
        elseif piece == MJGOOU
            bo.kingposB = idx
        end

        if piece > MJGOTE
            bo.Material -= MATERIAL[piece&0x0f]
        elseif piece >= MJFU
            bo.Material += MATERIAL[piece&0x0f]
        end
        bo.bb[piece] |= BitSet[idx]
    end

    for piece = MJFU:MJOU
        bo.Material -= bo.BlackPiecesInHands[piece]*MATERIAL_CAPTURE[piece]
        bo.Material += bo.WhitePiecesInHands[piece]*MATERIAL_CAPTURE[piece]
    end

    for piece = MJFU:MJRY
        bo.WhitePieces |=bo.bb[piece]
        bo.BlackPieces |=bo.bb[piece+MJGOTE]
    end

    bo.OccupiedSquares = bo.WhitePieces | bo.BlackPieces

#    for i = 1:(MJGONUM-1)
#        println(PIECENAMES[i+1])
#        DisplayBitBoard( bo.bb[i], false)
#    end
#    println("White Pieces:")
#    DisplayBitBoard( bo.WhitePieces, false)
#    println("Black Pieces:")
#    DisplayBitBoard( bo.BlackPieces, false)
#    println("Occupied Squares:")
#    DisplayBitBoard( bo.OccupiedSquares, false)

    bo
end

function InitSFEN(sfen::String, bo::Board)
    rank::Int = 9
    file::Int = 1
    promote::Int = 0
    sq::Array{Int,1}
    sengo::Int = SENTE

    fen, bw, mochi, num = split(sfen)
    sq = [0 for x = A9:I1]::Array{Int,1}

    for s in fen
        # println(s)
        if s == '/'
            rank -= 1
            file =  1
        elseif s == '+'
            promote = MJNARI
        elseif s >= '1' && s <= '9'
            file += int(s-'0')
        elseif s == 'P'
            sq[BOARDINDEX[file,rank]]=MJFU+promote
            file += 1
            promote = 0
        elseif s == 'L'
            sq[BOARDINDEX[file,rank]]=MJKY+promote
            file += 1
            promote = 0
        elseif s == 'N'
            sq[BOARDINDEX[file,rank]]=MJKE+promote
            file += 1
            promote = 0
        elseif s == 'S'
            sq[BOARDINDEX[file,rank]]=MJGI+promote
            file += 1
            promote = 0
        elseif s == 'G'
            sq[BOARDINDEX[file,rank]]=MJKI
            file += 1
            promote = 0
        elseif s == 'B'
            sq[BOARDINDEX[file,rank]]=MJKA+promote
            file += 1
            promote = 0
        elseif s == 'R'
            sq[BOARDINDEX[file,rank]]=MJHI+promote
            file += 1
            promote = 0
        elseif s == 'K'
            sq[BOARDINDEX[file,rank]]=MJOU
            file += 1
            promote = 0

            # ここから後手

        elseif s == 'p'
            sq[BOARDINDEX[file,rank]]=MJGOFU+promote
            file += 1
            promote = 0
        elseif s == 'l'
            sq[BOARDINDEX[file,rank]]=MJGOKY+promote
            file += 1
            promote = 0
        elseif s == 'n'
            sq[BOARDINDEX[file,rank]]=MJGOKE+promote
            file += 1
            promote = 0
        elseif s == 's'
            sq[BOARDINDEX[file,rank]]=MJGOGI+promote
            file += 1
            promote = 0
        elseif s == 'g'
            sq[BOARDINDEX[file,rank]]=MJGOKI
            file += 1
            promote = 0
        elseif s == 'b'
            sq[BOARDINDEX[file,rank]]=MJGOKA+promote
            file += 1
            promote = 0
        elseif s == 'r'
            sq[BOARDINDEX[file,rank]]=MJGOHI+promote
            file += 1
            promote = 0
        elseif s == 'k'
            sq[BOARDINDEX[file,rank]]=MJGOOU
            file += 1
            promote = 0
        end
    end

    if bw == "w"
        sengo = SENTE
    else # bw == "b"
        sengo = GOTE
    end

    mo_sente::Array{Int,1} = [0 for x in MJFU:MJOU]
    mo_gote::Array{Int,1}  = [0 for x in MJFU:MJOU]

    numberOfMochi::Int = 0
    for s in mochi
        if s == '-'
            break # 持ち駒なし
        elseif s == '1'
            if numberOfMochi == 1
                numberOfMochi = 11
            else
                numberOfMochi = 1
            end
        elseif s == '0'
            if numberOfMochi == 1
                numberOfMochi = 10
            end
        elseif '2' <= s <= '9'
            if numberOfMochi == 1
                numberOfMochi = 10 + int(s - '0')
            else
                numberOfMochi = int(s - '0')
            end
        elseif s == 'P'
            mo_sente[MJFU] = numberOfMochi == 0 ? 1: numberOfMochi
            numberOfMochi = 0
        elseif s == 'L'
            mo_sente[MJKY] = numberOfMochi == 0 ? 1: numberOfMochi
            numberOfMochi = 0
        elseif s == 'N'
            mo_sente[MJKE] = numberOfMochi == 0 ? 1: numberOfMochi
            numberOfMochi = 0
        elseif s == 'S'
            mo_sente[MJGI] = numberOfMochi == 0 ? 1: numberOfMochi
            numberOfMochi = 0
        elseif s == 'G'
            mo_sente[MJKI] = numberOfMochi == 0 ? 1: numberOfMochi
            numberOfMochi = 0
        elseif s == 'B'
            mo_sente[MJKA] = numberOfMochi == 0 ? 1: numberOfMochi
            numberOfMochi = 0
        elseif s == 'R'
            mo_sente[MJHI] = numberOfMochi == 0 ? 1: numberOfMochi
            numberOfMochi = 0
        elseif s == 'K' # 普通はない
            mo_sente[MJOU] = numberOfMochi == 0 ? 1: numberOfMochi
            numberOfMochi = 0

            # ここから後手

        elseif s == 'p'
            mo_gote[MJFU] = numberOfMochi == 0 ? 1: numberOfMochi
            numberOfMochi = 0
        elseif s == 'l'
            mo_gote[MJKY] = numberOfMochi == 0 ? 1: numberOfMochi
            numberOfMochi = 0
        elseif s == 'n'
            mo_gote[MJKE] = numberOfMochi == 0 ? 1: numberOfMochi
            numberOfMochi = 0
        elseif s == 's'
            mo_gote[MJGI] = numberOfMochi == 0 ? 1: numberOfMochi
            numberOfMochi = 0
        elseif s == 'g'
            mo_gote[MJKI] = numberOfMochi == 0 ? 1: numberOfMochi
            numberOfMochi = 0
        elseif s == 'b'
            mo_gote[MJKA] = numberOfMochi == 0 ? 1: numberOfMochi
            numberOfMochi = 0
        elseif s == 'r'
            mo_gote[MJHI] = numberOfMochi == 0 ? 1: numberOfMochi
            numberOfMochi = 0
        elseif s == 'k' # 普通はない
            mo_gote[MJOU] = numberOfMochi == 0 ? 1: numberOfMochi
            numberOfMochi = 0
        end
    end

    bo = SquareInit(sq,mo_sente,mo_gote,sengo,bo)
    bo
end

function DisplayMochiGoma( bo::Board, sengo::Int)
    senteNum::Int = 0
    goteNum::Int = 0
    if sengo == SENTE
        print("\e[34m", TEBANSTR[1], "\e[30m 持駒　")
        for i = MJOU:-1:MJFU
            mochi = bo.WhitePiecesInHands[i]
            if mochi > 0
                senteNum += 1

                print(MOCHISTR[i])
                if mochi > 1
                    print(mochi)
                end
            end
        end
        if senteNum == 0
            print("なし")
        end
        println()
    elseif sengo == GOTE
        print(TEBANSTR[2], " 持駒　")
        for i = MJOU:-1:MJFU
            mochi = bo.BlackPiecesInHands[i]
            if mochi > 0
                senteNum += 1

                print(MOCHISTR[i])
                if mochi > 1
                    print(mochi)
                end
            end
        end
        if senteNum == 0
            print("なし")
        end
        println()
    else
        println("Illegal Mochi Goma!")
    end
end

function DisplayBoard(bo::Board)
    println()
    DisplayMochiGoma(bo,GOTE)
    println(" ９　  ８　  ７　  ６　  ５　  ４　  ３　  ２　  １　 ")
    sq::Int = 0

    for i = 1:NumRank
        for j = 1:NumFile
            sq += 1
            koma = bo.square[sq]
            komamoji::UTF8String = doubleZenkakuSpace
            if koma == MJNONE
	        print("\e[36m ..  \e[30m ")
            else
                komamoji = KOMASTR[koma&0x0f]
                if koma < MJGOTE
                    print("\e[34m+")
                    print(KOMASTR[koma&0x0f],"\e[30m ")
                elseif koma > MJGOTE
                    print("\e[30m-")
                    print(KOMASTR[koma&0x0f],"\e[30m ")
                end
            end
        end
        println(DANSTR[i])
    end
    DisplayMochiGoma(bo,SENTE)
    println("Material = ", bo.Material)
    println()
end

function BBTest()
    #InitTables()
    bb::BitBoard = uint128(0)
    #println(bb)
    bb2 = setAttack( 5, 5, bb)
    #println(bb2)
    bo = Board()
    #DisplayBitBoard(bb2,false)
    for x=A9:I1
        b = BitSet[x]
        sqname = SQUARENAME[x]
        #println("SQ[$sqname] BitSet\[$x\] = ", b, " count one = ", count_ones(b),
        #        " leading_zeros = ", leading_zeros(b), " trailing_zeros = ", trailing_zeros(b))
        #DisplayBitBoard(b,false)
    end
    for ra = 1:NumRank, fi = 1:NumFile
        #println("BOARDINDEX[$fi,$ra] = $(BOARDINDEX[fi,ra])")
    end
    #board = SquareInit(HIRATE,[0,0,0,0,0,0,0,0],[0,0,0,0,0,0,0,0],SENTE,bo)
    sfenHirate = "lnsgkgsnl/1r5b1/ppppppppp/9/9/9/PPPPPPPPP/1B5R1/LNSGKGSNL w - 1"
    board = InitSFEN(sfenHirate, bo)
    #DisplayBoard(board)
    #sfen = "8l/1l+R2P3/p2pBG1pp/kps1p4/Nn1P2G2/P1P1P2PP/1PS6/1KSG3+r1/LN2+p3L w Sbgn3p 124"::String
    #bo2 = Board()
    #board2 = InitSFEN(sfen, bo2)
    #DisplayBoard(board2)
end

#BBTest()

max_stars_repo_name: kimrin/JapaneseChess

max_issues_repo_name: kimrin/JapaneseChess

max_forks_repo_name: kimrin/JapaneseChess


text:
# # Input and output formats
#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/guide/@__NAME__.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/guide/@__NAME__.ipynb)

# This section provides an overview of the input and output formats
# supported by DFTK, usually via integration with a third-party library.
#
# ## Reading input formats supported by ASE
# ASE is short for the
# [atomistic simulation environment](https://wiki.fysik.dtu.dk/ase/index.html),
# a Python package to simplify the process of setting up, running and
# analysing results from atomistic simulations across different programs.
# If it is installed it is automatically used by DFTK in order to read a wide range
# of input files, such as xyz, CIF, input files of various other codes
# (e.g. Quantum Espresso, VASP, ABINIT, CASTEP, ...).
# The full list of formats
# can be found in the [ASE IO documentation](https://wiki.fysik.dtu.dk/ase/ase/io/io.html).
#
# As an example we start the calculation of a simple antiferromagnetic iron crystal
# using a Quantum-Espresso input file, [Fe_afm.pwi](Fe_afm.pwi).
# From this file the lattice, atomic positions and the initial magnetisation are read.
# For more details about calculations on magnetic systems
# using collinear spin, see [Collinear spin and magnetic systems](@ref).
#
# First we parse the Quantum Espresso input file to DFTK datastructures:

using DFTK

qe_input  = "Fe_afm.pwi"
atoms     = load_atoms(qe_input)
positions = load_positions(qe_input)
lattice   = load_lattice(qe_input)
magnetic_moments = load_magnetic_moments(qe_input)

# At this point a file of any format supported by ASE could be passed instead,
# e.g. an `xyz` file or an ABINIT input file. Behind the scenes ASE takes care
# to select the right parser and extract the required structural information.

# Next we attach the pseudopotential information, since this information is currently
# not exposed inside the ASE datastructures.
# See [Creating slabs with ASE](@ref) for more details.

atoms = map(atoms) do el
    @assert el.symbol == :Fe
    ElementPsp(:Fe, psp=load_psp("hgh/pbe/fe-q16.hgh"))
end;

# Finally we run the calculation.

model = model_LDA(lattice, atoms, positions; magnetic_moments, temperature=0.01)
basis = PlaneWaveBasis(model; Ecut=10, kgrid=(2, 2, 2))
ρ0 = guess_density(basis, magnetic_moments)
scfres = self_consistent_field(basis, tol=1e-4, ρ=ρ0, mixing=KerkerMixing());

# !!! note "DFTK and ASE"
#     DFTK also supports using ASE to setup a calculation
#     and provides two-way conversion routines between key DFTK
#     and ASE datastructures. See [Creating slabs with ASE](@ref) for details.
#
# ## Writing VTK files for visualization
# For visualizing the density or the Kohn-Sham orbitals DFTK supports storing
# the result of an SCF calculations in the form of VTK files.
# These can afterwards be visualized using tools such
# as [paraview](https://www.paraview.org/).
# Using this feature requires
# the [WriteVTK.jl](https://github.com/jipolanco/WriteVTK.jl/) Julia package.

using WriteVTK
save_scfres("iron_afm.vts", scfres; save_ψ=true);

# This will save the iron calculation above into the file `iron_afm.vts`,
# using `save_ψ=true` to also include the KS orbitals.

# ## Writing and reading JLD2 files
# The full state of a DFTK self-consistent field calculation can be
# stored on disk in form of an [JLD2.jl](https://github.com/JuliaIO/JLD2.jl) file.
# This file can be read from other Julia scripts
# as well as other external codes supporting the HDF5 file format
# (since the JLD2 format is based on HDF5).

using JLD2
save_scfres("iron_afm.jld2", scfres);

# Since such JLD2 can also be read by DFTK to start or continue a calculation,
# these can also be used for checkpointing or for transferring results
# to a different computer.
# See [Saving SCF results on disk and SCF checkpoints](@ref) for details.

# (Cleanup files generated by this notebook)
rm("iron_afm.vts")
rm("iron_afm.jld2")


text:
"""
Given: A DNA string s of length at most 1000 bp.

Return: The reverse complement s^c of s.
"""

function reverse_complement(s::AbstractString)
	rules = Dict('A'=>'T','T'=>'A','G'=>'C','C'=>'G')
	return map(x->rules[x],reverse(s))
end

max_stars_repo_name: madisonestabrook/Julia

max_issues_repo_name: madisonestabrook/Julia

max_forks_repo_name: madisonestabrook/Julia


text:
my_integer = 5
my_float = 5.
my_range = 1:7
my_vector = [1, 3, 5]

# test
function g(x)
    if x>0 || x<-1
        return "nice"
    else
        return "WTF"
    end
end
g(-10)

max_stars_repo_name: paradoxuu/lalala

max_issues_repo_name: paradoxuu/lalala

max_forks_repo_name: paradoxuu/lalala


text:
module TensorCast

# This speeds up loading a bit, on Julia 1.5, about 1s in my test.
# https://github.com/JuliaPlots/Plots.jl/pull/2544/files
if isdefined(Base, :Experimental) && isdefined(Base.Experimental, Symbol("@optlevel"))
    @eval Base.Experimental.@optlevel 1
end

export @cast, @reduce, @matmul, @pretty

using LinearAlgebra, Random

using MacroTools, StaticArrays

using TransmuteDims, LazyStack
using LazyStack: stack_iter

include("macro.jl")
include("pretty.jl")

module Fast # shield non-macro code from @optlevel 1
    using ..TensorCast: pretty
    using LinearAlgebra, StaticArrays

    include("slice.jl")
    export sliceview, slicecopy, copy_glue, glue!, iscodesorted, countcolons

    include("view.jl")
    export diagview, mul!, rview, star

    include("static.jl")
    export static_slice, static_glue

end
using .Fast
const mul! = Fast.mul!

include("warm.jl") # saves 3s in my test

end # module


text:
include("../imports.jl")

@testset ExtendedTestSet "`DiscriminativeLRs`" begin
    model = Chain(Dense(3, 5), Dense(5, 3))
    pg = ParamGroups(IndexGrouper([1, 2]), model)
    o = Optimiser(
        DiscriminativeLRs(pg, Dict(1 => 0., 2 => 1.)),
        Descent(0.1)
    )
    x1 = model[1].weight
    x2 = model[2].weight
    # Weight of layer 1 has zeroed gradient
    @test apply!(o, x1, ones(size(x1))) == zeros(size(x1))
    # Weight of layer 2 has regular gradient
    @test apply!(o, x2, ones(size(x2))) != fill(0.1, size(x1))
end


text:
"""
    gauge(x)

Creates an `EChart` where a value is plotted as a needle indicating the value along some range.

## Methods
```julia
gauge(x::Union{Missing, Real})
```

## Arguments
* `breakpoints::AbstractVector = [0.2, 0.8, 1]` : percentage along gauge to switch indicator colors
* `colors::AbstractVector = ["#91c7ae", "#63869e", "#c23531"]` : colors for breakpoints
* `kwargs` : varargs to set any field of resulting `EChart` struct

## Notes

## Examples
```julia
gg = gauge(27.64)
```
"""
function gauge(x::Union{Missing, Real};
               breakpoints::AbstractVector = [0.2, 0.8, 1],
               colors::AbstractVector = ["#91c7ae", "#63869e", "#c23531"],
               kwargs...)

    length(breakpoints) != length(colors) ? error("breakpoints and colors vectors need to have the same length.") : nothing

    ec = newplot(kwargs, ec_charttype = "gauge")

    ec.series = [GaugeSeries(name = "Series 1",
                        _type = "gauge",
                        data = [Dict("value" => x)],
                        axisLine = AxisLine(lineStyle = LineStyle(width = 10, color = [[b,c] for (b,c) in zip(breakpoints, colors)]))
                        )
                ]

    return ec

end


text:
const ALL_TASK_OFFSETS =
    Dict{Tuple{String, VersionNumber}, Dict{Symbol, Int}}()

include("linux-x86-v1_3_1.jl")
include("linux-x86-v1_4_2.jl")
include("linux-x86-v1_5_3.jl")
include("linux-x86-v1_5_4.jl")
include("linux-x86-v1_6_1.jl")
include("linux-x86-v1_8_0.jl")

include("linux-x86_64-v1_3_1.jl")
include("linux-x86_64-v1_4_2.jl")
include("linux-x86_64-v1_5_3.jl")
include("linux-x86_64-v1_5_4.jl")
include("linux-x86_64-v1_6_1.jl")
include("linux-x86_64-v1_8_0.jl")

include("darwin-x86_64-v1_3_1.jl")
include("darwin-x86_64-v1_4_2.jl")
include("darwin-x86_64-v1_5_3.jl")
include("darwin-x86_64-v1_5_4.jl")
include("darwin-x86_64-v1_6_1.jl")

include("windows-x86-v1_3_1.jl")
include("windows-x86-v1_4_2.jl")
include("windows-x86-v1_5_3.jl")
include("windows-x86-v1_5_4.jl")
include("windows-x86-v1_6_1.jl")

include("windows-x86_64-v1_3_1.jl")
include("windows-x86_64-v1_4_2.jl")
include("windows-x86_64-v1_5_3.jl")
include("windows-x86_64-v1_5_4.jl")
include("windows-x86_64-v1_6_1.jl")

const ARCH = @static if string(Sys.ARCH)[1] == 'i'
    "x86"
else
    string(Sys.ARCH)
end

const PLATFORM = @static if Sys.islinux()
    "linux-" * ARCH
elseif Sys.iswindows()
    "windows-" * ARCH
elseif Sys.isapple()
    "darwin-" * ARCH
end

function find_offsets()
    candi = [v for (p, v) in keys(ALL_TASK_OFFSETS) if p == PLATFORM]
    sort!(candi)
    candi = [v for v in candi if v <= VERSION]
    isempty(candi) && error("No suitable offsets")
    ALL_TASK_OFFSETS[(PLATFORM, candi[end])]
end

const TASK_OFFSETS = find_offsets()

## utilities

# getter
function internal_getfield(task, field, ::Type{Ptr})
    ccall((:jl_getfield_ptr, libtask_julia),
          Ptr{Cvoid}, (Any, Csize_t),
          task, TASK_OFFSETS[field])
end

function internal_getfield(task, field, ::Type{Int32})
    ccall((:jl_getfield_int32_t, libtask_julia),
          Int32, (Any, Csize_t), task, TASK_OFFSETS[field])
end

# setter
function internal_setfield(task, field, ::Nothing)
    haskey(TASK_OFFSETS, field) && ccall(
        (:jl_setfield_nothing, libtask_julia),
        Any, (Any, Csize_t), task, TASK_OFFSETS[field])
end

function internal_setfield(task, field, p::Ptr{Nothing})
    p == C_NULL || error("this function is only for setting NULL value")
    haskey(TASK_OFFSETS, field) && ccall(
        (:jl_setfield_null, libtask_julia),
        Any, (Any, Csize_t), task, TASK_OFFSETS[field])
end

function internal_setfield(task, field, val::Csize_t)
    ccall((:jl_setfield_size_t, libtask_julia),
          Any, (Any, Csize_t, Csize_t),
          task, TASK_OFFSETS[field], val)
end

function internal_setfield(task, field, val::Int32)
    ccall((:jl_setfield_int32_t, libtask_julia),
          Any, (Any, Csize_t, Int32),
          task, TASK_OFFSETS[field], val)
end

function memset(task, val, offset_field; end_field=:END)
    haskey(TASK_OFFSETS, offset_field) && ccall(
        (:jl_memset, libtask_julia),
        Cvoid, (Any, Csize_t, Csize_t, Int),
        task, TASK_OFFSETS[offset_field], TASK_OFFSETS[end_field], val)
end


text:
module IO

#import GeoPhyInv.Data

mutable struct SegyHeader
	tracl::Int32
	tracr::Int32
	fldr::Int32
	tracf::Int32
	ep::Int32
	cdp::Int32
	cdpt::Int32
	trid::Int16
	nva::Int16
	nhs::Int16
	duse::Int16
	offset::Int32
	gelev::Int32
	selev::Int32
	sdepth::Int32
	gdel::Int32
	sdel::Int32
	swdep::Int32
	gwdep::Int32
	scalel::Int16
	scalco::Int16
	sx::Int32
	sy::Int32
	gx::Int32
	gy::Int32
	counit::Int16
	wevel::Int16
	swevel::Int16
	sut::Int16
	gut::Int16
	sstat::Int16
	gstat::Int16
	tstat::Int16
	laga::Int16
	lagb::Int16
	delrt::Int16
	muts::Int16
	mute::Int16
	ns::Int16
	dt::Int16
	gain::Int16
	igc::Int16
	igi::Int16
	corr::Int16
	sfs::Int16
	sfe::Int16
	slen::Int16
	styp::Int16
	stas::Int16
	stae::Int16
	tatyp::Int16
	afilf::Int16
	afils::Int16
	nofilf::Int16
	nofils::Int16
	lcf::Int16
	hcf::Int16
	lcs::Int16
	hcs::Int16
	year::Int16
	day::Int16
	hour::Int16
	minute::Int16
	sec::Int16
	timbas::Int16
	trwf::Int16
	grnors::Int16
	grnofr::Int16
	grnlof::Int16
	gaps::Int16
	otrav::Int16
	d1::Float32
	f1::Float32
	d2::Float32
	f2::Float32
	ungpow::Float32
	unscale::Float32
	ntr::Int32
	mark::Int16
	unass::Int16
end

segy_count = Dict{AbstractString,Int32}()
segy_count["tracl"]  = 0
segy_count["tracr"]  = 4
segy_count["fldr"]   = 8
segy_count["tracf"]  = 12
segy_count["ep"]     = 16
segy_count["cdp"]    = 20
segy_count["cdpt"]   = 24
segy_count["trid"]   = 28
segy_count["nva"]    = 30
segy_count["nhs"]    = 32
segy_count["duse"]   = 34
segy_count["offset"] = 36
segy_count["gelev"]  = 40
segy_count["selev"]  = 44
segy_count["sdepth"] = 48
segy_count["gdel"]   = 52
segy_count["sdel"]   = 56
segy_count["swdep"]  = 60
segy_count["gwdep"]  = 64
segy_count["scalel"] = 68
segy_count["scalco"] = 70
segy_count["sx"]     = 72
segy_count["sy"]     = 76
segy_count["gx"]     = 80
segy_count["gy"]     = 84
segy_count["counit"] = 88
segy_count["wevel"]  = 90
segy_count["swevel"] = 92
segy_count["sut"]    = 94
segy_count["gut"]    = 96
segy_count["sstat"]  = 98
segy_count["gstat"]  = 100
segy_count["tstat"]  = 102
segy_count["laga"]   = 104
segy_count["lagb"]   = 106
segy_count["delrt"]  = 108
segy_count["muts"]   = 110
segy_count["mute"]   = 112
segy_count["ns"]     = 114
segy_count["dt"]     = 116
segy_count["gain"]   = 118
segy_count["igc"]    = 120
segy_count["igi"]    = 122
segy_count["corr"]   = 124
segy_count["sfs"]    = 126
segy_count["sfe"]    = 128
segy_count["slen"]   = 130
segy_count["styp"]   = 132
segy_count["stas"]   = 134
segy_count["stae"]   = 136
segy_count["tatyp"]  = 138
segy_count["afilf"]  = 140
segy_count["afils"]  = 142
segy_count["nofilf"] = 144
segy_count["nofils"] = 146
segy_count["lcf"]    = 148
segy_count["hcf"]    = 150
segy_count["lcs"]    = 152
segy_count["hcs"]    = 154
segy_count["year"]   = 156
segy_count["day"]    = 158
segy_count["hour"]   = 160
segy_count["minute"] = 162
segy_count["sec"]    = 164
segy_count["timbas"] = 166
segy_count["trwf"]   = 168
segy_count["grnors"] = 170
segy_count["grnofr"] = 172
segy_count["grnlof"] = 174
segy_count["gaps"]   = 176
segy_count["otrav"]  = 178
segy_count["d1"]     = 180
segy_count["f1"]     = 184
segy_count["d2"]     = 188
segy_count["f2"]     = 192
segy_count["ungpow"] = 196
segy_count["unscale"]= 200
segy_count["ntr"]    = 204
segy_count["mark"]   = 208
segy_count["unass"]  = 210
segy_count["trace"]  = 240

function InitSegyHeader()
	h = SegyHeader(0,0,0,0,0,0,0,0,0,0,
	0,0,0,0,0,0,0,0,0,0,
	0,0,0,0,0,0,0,0,0,0,
	0,0,0,0,0,0,0,0,0,0,
	0,0,0,0,0,0,0,0,0,0,
	0,0,0,0,0,0,0,0,0,0,
	0,0,0,0,0,0,0,0,0,0,
	0,0,0,0,0,0,0,0,0,0)
	return h
end

function GrabSegyHeader(stream,swap_bytes,nt,file_header_size,j)
	position = file_header_size + (240+nt*4)*(j-1)
	seek(stream,position)
	h = InitSegyHeader()
	if (swap_bytes == false)
		for field in fieldnames(typeof(h))
			setfield!(h, field, read(stream,typeof(getfield(h,field))) )
		end
	else
		for field in fieldnames(typeof(h))
			setfield!(h, field, bswap(read(stream,typeof(getfield(h,field)))) )
		end
	end

	return h
end

function PutSegyHeader(stream,h,nt,file_header_size,j)
	position = file_header_size + (240+nt*4)*(j-1)
	seek(stream,position)
	for field in fieldnames(typeof(h))
		write(stream,getfield(h,field))
	end
end

"""
Read an input SeismicUnix file and output
the time-domain data as `Data.TD`.

# Arguments
* `fname::AbstractString` : name of the SeismicUnix file
TODO: conversion to `TD`, currently the ouput is just a data matrix.
"""
function readsu_data(;
		     fname::AbstractString="", verbose::Bool=false
		     )
end

"""
Read a SU file and return trace headers along with data matrix
TODO: use bswap accordingly?
      extend it to SEGY format as well
"""
function readsu(fname, verbose=false)

	file_hsize=0
	swap_bytes = false

	stream = open(fname)
	seek(stream, segy_count["ns"]+file_hsize)
	nt = read(stream,Int16)
	total = 60 + nt
	nrecords = round(Int,(filesize(stream)-file_hsize)/4/total)
	verbose && println("number of traces: ",nrecords)
	verbose && println("number of samples per trace: ",nt)

	h_segy = Array{SegyHeader}(undef, nrecords)
	seek(stream,file_hsize + segy_count["trace"])
	h_segy[1] = GrabSegyHeader(stream,swap_bytes,nt,file_hsize,1)
	dt = h_segy[1].dt/1000000 # convert to seconds

	d=zeros(nt, nrecords)
	# reading data and headers
	for j=1:nrecords
		position = file_hsize + total*(j-1)*4 + segy_count["trace"]
		# goto trace position
		seek(stream,position)
		# read trace
		for i in 1:nt
			d[i,j] = read(stream,Float32)
		end
		# read header, seek is inside this function
		h_segy[j] = GrabSegyHeader(stream,swap_bytes,nt,file_hsize,j)
	end
	close(stream)

	return reshape(d, (nt, nrecords)), h_segy

end

"""
Write a SU file using data and headers

# Arguments

* `fname` : filename
* `d` : data matrix, output of readsu, for example
* `h_segy` : header vector, see output of readsu
"""
function writesu(fname, d, h_segy::Vector{SegyHeader}=[InitSegyHeader() for irec in 1:size(d,2)])

	nt = size(d, 1)
	total = 60 + nt
	#nrecords = round(Int,(filesize(stream)-file_hsize)/4/total)
	nrecords = size(d, 2)

	file_hsize=0
	swap_bytes = false

	stream = open(fname, "w")

	# reading data and headers
	for j=1:nrecords
		# add nt in header if absent
		(h_segy[j].ns == 0) && (h_segy[j].ns = nt)

		PutSegyHeader(stream,h_segy[j],nt,file_hsize,j)

		# write down the trace
		dvecsp = convert(Vector{Float32}, d[:,j])
		position = file_hsize + total*(j-1)*4 + segy_count["trace"]
		seek(stream,position)
		write(stream, dvecsp)
	end
	close(stream)
end

#function TD(d::Data.TD, attrib=:su)
#
#	nrecords = sum(d.acqgeom.nr)
#	nss = d.acqgeom.nss
#
#	# create a data matrix of first field
#	dp = hcat(d.d[collect(1:nss),1][:,:]...);
#
#	# headers
#	h_segy = [InitSegyHeader() for irec in 1:nrecords]
#	irec = 0
#	for iss=1:nss
#		for ir=1:d.acqgeom.nr[iss]
#			irec += 1
#			h_segy[irec].ns = length(d.tgrid)
#			h_segy[irec].dt = step(d.tgrid)*1000000
#		end
#	end
#	writesu(fname, d, h_segy)
#end

end # module


text:
# Fermat's primality test: (Test for a number to be prime or not, based on Fermat's little theorem)
function fermatPrimalityTest(n, iterations = 5)
    if n < 4
        return n == 2 || n == 3
    end

    for i = 1:iterations
        a = rand(2:(n - 2), 1)[1]
        if ((BigInt(a) ^ (n - 1)) % n) != 1
            return false
        end
    end
    return true
end

# Fermat's factorization: (Representation of an odd integer as the difference of two squares)
function fermatFactorization(n::Int64)
    a = ceil(sqrt(n))
    b2 = a * a - n
    b = round(sqrt(b2))
    while (b * b != b2)
        a = a + 1
        b2 = a * a - n
        b = round(sqrt(b2))
    end
    return Int64(a - b)
end


text:
# This code is from https://github.com/JuliaAudio/SampledSignals.jl/blob/2d078e86489232f77af3696f0ea6d0e34016e7b0/test/support/util.jl
# Copyright (c) 2015: Spencer Russell. Released under the MIT "Expat" License.
mutable struct DummySampleSink{T} <: SampleSink
    samplerate::Float64
    buf::Array{T, 2}
end

DummySampleSink(eltype, samplerate::Number, channels::Int) =
    DummySampleSink{eltype}(samplerate, Array{eltype}(undef, 0, channels))
function DummySampleSink(eltype, samplerate::Union{typeof(1u"Hz"), typeof(1u"kHz"),
                                                   typeof(1.0u"Hz"), typeof(1.0u"kHz")},
                         nchannels::Int)
    samplerate = ustrip(uconvert(u"Hz", samplerate))
    DummySampleSink(eltype, samplerate, nchannels)
end

samplerate(sink::DummySampleSink) = sink.samplerate
nchannels(sink::DummySampleSink) = size(sink.buf, 2)
Base.eltype(sink::DummySampleSink{T}) where T = T

function SampledSignals.unsafe_write(sink::DummySampleSink, buf::Array,
                                     frameoffset, framecount)
    eltype(buf) == eltype(sink) || error("buffer type ($(eltype(buf))) doesn't match sink type ($(eltype(sink)))")
    nchannels(buf) == nchannels(sink) || error("buffer channel count ($(nchannels(buf))) doesn't match sink channel count ($(nchannels(sink)))")

    sink.buf = vcat(sink.buf, view(buf, (1:framecount) .+ frameoffset, :))

    framecount
end

# ######################
# Plotting
# ######################

PlotSpectroTemporal(x::AuditoryStimuli.DummySampleSink; kwargs...) = PlotSpectroTemporal(x.buf, x.samplerate; kwargs...)

function plot(x::AuditoryStimuli.DummySampleSink;
                xlab::String = "Time (s)",
                ylab::String = "Amplitude",
                kwargs...)

    t = 1:size(x.buf, 1)
    t = t ./ x.samplerate

    plot(t, x.buf, xlab = xlab, ylab = ylab; kwargs...)
end


text:
using SimulatedAnnealing
using Base.Test

tests = ["stable",
         "cooling",
         ]

println("Running tests:")

for t in tests
    println(" * $(t)")
    include("$(t).jl")
end


text:
"""
    allocate_data_dict( blocks::Array{String}, N_to_read::Integer,
                        snap_info::Array{InfoLine}, no_mass_block::Bool )

Helper function to allocate the data `Dict`.
"""
function allocate_data_dict(blocks::Array{String}, N_to_read::Integer,
                            snap_info::Array{InfoLine}, no_mass_block::Bool)

    # prepare dictionary for particle storage
    d = Dict{String, VecOrMat{T} where T}()

    for block ∈ blocks
        dim   = snap_info[getfield.(snap_info, :block_name) .== block][1].n_dim
        dtype = snap_info[getfield.(snap_info, :block_name) .== block][1].data_type
        if dim == 1
            d[block] = Vector{dtype}(undef, N_to_read)
        else
            d[block] = Matrix{dtype}(undef, dim, N_to_read)
        end
    end

    # allocate mass array, if it's not in a block
    if no_mass_block
        d["MASS"] = Vector{Float32}(undef, N_to_read)
    end

    return d
end


text:
using BGZFStreams
using Base.Test
using Compat

@testset "VirtualOffset" begin
    voff = VirtualOffset(0, 0)
    @test voff == VirtualOffset(0, 0)
    @test voff != VirtualOffset(0, 1)
    @test voff <= VirtualOffset(0, 1)
    @test voff != VirtualOffset(1, 0)
    @test voff <= VirtualOffset(1, 0)
    @test voff[1] == 0
    @test voff[2] == 0
    @test string(voff) == "BGZFStreams.VirtualOffset(0, 0)"

    voff += 1
    @test voff[1] == 0
    @test voff[2] == 1

    voff = VirtualOffset(1234, 555)
    @test voff[1] == 1234
    @test voff[2] == 555
    @test string(voff) == "BGZFStreams.VirtualOffset(1234, 555)"

    @test_throws ArgumentError VirtualOffset(1 << 48, 0)
    @test_throws ArgumentError VirtualOffset(0, 1 << 16)
end

@testset "BGZFStream" begin
    filename = joinpath(dirname(@__FILE__), "bar.bgz")
    stream = BGZFStream(filename, "r")
    @test read(stream, UInt8) === UInt8('b')
    @test read(stream, UInt8) === UInt8('a')
    @test read(stream, UInt8) === UInt8('r')
    @test eof(stream)
    @test flush(stream) === nothing
    @test close(stream) === nothing
    @test string(stream) == "BGZFStreams.BGZFStream{IOStream}(<mode=read>)"
    @test_throws ArgumentError read(stream, UInt8)

    stream = BGZFStream(filename, "r")
    data = zeros(Int8, 3)
    unsafe_read(stream, pointer(data), 3)
    @test data == UInt8['b', 'a', 'r']
    close(stream)
    @test_throws ArgumentError unsafe_read(stream, pointer(data), 3)

    stream = BGZFStream(filename, "r")
    @test virtualoffset(stream) === VirtualOffset(0, 0)
    read(stream, UInt8)
    read(stream, UInt8)
    @test virtualoffset(stream) === VirtualOffset(0, 2)
    seek(stream, VirtualOffset(0, 1))
    @test read(stream, UInt8) === UInt8('a')
    @test read(stream, UInt8) === UInt8('r')
    seekstart(stream)
    @test read(stream, UInt8) === UInt8('b')
    @test read(stream, UInt8) === UInt8('a')
    @test read(stream, UInt8) === UInt8('r')
    close(stream)

    # Empty data.
    empty_block = copy(BGZFStreams.EOF_BLOCK)
    stream = BGZFStream(IOBuffer(empty_block))
    @test_throws EOFError read(stream, UInt8)
    @test_throws EOFError read(stream, UInt8)
    stream = BGZFStream(IOBuffer(vcat(empty_block, empty_block)))
    @test_throws EOFError read(stream, UInt8)
    @test_throws EOFError read(stream, UInt8)
    stream = BGZFStream(IOBuffer(empty_block))
    @test isempty(read(stream))
    @test isempty(read(stream))
    stream = BGZFStream(IOBuffer(vcat(empty_block, empty_block)))
    @test isempty(read(stream))
    @test isempty(read(stream))

    # invalid data
    stream = BGZFStream(IOBuffer(vcat([0x1f, 0x00, 0x08, 0x04], rand(UInt8, 6))))
    @test_throws BGZFDataError read(stream)  # invalid gzip flag
    stream = BGZFStream(IOBuffer(vcat([0x1f, 0x8b, 0x00, 0x04], rand(UInt8, 6))))
    @test_throws BGZFDataError read(stream)  # invalid compression method
    stream = BGZFStream(IOBuffer(vcat([0x1f, 0x8b, 0x08, 0x01], rand(UInt8, 6))))
    @test_throws BGZFDataError read(stream)  # invalid flag

    filename = tempname()
    try
        stream = BGZFStream(filename, "w")
        @test virtualoffset(stream) == VirtualOffset(0, 0)
        @test write(stream, 'b') === 1
        @test write(stream, 'a') === 1
        @test write(stream, 'r') === 1
        @test virtualoffset(stream) == VirtualOffset(0, 3)
        @test write(stream, "baz") === 3
        @test eof(stream)
        @test flush(stream) === nothing
        @test close(stream) === nothing
        @test string(stream) == "BGZFStreams.BGZFStream{IOStream}(<mode=write>)"
    catch
        rethrow()
    finally
        if isfile(filename)
            rm(filename)
        end
    end

    # Pipe I/O
    filename = joinpath(dirname(@__FILE__), "bar.bgz")
    @test open(p -> read(BGZFStream(p)), `cat $(filename)`) == UInt8['b', 'a', 'r']
    @test_throws ArgumentError open(p -> virtualoffset(BGZFStream(p)), `cat $(filename)`)
    @test_throws ArgumentError open(p -> seek(BGZFStream(p), VirtualOffset(0, 0)), `cat $(filename)`)
    # TODO: test writing to a pipe?

    buffer = IOBuffer()
    stream = BGZFStream(buffer, "w")
    write(stream, "foo")
    close(stream)
    @test !isopen(buffer)
    @test_throws ArgumentError write(stream, 0x01)

    # Round trip
    for n in [0, 1, 2, 5, 10, 50, 100, 10_000, 100_000, 1_000_000]
        data = rand(0x00:0xf0, n)

        # bulk read/write
        buffer = IOBuffer()
        stream = BGZFStream(buffer, "w")
        # HACK: do not close the buffer after the stream is closed
        stream.onclose = io -> nothing
        write(stream, data)
        close(stream)
        seekstart(buffer)
        stream = BGZFStream(buffer)
        @test data == read(stream)
        close(stream)

        # read/write byte by byte
        buffer = IOBuffer()
        stream = BGZFStream(buffer, "w")
        # HACK: do not close the buffer after the stream is closed
        stream.onclose = io -> nothing
        for x in data
            write(stream, x)
        end
        close(stream)
        seekstart(buffer)
        stream = BGZFStream(buffer)
        data′ = UInt8[]
        while !eof(stream)
            push!(data′, read(stream, UInt8))
        end
        @test data == data′
        close(stream)
    end
end


text:
using Weave

weave("src/index.Rmd", informat="markdown",
      out_path = "src/", doctype = "github")


text:
using CUDAnative, CUDAdrv, CuArrays, Test

function vec_sum(a, b, res)
    tid = threadIdx().x
    i = (blockIdx().x-1) * blockDim().x + threadIdx().x

    tile = @cuDynamicSharedMem(Float64, res)

    tile[tid] = a[i]

    sync_threads()

    s = 1
    while s < blockDim().x
        if (tid % (2*s) == 1)
            tile[tid] += tile[tid+s]
        end
        sync_threads()
        s *= 2
    end

    if tid == 1
        b[blockIdx().x] = tile[1]
    end
    return nothing
end

function main()
    res = 32;

    a = [2.0 for i = 1:res]
    d_a = CuArray(a)
    d_a_out = similar(d_a)

    @cuda threads = (res, 1, 1) blocks = (1, 1, 1) shmem = sizeof(Float64)*res vec_sum(d_a, d_a_out, res)

    a_out = Array(d_a_out)
    a_check = Array(d_a)

    @test a_out[1] ≈ sum(a)
end

main()


text:
##############
#auxiliary functions to run the testes
##############

function Grammatrix_levenstein(Xᵢ,Xⱼ)
    K = Array{Float64,2}(undef,4,4)
    for (i,si) in enumerate(Xᵢ)
        for (j,sj) in  enumerate(Xⱼ)
            K[i,j] = exp(-evaluate(Levenshtein(),si,sj))
        end
    end
    return K
end

function Grammatrix_cossine(Xᵢ,Xⱼ,all_mod)

     K= Array{Float64,2}(undef,4,4)
    for (i,si) in enumerate(Xᵢ)
        for (j,sj) in  enumerate(Xⱼ)
            v1 = _word2vec(si,all_mod)
            v2 = _word2vec(sj,all_mod)
            K[i,j] = 1-cosine_dist(v1,v2)
        end
    end
    return K
end

max_stars_repo_name: SamDier/BOMoD

max_issues_repo_name: SamDier/BOMoD

max_forks_repo_name: SamDier/BOMoD


text:
module NamedColors
using ColorTypes
using REPL: levsort
import Colors: @colorant_str

export load_xkcd, load_resene, load_nbs, load_x11_style, load_winsor_newton, load_x11, load_css3, load_crayola, ALL_COLORS, @colorant_str, UnknownColorError, named_color

include("loading.jl")

"""
Loads a combined collection of all colors.
Name clashes are resolved byt the load order.
Which is defined, as from largest collection, to smalled collection.
This is the initialised for the globel const ALL_COLORS, used by "@colorant"
"""
function load_all_colors(load_functions = (
    load_resene,
    load_xkcd,
    load_nbs,
	load_winsor_newton,
    load_x11,
    load_crayola,
    load_css3,
    ))

    merge!(Dict{String, RGB24}(), map(x->x(), load_functions)...)
end

const ALL_COLORS = load_all_colors()

include("namelookup.jl")

end # module


text:
struct QboRef{Obj}
    value::Int
end
Base.convert(::Type{QboRef{Obj}}, x::Obj) where {Obj} = QboRef{Obj}(x.Id)
Base.convert(::Type{Maybe{QboRef{Obj}}}, x::Obj) where {Obj} = QboRef{Obj}(x.Id)
Base.convert(::Type{Optional{QboRef{Obj}}}, x::Obj) where {Obj} = QboRef{Obj}(x.Id)
Base.convert(::Type{Optional{QboRef{<:Any}}}, x::Missing) = missing

from_json(::Type{QboRef{T}}, data) where {T} = QboRef{T}(parse(Int, data["value"]))
from_json(::Type{QboRef{<:QBObject}}, data) = from_json(QboRef{key_to_type(data["type"])}, data)

function to_json(::Type{QboRef{<:QBObject}}, obj::QBObject)
    Dict(
        "type"=>String(typeof(obj).name.name),
        "value"=>obj.value
    )
end
to_json(::Type, obj) = to_json(obj)

function to_json(data::QboRef)
    Dict("value"=>data.value)
end