text:
# -*- mode: julia; -*-
Base.__precompile__(true)
"""
The WAV package is a pure Julia library for reading and writing the
[WAV audio file format](https://en.wikipedia.org/wiki/WAV).

It provides [`wavread`](@ref), [`wavwrite`](@ref) and
[`wavappend`](@ref) functions to read, write, and append to WAV files.
The function [`wavplay`](@ref) provides simple audio playback.

These functions behave similar to the former MATLAB functions of the
same name.

This module also provides `wavread` and `wavwrite` as `load` and
`save` methods for `format"WAV"` to the `FileIO` package.

To read and write `CUE` and `INFO` chunks, there are experimental
functions [`wav_cue_read`](@ref), [`wav_cue_write`](@ref),
[`wav_info_read`](@ref), [`wav_info_write`](@ref).

# Example
```julia
using WAV
fs = 8e3
t = 0.0:1/fs:prevfloat(1.0)
f = 1e3
y = sin.(2pi * f * t) * 0.1
wavwrite(y, "example.wav", Fs=fs)

y, fs = wavread("example.wav")
y = sin.(2pi * 2f * t) * 0.1
wavappend(y, "example.wav")

y, fs = wavread("example.wav")
wavplay(y, fs)
```
"""
module WAV
export wavread, wavwrite, wavappend, wavplay
export WAVChunk, WAVMarker, wav_cue_read, wav_cue_write, wav_info_write, wav_info_read
export WAVArray, WAVFormatExtension, WAVFormat
export isextensible, isformat, bits_per_sample
export WAVE_FORMAT_PCM, WAVE_FORMAT_IEEE_FLOAT, WAVE_FORMAT_ALAW, WAVE_FORMAT_MULAW
import Libdl
using FileIO
using Logging

include("AudioDisplay.jl")
include("WAVChunk.jl")

# The WAV specification states that numbers are written to disk in little endian form.
write_le(stream::IO, value) = write(stream, htol(value))
read_le(stream::IO, x::Type{T}) where {T} = ltoh(read(stream, T))

# used by WAVE_FORMAT_EXTENSIBLE
struct WAVFormatExtension
    nbits::UInt16 # overrides nbits in WAVFormat type
    channel_mask::UInt32
    sub_format::Array{UInt8, 1} # 16 byte GUID
    WAVFormatExtension() = new(0, 0, UInt8[])
    WAVFormatExtension(nb, cm, sb) = new(nb, cm, sb)
end

# Required WAV Chunk; The format chunk describes how the waveform data is stored
struct WAVFormat
    compression_code::UInt16
    nchannels::UInt16
    sample_rate::UInt32
    bytes_per_second::UInt32 # average bytes per second
    block_align::UInt16
    nbits::UInt16
    ext::WAVFormatExtension
    WAVFormat() = new(WAVE_FORMAT_PCM, 0, 0, 0, 16, 0, WAVFormatExtension())
    WAVFormat(cc, nchan, fs, bps, ba, nb, e) = new(cc, nchan, fs, bps, ba, nb, e)
end

const WAVE_FORMAT_PCM        = 0x0001 # PCM
const WAVE_FORMAT_IEEE_FLOAT = 0x0003 # IEEE float
const WAVE_FORMAT_ALAW       = 0x0006 # A-Law
const WAVE_FORMAT_MULAW      = 0x0007 # Mu-Law
const WAVE_FORMAT_EXTENSIBLE = 0xfffe # Extension!

isextensible(fmt::WAVFormat) = (fmt.compression_code == WAVE_FORMAT_EXTENSIBLE)
bits_per_sample(fmt::WAVFormat) = isextensible(fmt) ? fmt.ext.nbits : fmt.nbits

# DEFINE_GUIDSTRUCT("00000001-0000-0010-8000-00aa00389b71", KSDATAFORMAT_SUBTYPE_PCM);
const KSDATAFORMAT_SUBTYPE_PCM = [
0x01, 0x00, 0x00, 0x00,
0x00, 0x00,
0x10, 0x00,
0x80, 0x00,
0x00, 0xaa, 0x00, 0x38, 0x9b, 0x71
                                  ]
# DEFINE_GUIDSTRUCT("00000003-0000-0010-8000-00aa00389b71", KSDATAFORMAT_SUBTYPE_IEEE_FLOAT);
const KSDATAFORMAT_SUBTYPE_IEEE_FLOAT = [
0x03, 0x00, 0x00, 0x00,
0x00, 0x00,
0x10, 0x00,
0x80, 0x00,
0x00, 0xaa, 0x00, 0x38, 0x9b, 0x71
                                         ]
# DEFINE_GUID(KSDATAFORMAT_SUBTYPE_MULAW, 0x00000007, 0x0000, 0x0010, 0x80, 0x00, 0x00, 0xaa, 0x00, 0x38, 0x9b, 0x71);
const KSDATAFORMAT_SUBTYPE_MULAW = [
0x07, 0x00, 0x00, 0x00,
0x00, 0x00,
0x10, 0x00,
0x80, 0x00,
0x00, 0xaa, 0x00, 0x38, 0x9b, 0x71
                                    ]
#DEFINE_GUID(KSDATAFORMAT_SUBTYPE_ALAW, 0x00000006, 0x0000, 0x0010, 0x80, 0x00, 0x00, 0xaa, 0x00, 0x38, 0x9b, 0x71);
const KSDATAFORMAT_SUBTYPE_ALAW = [
0x06, 0x00, 0x00, 0x00,
0x00, 0x00,
0x10, 0x00,
0x80, 0x00,
0x00, 0xaa, 0x00, 0x38, 0x9b, 0x71
                                   ]

"""
    getformat(chunks::Vector{WAVChunk})::WAVFormat

Takes the vector of `WAVChunk` elements returned in the `opt` result
of function `wavread` and extracts format information, returned as an
instance of type `WAVFormat`.

The `WAVFormat` type is defined as:

```julia
struct WAVFormat
    compression_code::UInt16
    nchannels::UInt16
    sample_rate::UInt32
    bytes_per_second::UInt32 # average bytes per second
    block_align::UInt16
    nbits::UInt16
    ext::WAVFormatExtension
end
```

The `ext` field of type `WAVFormatExtension` is defined as:

```julia
struct WAVFormatExtension
    nbits::UInt16 # overrides nbits in WAVFormat type
    channel_mask::UInt32
    sub_format::Array{UInt8, 1} # 16 byte GUID
    WAVFormatExtension() = new(0, 0, Array(UInt8, 0))
    WAVFormatExtension(nb, cm, sb) = new(nb, cm, sb)
end
```

You can use the `isformat` function to test how the samples are
encoded, without worrying about the `WAVFormatExtension` type.
Extended WAV files were added to deal with some ambiguity in the
original specification.
"""
function getformat(chunks::Vector{WAVChunk})::WAVFormat
    data = chunks[findfirst(c -> c.id == Symbol("fmt "), chunks)].data
    buf = IOBuffer(data[5:end])
    return read_format(buf, convert(UInt32, length(data[5:end])))
end

"""
    isformat(fmt::WAVFormat, code)

The `isformat` function takes the `WAVFormat` object returned by
`WAV.getformat` and `code` can be one of the following constants:

* `WAVE_FORMAT_PCM`
* `WAVE_FORMAT_IEEE_FLOAT`
* `WAVE_FORMAT_ALAW`
* `WAVE_FORMAT_MULAW`
* `WAVE_FORMAT_EXTENSIBLE`

The function returns `true` when the samples are encoded in the
specified `code`.
"""
function isformat(fmt::WAVFormat, code)
    if code != WAVE_FORMAT_EXTENSIBLE && isextensible(fmt)
        subtype = UInt8[]
        if code == WAVE_FORMAT_PCM
            subtype = KSDATAFORMAT_SUBTYPE_PCM
        elseif code == WAVE_FORMAT_IEEE_FLOAT
            subtype = KSDATAFORMAT_SUBTYPE_IEEE_FLOAT
        elseif code == WAVE_FORMAT_ALAW
            subtype = KSDATAFORMAT_SUBTYPE_ALAW
        elseif code == WAVE_FORMAT_MULAW
            subtype = KSDATAFORMAT_SUBTYPE_MULAW
        else
            return false
        end
        return subtype == fmt.ext.sub_format
    end
    return fmt.compression_code == code
end

function WAVFormatExtension(bytes)
    if isempty(bytes)
        return WAVFormatExtension()
    end
    # split bytes into valid_bits_per_sample, channel_mask, and sub_format
    valid_bits_per_sample = (convert(UInt16, bytes[2]) << 8) | convert(UInt16, bytes[1])
    channel_mask = (convert(UInt32, bytes[6]) << 24) | (convert(UInt32, bytes[5]) << 16) | (convert(UInt32, bytes[4]) << 8) | convert(UInt32, bytes[3])
    sub_format = bytes[7:end]
    return WAVFormatExtension(valid_bits_per_sample, channel_mask, sub_format)
end

function read_header(io::IO)
    # check if the given file has a valid RIFF header
    riff = Vector{UInt8}(undef, 4)
    read!(io, riff)
    if riff !=  b"RIFF"
        error("Invalid WAV file: The RIFF header is invalid")
    end

    chunk_size = read_le(io, UInt32)

    # check if this is a WAV file
    format = Vector{UInt8}(undef, 4)
    read!(io, format)
    if format != b"WAVE"
        error("Invalid WAV file: the format is not WAVE")
    end
    return chunk_size
end

function write_header(io::IO, data_length::UInt32)
    write(io, b"RIFF") # RIFF header
    write_le(io, data_length) # chunk_size
    write(io, b"WAVE")
end
write_standard_header(io, data_length) = write_header(io, UInt32(data_length + 36))
write_extended_header(io, data_length) = write_header(io, UInt32(data_length + 60))

function read_format(io::IO, chunk_size::UInt32)
    # can I read in all of the fields at once?
    orig_chunk_size = convert(Int, chunk_size)
    if chunk_size < 16
        error("The WAVE Format chunk must be at least 16 bytes")
    end
    compression_code = read_le(io, UInt16)
    nchannels = read_le(io, UInt16)
    sample_rate = read_le(io, UInt32)
    bytes_per_second = read_le(io, UInt32)
    block_align = read_le(io, UInt16)
    nbits = read_le(io, UInt16)
    ext = UInt8[]
    chunk_size -= 16
    if chunk_size > 0
        extra_bytes_length = read_le(io, UInt16)
        if extra_bytes_length == 22
            ext = Vector{UInt8}(undef, extra_bytes_length)
            read!(io, ext)
        end
    end
    return WAVFormat(compression_code,
                     nchannels,
                     sample_rate,
                     bytes_per_second,
                     block_align,
                     nbits,
                     WAVFormatExtension(ext))
end

function format_length(fmt::WAVFormat)
    len = 16 # 16 is size of base format chunk
    if isextensible(fmt)
        len += 24 # 24 is the added length needed to encode the extension
    end
    return len
end

function write_format(io::IO, fmt::WAVFormat)
    len = format_length(fmt)
    # write the fmt subchunk header
    write(io, b"fmt ")
    write_le(io, convert(UInt32, len)) # subchunk length

    write_le(io, fmt.compression_code) # audio format (UInt16)
    write_le(io, fmt.nchannels) # number of channels (UInt16)
    write_le(io, fmt.sample_rate) # sample rate (UInt32)
    write_le(io, fmt.bytes_per_second) # byte rate (UInt32)
    write_le(io, fmt.block_align) # byte align (UInt16)
    write_le(io, fmt.nbits) # number of bits per sample (UInt16)

    if isextensible(fmt)
        write_le(io, convert(UInt16, 22))
        write_le(io, fmt.ext.nbits)
        write_le(io, fmt.ext.channel_mask)
        @assert length(fmt.ext.sub_format) == 16
        write(io, fmt.ext.sub_format)
    end
end

function WAVChunk(fmt::WAVFormat)
    io = IOBuffer()
    write_format(io, fmt)
    data = take!(io)[5:end]
    WAVChunk(Symbol("fmt "), data)
end

function pcm_container_type(nbits::Unsigned)
    if nbits > 32
        return Int64
    elseif nbits > 16
        return Int32
    elseif nbits > 8
        return Int16
    end
    return  UInt8
end

ieee_float_container_type(nbits) = (nbits == 32 ? Float32 : (nbits == 64 ? Float64 : error("$nbits bits is not supported for WAVE_FORMAT_IEEE_FLOAT.")))

function read_pcm_samples(io::IO, chunk_size, fmt::WAVFormat, subrange,
                          ::Type{sample_type}) where {sample_type}
    nbits = bits_per_sample(fmt)
    if isempty(subrange)
        return Array{sample_type, 2}(undef, 0, fmt.nchannels)
    end
    samples = Array{sample_type, 2}(undef, length(subrange), fmt.nchannels)
    nbytes = ceil(Integer, nbits / 8)
    bitshift = [0x0, 0x8, 0x10, 0x18, 0x20, 0x28, 0x30, 0x38, 0x40]
    mask = UInt64(0x1) << (nbits - 1)
    if nbits <= 8
        mask = UInt64(0)
    end
    skip(io, convert(UInt, (first(subrange) - 1) * nbytes * fmt.nchannels))
    raw_sample = Vector{UInt8}(undef, nbytes*length(samples))
    read!(io, raw_sample)
    skip(io, chunk_size - convert(UInt, last(subrange) * nbytes * fmt.nchannels))
    raw_sample = reshape(raw_sample, nbytes, size(samples, 2), size(samples, 1))
    for i = 1:size(samples, 1)
        for j = 1:size(samples, 2)
            my_sample = UInt64(0)
            for k = 1:nbytes
                my_sample |= convert(UInt64, raw_sample[k,j,i]) << bitshift[k]
            end
            my_sample >>= nbytes * 8 - nbits
            # sign extend negative values
            my_sample = xor(my_sample, mask) - mask
            samples[i, j] = convert(sample_type, signed(my_sample))
        end
    end
    samples
end

function read_ieee_float_samples(io::IO, chunk_size, fmt::WAVFormat, subrange, ::Type{floatType}) where {floatType}
    if isempty(subrange)
        return Array{floatType, 2}(undef, 0, fmt.nchannels)
    end
    nblocks = length(subrange)
    samples = Array{floatType, 1}(undef, fmt.nchannels* nblocks)
    nbits = bits_per_sample(fmt)
    skip(io, convert(UInt, (first(subrange) - 1) * (nbits / 8) * fmt.nchannels))
    samples = Vector{floatType}(undef, fmt.nchannels*nblocks)
    read!(io, samples) # read_le(stream::IO, x::Type{T}) where {T} = ltoh(fastread(stream, T))
    skip(io, chunk_size - convert(UInt, last(subrange) * (nbits / 8) * fmt.nchannels))
    copy(reshape(samples, Int(fmt.nchannels), Int(nblocks))')
end

# take the loop variable type out of the loop
function read_ieee_float_samples(io::IO, chunk_size, fmt::WAVFormat, subrange)
    floatType = ieee_float_container_type(bits_per_sample(fmt))
    read_ieee_float_samples(io, chunk_size, fmt, subrange, floatType)
end

function read_companded_samples(io::IO, chunk_size, fmt::WAVFormat, subrange, table)
    if isempty(subrange)
        return Array{eltype(table), 2}(undef, 0, fmt.nchannels)
    end
    nblocks = length(subrange)
    samples = Array{eltype(table), 2}(undef, nblocks, fmt.nchannels)
    skip(io, convert(UInt, (first(subrange) - 1) * fmt.nchannels))
    for i = 1:nblocks
        for j = 1:fmt.nchannels
            # add one to value from blocks because A-law stores values from 0 to 255.
            compressedByte::UInt8 = read(io, UInt8)
            # Julia indexing is 1-based; I need a value from 1 to 256
            samples[i, j] = table[compressedByte + 1]
        end
    end
    skip(io, convert(UInt, chunk_size - last(subrange) * fmt.nchannels))
    return samples
end

function read_mulaw_samples(io::IO, chunk_size, fmt::WAVFormat, subrange)
    # Quantized μ-law algorithm -- Use a look up table to convert
    # From Wikipedia, ITU-T Recommendation G.711 and G.191 specify the following intervals:
    #
    # ---------------------------------------+--------------------------------
    #  14 bit Binary Linear input code       | 8 bit Compressed code
    # ---------------------------------------+--------------------------------
    # +8158 to +4063 in 16 intervals of 256  |  0x80 + interval number
    # +4062 to +2015 in 16 intervals of 128  |  0x90 + interval number
    # +2014 to +991 in 16 intervals of 64    |  0xA0 + interval number
    # +990 to +479 in 16 intervals of 32     |  0xB0 + interval number
    # +478 to +223 in 16 intervals of 16     |  0xC0 + interval number
    # +222 to +95 in 16 intervals of 8       |  0xD0 + interval number
    # +94 to +31 in 16 intervals of 4        |  0xE0 + interval number
    # +30 to +1 in 15 intervals of 2         |  0xF0 + interval number
    # 0                                      |  0xFF
    # −1                                     |  0x7F
    # −31 to −2 in 15 intervals of 2         |  0x70 + interval number
    # −95 to −32 in 16 intervals of 4        |  0x60 + interval number
    # −223 to −96 in 16 intervals of 8       |  0x50 + interval number
    # −479 to −224 in 16 intervals of 16     |  0x40 + interval number
    # −991 to −480 in 16 intervals of 32     |  0x30 + interval number
    # −2015 to −992 in 16 intervals of 64    |  0x20 + interval number
    # −4063 to −2016 in 16 intervals of 128  |  0x10 + interval number
    # −8159 to −4064 in 16 intervals of 256  |  0x00 + interval number
    # ---------------------------------------+--------------------------------
    MuLawDecompressTable =
    [
    -32124,-31100,-30076,-29052,-28028,-27004,-25980,-24956,
    -23932,-22908,-21884,-20860,-19836,-18812,-17788,-16764,
    -15996,-15484,-14972,-14460,-13948,-13436,-12924,-12412,
    -11900,-11388,-10876,-10364, -9852, -9340, -8828, -8316,
    -7932, -7676, -7420, -7164, -6908, -6652, -6396, -6140,
    -5884, -5628, -5372, -5116, -4860, -4604, -4348, -4092,
    -3900, -3772, -3644, -3516, -3388, -3260, -3132, -3004,
    -2876, -2748, -2620, -2492, -2364, -2236, -2108, -1980,
    -1884, -1820, -1756, -1692, -1628, -1564, -1500, -1436,
    -1372, -1308, -1244, -1180, -1116, -1052,  -988,  -924,
    -876,  -844,  -812,  -780,  -748,  -716,  -684,  -652,
    -620,  -588,  -556,  -524,  -492,  -460,  -428,  -396,
    -372,  -356,  -340,  -324,  -308,  -292,  -276,  -260,
    -244,  -228,  -212,  -196,  -180,  -164,  -148,  -132,
    -120,  -112,  -104,   -96,   -88,   -80,   -72,   -64,
    -56,   -48,   -40,   -32,   -24,   -16,    -8,     -1,
    32124, 31100, 30076, 29052, 28028, 27004, 25980, 24956,
    23932, 22908, 21884, 20860, 19836, 18812, 17788, 16764,
    15996, 15484, 14972, 14460, 13948, 13436, 12924, 12412,
    11900, 11388, 10876, 10364,  9852,  9340,  8828,  8316,
    7932,  7676,  7420,  7164,  6908,  6652,  6396,  6140,
    5884,  5628,  5372,  5116,  4860,  4604,  4348,  4092,
    3900,  3772,  3644,  3516,  3388,  3260,  3132,  3004,
    2876,  2748,  2620,  2492,  2364,  2236,  2108,  1980,
    1884,  1820,  1756,  1692,  1628,  1564,  1500,  1436,
    1372,  1308,  1244,  1180,  1116,  1052,   988,   924,
    876,   844,   812,   780,   748,   716,   684,   652,
    620,   588,   556,   524,   492,   460,   428,   396,
    372,   356,   340,   324,   308,   292,   276,   260,
    244,   228,   212,   196,   180,   164,   148,   132,
    120,   112,   104,    96,    88,    80,    72,    64,
    56,    48,    40,    32,    24,    16,     8,     0
     ]
    @assert length(MuLawDecompressTable) == 256
    return read_companded_samples(io, chunk_size, fmt, subrange, MuLawDecompressTable)
end

function read_alaw_samples(io::IO, chunk_size, fmt::WAVFormat, subrange)
    # Quantized A-law algorithm -- Use a look up table to convert
    ALawDecompressTable =
    [
    -5504, -5248, -6016, -5760, -4480, -4224, -4992, -4736,
    -7552, -7296, -8064, -7808, -6528, -6272, -7040, -6784,
    -2752, -2624, -3008, -2880, -2240, -2112, -2496, -2368,
    -3776, -3648, -4032, -3904, -3264, -3136, -3520, -3392,
    -22016,-20992,-24064,-23040,-17920,-16896,-19968,-18944,
    -30208,-29184,-32256,-31232,-26112,-25088,-28160,-27136,
    -11008,-10496,-12032,-11520,-8960, -8448, -9984, -9472,
    -15104,-14592,-16128,-15616,-13056,-12544,-14080,-13568,
    -344,  -328,  -376,  -360,  -280,  -264,  -312,  -296,
    -472,  -456,  -504,  -488,  -408,  -392,  -440,  -424,
    -88,   -72,   -120,  -104,  -24,   -8,    -56,   -40,
    -216,  -200,  -248,  -232,  -152,  -136,  -184,  -168,
    -1376, -1312, -1504, -1440, -1120, -1056, -1248, -1184,
    -1888, -1824, -2016, -1952, -1632, -1568, -1760, -1696,
    -688,  -656,  -752,  -720,  -560,  -528,  -624,  -592,
    -944,  -912,  -1008, -976,  -816,  -784,  -880,  -848,
    5504,  5248,  6016,  5760,  4480,  4224,  4992,  4736,
    7552,  7296,  8064,  7808,  6528,  6272,  7040,  6784,
    2752,  2624,  3008,  2880,  2240,  2112,  2496,  2368,
    3776,  3648,  4032,  3904,  3264,  3136,  3520,  3392,
    22016, 20992, 24064, 23040, 17920, 16896, 19968, 18944,
    30208, 29184, 32256, 31232, 26112, 25088, 28160, 27136,
    11008, 10496, 12032, 11520, 8960,  8448,  9984,  9472,
    15104, 14592, 16128, 15616, 13056, 12544, 14080, 13568,
    344,   328,   376,   360,   280,   264,   312,   296,
    472,   456,   504,   488,   408,   392,   440,   424,
    88,    72,   120,   104,    24,     8,    56,    40,
    216,   200,   248,   232,   152,   136,   184,   168,
    1376,  1312,  1504,  1440,  1120,  1056,  1248,  1184,
    1888,  1824,  2016,  1952,  1632,  1568,  1760,  1696,
    688,   656,   752,   720,   560,   528,   624,   592,
    944,   912,  1008,   976,   816,   784,   880,   848
     ]
    @assert length(ALawDecompressTable) == 256
    return read_companded_samples(io, chunk_size, fmt, subrange, ALawDecompressTable)
end

function compress_sample_mulaw(sample)
    MuLawCompressTable =
    [
    0,0,1,1,2,2,2,2,3,3,3,3,3,3,3,3,
    4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,
    5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,
    5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,
    6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,
    6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,
    6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,
    6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,
    7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,
    7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,
    7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,
    7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,
    7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,
    7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,
    7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,
    7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7
     ]
    @assert length(MuLawCompressTable) == 256
    cBias = 0x84
    cClip = 32635

    sampleSign = (sample >>> 8) & 0x80
    if sampleSign != 0
        sample = -sample
    end
    if sample > cClip
        sample = cClip
    end
    sample = sample + cBias
    sampleExponent = MuLawCompressTable[(sample >>> 7) + 1]
    mantissa = (sample >> (sampleExponent+3)) & 0x0F
    UInt8((~ (sampleSign | (sampleExponent << 4) | mantissa)) & 0xff)
end

function compress_sample_alaw(sample)
    ALawCompressTable =
    [
    1,1,2,2,3,3,3,3,
    4,4,4,4,4,4,4,4,
    5,5,5,5,5,5,5,5,
    5,5,5,5,5,5,5,5,
    6,6,6,6,6,6,6,6,
    6,6,6,6,6,6,6,6,
    6,6,6,6,6,6,6,6,
    6,6,6,6,6,6,6,6,
    7,7,7,7,7,7,7,7,
    7,7,7,7,7,7,7,7,
    7,7,7,7,7,7,7,7,
    7,7,7,7,7,7,7,7,
    7,7,7,7,7,7,7,7,
    7,7,7,7,7,7,7,7,
    7,7,7,7,7,7,7,7,
    7,7,7,7,7,7,7,7
     ]
    @assert length(ALawCompressTable) == 128
    cBias = 0x84
    cClip = 32635
    sampleSign = ((~sample >>> 8) & 0x80)
    if sampleSign == 0
        sample = -sample
    end
    if sample > cClip
        sample = cClip
    end
    compressedByte = 0
    if sample >= 256
        sampleExponent = ALawCompressTable[((sample >>> 8) & 0x7f) + 1]
        mantissa = (sample >>> (sampleExponent + 3) ) & 0x0f
        compressedByte = ((sampleExponent << 4) | mantissa) & 0xff
    else
        compressedByte = (sample >>> 4) & 0xff
    end
    compressedByte = xor(xor(sampleSign, 0x55), compressedByte)
    UInt8(compressedByte & 0xff)
end

function write_companded_samples(io::IO, samples::AbstractArray{T}, compander::Function) where T <: Integer
    for i = 1:size(samples, 1)
        for j = 1:size(samples, 2)
            write_le(io, compander(samples[i, j]))
        end
    end
end

function write_companded_samples(io::IO, samples::AbstractArray{T}, compander::Function) where T <: AbstractFloat
    samples = convert(Array{Int16}, round.(samples * typemax(Int16)))
    write_companded_samples(io, samples, compander)
end

# PCM data is two's-complement except for resolutions of 1-8 bits, which are represented as offset binary.

# support every bit width from 1 to 8 bits
convert_pcm_to_double(samples::AbstractArray{UInt8}, nbits::Integer) = convert(Array{Float64}, samples) ./ (2.0^nbits - 1) .* 2.0 .- 1.0
convert_pcm_to_double(::AbstractArray{Int8}, ::Integer) = error("WAV files use offset binary for less than 9 bits")
# support every bit width from 9 to 64 bits
convert_pcm_to_double(samples::AbstractArray{T}, nbits::Integer) where T <: Signed = convert(Array{Float64}, samples) / (2.0^(nbits - 1) - 1)

function read_data(io::IO, chunk_size, fmt::WAVFormat, format, subrange)
    # "format" is the format of values, while "fmt" is the WAV file level format
    convert_to_double = x -> convert(Array{Float64}, x)

    if subrange === Nothing
        Base.depwarn("`wavread(..., subrange=Nothing)` is deprecated, use `wavread(..., subrange=:)` instead.", :read_data)
        subrange = (:)
    end
    if subrange === (:)
        # each block stores fmt.nchannels channels
        subrange = 1:convert(UInt, chunk_size / fmt.block_align)
    end
    if isformat(fmt, WAVE_FORMAT_PCM)
        samples = read_pcm_samples(io, chunk_size, fmt, subrange,
                                   pcm_container_type(bits_per_sample(fmt)))
        convert_to_double = x -> convert_pcm_to_double(x, bits_per_sample(fmt))
    elseif isformat(fmt, WAVE_FORMAT_IEEE_FLOAT)
        samples = read_ieee_float_samples(io, chunk_size, fmt, subrange)
    elseif isformat(fmt, WAVE_FORMAT_MULAW)
        samples = read_mulaw_samples(io, chunk_size, fmt, subrange)
        convert_to_double = x -> convert_pcm_to_double(x, 16)
    elseif isformat(fmt, WAVE_FORMAT_ALAW)
        samples = read_alaw_samples(io, chunk_size, fmt, subrange)
        convert_to_double = x -> convert_pcm_to_double(x, 16)
    else
        error("$(fmt.compression_code) is an unsupported compression code!")
    end
    if format == "double"
        samples = convert_to_double(samples)
    end
    samples
end

function write_pcm_samples(io::IO, fmt::WAVFormat, samples::AbstractArray{T}) where T <: Integer
    nbits = bits_per_sample(fmt)
    # number of bytes per sample
    nbytes = ceil(Integer, nbits / 8)
    for i = 1:size(samples, 1)
        for j = 1:size(samples, 2)
            my_sample = samples[i, j]
            # shift my_sample into the N most significant bits
            my_sample <<= nbytes * 8 - nbits
            for k = 1:nbytes
                write_le(io, convert(UInt8, my_sample & 0xff))
                my_sample = my_sample >> 8
            end
        end
    end
end

function write_pcm_samples(io::IO, fmt::WAVFormat, samples::AbstractArray{T}) where T <: AbstractFloat
    nbits = bits_per_sample(fmt)
    # Scale the floating point values to the PCM range
    if nbits > 8
        # two's complement
        samples = convert(Array{pcm_container_type(nbits)}, round.(samples * (2.0^(nbits - 1) - 1)))
    else
        # offset binary
        samples = convert(Array{UInt8}, round.((samples .+ 1.0) / 2.0 * (2.0^nbits - 1)))
    end
    return write_pcm_samples(io, fmt, samples)
end

function write_ieee_float_samples(io::IO, samples)
    # Interleave the channel samples before writing to the stream.
    for i = 1:size(samples, 1) # for each sample
        for j = 1:size(samples, 2) # for each channel
            write_le(io, samples[i, j])
        end
    end
end

# take the loop variable type out of the loop
function write_ieee_float_samples(io::IO, fmt::WAVFormat, samples)
    floatType = ieee_float_container_type(bits_per_sample(fmt))
    write_ieee_float_samples(io, convert(Array{floatType}, samples))
end

function write_data(io::IO, fmt::WAVFormat, samples::AbstractArray)
    if isformat(fmt, WAVE_FORMAT_PCM)
        return write_pcm_samples(io, fmt, samples)
    elseif isformat(fmt, WAVE_FORMAT_IEEE_FLOAT)
        return write_ieee_float_samples(io, fmt, samples)
    elseif isformat(fmt, WAVE_FORMAT_MULAW)
        return write_companded_samples(io, samples, compress_sample_mulaw)
    elseif isformat(fmt, WAVE_FORMAT_ALAW)
        return write_companded_samples(io, samples, compress_sample_alaw)
    else
        error("$(fmt.compression_code) is an unsupported compression code.")
    end
end

make_range(subrange) = subrange
make_range(subrange::Number) = 1:convert(Int, subrange)

"""
    wavread(io::IO; subrange=:, format="double")
    wavread(filename::AbstractString; subrange=:, format="double")

Reads the samples from a WAV file. The samples are converted to
floating point values in the range −1.0 to 1.0 by default.

The available options, and the default values, are:

* `format` selects the form of data returned:
  * `format="double"` (default) returns double-precision floating point
    (`Float64`) values in the range −1.0 to 1.0.
  * `format="native"` returns the values as encoded in the file.
  * `format="size"` returns a 2-tuple (`n`, `m`) containing the number of
    samples 'n' and the number of channels 'm' in the file (like `size(y)`),
    rather than the regular 4-tuple `(y, Fs, nbits, opt)` with the
    actual samples in `y`.

* `subrange` controls which samples are returned.
  The default (`:`) returns all samples in the file.
  Passing an integer `N` (or equivalently the range `1:N`) returns
  the first `N` samples of each channel.
  Passing a unitrange `I:J` returns `length(I:J)` consecutive
  samples from each channel, starting with the `I`-th sample.

The function returns a 4-tuple with elements

* `y`: A 2-dimensional array containing the waveform samples,
  where the row index represents the time axis
  and the colum index the channel number
* `Fs`: The sampling frequency in hertz
* `nbits`: the number of bits used to encode each sample
* `opt`: A vector of [`WAVChunk`](@ref) elements representing
  optional chunksfound in the WAV file.

The elements in the `opt` vector depend on the contents
of the WAV file. A `WAVChunk` is defined as

```julia
struct WAVChunk
    id::Symbol
    data::Vector{UInt8}
end
```

where `id` is the four-character chunk ID. All valid WAV files will
contain a `fmt` chunk, with `id==Symbol("fmt ")` (note the trailing
space).

In order to obtain the contents of the format chunk, call
`WAV.getformat(opt)`.

The following methods are also defined to make this function
compatible with MATLAB’s former `wavread` function:

    wavread(filename::AbstractString, fmt::AbstractString) = wavread(filename, format=fmt)
    wavread(filename::AbstractString, n) = wavread(filename, subrange=n)
    wavread(filename::AbstractString, n, fmt) = wavread(filename, subrange=n, format=fmt)

# Example

    y, fs, nbits, opt = wavread("example.wav")

See also: [`wavwrite`](@ref)
"""
function wavread(io::IO; subrange=(:), format="double")
    chunk_size = read_header(io)
    samples = Array{Float64, 1}()
    nbits = 0
    sample_rate = Float32(0.0)
    opt = WAVChunk[]

    # Subtract the size of the format field from chunk_size; now it holds the size
    # of all the sub-chunks
    chunk_size -= 4
    # GitHub Issue #18: Check if there is enough data to read another chunk
    subchunk_header_size = 4 + sizeof(UInt32)
    fmt = WAVFormat()
    data_position = 0
    data_size = 0
    while chunk_size >= subchunk_header_size
        # Read subchunk ID and size
        subchunk_id = Vector{UInt8}(undef, 4)
        read!(io, subchunk_id)
        subchunk_size = read_le(io, UInt32)
        nextchunk_start = position(io) + subchunk_size
        if subchunk_size > chunk_size
            chunk_size = 0
            break
        end
        chunk_size -= subchunk_header_size + subchunk_size
        # check the subchunk ID
        if subchunk_id == b"fmt "
            fmt = read_format(io, subchunk_size)
            sample_rate = Float32(fmt.sample_rate)
            nbits = bits_per_sample(fmt)
            push!(opt, WAVChunk(fmt))
        elseif subchunk_id == b"data"
            data_position = position(io)
            data_size = subchunk_size
        else
            subchunk_data = Vector{UInt8}(undef, subchunk_size)
            read!(io, subchunk_data)
            push!(opt, WAVChunk(Symbol(subchunk_id), subchunk_data))
        end
        seek(io, nextchunk_start)
    end
    if data_position != 0
        seek(io, data_position)
        if format == "size"
            return convert(Int, data_size / fmt.block_align), convert(Int, fmt.nchannels)
        end
        samples = read_data(io, data_size, fmt, format, make_range(subrange))
    end
    return samples, sample_rate, nbits, opt
end

function wavread(filename::AbstractString; subrange=(:), format="double")
    open(filename, "r") do io
        wavread(io, subrange=subrange, format=format)
    end
end

# These are the MATLAB compatible signatures
wavread(filename::AbstractString, fmt::AbstractString) = wavread(filename, format=fmt)
wavread(filename::AbstractString, n) = wavread(filename, subrange=n)
wavread(filename::AbstractString, n, fmt) = wavread(filename, subrange=n, format=fmt)

get_default_compression(::AbstractArray{T}) where T <: Integer = WAVE_FORMAT_PCM
get_default_compression(::AbstractArray{T}) where T<:AbstractFloat = WAVE_FORMAT_IEEE_FLOAT
get_default_pcm_precision(::AbstractArray{UInt8}) = 8
get_default_pcm_precision(::AbstractArray{Int16}) = 16
get_default_pcm_precision(::Any) = 24

function get_default_precision(samples, compression)
    if compression == WAVE_FORMAT_ALAW || compression == WAVE_FORMAT_MULAW
        return 8
    elseif compression == WAVE_FORMAT_IEEE_FLOAT
        return 32
    end
    get_default_pcm_precision(samples)
end

"""
    wavwrite(y::AbstractArray, io::IO;
             Fs=8000, nbits=0, compression=0, chunks::Vector{WAVChunk}=WAVChunk[])
    wavwrite(y::AbstractArray, filename::String;
             Fs=8000, nbits=0, compression=0, chunks::Vector{WAVChunk}=WAVChunk[])

Writes sample matrix `y` in RIFF/WAVE format to a file. Each column of
the data represents a different channel. Stereo files contain two
columns (left and right).

The second argument accepts either an `IO` object or a filename
(`String`).

The function choses by default a sample rate of 8 kHz and an output
data type and bits-per-sample number based on `eltype(y)`. These
defaults can be changed using optional keyword arguments:

   * `Fs`: sampling frequency in hertz
   * `nbits`: specify the number of bits to be used to encode each
     sample; the default (0) is an automatic choice based on
     the values of `compression` and `eltype(samples)`
   * `compression` controls the type of encoding used in the file,
     and can be one of
     * `WAVE_FORMAT_PCM`: `UInt8`, `Int16`, etc. integer encoding
     * `WAVE_FORMAT_IEEE_FLOAT`: `Float32` or `Float64` encoding
     * `WAVE_FORMAT_ALAW`: ITU-T G.711 A-law encoding (8-bit log-scale,
       used in European telephone networks)
     * `WAVE_FORMAT_MULAW`: ITU-T G.711 µ-law encoding (8-bit log-scale,
       used in American telephone networks)
     The default (`0`) is to pick an encoding automatically based on
     `typeof(samples)` (see below).
   * `chunks` (default = `WAVChunk[]`): a vector of `WAVChunk` objects to be written to the file (in addition to the format chunk). See below for some utilities for creating `CUE` and `INFO` chunks.

Unless otherwise specified via `nbits` and `compression`, the type of
the input array `y` determines the data format used in the generated
file. The function attempts to picks among the encodings commonly
supported by other WAV-reading audio software the one that best
preserves the provided input type.

For `Integer` input arrays, and `compression=WAVE_FORMAT_PCM` or
`compression=0`, the permitted sample-value ranges are:

|`nbits`| `eltype(y)` | supported range     | output type |
|-------|-------------|---------------------|-------------|
| 0     | UInt8,Int16 |    full range       | UInt8,Int16 |
| 0     | Int32       |  –2^23 ≤ y < 2^23   | Int32       |
| 8     | <: Integer  |      0 ≤ y ≤ 255    | UInt8       |
| 16    | <: Integer  | –32768 ≤ y ≤ +32767 | Int16       |
| 24    | <: Integer  |  –2^23 ≤ y < 2^23   | Int32       |
| 32    | <: Integer  |  –2^31 ≤ y < 2^31   | Int32       |

If `y` is a floating-point array, and
`compression=WAVE_FORMAT_IEEE_FLOAT` or `compression=0`, the full
range of Float32 values are supported:

|`nbits`| `eltype(y)`     | supported range     | output type |
|-------|-----------------|---------------------|-------------|
| 0     | Float32,Float64 | –Inf32 ≤ y ≤ +Inf32 | Float32     |
| 32    | Float32,Float64 | –Inf32 ≤ y ≤ +Inf32 | Float32     |
| 64    | Float32,Float64 | –Inf64 ≤ y ≤ +Inf64 | Float64     |

If `y` is a floating-point array, and `compression=WAVE_FORMAT_PCM`,
the input data ranges are:

|`nbits`| `eltype(y)`      | supported range | output type |
|-------|------------------|-----------------|-------------|
| 8     | Float32, Float64 | –1.0 ≤ y ≤ +1.0 | UInt8       |
| 16    | Float32, Float64 | –1.0 ≤ y ≤ +1.0 | Int16       |
| 24    | Float32, Float64 | –1.0 ≤ y ≤ +1.0 | Int32       |
| 32    | Float32, Float64 | –1.0 ≤ y ≤ +1.0 | Int32       |

The output format column shows the element type returned by
`wavread(..., format=-"native")`.

The following methods are also defined to make this function
compatible with MATLAB’s former `wavwrite` function:

```julia
wavwrite(y::Array, f::Real, filename::String) = wavwrite(y, filename, Fs=f)
wavwrite(y::Array, f::Real, N::Real, filename::String) = wavwrite(y, filename, Fs=f, nbits=N)
wavwrite(y::Array{T}, io::IO) where {T<:Integer} = wavwrite(y, io, nbits=sizeof(T)*8)
wavwrite(y::Array{T}, filename::String) where {T<:Integer} = wavwrite(y, filename, nbits=sizeof(T)*8)
wavwrite(y::Array{Int32}, io::IO) = wavwrite(y, io, nbits=24)
wavwrite(y::Array{Int32}, filename::String) = wavwrite(y, filename, nbits=24)
wavwrite(y::Array{T}, io::IO) where {T<:FloatingPoint} = wavwrite(y, io, nbits=sizeof(T)*8, compression=WAVE_FORMAT_IEEE_FLOAT)
wavwrite(y::Array{T}, filename::String) where {T<:FloatingPoint} = wavwrite(y, filename, nbits=sizeof(T)*8, compression=WAVE_FORMAT_IEEE_FLOAT)
```

See also: [`wavread`](@ref), [`wavappend`](@ref)
"""
function wavwrite(samples::AbstractArray, io::IO; Fs=8000, nbits=0, compression=0,
                  chunks::Vector{WAVChunk}=WAVChunk[])
    if compression == 0
        compression = get_default_compression(samples)
    elseif compression == WAVE_FORMAT_ALAW || compression == WAVE_FORMAT_MULAW
        nbits = 8
    end
    if nbits == 0
        nbits = get_default_precision(samples, compression)
    end
    compression_code = compression
    nchannels = size(samples, 2)
    sample_rate = Fs
    my_nbits = ceil(Integer, nbits / 8) * 8
    block_align = my_nbits / 8 * nchannels
    bps = sample_rate * block_align
    data_length::UInt32 = size(samples, 1) * block_align
    extra_length::UInt32 = 8 * length(chunks) + sum([length(c.data) for c in chunks])
    ext = WAVFormatExtension()

    if nchannels > 2 || my_nbits > 16 || my_nbits != nbits
        compression_code = WAVE_FORMAT_EXTENSIBLE
        valid_bits_per_sample = nbits
        channel_mask = 0
        sub_format = UInt8[]
        if compression == WAVE_FORMAT_PCM
            sub_format = KSDATAFORMAT_SUBTYPE_PCM
        elseif compression == WAVE_FORMAT_IEEE_FLOAT
            sub_format = KSDATAFORMAT_SUBTYPE_IEEE_FLOAT
        elseif compression == WAVE_FORMAT_ALAW
            sub_format = KSDATAFORMAT_SUBTYPE_ALAW
        elseif compression == WAVE_FORMAT_MULAW
            sub_format = KSDATAFORMAT_SUBTYPE_MULAW
        else
            error("Unsupported extension sub format: $compression")
        end
        ext = WAVFormatExtension(valid_bits_per_sample, channel_mask, sub_format)
        write_extended_header(io, data_length + extra_length)
    else
        write_standard_header(io, data_length + extra_length)
    end
    fmt = WAVFormat(compression_code,
                    nchannels,
                    sample_rate,
                    bps,
                    block_align,
                    my_nbits,
                    ext)
    write_format(io, fmt)

    for eachchunk in chunks
        write(io, eachchunk.id)
        write_le(io, UInt32(length(eachchunk.data)))
        for eachbyte in eachchunk.data
            write(io, eachbyte)
        end
    end

    # write the data subchunk header
    write(io, b"data")
    write_le(io, data_length) # UInt32
    write_data(io, fmt, samples)
end

function wavwrite(samples::AbstractArray, filename::AbstractString; Fs=8000, nbits=0, compression=0,
                  chunks::Vector{WAVChunk}=WAVChunk[])
    open(filename, "w") do io
        wavwrite(samples, io, Fs=Fs, nbits=nbits, compression=compression, chunks=chunks)
    end
end

"""
    wavappend(samples::AbstractArray, io::IO)
    wavappend(samples::AbstractArray, filename::AbstractString)

Append samples to an existing WAV file. All parameters (data type and
range, output format, number of bits, number of channels, etc.) are
assumed to match.
"""
function wavappend(samples::AbstractArray, io::IO)
    seekstart(io)
    chunk_size = read_header(io)
    subchunk_id = Vector{UInt8}(undef, 4)
    read!(io, subchunk_id)
    subchunk_size = read_le(io, UInt32)
    if subchunk_id != b"fmt "
        error("First chunk is not the format")
    end
    fmt = read_format(io, subchunk_size)

    if fmt.nchannels != size(samples,2)
        error("Number of channels do not match")
    end

    # Compute data length of current chunk to-be-appended.
    data_length = size(samples, 1) * fmt.block_align
    # Update `chunksize`: add length of new data.
    seek(io,4)
    write_le(io, convert(UInt32, chunk_size + data_length))
    # Get `subchunk2size`.
    seek(io, 24 + subchunk_size)
    data_length_old = read_le(io, UInt32)
    # Update `subchunk2size`: add length of new data.
    seek(io, 24 + subchunk_size)
    write_le(io, convert(UInt32, data_length_old + data_length))

    seekend(io)
    write_data(io, fmt, samples)
end

function wavappend(samples::AbstractArray, filename::AbstractString)
    open(filename, "a+") do io
        wavappend(samples,io)
    end
end

wavwrite(y::AbstractArray, f::Real, filename::AbstractString) = wavwrite(y, filename, Fs=f)
wavwrite(y::AbstractArray, f::Real, n::Real, filename::AbstractString) = wavwrite(y, filename, Fs=f, nbits=n)

# support for writing native arrays...
wavwrite(y::AbstractArray{T}, io::IO) where T <: Integer = wavwrite(y, io, nbits=sizeof(T)*8)
wavwrite(y::AbstractArray{T}, filename::AbstractString) where T <: Integer = wavwrite(y, filename, nbits=sizeof(T)*8)
wavwrite(y::AbstractArray{Int32}, io::IO) = wavwrite(y, io, nbits=24)
wavwrite(y::AbstractArray{Int32}, filename::AbstractString) = wavwrite(y, filename, nbits=24)
wavwrite(y::AbstractArray{T}, io::IO) where T <: AbstractFloat = wavwrite(y, io, nbits=sizeof(T)*8, compression=WAVE_FORMAT_IEEE_FLOAT)
wavwrite(y::AbstractArray{T}, filename::AbstractString) where T <: AbstractFloat = wavwrite(y, filename, nbits=sizeof(T)*8, compression=WAVE_FORMAT_IEEE_FLOAT)

# FileIO integration support
load(s::Stream{format"WAV"}; kwargs...) = wavread(s.io; kwargs...)
save(s::Stream{format"WAV"}, data; kwargs...) = wavwrite(data, s.io; kwargs...)

load(f::File{format"WAV"}; kwargs...) = wavread(f.filename; kwargs...)
save(f::File{format"WAV"}, data; kwargs...) = wavwrite(data, f.filename; kwargs...)

"""
    wavplay(data, fs)
    wavplay(filename::AbstractString)

Plays the audio waveform `data` at sampling frequency `fs` (in hertz).
To play a stereo signal, provide two columns in array `data` (left and
right channel), as in [`wavwrite`](@ref).

The `filename` form reads both waveform data and sampling frequency
from the named WAV file to play it.

The supported backends are:
* AudioQueue (macOS)
* PulseAudio (Linux, libpulse-simple)
* PlaySound  (Windows)

See also: [`wavwrite`](@ref)
"""
function wavplay end
wavplay(fname::AbstractString) = wavplay(wavread(fname)[1:2]...)
@static if Sys.islinux()
    include("wavplay-pulse.jl")
elseif Sys.isapple()
    include("wavplay-audioqueue.jl")
elseif Sys.iswindows()
    include("wavplay-win32.jl")
else
    wavplay(data, fs) = @warn "wavplay is not currently implemented on $(Sys.KERNEL)"
end

end # module


text:
module DefaultScheduler

import ..DagScheduler: Dagger, RunEnv, NodeEnv, Thunk, Chunk, DRef, FileRef, ComputeCost, TaskIdType, istask,
        cost_function, taskid, walk_dag, fval, SPLIT_CAPACITY_RATIO, SPLIT_COST_RATIO,
        find_node

#-------------------------------------
# Execution stages based on dependency
#-------------------------------------
function execution_stages(dag_node::Thunk, deps::Dict{Thunk, Set{Thunk}}=Dagger.dependents(dag_node), root::Thunk=dag_node)
    newinps = Set{Thunk}()
    for inp in dag_node.inputs
        if istask(inp)
            for cinp in execution_stages(inp, deps, root)
                push!(newinps, cinp)
            end
        end
    end
    ((dag_node == root) || (length(deps[dag_node]) > 1)) ? [copy_thunk(dag_node, newinps)] : [newinps...]
end

#----------------------------------------------
# Execution stages based on costs and affinity
#----------------------------------------------
cost_for_nodes(runenv::RunEnv, root_t::Thunk) = [node.brokerid => cost_for_node(node, root_t) for node in runenv.nodes]
function cost_for_node(node::NodeEnv, root_t::Thunk)
    nc = Dict{TaskIdType,ComputeCost}()
    compute_cost(root_t, [node.brokerid, node.executorids...], node.host, nc)
    nc
end

function extend_stages_by_affinity(runenv::RunEnv, stages::Vector{Thunk}, root_t::Thunk, costs::Vector{Pair{UInt64,Dict{TaskIdType,ComputeCost}}}, nodecap::Float64, deps::Dict{Thunk, Set{Thunk}}=Dagger.dependents(root_t))
    for stage in stages
        #add_stages_by_compute_cost(stage, root_t, costs, nodecap)

        walk_dag(stage, false) do n,d
            add_stages_by_compute_cost(n, root_t, costs, nodecap, deps)
        end
    end
    nothing
end

compute_cost(dag_node, processorids::Vector{UInt64}, host::IPAddr, costs::Dict{TaskIdType,ComputeCost}) = ComputeCost(0.0, 0.0)
function compute_cost(dag_node::Thunk, processorids::Vector{UInt64}, host::IPAddr, costs::Dict{TaskIdType,ComputeCost})
    inps = dag_node.inputs
    costs[taskid(dag_node)] = isempty(inps) ? ComputeCost(1.0, 0.0) : cost_function(dag_node)(map((inp)->compute_cost(inp, processorids, host, costs), inps))
end
function compute_cost(dag_node::Chunk, processorids::Vector{UInt64}, host::IPAddr, costs::Dict{TaskIdType,ComputeCost})
    # network transfer cost (arbitrarily set as 1 per 1MB - needs to be parameterized)
    netcost = 0.0

    if (isa(dag_node.handle, DRef) && !(dag_node.handle.owner in processorids)) || (isa(dag_node.handle, FileRef) && !(dag_node.handle.host == host))
        netcost = ceil(Int, dag_node.handle.size/10^6)
    end

    # cpu cost arbitrarily set to 1 for all functions - needs to be provided in configuration/function annotation
    ComputeCost(1.0, netcost)
end

function schedule_target(id::TaskIdType, costs::Vector{Pair{UInt64,Dict{TaskIdType,ComputeCost}}})
    target = UInt64(0)
    target_cost = ComputeCost(0.0, 0.0)
    for (brokerid,brokercosts) in costs
        brokercost = get(brokercosts, id, target_cost)
        if (target === UInt64(0)) || (brokercost < target_cost)
            target = brokerid
            target_cost = brokercost
        end
    end

    cost_comparison = Vector{Pair{UInt64,Float64}}()
    total_costs = Vector{Pair{UInt64,Float64}}()
    for (brokerid,brokercosts) in costs
        push!(total_costs, brokerid => fval(get(brokercosts, id, ComputeCost(0.0, 0.0))))
        push!(cost_comparison, brokerid => (get(brokercosts, id, target_cost) / target_cost))
    end
    sort!(total_costs, lt=(x,y)->x[2]<y[2])
    sort!(cost_comparison, lt=(x,y)->x[2]<y[2])

    total_costs, cost_comparison
end

function optimum_affinity(tid::TaskIdType, costs::Vector{Pair{UInt64,Dict{TaskIdType,ComputeCost}}}, avg_node_capacity::Float64)
    total_costs, cost_comparison = schedule_target(tid, costs)
    insignificant_cost(total_costs, avg_node_capacity) || strong_affinity(cost_comparison)
end

function strong_affinity(cost_comparison)
    (length(cost_comparison) < 2) ||                                    # affinity with just one node
    (cost_comparison[2][2] >= SPLIT_COST_RATIO)                         # cost on one node substantially less than others
end

function insignificant_cost(total_costs, avg_node_capacity)
    (total_costs[1][2] <= (SPLIT_CAPACITY_RATIO * avg_node_capacity))   # total cost is not significant enough to split
end

function add_stages_by_compute_cost(staged_node::Thunk, full_dag::Thunk, costs::Vector{Pair{UInt64,Dict{TaskIdType,ComputeCost}}}, avg_node_capacity::Float64, deps::Dict{Thunk, Set{Thunk}})
    sn_tc, sn_cc = schedule_target(taskid(staged_node), costs)
    (insignificant_cost(sn_tc, avg_node_capacity) || strong_affinity(sn_cc)) && return

    # if no clear affinity, search till we find a node with one
    stop_nodes = Vector{TaskIdType}()
    for inp in staged_node.inputs
        istask(inp) && push!(stop_nodes, taskid(inp))
    end

    let chld=Set{Thunk}(), start_node=find_node(taskid(staged_node))
        walk_dag(start_node, false, stop_nodes) do n,d
            istask(n) || return

            total_costs, cost_comparison = schedule_target(taskid(n), costs)
            ((total_costs[1][2]/sn_tc[1][2] <= 0.5) ||  insignificant_cost(total_costs, avg_node_capacity) || !strong_affinity(cost_comparison)) && return

            # find children of this node if there are any
            nchld = filter(c->ischild(n, taskid(c)), chld)
            setdiff!(chld, nchld)

            # keep children only if substantial cost difference or multi-dependency node
            filter!(nchld) do c
                (length(deps[c]) > 1) && (return true)
                c_tc, c_cc = schedule_target(taskid(c), costs)
                total_costs[1][2]/c_tc[1][2] <= 0.5
            end

            # create a new staged node (with its children) and insert
            staged_n = copy_thunk(n, nchld)
            push!(chld, staged_n)
        end

        # insert new split points at staged node
        if !isempty(chld)
            #println("found nodes!! $(chld)")
            for inp in staged_node.inputs
                push!(chld, inp)
            end
            staged_node.inputs = tuple(chld...)
        end
    end
end

copy_thunk(dag_node, newinps) = Dagger.Thunk(dag_node.f, newinps...;
            id = dag_node.id,
            get_result = dag_node.get_result,
            meta = dag_node.meta,
            persist = dag_node.persist,
            cache = dag_node.cache,
            cache_ref = dag_node.cache_ref,
            affinity = dag_node.affinity)

end # module DefaultScheduler

function schedule(runenv::RunEnv, root_t::Thunk)
    execstages = DefaultScheduler.execution_stages(root_t)

    if length(runenv.nodes) > 1
        costs = DefaultScheduler.cost_for_nodes(runenv, root_t)
        nodecap = mean_node_capacity(runenv)
        DefaultScheduler.extend_stages_by_affinity(runenv, execstages, root_t, costs, nodecap)
    else
        costs = []
    end

    execstages, costs
end

function affinity(tid::TaskIdType, brokerid::UInt64, costs::Vector{Pair{UInt64,Dict{TaskIdType,ComputeCost}}})
    tc, cc = DefaultScheduler.schedule_target(tid, costs)

    if !isempty(cc)
        for (b,c) in cc
            (b === brokerid) && (return c)
        end
    end

    1.0
end

ext: jl


text:
module TransformExport

using Codex
using CSV
using DataFrames
using DataFrames: ColumnIndex
using Dates

export
    read_csv,
    responses,
    simplify,
    substitute_names,
    names2usernames,
    process,
    rm_descriptions,
    split_datetime

"""
    read_csv(path; delim)::DataFrame

Copies CSV at `path` into memory.
"""
function read_csv(path; delim=',')::DataFrame
    try
        CSV.read(path, DataFrame; delim, stringtype=String)
    catch
        # CSV version 0.8.
        CSV.read(path, DataFrame; delim)
    end
end

"""
    responses(dir::String)::Dict{String,DataFrame}

Return responses for an export folder such as "2020-08".
"""
function responses(dir::String)::Dict{String,DataFrame}
    dir = joinpath(dir, "responses")
    files = filter(endswith(".csv"), readdir(dir))
    names = Codex.rmextension.(files)
    paths = joinpath.(dir, files)
    dfs = map(path -> read_csv(path, delim=';'), paths)
    Dict(zip(names, dfs))
end

rm_timing!(df::DataFrame) = select!(df, Not(r".+\_timing"))
rm_boring_timestamps!(df::DataFrame) = select!(df, Not(r".+\_(?<!completed_)(from|at)"))
rm_boring_foreign_keys!(df::DataFrame) = select!(df, Not(r".+\_(?<!by_)id"))
rm_locale!(df::DataFrame) = "locale" in names(df) ? select!(df, Not("locale")) : df
rm_empty_rows!(df::DataFrame) = dropmissing!(df, :filled_out_by_id)
rename_id_col!(df::DataFrame) = rename!(df, Dict(:filled_out_by_id => :id))

"""
    simplify(df)::DataFrame

Renames id column after removing extraneous rows and columns, that is, removes empty rows and
removes columns such as `protocol_subscription_id`, `open_from` and `v2_1_timing`.
"""
simplify(df::DataFrame)::DataFrame = DataFrame(
    apply([rm_timing!, rm_boring_timestamps!, rm_boring_foreign_keys!, rm_locale!, rm_empty_rows!, rename_id_col!], df)
)

"""
    _description_regex()

Return regex for matching a description such as `1 (lorem)` or `2 <br /> (ipsum)`.
"""
_description_regex() = r"([0-9]+) (<br\w*\/> )?\([^\)]+\)"

"""
    _rm_description(e::String)::String

Apply regex replace on element `e`.
"""
_rm_description(e) = e === missing ? e : parse(Int, replace(e, _description_regex() => s"\1"))

"""
    _rm_descriptions(col)::Array{Int,1}

Apply a regex replace and type conversion to all elements of the column `col`.
"""
_rm_descriptions(col) = map(_rm_description, col)

"""
    _contains_description(col)::Bool

Return whether the column `col` contains descriptions.
"""
function _contains_description(col)::Bool
    function is_match(value)
        if value === missing
            false
        else
            contains(value, _description_regex())
        end
    end
    any(is_match, col)
end

"""
    rm_descriptions(df)::DataFrame

Find responses containing a description, for example `6 (heel erg)`, and remove the description.
"""
function rm_descriptions(df)::DataFrame
    df = DataFrame(df)
    # Not using in-place replacement since the type has to change.
    function map_col(i::Number)::Array
        col = df[!, 1]
        _rm_descriptions(col)
    end
    for i in 1:ncol(df)
        col = df[!, i]
        if contains(string(typeof(col)), "String") && _contains_description(col)
            df[!, i] = _rm_descriptions(col)
        end
    end
    df
end

"""
    substitute_names(df, with::DataFrame)::DataFrame
    substitute_names(with)::Function

Replaces `person_id`s by the first name as listed in `with`.
"""
function substitute_names(df, with::DataFrame)::DataFrame
    mapping = Dict(zip(with.person_id, with.first_name))
    select!(df, :id => ByRow(id -> mapping[id]), Not(:id))
    rename!(df, Dict(:id_function => :id))
end
substitute_names(with) = df -> substitute_names(df, with)

"""
    parsedatetime(str)::DateTime

Parse a date and time string from the export to a Julia DateTime object.
"""
parsedatetime(s) = DateTime(s, "d-m-y H:M:S")

"""
    process(in_dir, out_dir; fns)

Processes the responses from the export folder, applies the functions `fns` and places the files at `out_dir`.
"""
function process(in_dir, out_dir; fns=nothing)
    # mkpath(out_dir)
    out_dir
end

function names2usernames(df, id_username)::DataFrame
    df.id = Codex.map_by_df(df.id, id_username, :id, :username)
    df
end

"""
    split_datetime(df::DataFrame, datetime_col::ColumnIndex)::DataFrame

Split the datetime column `datetime_col` into two columns, namely one for date and one for time.
"""
function split_datetime(df::DataFrame, datetime_col::T)::DataFrame where {T<:ColumnIndex}
    df = DataFrame(df)
    dates = map(x -> first(split(x, " ")), df[!, datetime_col])
    dates = Date.(dates, dateformat"dd-mm-yyyy")
    df.date = dates
    df.time = map(x -> split(x, " ")[2], df[!, datetime_col])
    select!(df, Not(datetime_col))
    df
end

end # module


text:
struct Marks
    # Union to avoid creating a vector in most cases
    soft::Dict{String, Union{Symbol, Vector{Symbol}}} # TODO: should be a MultiDict
    hard::Vector{Symbol} # attached to all instances of a @testset
end

Marks() = Marks(Dict{String, Union{Symbol, Vector{Symbol}}}(), Symbol[])

const _pass = :__pass__
const _fail = :__fail__

# marks contains at most one of _pass, _fail

# return true (pass), false (fail), or nothing (unrun)
function pastresult(marks::Marks, subject)
    ms = get(marks.soft, subject, Symbol())
    sym::Symbol = ms isa Symbol ? ms : isempty(ms) ? Symbol() : ms[1]
    sym === _pass ? true : sym === _fail ? false : nothing
end

function setresult!(marks::Marks, subject, success::Bool)
    soft = marks.soft
    ms = get(soft, subject, Symbol())
    res = success ? _pass : _fail
    if ms isa Symbol
        # ms could be Symbol() from get's default or delmark!
        if ms ∈ (_pass, _fail, Symbol())
            soft[subject] = res
        else
            # res always in first position
            soft[subject] = [res, ms]
        end
    else # ms isa Vector
        if !isempty(ms) && ms[1] ∈ (_pass, _fail)
            ms[1] = res
        else
            pushfirst!(ms, res)
        end
    end
end

function markiter(marks::Marks, subject, skipres::Bool)
    ms = get(marks.soft, subject, Symbol())
    if ms isa Symbol
        if ms === Symbol() || skipres && ms ∈ (_pass, _fail)
            marks.hard
        else
            Iterators.flatten((marks.hard, (ms,)))
        end
    else
        Iterators.flatten((marks.hard,
                           if skipres
                               Iterators.filter(m -> m ∉ (_pass, _fail), ms)
                           else
                               ms
                           end))
    end
end

function addmark!(marks::Marks, subject, m::Symbol)
    m ∈ marks.hard && return false
    soft = marks.soft
    ms = get(soft, subject, Symbol())
    if ms isa Symbol
        if ms === m
            false
        elseif ms === Symbol()
            soft[subject] = m
            true
        else
            soft[subject] = [ms, m]
            true
        end
    elseif m ∉ ms
        push!(ms, m)
        true
    else
        false
    end
end

# return whether warning is issued
function delmark!(marks::Marks, subject, m::Symbol, warned::Bool=false)
    if m ∈ marks.hard && !warned
        @warn "cannot remove statically attached label (@testset :$m ...)"
        return true
    end
    soft = marks.soft
    ms = get(soft, subject, Symbol())
    if ms isa Symbol
        if ms === m
            soft[subject] = Symbol()
        end
    else
        p = findfirst(==(m), ms)
        if p !== nothing
            deleteat!(ms, p)
        end
    end
    false
end

function hasmark(marks::Marks, subject, m::Symbol)
    m ∈ marks.hard && return true
    ms = get(marks.soft, subject, Symbol())
    if ms isa Symbol
        ms === m
    else
        m ∈ ms
    end
end

max_issues_repo_name: rfourquet/InlineTest

max_forks_repo_name: rfourquet/InlineTest


text:
using InteractiveChaos, Makie

i = 1
p_index = 1

systems = [(Systems.logistic(), 2.5, 4.0, "r"),
           (Systems.henon(), 0.8, 1.4, "a"),
           (Systems.standardmap(), 0.6, 1.2, "k")]

ds, p_min, p_max, parname = systems[1]

oddata = interactive_orbitdiagram(
                             ds, i, p_index, p_min, p_max;
                             parname = parname, n = 100)

ps, us = scaleod(oddata...)


text:
# allows for registering listeners that handle functions of the form:
#		somefunction(listener, args...)

# a function can be broadcast to all listeners using the form:
#		broadcastto(somefunction, (uid, exch), args...)

# to register a listener:
#		listenfor(somefunction, listener, vectorOfUIDsToFilter, vectorOfExchangesToFilter)

# Example:

# julia> callback(x, y) = println(x, " ", y)
#   callback (generic function with 1 method)
# julia> listenfor(callback, "someobject", [5,6])
# julia> listenfor(callback, "someotherobject", [6,7])
# julia> broadcastto(callback, 5, "payload")
#   someobject payload
# julia> broadcastto(callback, 6, "payload")
#   someobject payload
#   someotherobject payload
# julia> broadcastto(callback, 7, "payload")
#   someotherobject payload

##############################################
# matrix version
##############################################

type ListenerFilter
	# each vector in the matrix is a list of listeners for that (uid, exch) pair
	# note: the final row (maxuid+1) is for those listeners that don't filter on UID
	# note: the final col (NUMEXCH+1) is for those listeners that don't filter on exch
	# note: the bottom right corner (maxuid+1, NUMEXCH+1) is for listeners with no filter
	maxuid::Int
	filterMatrix::Matrix{Vector{Any}}
	nrows::Int
	ncols::Int

	function ListenerFilter(maxuid::Integer)
		nr = maxuid + 1
		nc = typemax(Exchange).val + 1
		a = Array(Vector{Any}, nr, nc)
		for r in 1:nr
			for c in 1:nc
				a[r,c] = Any[]
			end
		end
		new(maxuid, a, nr, nc)
	end
end

function addListener(lf::ListenerFilter, listener, items1::VecI = Int[], items2::VecI = Int[])
	rows = (length(items1) == 0 ? [lf.maxuid+1] :  items1)
	cols = (length(items2) == 0 ? [typemax(Exchange).val+1] : items2)
	for r in rows
		for c in cols
			# println("$r $c $rows $cols")
			push!(lf.filterMatrix[r,c], listener)
		end
	end
	# println("ADDED LISTENER: ", lf, " ", listener, " ", items1, " ", items2)
end

type Broadcaster
	maxuid::Int
	registry::Dict{Function, ListenerFilter}

	Broadcaster() = new()
	# Broadcaster(maxuid::Int) = new(maxuid, Dict{Function, ListenerFilter}())
end

Base.string(b::Broadcaster) = "Broadcaster{maxuid=$(b.maxuid), functions=$(keys(b.registry))}"
Base.print(io::IO, b::Broadcaster) = print(io, string(b))
Base.show(io::IO, b::Broadcaster) = print(io, string(b))

# set up a const global (singleton) BROADCASTER
const BROADCASTER = Broadcaster()

function initBroadcaster(maxuid::Integer)
	BROADCASTER.maxuid = maxuid
	BROADCASTER.registry = Dict{Function, ListenerFilter}()
end

# function listen(bc::Broadcaster, f::Function, listener, items1::Vector{UID} = UID[], items2::Vector{Exchange} = Exchange[])
function listenfor(f::Function, listener, items1::Vector{UID} = UID[], items2::Vector{Exchange} = Exchange[])
	DEBUG("listenfor: $f $listener")
	lf = get!(BROADCASTER.registry, f, ListenerFilter(BROADCASTER.maxuid))
	addListener(lf, listener, round(Int, items1), [Int(exch.val) for exch in items2])
end

listenfor(f::Function, listener, uid::UID) = listenfor(f, listener, [uid])
listenfor(f::Function, listener, exch::Exchange) = listenfor(f, listener, UID[], [exch])
listenfor(f::Function, listener, uid::UID, exch::Exchange) = listenfor(f, listener, [uid], [exch])

# call this function in order to call function "f" for all listeners, with that listener as the first arg
# note we can set UID to 0, or Exchange to NONE, and it will not filter on that.
#		this allows us the flexibility to send/listen to exchange-only messages, or symbol-only messages

# function broadcastto(bc::Broadcaster, f::Function, vals::Tuple{UID, Exchange}, args...)
function broadcastto(f::Function, vals::Tuple{UID, Exchange}, args...)
	if !haskey(BROADCASTER.registry, f)
		return
	end

	lf = BROADCASTER.registry[f]
	r = Int(vals[1])
	c = Int(vals[2])

	for listener in lf.filterMatrix[lf.nrows, lf.ncols]  # no filter
		f(listener, args...)
	end

	if r != 0
		for listener in lf.filterMatrix[r, lf.ncols] # filter on UID
			f(listener, args...)
		end
	end

	if c != 0
		for listener in lf.filterMatrix[lf.nrows, c] # filter on Exchange
			f(listener, args...)
		end
	end

	if r != 0 && c != 0
		for listener in lf.filterMatrix[r, c] # filter on both
			f(listener, args...)
		end
	end
end


text:
function __init__()::Nothing
    _import_all_on_init()
    _print_welcome_message()
    return nothing
end


text:
@doc doc"""
    latexalign()
Generate a ``LaTeX`` align environment from an input.

# Examples
## use with arrays

```julia
lhs = [:(dx/dt), :(dy/dt), :(dz/dt)]
rhs = [:(y-x), :(x*z-y), :(-z)]
latexalign(lhs, rhs)
```

```LaTeX
\begin{align}
\frac{dx}{dt} =& y - x \\\\
\frac{dy}{dt} =& x \cdot z - y \\\\
\frac{dz}{dt} =& - z \\\\
\end{align}
```

## use with ParameterizedFunction

```julia-repl
julia> using DifferentialEquations
julia> ode = @ode_def foldChangeDetection begin
    dm = r_m * (i - m)
    dy = r_y * (p_y * i/m - y)
end i r_m r_y p_y

julia> latexalign(ode)
```
```LaTeX
\begin{align}
\frac{dm}{dt} =& r_{m} \cdot \left( i - m \right) \\\\
\frac{dy}{dt} =& r_{y} \cdot \left( \frac{p_{y} \cdot i}{m} - y \right) \\\\
\end{align}
```

"""
latexalign(args...; kwargs...) = latexify(args...; kwargs..., env=:align)

function _latexalign(arr::AbstractMatrix; separator=" =& ", double_linebreak=false, starred=false, rows=:all, aligned=false, kwargs...)
    eol = double_linebreak ? " \\\\\\\\\n" : " \\\\\n"
    arr = latexraw.(arr; kwargs...)
    separator isa String && (separator = fill(separator, size(arr)[1]))
    str = ""
    if aligned
        str *= "\\begin{aligned}\n"
    else
        str *= "\\begin{align$(starred ? "*" : "")}\n"
    end
    if rows == :all
        iterate_rows = 1:(size(arr)[1])
    else
        iterate_rows = rows
    end

    for i in iterate_rows
        if i != last(iterate_rows)
            str *= join(arr[i,:], separator[i]) * eol
        else
            str *= join(arr[i,:], separator[i]) * "\n"
        end
    end
    if aligned
        str *= "\\end{aligned}\n"
    else
        str *= "\\end{align$(starred ? "*" : "")}\n"
    end
    latexstr = LaTeXString(str)
    COPY_TO_CLIPBOARD && clipboard(latexstr)
    return latexstr
end

function _latexalign(lhs::AbstractArray, rhs::AbstractArray; kwargs...)
    return latexalign(hcat(lhs, rhs); kwargs...)
end

function _latexalign(lhs::Tuple, rhs::Tuple; kwargs...)
    return latexalign(hcat(collect(lhs), collect(rhs)); kwargs...)
end

_latexalign(args::Tuple...; kwargs...) = latexalign(reduce(hcat, [collect(i) for i in args]); kwargs...)

_latexalign(arg::Tuple; kwargs...) = latexalign(reduce(hcat, [collect(i) for i in arg]); kwargs...)

function _latexalign(nested::AbstractVector{AbstractVector}; kwargs...)
    return latexalign(reduce(hcat, nested); kwargs...)
end

function _latexalign(d::AbstractDict; kwargs...)
    latexalign(collect(keys(d)), collect(values(d)); kwargs...)
end

"""
    _latexalign(vec::AbstractVector)

Go through the elements, split at any = sign, pass on as a matrix.
"""
function _latexalign(vec::AbstractVector; kwargs...)
    lvec = _latexraw.(vec; kwargs...)
    ## turn the array into a matrix
    lmat = reduce(hcat, split.(lvec, " = "))
    ## turn the matrix ito arrays of left-hand-side, right-hand-side.
    larr = [lmat[i,:] for i in 1:size(lmat, 1)]
    length(larr) < 2 && error("Invalid intput to _latexalign().")
    return latexalign( reduce(hcat, larr) ; kwargs...)
end


text:
"""
    struct SymMatrix{T} <: AbstractMatrix{T}
        Q::Vector{T}
        n::Int
    end

Symmetric ``n \\times n`` matrix storing the vectorized upper triangular
part of the matrix in the `Q` vector (this corresponds to the vectorized
form of `MathOptInterface.PositiveSemidefiniteConeTriangle`).
It implement the `AbstractMatrix` interface except for `setindex!` as it might
break its symmetry. The [`symmetric_setindex!`](@ref) function should be used
instead.
"""
struct SymMatrix{T} <: AbstractMatrix{T}
    Q::Vector{T}
    n::Int
end

# i <= j
#function trimap(i, j, n) # It was used when Q was the lower triangular part
#    div(n*(n+1), 2) - div((n-i+1)*(n-i+2), 2) + j-i+1
#end

# j <= i
function trimap(i, j)
    div((i-1)*i, 2) + j
end

function trimat(::Type{T}, f, n, σ) where {T}
    Q = Vector{T}(undef, trimap(n, n))
    for i in 1:n
        for j in 1:i
            Q[trimap(i, j)] = f(σ[i], σ[j])
        end
    end
    SymMatrix{T}(Q, n)
end

Base.size(Q::SymMatrix) = (Q.n, Q.n)

"""
    symmetric_setindex!(Q::SymMatrix, value, i::Integer, j::Integer)

Set `Q[i, j]` and `Q[j, i]` to the value `value`.
"""
function symmetric_setindex!(Q::SymMatrix, value, i::Integer, j::Integer)
    Q.Q[trimap(max(i, j), min(i, j))] = value
end

function Base.getindex(Q::SymMatrix, i::Integer, j::Integer)
    return Q.Q[trimap(max(i, j), min(i, j))]
end
Base.getindex(Q::SymMatrix, I::Tuple) = Q[I...]
Base.getindex(Q::SymMatrix, I::CartesianIndex) = Q[I.I]


text:
function drelu_dZ(NN::NeuralNetwork, l::Int)
    return @. NN.cache.Z[l] > 0
end

function dlogistic_dZ(NN::NeuralNetwork, l::Int)
    return @. NN.cache.A[l] * (1 - NN.cache.A[l])
end

function dtanh_dZ(NN::NeuralNetwork, l::Int)
    return @. 1 - NN.cache.A[l]^2
end

function dsoftmax_dZ(NN::NeuralNetwork, l::Int)
    return @. NN.cache.A[l] * (1 - NN.cache.A[l])
end

function propagate_back!(NN::NeuralNetwork, mini_batch::MiniBatch)
    for l = NN.L:-1:1
        compute_dJ_dA!(NN, l, mini_batch)
        compute_dJ_dZ!(NN, l)
        compute_dJ_dW!(NN, l, mini_batch)
        compute_dJ_db!(NN, l, mini_batch)
    end
end

function compute_dJ_dA!(NN::NeuralNetwork, l::Int, mini_batch::MiniBatch)
    if l ≠ NN.L
        NN.cache.dA[l] = NN.params.W[l+1]' * NN.cache.dZ[l+1]
    else
        if NN.params.loss_function == cross_entropy_binary
            NN.cache.dA[l] = @. -mini_batch.Y / NN.cache.A[NN.L] +
               (1 - mini_batch.Y) / (1 - NN.cache.A[NN.L])
        elseif NN.params.loss_function == cross_entropy_multi
            NN.cache.dA[l] = @. -mini_batch.Y / NN.cache.A[NN.L]
        elseif NN.params.loss_function == mean_squared_error
            NN.cache.dA[l] = @. NN.cache.A[NN.L] - mini_batch.Y
        else
            error("unknown loss function provided in params")
        end
    end
end

function compute_dJ_dZ!(NN::NeuralNetwork, l::Int)
    if NN.params.activation_functions[l] == tanh
        NN.cache.dZ[l] = NN.cache.dA[l] .* drelu_dZ(NN, l)
    elseif NN.params.activation_functions[l] == relu
        NN.cache.dZ[l] = NN.cache.dA[l] .* dtanh_dZ(NN, l)
    elseif NN.params.activation_functions[l] == logistic
        NN.cache.dZ[l] = NN.cache.dA[l] .* dlogistic_dZ(NN, l)
    elseif NN.params.activation_functions[l] == softmax
        NN.cache.dZ[l] = NN.cache.dA[l] .* dsoftmax_dZ(NN, l)
    else
        error("unknown activation function(s) provided in params")
    end
end

function compute_dJ_dW!(NN::NeuralNetwork, l::Int, mini_batch::MiniBatch)
    if l ≠ 1
        NN.cache.dW[l] = (NN.cache.dZ[l] ./ size(mini_batch.X, 2)) * NN.cache.A[l-1]'
    else
        NN.cache.dW[l] = (NN.cache.dZ[l] ./ size(mini_batch.X, 2)) * mini_batch.X'
    end

    # Add in regularization if using gradient descent
    if NN.hparams.optimization == :gd
        NN.cache.dW[l] += (NN.hparams.λ / size(mini_batch.X, 2)) .* NN.params.W[l]
    end
end

function compute_dJ_db!(NN::NeuralNetwork, l::Int, mini_batch::MiniBatch)
    NN.cache.db[l] = sum(NN.cache.dZ[l], dims = 2) ./ size(mini_batch.X, 2)
end


text:
module JlrsTests
mutable struct MutableStruct
    x
    y::UInt64
end

struct ParameterStruct{T}
    a::T
end

struct ValueTypeStruct{T}
end

valuedispatch(::ValueTypeStruct{1})::Int = 3
valuedispatch(::ValueTypeStruct{2})::Float64 = 3.0

function inlinetuple()::Tuple{UInt32, UInt16, Int64}
    (1, 2, 3)
end

function symbol()::Symbol
    :x
end

function base()::Module
    Base
end

function datatype()::DataType
    Bool
end

function callrust(ptr::Ptr)::Bool
    ccall(ptr, Bool, ())
end

function callrustwitharr(ptr::Ptr, arr::Array{Float64, 1})::Bool
    ccall(ptr, Bool, (Array,), arr)
end

function vecofmodules()::Vector{Module}
    [Base; Core; Main]
end

function anothervecofmodules()::Vector{Module}
    [Base; Core; Main]
end

function funcwithkw(a::Int; b::Int=1)
    a + b
end

function funcwithkw(a::Int, rest...; b::Int=1)
    a + sum(rest) + b
end

function funcwithabstractkw(a::Float32; b::Real=1.0f0)
    a + b
end
end

max_stars_repo_name: gobuki/jlrs

max_issues_repo_name: gobuki/jlrs

max_forks_repo_name: gobuki/jlrs


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct OneDOneQMachine <: Machine
        R::Float64
        Xd::Float64
        Xq::Float64
        Xd_p::Float64
        Xq_p::Float64
        Td0_p::Float64
        Tq0_p::Float64
        MVABase::Float64
        ext::Dict{String, Any}
        states::Vector{Symbol}
        n_states::Int64
        internal::InfrastructureSystemsInternal
    end

Parameters of 4-states synchronous machine: Simplified Marconato model
 The derivative of stator fluxes (ψd and ψq) is neglected and ωψd = ψd and
 ωψq = ψq is assumed (i.e. ω=1.0). This is standard when
 transmission network dynamics is neglected.

# Arguments
- `R::Float64`: Resistance after EMF in machine per unit, validation range: (0, nothing)
- `Xd::Float64`: Reactance after EMF in d-axis per unit, validation range: (0, nothing)
- `Xq::Float64`: Reactance after EMF in q-axis per unit, validation range: (0, nothing)
- `Xd_p::Float64`: Transient reactance after EMF in d-axis per unit, validation range: (0, nothing)
- `Xq_p::Float64`: Transient reactance after EMF in q-axis per unit, validation range: (0, nothing)
- `Td0_p::Float64`: Time constant of transient d-axis voltage, validation range: (0, nothing)
- `Tq0_p::Float64`: Time constant of transient q-axis voltage, validation range: (0, nothing)
- `MVABase::Float64`: Nominal Capacity in MVA, validation range: (0, nothing)
- `ext::Dict{String, Any}`
- `states::Vector{Symbol}`
- `n_states::Int64`
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct OneDOneQMachine <: Machine
    "Resistance after EMF in machine per unit"
    R::Float64
    "Reactance after EMF in d-axis per unit"
    Xd::Float64
    "Reactance after EMF in q-axis per unit"
    Xq::Float64
    "Transient reactance after EMF in d-axis per unit"
    Xd_p::Float64
    "Transient reactance after EMF in q-axis per unit"
    Xq_p::Float64
    "Time constant of transient d-axis voltage"
    Td0_p::Float64
    "Time constant of transient q-axis voltage"
    Tq0_p::Float64
    "Nominal Capacity in MVA"
    MVABase::Float64
    ext::Dict{String, Any}
    states::Vector{Symbol}
    n_states::Int64
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function OneDOneQMachine(R, Xd, Xq, Xd_p, Xq_p, Td0_p, Tq0_p, MVABase, ext=Dict{String, Any}(), )
    OneDOneQMachine(R, Xd, Xq, Xd_p, Xq_p, Td0_p, Tq0_p, MVABase, ext, [:eq_p, :ed_p], 2, InfrastructureSystemsInternal(), )
end

function OneDOneQMachine(; R, Xd, Xq, Xd_p, Xq_p, Td0_p, Tq0_p, MVABase, ext=Dict{String, Any}(), )
    OneDOneQMachine(R, Xd, Xq, Xd_p, Xq_p, Td0_p, Tq0_p, MVABase, ext, )
end

# Constructor for demo purposes; non-functional.
function OneDOneQMachine(::Nothing)
    OneDOneQMachine(;
        R=0,
        Xd=0,
        Xq=0,
        Xd_p=0,
        Xq_p=0,
        Td0_p=0,
        Tq0_p=0,
        MVABase=0,
        ext=Dict{String, Any}(),
    )
end

"""Get OneDOneQMachine R."""
get_R(value::OneDOneQMachine) = value.R
"""Get OneDOneQMachine Xd."""
get_Xd(value::OneDOneQMachine) = value.Xd
"""Get OneDOneQMachine Xq."""
get_Xq(value::OneDOneQMachine) = value.Xq
"""Get OneDOneQMachine Xd_p."""
get_Xd_p(value::OneDOneQMachine) = value.Xd_p
"""Get OneDOneQMachine Xq_p."""
get_Xq_p(value::OneDOneQMachine) = value.Xq_p
"""Get OneDOneQMachine Td0_p."""
get_Td0_p(value::OneDOneQMachine) = value.Td0_p
"""Get OneDOneQMachine Tq0_p."""
get_Tq0_p(value::OneDOneQMachine) = value.Tq0_p
"""Get OneDOneQMachine MVABase."""
get_MVABase(value::OneDOneQMachine) = value.MVABase
"""Get OneDOneQMachine ext."""
get_ext(value::OneDOneQMachine) = value.ext
"""Get OneDOneQMachine states."""
get_states(value::OneDOneQMachine) = value.states
"""Get OneDOneQMachine n_states."""
get_n_states(value::OneDOneQMachine) = value.n_states
"""Get OneDOneQMachine internal."""
get_internal(value::OneDOneQMachine) = value.internal


text:
module Test2D

using Test
using Trixi2Vtk

include("test_trixi2vtk.jl")

# Start with a clean environment: remove Trixi output directory if it exists
outdir = "out"
isdir(outdir) && rm(outdir, recursive=true)

# Create empty artifacts directory where all files that should be preserved will be stored
artifacts_dir = joinpath(pathof(Trixi2Vtk) |> dirname |> dirname, "artifacts")
isdir(artifacts_dir) && rm(artifacts_dir, recursive=true)
mkdir(artifacts_dir)

@testset "2D" begin
  run_trixi(joinpath("2d", "elixir_advection_extended.jl"), maxiters=1)

  @testset "uniform mesh" begin
    test_trixi2vtk("solution_000000.h5", outdir,
        hashes=[("solution_000000.vtu", "1ec2c93c0c9c4f4992dea54afaf2a348ece0160e"),
                ("solution_000000_celldata.vtu", "e396c3ba63276347966d4264cf0f52d592221830")])

    # Store output files as artifacts to facilitate debugging of failing tests
    outfiles = ("solution_000000.vtu", "solution_000000_celldata.vtu")
    testname = "uniform-mesh"
    for outfile in outfiles
      println("Copying '", abspath(joinpath(outdir, outfile)),
              "' to '", abspath(joinpath(artifacts_dir, testname * "-" * outfile)),
              "'...")
      cp(joinpath(outdir, outfile), joinpath(artifacts_dir, testname * "-" * outfile))
    end
  end

  @testset "uniform mesh with vti output" begin
    if Sys.isapple()
      # This test fails on MacOS due to differing binary VTI files (even though the contents match)
      test_trixi2vtk("restart_000001.h5", outdir,
          format=:vti)
    else
      test_trixi2vtk("restart_000001.h5", outdir,
          hashes=[("restart_000001.vti", "664f25ab018a373774b5aad69ad3f2f5a3b21649"),
                  ("restart_000001_celldata.vtu", "e396c3ba63276347966d4264cf0f52d592221830")],
          format=:vti)
    end

    # Store output files as artifacts to facilitate debugging of failing tests
    outfiles = ("restart_000001.vti", "restart_000001_celldata.vtu")
    testname = "uniform-mesh-with-vti-output"
    for outfile in outfiles
      println("Copying '", abspath(joinpath(outdir, outfile)),
              "' to '", abspath(joinpath(artifacts_dir, testname * "-" * outfile)),
              "'...")
      cp(joinpath(outdir, outfile), joinpath(artifacts_dir, testname * "-" * outfile))
    end
  end
end

# Clean up afterwards: delete Trixi output directory
@test_nowarn rm(outdir, recursive=true)

end


text:
@pure function _find(names::Tuple{Vararg{Symbol}}, name::Symbol)
    for i in 1:length(names)
        if names[i] === name
            return i
        end
    end
    error("Cannot find $name in $names")
end

@pure function _issubset(a::Tuple{Vararg{Symbol}}, b::Tuple{Vararg{Symbol}})
    for sa ∈ a
        if sa ∉ b
            return false
        end
    end
    return true
end

@pure function _issetequal(a::NTuple{N, Symbol}, b::NTuple{N, Symbol}) where {N}
    for sa ∈ a
        if sa ∉ b
            return false
        end
    end
    return true
end
@pure _issetequal(::Tuple{Vararg{Symbol}}, ::Tuple{Vararg{Symbol}}) = false

@pure function _headsubset(a::Tuple{Vararg{Symbol}}, b::Tuple{Vararg{Symbol}})
    out = ()
    for s ∈ a
        if s ∈ b
            out = (out..., s)
        else
            return out
        end
    end
    return out
end

@generated function _headsubset(::Val{a}, ::Val{b}) where {a, b}
    out = ()
    for s ∈ a
        if s ∈ b
            out = (out..., s)
        else
            return out
        end
    end
    return out
end

@pure function _intersect(a::Tuple{Vararg{Symbol}}, b::Tuple{Vararg{Symbol}})
    out = ()
    for s ∈ a
        if s ∈ b
            out = (out..., s)
        end
    end
    return out
end

@generated function _intersect(::Val{a}, ::Val{b}) where {a, b}
    out = ()
    for s ∈ a
        if s ∈ b
            out = (out..., s)
        end
    end
    return out
end

@pure function _setdiff(a::Tuple{Vararg{Symbol}}, b::Tuple{Vararg{Symbol}})
    out = ()
    for s ∈ a
        if s ∉ b
            out = (out..., s)
        end
    end
    return out
end

@generated function _eltypes(a::Tuple{Vararg{AbstractVector}})
    Ts = []
    for V in a.parameters
        push!(Ts, eltype(V))
    end
    return Tuple{Ts...}
end

@generated function _makevectors(::Type{Ts}, dims::Tuple{Int}) where {Ts <: Tuple}
    exprs = [:(Vector{$T}(undef, dims)) for T ∈ Ts.parameters]
    return quote
        @_inline_meta
        return tuple($(exprs...))
    end
end

# index of the last value of vector a that is less than but not equal to x;
# returns 0 if x is less than all values of v.
function searchsortedlastless(v::AbstractVector, x, lo::Int = first(keys(v)), hi::Int = last(keys(v)))
    lo = lo-1
    hi = hi+1
    @inbounds while lo < hi-1
        m = (lo+hi)>>>1
        y = v[m]
        if isless(y, x)
            lo = m
        else
            hi = m
        end
    end
    return lo
end

function searchsortedfirstgreater(v::AbstractVector, x, lo::Int = first(keys(v)), hi::Int = last(keys(v)))
    lo = lo-1
    hi = hi+1
    @inbounds while lo < hi-1
        m = (lo+hi)>>>1
        if isless(x, v[m])
            hi = m
        else
            lo = m
        end
    end
    return hi
end

@pure _add(n::Int, m::Int) = n + m

_all(f, ::Tuple{}, ::Tuple{}) = true
@generated function _all(f, t1::NTuple{n, Any}, t2::NTuple{n, Any}) where {n}
    exprs = [:(f(t1[$i], t2[$i])) for i = 1:n]
    out = reduce((final, expr) -> :($final && $expr), exprs) # short-circuiting...
    return quote
        @_inline_meta
        return $out
    end
end

_valuetype(::Type{NamedTuple{names, T}}) where {names, T} = T
_valuetype(::NamedTuple{names, T}) where {names, T} = T

@generated function _cat_types(::Type{T1}, ::Type{T2}) where {T1 <: Tuple, T2 <: Tuple}
    return Tuple{T1.parameters..., T2.parameters...}
end


text:
module NLP1

using JuMP, BARON, Test

m = Model(BARON.Optimizer)
ub = [6,4]
@variable(m, 0 ≤ x[i=1:2] ≤ ub[i])

@constraint(m, x[1]*x[2] ≤ 4)

@objective(m, Min, -x[1] - x[2])

optimize!(m)

@test isapprox(value(x[1]), 6, rtol=1e-6)
@test isapprox(value(x[2]), 2/3, rtol=2e-6)
@test isapprox(objective_value(m), -20/3, rtol=1e-6)

end # module


text:
saivdr_permdctmtx(sz::Integer...) = saivdr_permdctmtx(sz)

saivdr_permdctmtx(sz::NTuple{1}) = permdctmtx(sz)

function saivdr_permdctmtx(sz::NTuple{2})
    T = Float64
    mtx = MDCDL.representationmatrix(x->dct(T.(x)), sz)

    # cls[1,1] = 1
    # cls[2,1] = 3
    # cls[1,2] = 4
    # cls[2,2] = 2
    cls = [ 1 4; 3 2 ]
    indexorders = map(CartesianIndices(sz)) do ci
        cls[(mod.(ci.I .- 1, 2) .+ 1)...]
    end
    permids = sortperm(vec(indexorders); alg=Base.DEFAULT_STABLE)

    @views vcat([ transpose(mtx[pi,:]) for pi in permids ]...)
end

# function saivdr_permdctmtx(sz::NTuple{3})
#     T = Float64
#     mtx = MDCDL.representationmatrix(x->dct(T.(x)), sz)
#
#     # cls[1,1,1] = 1
#     # cls[2,1,1] = 8
#     # cls[1,2,1] = 6
#     # cls[2,2,1] = 3
#     # cls[1,1,2] = 5
#     # cls[2,1,2] = 4
#     # cls[1,2,2] = 2
#     # cls[2,2,2] = 7
#     cls = cat([1 6; 8 3], [5 2; 4 7], dims=3)
#     indexorders = map(CartesianIndices(sz)) do ci
#         cls[(mod.(ci.I .- 1, 2) .+ 1)...]
#     end
#     permids = sortperm(vec(indexorders); alg=Base.DEFAULT_STABLE)
#
#     @views vcat([ transpose(mtx[pi,:]) for pi in permids ]...)
# end

function analyze_compatible(nsolt::AbstractNsolt{T,2}, px::AbstractMatrix, nBlocks::NTuple; kwargs...) where {T}
    ty = initialStep(nsolt, px; kwargs...)

    nShifts = fld.(size(px, 2), nBlocks)
    rotatedimsfcns = (t->rotatedimspv(t, nBlocks[1]), copy)
    exty = extendAtoms(nsolt, ty, nShifts, rotatedimsfcns; kwargs...)
    exty = rotatedimspv(exty, nBlocks[2])
    return shiftFilterSymmetry(nsolt, exty)
end

# function analyze_compatible(nsolt::AbstractNsolt{T,3}, px::AbstractMatrix, nBlocks::NTuple; kwargs...) where {T}
#     ty = initialStep(nsolt, px; kwargs...)
#
#     nShifts = fld.(size(px, 2), nBlocks)
#     rotatedimsfcns = (t->rotatedimspv(t, nShifts[1]), t->rotatedimspv(t, nShifts[2], copy)
#     exty = extendAtoms(nsolt, ty, nShifts, rotatedimsfcns; kwargs...)
#     exty = rotatedimspv(exty, nBlocks[3])
#     return shiftFilterSymmetry(nsolt, exty)
# end

max_stars_repo_name: nixir/MDCDL

max_issues_repo_name: nixir/MDCDL

max_forks_repo_name: nixir/MDCDL


text:
const JLPATH_PREFIX = "jlpkg:"

function c_strdup(str)
	ptr = @ccall libngspice.tmalloc((sizeof(str)+1)::Cint)::Ptr{UInt8}
	@GC.preserve str Base.unsafe_copyto!(ptr, pointer(str), sizeof(str))
	unsafe_store!(ptr, UInt8(0), sizeof(str) + 1)
	ptr
end

function path_resolve(name::Cstring, dir::Cstring)::Cstring
    name = unsafe_string(name)
    dir = unsafe_string(dir)

    if startswith(name, JLPATH_PREFIX)
        path = name[sizeof(JLPATH_PREFIX)+1:end]
        components = splitpath(path)

        @assert components[1] != "/"

        pkg_path = Base.locate_package(Base.PkgId(components[1]))
        if pkg_path !== nothing
            pkg_path = realpath(joinpath(dirname(pkg_path), ".."))
            file_path = joinpath(pkg_path, components[2:end]...)
            return c_strdup(file_path)
        end

        return Ptr{Cchar}(-1)
    end
    return C_NULL
end

max_stars_repo_name: JuliaComputing/Ngspice

max_issues_repo_name: JuliaComputing/Ngspice

max_forks_repo_name: JuliaComputing/Ngspice


text:
# # 210: 2D Nonlinear Poisson with reaction
# ([source code](SOURCE_URL))

module Example210_NonlinearPoisson2D_Reaction

using Printf
using VoronoiFVM
using ExtendableGrids
using GridVisualize

function main(;n=10,Plotter=nothing,verbose=false, unknown_storage=:sparse)

    h=1.0/convert(Float64,n)
    X=collect(0.0:h:1.0)
    Y=collect(0.0:h:1.0)

    grid=VoronoiFVM.Grid(X,Y)
    data = (eps=1.0e-2, k=1.0)

    function reaction!(f,u,node,data)
        f[1]=data.k*(u[1]-u[2])
        f[2]=data.k*(u[2]-u[1])
    end

    function flux!(f,u0,edge,data)
        u=unknowns(edge,u0)
        f[1]=data.eps*(u[1,1]-u[1,2])
        f[2]=data.eps*(u[2,1]-u[2,2])
    end

    function source!(f,node,data)
        x1=node[1]-0.5
        x2=node[2]-0.5
        f[1]=exp(-20*(x1^2+x2^2))
    end

    function storage!(f,u,node,data)
        f[1]=u[1]
        f[2]=u[2]
    end

    physics=VoronoiFVM.Physics(num_species=2,
                               data=data,
                               flux=flux!,
                               storage=storage!,
                               reaction=reaction!,
                               source=source!)

    sys=VoronoiFVM.System(grid,physics,unknown_storage=unknown_storage)

    enable_species!(sys,1,[1])
    enable_species!(sys,2,[1])

    inival=unknowns(sys)
    U=unknowns(sys)
    inival.=0.0

    control=VoronoiFVM.NewtonControl()
    control.verbose=verbose
    control.tol_linear=1.0e-5
    control.max_lureuse=0
    tstep=0.01
    time=0.0
    istep=0
    u15=0
    p=GridVisualizer(Plotter=Plotter,layout=(2,1))
    while time<1
        time=time+tstep
        solve!(U,inival,sys,control=control,tstep=tstep)
        inival.=U
        if verbose
            @printf("time=%g\n",time)
        end
        u15=U[15]
        tstep*=1.0
        istep=istep+1
        scalarplot!(p[1,1],grid,U[1,:],clear=true)
        scalarplot!(p[2,1],grid,U[2,:],show=true)
    end
    return u15
end

function test()
    main(unknown_storage=:sparse) ≈ 0.014566189535134827 &&
        main(unknown_storage=:dense) ≈ 0.014566189535134827
end
end


text:
using GZip

import Base: eof, close, position
export BamReader, close, value, eof, advance!, eachposition

mutable struct BamReader
    bamStream
    readOrientation #useReverseReads::Bool
    done::Bool
    position::Int64
    contigs::ReferenceContigs
end

function BamReader(bamFileName::String, readOrientation, contigs)
    f = GZip.open(bamFileName)

    # make sure this is a BAM file
    code = read!(f, Array{UInt8}(undef, 4))
    @assert code == b"BAM\1"

    # get through the header data
    l_text = read(f, Int32)
    skip(f, l_text)

    # make sure the contigs match our reference
    n_ref = read(f, Int32)
    if !(n_ref == contigs.count)
        println("Your bam file is not aligned to the UCSC hg38 genome.")
        println("See the step 3.1 at https://github.com/hiranumn/AIControl.jl to realign your genome")
        println("to the specific version of hg38 using bowtie2.")
        exit()
    end

    for j in 1:n_ref
        l_name = read(f, Int32)
        refName = String(read(f, Array{UInt8}(undef, l_name))[1:end-1]) # ignore the null terminator
        l_ref = read(f, Int32)
        if !(l_ref == contigs.sizes[j]) || !(refName == contigs.names[j])
            println("Your bam file is not aligned to the UCSC hg38 genome.")
            println("See the step 3.1 at https://github.com/hiranumn/AIControl.jl to realign your genome")
            println("to the specific version of hg38 using bowtie2.")
            exit()
        end
    end

    r = BamReader(f, readOrientation, false, 1, contigs)
    advance!(r)
    r
end
close(reader::BamReader) = GZip.close(reader.bamStream)
value(reader::BamReader) = 1
position(reader::BamReader) = reader.position
eof(reader::BamReader) = reader.position == -1

function advance!(r::BamReader)
    f = r.bamStream
    while !r.done
        if peek(f) == -1 # eof does not work on the BAM files either in C++ or here (BAM vs. gzip issue?)
            r.done = true
            r.position = -1
            return
        end

        buf = Array{Int32}(undef, 5) # [block_size, refID, pos, bin_mq_nl, flag_nc]
        gzread(f, pointer(buf), 20)
        block_size = buf[1]
        refID = buf[2] + 1 # the reference contig this read maps to

        # get the read position
        if refID != 0
            r.position = buf[3] + r.contigs.offsets[refID] + 1; # convert to 1 based indexing
        end

        forward = (buf[5] & 1048576) == 0 # see if we are reverse complemented
        skip(f, block_size-16) # skip the rest of the entry

        # break if we found a read in the right direction
        if refID != 0 && (r.readOrientation == :any || (forward && r.readOrientation == :forward) || (!forward && r.readOrientation == :reverse))
            return
        end
    end
end

# here we want to update the reader
#eachposition(r::BamReader) = BamReaderIterator(r)
#struct BamReaderIterator
#	reader::BamReader
#end
#Base.start(it::BamReaderIterator) = it.reader.position
#Base.done(it::BamReaderIterator, position) = position == -1
#function Base.next(it::BamReaderIterator, position)
#	pos = it.reader.position
#	advance!(it.reader)
#	pos,it.reader.position
#end