text:
#=
UTF16Str and UCS2Str types (UTF-16 encoding and pure BMP UCS-2)

Copyright 2017-2018 Gandalf Software, Inc., Scott P. Jones,
and other contributors to the Julia language
Licensed under MIT License, see LICENSE.md
Based in (small) part on code for UTF16String that used to be in Julia
=#

const _trail_mask = 0xdc00_dc00_dc00_dc00
const _hi_bit_16  = 0x8000_8000_8000_8000

@inline _mask_surr(v)  = xor((v | v<<1 | v<<2 | v<<3 | v<<4 | v<<5) & _hi_bit_16, _hi_bit_16)
@inline _get_masked(v::UInt64) = _mask_surr(xor(v, _trail_mask))
@inline _get_masked(qpnt::Ptr{UInt64}) = _get_masked(unsafe_load(qpnt))
@inline _get_lead(qpnt) = xor(_get_masked(qpnt), _hi_bit_16)

@inline get_utf16(ch) = (0xd7c0 + (ch >> 10))%UInt16, (0xdc00 + (ch & 0x3ff))%UInt16

@inline function _align_len_utf16(pnt, cnt, v)
    len = 0
    fin = pnt + cnt
    while (pnt += CHUNKSZ) < fin
        len += count_ones(v)
        v = _get_lead(pnt)
    end
    len + count_ones((cnt & CHUNKMSK) == 0 ? v : (v & _mask_bytes(cnt)))
end

_length_al(::CodeUnitMulti, ::Type{UTF16CSE}, beg::Ptr{UInt16}, cnt::Int) =
    (pnt = reinterpret(Ptr{UInt64}, beg); _align_len_utf16(pnt, cnt<<1, _get_lead(pnt)))

function _length(::CodeUnitMulti, ::Type{UTF16CSE}, beg::Ptr{UInt16}, cnt::Int)
    align = reinterpret(UInt, beg)
    pnt = reinterpret(Ptr{UInt64}, align & ~CHUNKMSK)
    v = _get_lead(pnt)
    if (align &= CHUNKMSK) != 0
        msk = _mask_bytes(align)
        v = (v & ~msk) | (msk & _trail_mask)
        cnt += (align>>>1)
    end
    _align_len_utf16(pnt, cnt<<1, v)
end

function _nextind(::CodeUnitMulti, str::MS_UTF16, pos::Int, nchar::Int)
    nchar < 0 && ncharerr(nchar)
    siz = ncodeunits(str)
    @boundscheck 0 <= pos <= siz || boundserr(str, pos)
    @preserve str begin
        beg = pointer(str)
        pnt = bytoff(beg, pos - 1)
        fin = bytoff(beg, siz)
        cu = get_codeunit(pnt)
        nchar == 0 && (is_surrogate_trail(cu) ? index_error(str, pos) : return pos)
        is_surrogate_trail(cu) && (pnt += 2)
        # pnt should now point to a valid start of a character
        # This could be sped up, by looking at chunks, and if all ASCII (common case),
        # simply move forward 8
        while (nchar -= 1) > 0 && pnt < fin
            pnt += ifelse(is_surrogate_lead(get_codeunit(pnt)), 4, 2)
        end
        chrdiff(pnt, beg)
    end
end

function _prevind(::CodeUnitMulti, str::MS_UTF16, pos::Int, nchar::Int)
    nchar < 0 && ncharerr(nchar)
    @boundscheck 0 < pos <= ncodeunits(str)+1 || boundserr(str, pos)
    @preserve str begin
        beg = pointer(str)
        pnt = bytoff(beg, pos - 1)
        nchar == 0 && (is_surrogate_trail(get_codeunit(pnt)) ? index_error(str, pos) : return pos)
        # This could be sped up, by looking at chunks, and if all ASCII (common case),
        # simply move back 8
        while pnt >= beg
            pnt -= ifelse(is_surrogate_trail(get_codeunit(pnt)), 4, 2)
            (nchar -= 1) > 0 || break
        end
        chrdiff(pnt + 2, beg)
    end
end

# Check for any surrogate characters
function is_bmp(str::MS_UTF16)
    (siz = sizeof(str)) == 0 && return true
    # Todo: handle unaligned for ARM32
    @preserve str begin
        siz < CHUNKSZ && return (_get_masked(_pnt64(str)) & _mask_bytes(siz)) == 0

        pnt, fin = _calcpnt(str, siz)
        while (pnt += CHUNKSZ) <= fin
            _get_masked(pnt) == 0 || return false
        end
        pnt - CHUNKSZ == fin || (_get_masked(pnt) & _mask_bytes(siz)) == 0
    end
end

@inline function _check_bmp_utf16_al(pnt, cnt, v)
    fin = pnt + cnt
    v = _get_masked(v)
    while (pnt += CHUNKSZ) < fin
        v == 0 || return false
        v = _get_masked(pnt)
    end
    ((cnt & CHUNKMSK) == 0 ? v : (v & _mask_bytes(cnt))) == 0
end
@inline _check_bmp_utf16_al(pnt, cnt) = _check_bmp_utf16_al(pnt, cnt, unsafe_load(pnt))

@inline function _check_bmp_utf16_ul(beg, cnt)
    align = reinterpret(UInt, beg)
    pnt = reinterpret(Ptr{UInt64}, align & ~CHUNKMSK)
    v = unsafe_load(pnt)
    if (align &= CHUNKMSK) != 0
        v &= ~_mask_bytes(align)
        cnt += align
    end
    _check_bmp_utf16_al(pnt, cnt, v)
end

is_bmp(str::Str{UTF16CSE}) =
    (cnt = sizeof(str)) == 0 ? true :
    @preserve str _check_bmp_utf16_al(reinterpret(Ptr{UInt64}, pointer(str)), cnt)

is_bmp(str::SubString{<:Str{UTF16CSE}}) =
    (cnt = sizeof(str)) == 0 ? true : @preserve str _check_bmp_utf16_ul(pointer(str), cnt)

is_bmp(str::MaybeSub{<:Str{<:UCS2_CSEs}}) = true

is_ascii(str::Str{_UCS2CSE}) = false
is_latin(str::Str{_UCS2CSE}) = false

is_ascii(str::SubString{<:Str{_UCS2CSE}}) = is_ascii(Str{UCS2CSE}(str.data))
is_latin(str::SubString{<:Str{_UCS2CSE}}) = is_latin(Str{UCS2CSE}(str.data))

# Speed this up accessing 64 bits at a time
@propagate_inbounds function _cnt_non_bmp(len, pnt::Ptr{UInt16})
    cnt = 0
    fin = bytoff(pnt, len)
    while pnt < fin
        cnt += is_surrogate_lead(get_codeunit(pnt))
        pnt += 2
    end
    cnt
end

@inline _lastindex(::CodeUnitMulti, str::MS_UTF16) =
    ((len = ncodeunits(str)) != 0
     ? (is_surrogate_codeunit(get_codeunit(pointer(str), len)) ? len-1 : len) : 0)

get_supplementary(lead::Unsigned, trail::Unsigned) = (UInt32(lead-0xd7f7)<<10 + trail)

function _nextcpfun(::CodeUnitMulti, ::Type{UTF16CSE}, pnt)
    ch = get_codeunit(pnt)
    (is_surrogate_lead(ch)
     ? (get_supplementary(ch, get_codeunit(pnt + 2)), pnt + 4)
     : (ch%UInt32, pnt + 2))
end

@propagate_inbounds function _next(::CodeUnitMulti, T, str::MS_UTF16, pos::Int)
    @boundscheck pos <= ncodeunits(str) || boundserr(str, pos)
    @preserve str begin
        pnt = bytoff(pointer(str), pos)
        ch = get_codeunit(pnt - 2)
        (is_surrogate_lead(ch)
         ? (T(get_supplementary(ch, get_codeunit(pnt))), pos + 2)
         : (T(ch), pos + 1))
    end
end

@inline _thisind(::CodeUnitMulti, str::MS_UTF16, len, pnt, pos) =
    Int(pos) - is_surrogate_trail(get_codeunit(pnt, pos))

@propagate_inbounds @inline function _nextind(::CodeUnitMulti, str::MS_UTF16, pos::Int)
    pos == 0 && return 1
    @boundscheck 1 <= pos <= ncodeunits(str) || boundserr(str, pos)
    @preserve str pos + 1 + is_surrogate_lead(get_codeunit(pointer(str), pos))
end

@propagate_inbounds @inline function _prevind(::CodeUnitMulti, str::MS_UTF16, pos::Int)
    (pos -= 1) == 0 && return 0
    numcu = ncodeunits(str)
    @boundscheck 0 < pos <= numcu || boundserr(str, pos + 1)
    @preserve str pos - is_surrogate_trail(get_codeunit(pointer(str), pos))
end

function _reverse(::CodeUnitMulti, ::Type{UTF16CSE}, len, pnt::Ptr{T}) where {T<:CodeUnitTypes}
    buf, beg = _allocate(UInt16, len)
    out = bytoff(beg, len)
    while out > beg
        ch = get_codeunit(pnt)
        is_surrogate_lead(ch) && set_codeunit!(out -= 2, get_codeunit(pnt += 2))
        set_codeunit!(out -= 2, ch)
        pnt += 2
    end
    Str(UTF16CSE, buf)
end

@inline _isvalid_char_pos(::CodeUnitMulti, ::Type{UTF16CSE}, str, pos::Integer) =
    @preserve str !is_surrogate_trail(get_codeunit(pointer(str), pos))

function is_valid(::Type{<:Str{<:UCS2_CSEs}}, data::AbstractArray{UInt16})
    @inbounds for ch in data
        is_surrogate_codeunit(ch) && return false
    end
    true
end

function is_valid(::Type{<:Str{<:UCS2_CSEs}}, pnt::Ptr{UInt16}, len)
    fin = bytoff(pnt, len)
    while pnt < fin # check for surrogates
        is_surrogate_codeunit(get_codeunit(pnt)) && return false
        pnt += 2
    end
    true
end

# These need to change to output directly, not converted to UTF-8

# Single character conversion

@inline function _convert_utf_n(::Type{C}, ch::UInt32) where {C<:UTF16CSE}
    buf, pnt = _allocate(UInt16, 2)
    # output surrogate pair
    c1, c2 = get_utf16(ch)
    set_codeunit!(pnt,     c1)
    set_codeunit!(pnt + 1, c2)
    Str(C, buf)
end

convert(::Type{<:Str{C}}, ch::Unsigned) where {C<:UTF16CSE} =
    (is_unicode(ch)
     ? (ch <= 0xffff ? _convert(UTF16CSE, ch%UInt16) : _convert_utf_n(C, ch%UInt32))
     : unierror(UTF_ERR_INVALID, 0, ch))

# Type _UCS2CSE must have at least 1 character > 0xff, so use other types as well
convert(::Type{<:Str{_UCS2CSE}}, ch::Unsigned) =
    (ch > 0xff
     ? (is_bmp(ch) ? _convert(_UCS2CSE, ch%UInt16) : unierror(UTF_ERR_INVALID, 0, ch))
     : _convert(_LatinCSE, ch%UInt8))

function convert(::Type{<:Str{C}}, str::AbstractString) where {C<:UCS2_CSEs}
    is_empty(str) && return empty_str(C)
    # Might want to have an invalids_as argument
    len, flags, num4byte, num3byte = unsafe_check_string(str)
    num4byte == 0 || unierror(UTF_ERR_INVALID_UCS2)
    buf, pnt = _allocate(UInt16, len)
    @inbounds for ch in str
        set_codeunit!(pnt, ch%UInt16)
        pnt += 2
    end
    Str((C === _UCS2CSE && num3byte != 0) ? _UCS2CSE : UCS2CSE, buf)
end

function convert(::Type{<:Str{C}}, str::MS_ByteStr) where {C<:UCS2_CSEs}
    # Might want to have an invalids_as argument
    # handle zero length string quickly
    (siz = sizeof(str)) == 0 && return empty_str(C)
    @preserve str begin
        pnt = pointer(str)
        # Check that is correct UTF-8 encoding and get number of words needed
        len, flags, num4byte, num3byte = fast_check_string(pnt, siz)
        num4byte == 0 || unierror(UTF_ERR_INVALID_UCS2)
        # Optimize case where no characters > 0x7f
        Str((C === _UCS2CSE && num3byte != 0) ? _UCS2CSE : UCS2CSE,
            flags == 0 ? _cvtsize(UInt16, pnt, len) : _encode_utf16(pnt, len))
    end
end

# handle zero length string quickly, just widen these
convert(::Type{<:Str{C}},
        str::MaybeSub{<:Str{Union{ASCIICSE, Latin_CSEs}}}) where {C<:UCS2_CSEs} =
    (siz = sizeof(str)) == 0 ? empty_str(C) : Str(C, _cvtsize(UInt16, pointer(str), siz))

function convert(::Type{<:Str{C}}, str::MS_UTF16) where {C<:UCS2_CSEs}
    # Might want to have an invalids_as argument
    # handle zero length string quickly
    (siz = sizeof(str)) == 0 && return empty_str(C)
    # Check if conversion is valid
    is_bmp(str) || unierror(UTF_ERR_INVALID_UCS2)
    @preserve str Str(C, _cvtsize(UInt16, pointer(str), len))
end

function is_valid(::Type{<:Str{UTF16CSE}}, data::AbstractArray{UInt16})
    (len = length(data) - 1) < 0 && return true
    pos = 0
    @inbounds while pos < len # check for unpaired surrogates
        ch = data[pos += 1]
        is_surrogate_codeunit(ch) && is_surrogate_lead(ch) &&
            is_surrogate_trail(data[pos += 1]) ||
            return false
    end
    @inbounds return pos > len || !is_surrogate_codeunit(get_codeunit(data[pos + 1]))
end

# This can be sped up, to check 4 words at a time, only checking for unpaired
# or out of order surrogates when one is found in the UInt64
function is_valid(::Type{<:Str{UTF16CSE}}, pnt::Ptr{UInt16}, len)
    len == 0 && return true
    fin = bytoff(pnt, len - 1)
    while pnt < fin # check for unpaired surrogates
        ch = get_codeunit(pnt)
        is_surrogate_codeunit(ch) && is_surrogate_lead(ch) &&
            is_surrogate_trail(get_codeunit(pnt += 2)) ||
            return false
        pnt += 2
    end
    pnt > fin || !is_surrogate_codeunit(get_codeunit(pnt))
end

function convert(::Type{<:Str{UTF16CSE}}, str::AbstractString)
    is_empty(str) && return empty_utf16
    len, flags, num4byte = unsafe_check_string(str)
    buf, pnt = _allocate(UInt16, len + num4byte)
    @inbounds for ch in str
        c = ch%UInt32
        if c > 0x0ffff
            # output surrogate pair
            w1, w2 = get_utf16(c)
            set_codeunit!(pnt, w1)
            set_codeunit!(pnt+2, w2)
            pnt += 4
        else
            set_codeunit!(pnt, c%UInt16)
            pnt += 2
        end
    end
    Str(UTF16CSE, buf)
end

function convert(::Type{<:Str{UTF16CSE}}, str::MS_ByteStr)
    # handle zero length string quickly
    (siz = sizeof(str)) == 0 && return empty_utf16
    @preserve str begin
        pnt = pointer(str)
        # Check that is correct UTF-8 encoding and get number of words needed
        len, flags, num4byte = fast_check_string(pnt, siz)
        # Optimize case where no characters > 0x7f
        Str(UTF16CSE, flags == 0 ? _cvtsize(UInt16, pnt, len) : _encode_utf16(pnt, len + num4byte))
    end
end

function convert(::Type{<:Str{UTF16CSE}}, str::MS_UTF8)
    # handle zero length string quickly
    (siz = sizeof(str)) == 0 && return empty_utf16
    @preserve str begin
        pnt = pointer(str)
        len, flags, num4byte = count_chars(UTF8Str, pnt, siz)
        # Optimize case where no characters > 0x7f
        Str(UTF16CSE, flags == 0 ? _cvtsize(UInt16, pnt, len) : _encode_utf16(pnt, len + num4byte))
    end
end

"""
Converts an already validated UTF-8 encoded vector of `UInt8` to a `UTF16Str`

Input Arguments:

*   `pnt` `Ptr{UInt8}` of UTF-8 encoded data
*   `len` length of output in 16-bit words

Returns:

*   `UTF16Str`
"""
function _encode_utf16(pnt::Ptr{UInt8}, len)
    buf, out = _allocate(UInt16, len)
    fin = out + (len<<1)
    @inbounds while out < fin
        ch = get_codeunit(pnt)%UInt16
        # Handle ASCII characters
        if ch <= 0x7f
        # Handle range 0x80-0x7ff
        elseif ch < 0xe0
            ch = ((ch & 0x1f) << 6) | (get_codeunit(pnt += 1) & 0x3f)
        # Handle range 0x800-0xffff
        elseif ch < 0xf0
            ch = get_utf8_3byte(pnt += 2, ch)
        # Handle range 0x10000-0x10ffff
        else
            ch32 = get_utf8_4byte(pnt += 3, ch)
            # output surrogate pair
            w1, ch = get_utf16(ch32)
            set_codeunit!(out, w1)
            out += 2
        end
        set_codeunit!(out, ch)
        out += 2
        pnt += 1
    end
    buf
end

_encode_utf16(dat::Vector{UInt8}, len) = @preserve dat _encode_utf16(pointer(dat), len)
_encode_utf16(str::MS_ByteStr, len)    = @preserve str _encode_utf16(pointer(str), len)

@inline _cvt_16_to_utf8(::Type{<:Str{UTF16CSE}}, pnt, len)     = _transcode_utf8(pnt, len)
@inline _cvt_16_to_utf8(::Type{<:Str{<:UCS2_CSEs}}, pnt, len)  = _encode_utf8(pnt, len)
@inline _cvt_16_to_utf8(::Type{<:Str{<:UTF32_CSEs}}, pnt, len) = _encode_utf8(pnt, len)

function _cvt_utf8(::Type{T}, str::S) where {T<:Union{String, Str{UTF8CSE}}, S}
    # handle zero length string quickly
    (len = ncodeunits(str)) == 0 && return empty_str(T)
    @preserve str begin
        # get number of bytes to allocate (use faster count for validated strings)
        pnt = pointer(str)
        len, flags, num4byte, num3byte, num2byte, latin1 = count_chars(S, pnt, len)
        Str(cse(T),
            (flags == 0
             ? _cvtsize(UInt8, pnt, len)
             : _cvt_16_to_utf8(S, pnt, len + latin1 + num2byte + num3byte*2 + num4byte*3)))
    end
end

# Split this way to avoid ambiguity errors
convert(::Type{String},  str::MaybeSub{<:Str{Union{Word_CSEs,Quad_CSEs}}}) =
    _cvt_utf8(String, str)
convert(::Type{<:Str{UTF8CSE}}, str::MaybeSub{<:Str{Union{Word_CSEs,Quad_CSEs}}}) =
    _cvt_utf8(UTF8Str, str)

"""
Converts an already validated UTF-32 encoded vector of `UInt32` to a `UTF16Str`

Input Arguments:

*   `dat` `Vector{UInt32}` of UTF-32 encoded data
*   `len` length of output in 16-bit words

Returns:

*   `Vector{UInt8}`
"""
function _encode_utf16(dat::Ptr{UInt32}, len)
    buf, pnt = _allocate(UInt16, len)
    fin = pnt + (len<<1)
    @inbounds while pnt < fin
        ch = get_codeunit(dat)
        dat += 4
        if ch > 0x0ffff
            # Output surrogate pair for 0x10000-0x10ffff
            w1, w2 = get_utf16(ch)
            set_codeunit!(pnt,   w1)
            set_codeunit!(pnt+2, w2)
            pnt += 4
        else
            set_codeunit!(pnt, ch%UInt16)
            pnt += 2
        end
    end
    buf
end

# Copies because not safe to expose the internal array (would allow mutation)
function convert(::Type{Vector{UInt16}}, str::Str{<:Word_CSEs})
    len = ncodeunits(str)
    vec = create_vector(UInt16, len)
    @inbounds unsafe_copyto!(pointer(vec), pointer(str), len)
    vec
end

# Todo: Some of these need to be fixed to account for SubStr, when that is added
convert(::Type{T},  str::MaybeSub{T}) where {T<:Str{<:Union{UCS2_CSEs, UTF32_CSEs}}} = str
convert(::Type{<:Str{UTF16CSE}}, str::MaybeSub{<:Str{<:UCS2_CSEs}}) = Str(UTF16CSE, str.data)

unsafe_convert(::Type{Ptr{UInt16}}, s::MS_UTF16) = pointer(s)

function convert(::Type{<:Str{UTF16CSE}}, dat::AbstractArray{UInt16})
    is_empty(dat) && return empty_utf16
    len, flags, num4byte = unsafe_check_string(dat, 1, lastindex(dat))
    # Optimize case where no surrogate characters
    Str(UTF16CSE, flags == 0 ? _cvtsize(UInt16, dat, len) : _encode_utf16(dat, len + num4byte))
end

function convert(::Type{<:Str{UTF16CSE}}, dat::Vector{UInt16})
    (len = length(dat)) == 0 && return empty_utf16
    @preserve dat fast_check_string(pointer(dat), len)
    Str(UTF16CSE, _copysub(dat))
end

function convert(::Type{<:Str{UTF16CSE}}, str::MaybeSub{<:Str{Text2CSE}})
    (len = ncodeunits(str)) == 0 && return empty_utf16
    @preserve str fast_check_string(pointer(str), len)
    Str(UTF16CSE, _copysub(str))
end

convert(::Type{<:Str{Text2CSE}},
        str::MaybeSub{<:Str{C}}) where {C<:Union{ASCIICSE,Text1CSE,BinaryCSE,Latin_CSEs}} =
            @preserve str Str(Text2CSE, _cvtsize(UInt16, pointer(str), ncodeunits(str)))

function convert(::Type{<:Str{C}},
                 str::MaybeSub{<:Str{T}}) where {C<:UCS2_CSEs,T<:Union{Text2CSE,Text4CSE}}
    (len = ncodeunits(str)) == 0 && return C === _UCS2CSE ? empty_ascii : empty_ucs2
    @preserve str begin
        pnt = pointer(str)
        len, flags, num4byte, num3byte, num2byte, latin1byte = fast_check_string(pnt, len)
        num4byte == 0 || unierror(UTF_ERR_INVALID_UCS2)
        (C === UCS2CSE || (num2byte+num3byte) != 0) &&
            return Str(C, T===Text2CSE ? _copysub(str) : _cvtsize(UInt16, pnt, len))
        Str(latin1byte == 0 ? ASCIICSE : _LatinCSE, _cvtsize(UInt8, pnt, len))
    end
end

_convert(pnt::Ptr{T}, len, T1) where {T<:Union{UInt16,UInt16_U,UInt16_S,UInt16_US}} =
    ((ch = unsafe_load(pnt)) == 0xfffe
     ? _convert(reinterpret(Ptr{T1}, pnt + 2), len - 1)
     : (ch == 0xfeff ? _convert(pnt + 2, len - 1) : _convert(pnt, len)))

function _convert(pnt::Ptr{T}, len) where {T}
    BT = basetype(T)
    buf, out = _allocate(BT, len)
    fin = bytoff(out, len)
    while out < fin
        set_codeunit!(out, get_codeunit(pnt))
        pnt += sizeof(T)
        out += sizeof(BT)
    end
    buf, out
end

function convert(::Type{<:Str{UTF16CSE}}, bytes::AbstractArray{UInt8})
    is_empty(bytes) && return empty_utf16
    # Note: there are much better ways of detecting what the likely encoding is here,
    # this only deals with big or little-ending UTF-32
    # It really should detect at a minimum UTF-8, UTF-16 big and little
    len = length(bytes)
    isodd(len) && unierror(UTF_ERR_ODD_BYTES_16, len, 0)
    len >>>= 1
    @preserve bytes begin
        pnt = pointer(bytes)
        if isodd(reinterpret(UInt, pnt))
            buf, out = _convert(reinterpret(Ptr{UInt16_U}, pnt), len, swappedtype(UInt16_U))
        else
            buf, out = _convert(reinterpret(Ptr{UInt16}, pnt), len, swappedtype(UInt16))
        end
        is_valid(UTF16Str, out, len) || unierror(UTF_ERR_INVALID, 0, 0)
        Str(UTF16CSE, buf)
    end
end

@inline function pushchar!(rst, uc)
    if uc <= 0x0ffff
        push!(rst, uc%UInt16)
    else
        w1, w2 = get_utf16(uc)
        push!(rst, w1, w2)
    end
end

"""Handle case where result vector is longer"""
function _maprest(fun, str, len, pnt, fin, buf, out, uc)
    rst = Vector{UInt16}()
    pushchar!(rst, uc)
    while pnt < fin
        ch = get_codeunit(pnt)%UInt32
        # check for surrogate pair
        is_surrogate_lead(ch) &&
            (ch = get_supplementary(ch, get_codeunit(pnt += 2)))
        pushchar!(rst, check_valid(UInt32(fun(ch%UTF32Chr)), chrdiff(pnt, pointer(str))))
        pnt += 2
    end
    # We now have a vector to add to the end of buf
    lenrst = length(rst)
    totbuf, totpnt = _allocate(UInt16, len + lenrst)
    unsafe_copyto!(totpnt, 1, pnt, 1, (out - pointer(buf))>>>1)
    unsafe_copyto!(totpnt, out + 1, pointer(rst), 1, lenrst)
    Str(UTF16CSE, totbuf)
end

function _map(::Type{C}, ::Type{T}, fun, len, str) where {C<:CSE, T<:Str}
    pnt = pointer(str)
    buf = Base.StringVector(len * sizeof(UInt16))
    beg = out = reinterpret(Ptr{UInt16}, pointer(buf))
    surrflag = false
    fin = pnt + sizeof(str)
    outend = out + sizeof(str)
    while pnt < fin
        ch = get_codeunit(pnt)%UInt32
        # check for surrogate pair
        T == UTF16Str && is_surrogate_lead(ch) &&
            (ch = get_supplementary(ch, get_codeunit(pnt += 2)))
        # Note: no checking for invalid here, UTF16Str is always valid
        uc = check_valid(UInt32(fun(ch%UTF32Chr)), chrdiff(pnt, pointer(str)))
        if uc < 0x10000
            out < outend || return _maprest(fun, str, len, buf, pnt, fin, out, uc)
            set_codeunit!(out, uc%UInt16)
            out += 2
        else
            out + 2 < outend || return _maprest(fun, str, len, buf, pnt, fin, out, uc)
            w1, w2 = get_utf16(uc)
            set_codeunit!(out, w1)
            set_codeunit!(out + 2, w2)
            surrflag = true
            out += 4
        end
        pnt += 2
    end
    out < outend && resize!(buf, chrdiff(UInt16, out - beg))
    s = String(buf)
    if !surrflag
        Str(C, s)
    elseif T == _UCS2Str
        # Convert to 32-bit, to keep result in UniStr type union
        # TODO: check this
        convert(_UTF32Str, Str(UTF16CSE, s))
    else
        Str(UTF16CSE, s)
    end
end

map(fun, str::MaybeSub{T}) where {C<:Union{UCS2_CSEs, UTF16CSE},T<:Str{C}} =
    (len = ncodeunits(str)) == 0 ? empty_str(T) : @preserve str _map(C, T, fun, len, str)


text:
using Documenter
import Match

makedocs(modules=[Match],
    clean=false,
    format=:html,
    sitename="Match.jl",
    authors="Kevin Squire and contributors.",
    pages=[
        "Home" => "index.md"
    ],
    # Use clean URLs, unless built as a "local" build
    html_prettyurls=!("local" in ARGS),)

deploydocs(repo="github.com/kmsquire/Match.jl.git",
    target="build",
    julia="0.6",
    deps=nothing,
    make=nothing,)

ext: jl


text:
@inline (bc::BoundaryFunction{:x, Y, Z})(j, k, grid, time, args...) where {Y, Z} =
    bc.func(ynode(Y, j, grid), znode(Z, k, grid), time)

@inline (bc::BoundaryFunction{:y, X, Z})(i, k, grid, time, args...) where {X, Z} =
    bc.func(xnode(X, i, grid), znode(Z, k, grid), time)

@inline (bc::BoundaryFunction{:z, X, Y})(i, j, grid, time, args...) where {X, Y} =
    bc.func(xnode(X, i, grid), ynode(Y, j, grid), time)


text:
function printpkgstyle(io::IO, cmd::Symbol, text::String, ignore_indent::Bool=false; color=:green)
    indent = textwidth(string(:Precompiling)) # "Precompiling" is the longest operation
    ignore_indent && (indent = 0)
    printstyled(io, lpad(string(cmd), indent), color=color, bold=true)
    println(io, " ", text)
end

function linewrap(str::String; io = stdout_f(), padding = 0, width = Base.displaysize(io)[2])
    text_chunks = split(str, ' ')
    lines = String[""]
    for chunk in text_chunks
        new_line_attempt = string(last(lines), chunk, " ")
        if length(strip(new_line_attempt)) > width - padding
            lines[end] = strip(last(lines))
            push!(lines, string(chunk, " "))
        else
            lines[end] = new_line_attempt
        end
    end
    return lines
end

const URL_regex = r"((file|git|ssh|http(s)?)|(git@[\w\-\.]+))(:(//)?)([\w\.@\:/\-~]+)(\.git)?(/)?"x
isurl(r::String) = occursin(URL_regex, r)

stdlib_dir() = normpath(joinpath(Sys.BINDIR::String, "..", "share", "julia", "stdlib", "v$(VERSION.major).$(VERSION.minor)"))
stdlib_path(stdlib::String) = joinpath(stdlib_dir(), stdlib)

function pathrepr(path::String)
    # print stdlib paths as @stdlib/Name
    if startswith(path, stdlib_dir())
        path = "@stdlib/" * basename(path)
    end
    return "`" * Base.contractuser(path) * "`"
end

function set_readonly(path)
    for (root, dirs, files) in walkdir(path)
        for file in files
            filepath = joinpath(root, file)
            # `chmod` on a link would change the permissions of the target.  If
            # the link points to a file within the same root, it will be
            # chmod'ed anyway, but we don't want to make directories read-only.
            # It's better not to mess with the other cases (links to files
            # outside of the root, links to non-file/non-directories, etc...)
            islink(filepath) && continue
            fmode = filemode(filepath)
            try
                chmod(filepath, fmode & (typemax(fmode) ⊻ 0o222))
            catch
            end
        end
    end
    return nothing
end
set_readonly(::Nothing) = nothing

# try to call realpath on as much as possible
function safe_realpath(path)
    if ispath(path)
        try
            return realpath(path)
        catch
            return path
        end
    end
    a, b = splitdir(path)
    return joinpath(safe_realpath(a), b)
end

# Windows sometimes throw on `isdir`...
function isdir_nothrow(path::String)
    try isdir(path)
    catch e
        false
    end
end

function isfile_nothrow(path::String)
    try isfile(path)
    catch e
        false
    end
end

function casesensitive_isdir(dir::String)
    dir = abspath(dir)
    lastdir = splitpath(dir)[end]
    isdir_nothrow(dir) && lastdir in readdir(joinpath(dir, ".."))
end

get_bool_env(name::String; default::String="false") =
    lowercase(get(ENV, name, default)) in ("t", "true", "y", "yes", "1")

## ordering of UUIDs ##
if VERSION < v"1.2.0-DEV.269"  # Defined in Base as of #30947
    Base.isless(a::UUID, b::UUID) = a.value < b.value
end


text:
# Map projections

export UNSEEN,
    lat2colat, colat2lat,
    project,
    equiprojinv, mollweideprojinv, orthoinv,
    equirectangular, mollweide, orthographic,
    orthographic2, gnomonic

import RecipesBase

const UNSEEN = -1.6375e+30

"""
    project(invprojfn, m::Map{T, O}, bmpwidth, bmpheight; kwargs...) where {T <: Number, O <: Order}

Return a 2D bitmap (array) containing a cartographic projection of the
map and a 2D bitmap containing a boolean mask. The size of the bitmap
is `bmpwidth`×`bmpheight` pixels. The function `projfn` must be a
function which accepts as input two parameters `x` and `y` (numbers
between -1 and 1).

The following keywords can be used in the call:

- `center`: 2-tuple specifying the location (colatitude, longitude) of the sky
  point that is to be placed in the middle of the image (in radians)
- `unseen`: by default, Healpix maps use the value -1.6375e+30 to mark
  unseen pixels. You can specify a different value using this
  keyword. This should not be used in common applications.

Return a `Array{Union{Missing, Float32}}` containing the intensity of
each pixel. Pixels falling outside the projection are marked as NaN,
and unseen pixels are marked as `missing`.
"""
function project(invprojfn, m::Map{T,O}, bmpwidth, bmpheight,
                 projparams = Dict()) where {T <: Number,O <: Healpix.Order}

    center = get(projparams, :center, (0, 0))
    unseen = get(projparams, :unseen, UNSEEN)
    desttype = get(projparams, :desttype, Float32)

    img = Array{desttype}(undef, bmpheight, bmpwidth)
    masked = zeros(Bool, bmpheight, bmpwidth)

    anymasked = false
    for j in 1:bmpheight
        y = 2 * (j - 1) / (bmpheight - 1) - 1
        for i in 1:bmpwidth
            x = 2 * (i - 1) / (bmpwidth - 1) - 1
            visible, lat, long = invprojfn(x, y)
            if visible
                value = m.pixels[Healpix.ang2pix(m, lat2colat(lat), long)]
                if ismissing(value) || isnan(value) || (
                    !ismissing(unseen) && unseen == value)
                    img[j, i] = NaN
                    masked[j, i] = true
                    anymasked = true
                else
                    img[j, i] = value
                end
            else
                img[j, i] = NaN
            end
        end
    end

    img, masked, anymasked
end

################################################################################

lat2colat(x) = π / 2 - x
colat2lat(x) = π / 2 - x

@doc raw"""
   lat2colat(x)
   colat2lat(x)

Convert colatitude into latitude and vice versa. Both `x` and the
result are expressed in radians.
"""
lat2colat, colat2lat

################################################################################

"""
    equiproj(lat, lon)

Equirectangular projection. Given the latitude `lat` (in radians) and the
longitude (in radians), return a tuple (Bool, Number, Number) where the
first Boolean is a flag telling if the point falls within the projection (true)
or not (false), and the two numbers are the x and y coordinates of the point
on the projection plane (both are in the range [−1, 1]).
"""
function equiproj(lat, lon)
    # We use `rem2pi` because we need angles in the range [-π, +π]
    x, y = ((rem2pi(lon + π, RoundNearest) - π) / π, 2 * (rem2pi(lat + π, RoundNearest) - π) / π)
    (true, x, y)
end

"""
    equiprojinv(x, y)

Inverse equirectangular projection. Given a point (x, y)
on the plane [-1, 1] × [-1, 1], return a tuple (Bool, Number, Number)
where the first Boolean is a flag telling if the point falls
within the projection (true) or not (false), and the two numbers
are the latitude and longitude in radians.
"""
function equiprojinv(x, y)
    ((-1 ≤ x ≤ 1) && (-1 ≤ y ≤ 1)) || return (false, 0, 0)

    (true, π / 2 * y, π * x)
end

"""
    mollweideproj(lat, lon)

Mollweide projection. Given the latitude `lat` (in radians) and the
longitude (in radians), return a tuple (Bool, Number, Number) where the
first Boolean is a flag telling if the point falls within the projection (true)
or not (false), and the two numbers are the x and y coordinates of the point
on the projection plane (both are in the range [−1, 1]).
"""
function mollweideproj(lat, lon)
    sinlat, coslat = sincos(lat)
    (true, sinlat, 2 / π * lon * coslat)
end

"""
    mollweideprojinv(x, y)

Inverse Mollweide projection. Given a point (x, y) on the plane,
with x ∈ [-1, 1], y ∈ [-1, 1], return a 3-tuple of type
(Bool, Number, Number). The boolean specifies if (x, y) falls within
the map (true) or not (false), the second and third arguments are
the latitude and longitude in radians.
"""
function mollweideprojinv(x, y)
    # See https://en.wikipedia.org/wiki/Mollweide_projection, we set
    #
    #     R = 1/√2
    #
    # x ∈ [-1, 1], y ∈ [-1, 1]

    x^2 + y^2 ≥ 1 && return (false, 0.0, 0.0)

    sinθ = y
    cosθ = sqrt(1 - sinθ^2)
    θ = asin(sinθ)

    lat = asin((2θ + 2 * sinθ * cosθ) / π)
    long = -2 * π * x / (2cosθ)
    (true, lat, long)

end

"""
    orthoinv(x, y, ϕ1, λ0)

Inverse orthographic projection centered on (ϕ1, λ0). Given a
point (x, y) on the plane, with x ∈ [-1, 1], y ∈ [-1, 1], return
a 3-tuple of type (Bool, Number, Number). The boolean specifies
if (x, y) falls within the map (true) or not (false), the second
and third arguments are the latitude and longitude in radians.
"""
function orthoinv(x, y, ϕ1, λ0)
    # Assume R = 1/√2. The notation ϕ1, λ0 closely follows
    # the book "Map projections — A working manual" by
    # John P. Snyder (page 145 and ff.)

    R = 1
    ρ = √(x^2 + y^2)
        if ρ > R
        return (false, zero(ϕ1), zero(λ0))
    end

    c = asin(ρ / R)
    sinc, cosc = sin(c), cos(c)
    if cosc < 0
        return (false, zero(ϕ1), zero(λ0))
    end

    if ρ ≈ 0
        return (true, ϕ1, λ0)
    end

    ϕ = asin(cosc * sin(ϕ1) + y * sinc * cos(ϕ1) / ρ)
    if ϕ1 ≈ π / 2
        λ = λ0 + atan(x, -y)
    elseif ϕ1 ≈ -π / 2
        λ = λ0 + atan(x, y)
    else
        λ = λ0 + atan(x * sinc, (ρ * cos(ϕ1) * cosc - y * sin(ϕ1) * sinc))
    end

    (true, ϕ, λ)
end

"""
    function ortho2inv(x, y, ϕ1, λ0)

Inverse stereo orthographic projection centered on (ϕ1, λ0). Given
a point (x, y) on the plane, with x ∈ [-1, 1], y ∈ [-1, 1], return
a 3-tuple of type (Bool, Number, Number). The boolean specifies
if (x, y) falls within the map (true) or not (false), the second
and third arguments are the latitude and longitude in radians.
"""
function ortho2inv(x, y, ϕ1, λ0)
    x ≤ 0 && return orthoinv(2x + 1, y, ϕ1, λ0)
    orthoinv(2x - 1, y, ϕ1, λ0 + π)
end

"""
    function gnomonicinv(x, y, ϕ1, λ0, fov_rad)

Gnomonic projection centered on (ϕ1, λ0), with a field of view
equal to `fov_rad` (in radians).  Given a point (x, y) on the plane,
with x ∈ [-1, 1], y ∈ [-1, 1], return a 3-tuple of type (Bool,
Number, Number). The boolean specifies if (x, y) falls within
the map (true) or not (false), the second and third arguments
are the latitude and longitude in radians.
"""
function gnominv(x, y, ϕ1, λ0, ψ0, fov_rad)
    # We fix a Earth radius such that the field of view of the
    # projection is the one expected. Note that the formula
    # diverges if `fov_rad` is 90° (as expected).
    R = 1. / tan(fov_rad)

    # We use basic geometry to cast a 3D ray from the center of
    # the sphere to the tangent plane, and then we compute the
    # intersection between the ray and the sphere.
    gamma = 1 / sqrt(1 + (x^2 + y^2) / R^2)
    vecx, vecy, vecz = gamma * R, -gamma * x, -gamma * y

    # We implement the rotations in the following order:
    # 1. Rotation by -ψ0 (orientation) around the x axis
    # 2. Rotation by λ0 (longitude) around the z axis
    # 3. Rotation by ϕ1 (latitude) around the y axis

    sin_ϕ1, cos_ϕ1 = sincos(ϕ1)
    sin_λ0, cos_λ0 = sincos(λ0)
    sin_ψ0, cos_ψ0 = sincos(-ψ0) # Change the sign to match Healpy conventions

    vecxrot = (cos_ϕ1 * cos_λ0 * vecx +
               (-cos_ψ0 * cos_ϕ1 * sin_λ0 + sin_ψ0 * sin_ϕ1) * vecy +
               (-sin_ψ0 * cos_ϕ1 * sin_λ0 + sin_ϕ1 * cos_ψ0) * vecz)
    vecyrot = ((sin_λ0) * vecx +
               (cos_ψ0 * cos_λ0 + sin_ψ0 * sin_ϕ1 * sin_λ0) * vecy +
               (-sin_ψ0 * cos_λ0) * vecz)
    veczrot = ((-sin_ϕ1 * cos_λ0) * vecx +
               (sin_ψ0 * sin_ϕ1 * sin_λ0 + cos_ϕ1 * sin_ψ0) * vecy +
               (-sin_ψ0 * sin_ϕ1 * sin_λ0 + cos_ϕ1 * cos_ψ0) * vecz)

    theta, phi = vec2ang(vecxrot, vecyrot, veczrot)
    (true, colat2lat(theta), phi)
end

################################################################################

"""
    equirectangular(m::Map{T,O}; kwargs...) where {T <: Number, O <: Order}

High-level wrapper around `project` for equirectangular projections.
"""
function equirectangular(m::Map{T,O}, projparams = Dict()) where {T <: Number,O <: Order}
    width = get(projparams, :width, 720)
    height = get(projparams, :height, width)
    project(equiprojinv, m, width, height, projparams)
end

"""
    mollweide(m::Map{T,O}, projparams = Dict()) where {T <: Number,O <: Order}

High-level wrapper around `project` for Mollweide projections.

The following parameters can be set in the `projparams` dictionary:

- `width`: width of the image, in pixels (default: 720 pixels)
- `height`: height of the image, in pixels; if not specified, it will be assumed
  to be equal to `width`
"""
function mollweide(m::Map{T,O}, projparams = Dict()) where {T <: Number,O <: Order}
    width = get(projparams, :width, 720)
    height = get(projparams, :height, width ÷ 2)
    project(mollweideprojinv, m, width, height, projparams)
end

"""
    orthographic(m::Map{T,O}, projparams = Dict()) where {T <: Number,O <: Order}

High-level wrapper around `project` for orthographic projections.

The following parameters can be set in the `projparams` dictionary:

- `width`: width of the image, in pixels (default: 720 pixels)
- `height`: height of the image, in pixels; if not specified, it will be assumed
  to be equal to `width`
- `center`: position of the pixel in the middle of the left globe (*latitude* and
  longitude).
"""
function orthographic(m::Map{T,O}, projparams = Dict()) where {T <: Number,O <: Order}
    width = get(projparams, :width, 720)
    height = get(projparams, :height, width)
    ϕ0, λ0 = get(projparams, :center, (0, 0))
    project(m, width, height, projparams) do x, y
        orthoinv(x, y, ϕ0, λ0)
    end
end

"""
    orthographic2(m::Map{T,O}, projparams = Dict()) where {T <: Number,O <: Order}

High-level wrapper around `project` for stereo orthographic projections.

The following parameters can be set in the `projparams` dictionary:

- `width`: width of the image, in pixels (default: 720 pixels)
- `height`: height of the image, in pixels; if not specified, it will be assumed
  to be equal to `width`
- `center`: position of the pixel in the middle of the left globe (*latitude* and
  longitude). Default is (0, 0).
"""
function orthographic2(m::Map{T,O}, projparams = Dict()) where {T <: Number,O <: Order}
    width = get(projparams, :width, 720)
    height = get(projparams, :height, width ÷ 2)
    ϕ0, λ0 = get(projparams, :center, (0, 0))
    project(m, width, height, projparams) do x, y
ortho2inv(x, y, ϕ0, λ0)
    end
end

"""
    gnomonic(m::Map{T,O}, projparams = Dict()) where {T <: Number,O <: Order}

High-level wrapper around `project` for gnomonic projections.

The following parameters can be set in the `projparams` dictionary:

- `width`: width of the image, in pixels (default: 720 pixels)
- `height`: height of the image, in pixels; if not specified, it will be assumed
  to be equal to `width`
- `center`: position and orientation of the pixel in the middle. It is a 3-element
  tuple containing:

  1. The *latitude* of the pixel, in radians
  2. The longitude of the pixel, in radians
  3. The rotation to be applied to the image, in radians

- `fov_rad`: size of the image along the x and y axes, in radians (default: 15°)

# Example

````julia
plot(m, gnomonic, Dict(:fov_rad = deg2rad(1.5), :center = (0, 0, deg2rad(45))))
````
"""
function gnomonic(m::Map{T,O}, projparams = Dict()) where {T <: Number,O <: Order}
    width = get(projparams, :width, 720)
    height = get(projparams, :height, width)
    ϕ0, λ0, ψ0 = get(projparams, :center, (0, 0, 0))
    fov_rad = get(projparams, :fov_rad, deg2rad(15 / 60))
    project(m, width, height, projparams) do x, y
        gnominv(x, y, ϕ0, λ0, ψ0, fov_rad)
    end
end

################################################################################

RecipesBase.@recipe function plot(m::Map{T,O},
                                  projection = mollweide,
                                  projparams = Dict()) where {T <: Number,O <: Order}

    img, mask, anymasked = projection(m, projparams)

    if anymasked
        RecipesBase.@series begin
            seriestype --> :shape
            primary --> false
            c --> :grey
            line --> :grey

            width, height = size(mask)
            xm = Float64[]
            ym = Float64[]
            for y in 1:height
                curx = 1
                # Instead of drawing each single masked pixel as a
                # square, squash together long runs of masked pixels
                # into one rectangle
                while curx ≤ width
                    if mask[curx, y]
                        startx = curx
                        while curx ≤ width && mask[curx, y]
                            curx += 1
                        end

                        append!(xm, [
                            startx - 0.5,
                            startx - 0.5,
                            curx + 0.5,
                            curx + 0.5, NaN
                        ])
                        append!(ym, [
                            y - 0.5,
                            y + 0.5,
                            y + 0.5,
                    y - 0.5, NaN])
                    else
                    curx += 1
                    end
                end
            end

            ym, xm
        end
    end

    seriestype --> :heatmap
    aspect_ratio --> 1
    colorbar --> :bottom
    framestyle --> :none

    img
end

@doc raw"""
    plot(m::Map{T,O}, projection = mollweide, projparams = Dict())

Draw a representation of the map, using some specific projection. The
parameter `projection` must be a function returning the
bitmap. Possible values for `projection` are the following:

- `equirectangular`
- `mollweide`
- `orthographic`
- `orthographic2`
- `gnomonic`

You can define your own projections.

The dictionary `projparams` allows to hack a number of parameters used
in the projection.

# References

See also [`equirectangular`](@ref), [`mollweide`](@ref),
[`orthographic`](@ref), [`orthographic2`](@ref), and
[`gnomonic`](@ref).
"""
plot


text:
@testset "1380.lucky-numbers-in-a-matrix.jl" begin
    @test lucky_numbers(hcat([[3,7,8],[9,11,13],[15,16,17]]...)) == [15]
    @test lucky_numbers(hcat([[1,10,4,2],[9,3,8,7],[15,16,17,12]]...)) == [12]
    @test lucky_numbers(hcat([[7,8],[1,2]]...)) == [7]
end


text:
#
# ```@setup rawspectra
# # the example command line input for this script
# ARGS = ["example.toml", "P143hm1", "P143hm2", "--plot"]
# ```

# # [Raw Spectra (rawspectra.jl)](@id rawspectra)
# The first step in the pipeline is simply to compute the pseudo-spectrum between the
# maps ``X`` and ``Y``.
# We define the pseudo-spectrum ``\widetilde{C}_{\ell}`` as the
# result of the estimator on spherical harmonic coefficients of the *masked* sky,
# ```math
# \widetilde{C}_{\ell} = \frac{1}{2\ell+1} \sum_m a^{i,X}_{\ell m} a^{j,Y}_{\ell m}.
# ```
# Since we mask the galaxy and point sources, this is a biased estimate of the underlying
# power spectrum. The mask couples modes together, and also removes power from parts of the
# sky. This coupling is described by a linear operator ``\mathbf{M}``, the mode-coupling
# matrix.
# For more details on spectra and mode-coupling, please refer to the [documentation for
# PowerSpectra.jl](https://xzackli.github.io/PowerSpectra.jl/dev/spectra/).
# If this matrix is known, then one can perform a linear solve to obtain an unbiased
# estimate of the underlying power spectrum ``C_{\ell}``,
# ```math
# \langle\widetilde{C}_{\ell}\rangle =
# \mathbf{M}^{XY}(i,j)_{\ell_1 \ell_2} \langle {C}_{\ell} \rangle,
# ```
#
# This script performs this linear solve, *without accounting for beams*. The noise spectra
# are estimated from the difference of auto- and cross-spectra,
# The command-line syntax for using this component to compute mode-coupled spectra is
#
# ```@raw html
# <pre class="shell">
# <code class="language-shell hljs">$ julia rawspectra.jl global.toml [map1] [map2]</code></pre>
# ```
# `[map1]` and `[map2]` must be names of maps described in global.toml.
#
# This page shows the results of running the command
# ```@raw html
# <pre class="shell">
# <code class="language-shell hljs">$ julia src/rawspectra.jl example.toml P143hm1 P143hm2</code></pre>
# ```
# The first step is just to unpack the command-line arguments, which consist of the
# TOML config file and the map names, which we term channels 1 and 2.

configfile, mapid1, mapid2 = ARGS

#
# ## File Loading and Cleanup
# We start by loading the necessary packages.

using TOML
using Healpix
using PowerSpectra
using Plots

#
config = TOML.parsefile(configfile)

# For both input channels, we need to do some pre-processing steps.
# 1. We first load in the ``I``, ``Q``, and ``U`` Stokes vectors, which are FITS
#    columns 1, 2, and 3 in the Planck map files. These must be converted from nested to ring
#    ordered, to perform SHTs.
# 2. We read in the corresponding masks in temperature and polarization.
# 3. We zero the missing pixels in the maps, and also zero the corresponding pixels in the
#    masks.
# 4. We convert from ``\mathrm{K}_{\mathrm{CMB}}`` to ``\mu \mathrm{K}_{\mathrm{CMB}}``.
# 5. Apply a small polarization amplitude adjustment, listed as `poleff` in the config.
# 6. We also remove some noisy pixels with ``\mathrm{\sigma}(Q) > 10^6\,\mu\mathrm{K}^2``
#    or ``\mathrm{\sigma}(U) > 10^6\,\mu\mathrm{K}^2``, or if they are negative. This
#    removes a handful of pixels in
#    the 2018 maps at 100 GHz which interfere with covariance estimation.
# 7. Estimate and subtract the pseudo-monopole and pseudo-dipole.

"""Return (polarized map, maskT, maskP) given a config and map identifier"""
function load_maps_and_masks(config, mapid, maptype=Float64)
    ## read map file
    println("Reading ", config["map"][mapid])
    mapfile = joinpath(config["scratch"], "maps", config["map"][mapid])
    polmap = PolarizedHealpixMap(
        nest2ring(readMapFromFITS(mapfile, 1, maptype)),  # I
        nest2ring(readMapFromFITS(mapfile, 2, maptype)),  # Q
        nest2ring(readMapFromFITS(mapfile, 3, maptype)))  # U

    ## read maskT and maskP
    maskfileT = joinpath(config["scratch"], "masks", config["maskT"][mapid])
    maskfileP = joinpath(config["scratch"], "masks", config["maskP"][mapid])
    maskT = readMapFromFITS(maskfileT, 1, maptype)
    maskP = readMapFromFITS(maskfileP, 1, maptype)

    ## read Q and U pixel variances, and convert to μK
    covQQ = nest2ring(readMapFromFITS(mapfile, 8, maptype)) .* 1e12
    covUU = nest2ring(readMapFromFITS(mapfile, 10, maptype)) .* 1e12

    ## identify missing pixels and also pixels with crazy variances
    for p in eachindex(maskT)
        if (polmap.i[p] < -1.6e30) | (covQQ[p] > 1e6) | (covUU[p] > 1e6) |
                (covQQ[p] < 0) | (covUU[p] < 0)
            maskT[p] = 0.
            maskP[p] = 0.
        end
    end

    ## go from KCMB to μKCMB, and apply polarization factor
    poleff = config["poleff"][mapid]
    scale!(polmap, 1e6, 1e6 * poleff)  # apply 1e6 to (I) and 1e6 * poleff to (Q,U)

    ## fit and remove pseudo-monopole/dipole in I
    monopole, dipole = fitdipole(polmap.i * maskT)
    subtract_monopole_dipole!(polmap.i, monopole, dipole)

    return polmap, maskT, maskP
end

# We'll use this function for the half-missions involved here.

m₁, maskT₁, maskP₁ = load_maps_and_masks(config, mapid1)
m₂, maskT₂, maskP₂ = load_maps_and_masks(config, mapid2)
plot(m₁.i, clim=(-200,200))  # plot the intensity map

#
if "--plot" ∈ ARGS
    plot(maskT₁)  # show the temperature mask
end

#
run_name = config["general"]["name"]
function save_if_needed(mapobj, mapfile)
    if isfile(mapfile) == false
        saveToFITS(mapobj, mapfile)
    end
end

save_if_needed(maskT₁, joinpath(config["scratch"], "masks", "$(run_name)_$(mapid1)_maskT.fits"))
save_if_needed(maskP₁, joinpath(config["scratch"], "masks", "$(run_name)_$(mapid1)_maskP.fits"))
save_if_needed(maskT₂, joinpath(config["scratch"], "masks", "$(run_name)_$(mapid2)_maskT.fits"))
save_if_needed(maskP₂, joinpath(config["scratch"], "masks", "$(run_name)_$(mapid2)_maskP.fits"))

#
# # Computing Spectra and Saving
#
# Once you have cleaned up maps and masks, you compute the
# calculation is described in [PowerSpectra - Mode Coupling](https://xzackli.github.io/PowerSpectra.jl/dev/spectra/#Mode-Coupling-for-EE,-EB,-BB).
# That package has a utility function [`master`](https://xzackli.github.io/PowerSpectra.jl/dev/module_index/#PowerSpectra.master-Tuple{Healpix.PolarizedHealpixMap,%20Healpix.Map,%20Healpix.Map,%20Healpix.PolarizedHealpixMap,%20Healpix.Map,%20Healpix.Map})
# that performs the full MASTER calculation on two ``IQU`` maps with associated masks.
#

## do the mode coupling on all T, E, and B spectra
Cl = master(m₁, maskT₁, maskP₁,
            m₂, maskT₂, maskP₂)
nside = maskT₁.resolution.nside  # get the resolution from any of the maps
lmax = nside2lmax(nside)
println(keys(Cl))  # check what spectra were computed

# The PowerSpectra.jl package has the Planck bestfit theory and beams as utility functions,
# for demo and testing purposes. We can use it that for plotting here.

if "--plot" ∈ ARGS
    spec = :TT
    Wl = PowerSpectra.planck_beam_Wl("143", "hm1", "143", "hm2", spec, spec; lmax=lmax)
    pixwinT = SpectralVector(pixwin(nside)[1:(lmax+1)])
    ell = eachindex(Wl)
    prefactor = ell .* (ell .+ 1) ./ (2π)
    plot( prefactor .*  Cl[spec] ./ (Wl .* pixwinT.^2),
        label="\$D_{\\ell}\$", xlim=(0,2nside))
    theory = PowerSpectra.planck_theory_Dl()
    plot!(theory[spec], label="theory $(spec)")
end

# Now we save our spectra.
## set up spectra path
using CSV, DataFrames
spectrapath = joinpath(config["scratch"], "rawspectra")
mkpath(spectrapath)

## assemble a table with the ells and spectra
df = DataFrame()
df[!,:ell] = eachindex(Cl[first(keys(Cl))])
for spec in keys(Cl)
    df[!,spec] = parent(Cl[spec])
end

CSV.write(joinpath(spectrapath, "$(run_name)_$(mapid1)x$(mapid2).csv"), df)

max_forks_repo_name: sgiardie/PSpipe


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

#
# package typedb.client.concept.answer;
#
# import typedb.client.api.answer.Numeric;
# import typedb.client.common.exception.TypeDBClientException;
# import typedb.protocol.AnswerProto;
#
# import javax.annotation.Nullable;
#
# import static typedb.client.common.exception.ErrorMessage.Internal.ILLEGAL_CAST;
# import static typedb.client.common.exception.ErrorMessage.Query.BAD_ANSWER_TYPE;
#
# public class NumericImpl implements Numeric {
#
#     @Nullable
#     private final Long longValue;
#     @Nullable
#     private final Double doubleValue;
#
#     private NumericImpl(@Nullable Long longValue, @Nullable Double doubleValue) {
#         this.longValue = longValue;
#         this.doubleValue = doubleValue;
#     }
#
#     public static NumericImpl of(AnswerProto.Numeric numeric) {
#         switch (numeric.getValueCase()) {
#             case LONG_VALUE:
#                 return NumericImpl.ofLong(numeric.getLongValue());
#             case DOUBLE_VALUE:
#                 return NumericImpl.ofDouble(numeric.getDoubleValue());
#             case NAN:
#                 return NumericImpl.ofNaN();
#             default:
#                 throw new TypeDBClientException(BAD_ANSWER_TYPE, numeric.getValueCase());
#         }
#     }
#
#     private static NumericImpl ofLong(long value) {
#         return new NumericImpl(value, null);
#     }
#
#     private static NumericImpl ofDouble(double value) {
#         return new NumericImpl(null, value);
#     }
#
#     private static NumericImpl ofNaN() {
#         return new NumericImpl(null, null);
#     }
#
#     @Override
#     public boolean isLong() {
#         return longValue != null;
#     }
#
#     @Override
#     public boolean isDouble() {
#         return doubleValue != null;
#     }
#
#     @Override
#     public boolean isNaN() {
#         return !isLong() && !isDouble();
#     }
#
#     @Override
#     public long asLong() {
#         if (isLong()) return longValue;
#         else throw new TypeDBClientException(ILLEGAL_CAST, Long.class);
#     }
#
#     @Override
#     public Double asDouble() {
#         if (isDouble()) return doubleValue;
#         else throw new TypeDBClientException(ILLEGAL_CAST, Double.class);
#     }
#
#     @Override
#     public Number asNumber() {
#         if (isLong()) return longValue;
#         else if (isDouble()) return doubleValue;
#         else throw new TypeDBClientException(ILLEGAL_CAST, Number.class);
#     }
# }


text:
"""
    MRCData{T<:Number,N,EH<:MRCExtendedHeader,D<:AbstractArray{T,N}} <: AbstractArray{T,N}

Container for electron density data loaded from an MRC file.

Each instance of `MRCData` carries with it a header and extended header, respectively
accessed with [`header`](@ref) and [`extendedheader`](@ref).

If changes are made to the header, extended header, or the data, call [`updateheader!`](@ref)
before writing to file to ensure the header is stil consistent with the data.

    MRCData()

Create an empty array.

    MRCData(header[, extendedheader])

Create an array whose size is indicated by the entries in `header`.

    MRCData(size)
    MRCData(size...)

Create an array of the specified size.
"""
struct MRCData{T<:Number,N,EH,D} <: AbstractArray{T,N}
    header::MRCHeader
    extendedheader::EH
    data::D
end
function MRCData(header, extendedheader, data::AbstractArray{T,N}) where {T,N}
    return MRCData{T,N,typeof(extendedheader),typeof(data)}(header, extendedheader, data)
end
function MRCData(header=MRCHeader(), extendedheader=MRCExtendedHeader())
    data_size = size(header)
    data_length = prod(data_size)
    dtype = datatype(header)
    dims = ndims(header)
    s = ntuple(i -> data_size[i], dims)
    data = Array{dtype,dims}(undef, s)
    return MRCData(header, extendedheader, data)
end
function MRCData(size::NTuple{3,<:Integer})
    header = MRCHeader()
    for i in 1:3
        setproperty!(header, (:nx, :ny, :nz)[i], size[i])
    end
    return MRCData(header)
end
MRCData(size::Integer...) = MRCData(size)

"""
    header(data::MRCData) -> MRCHeader

Get header.
"""
header(d::MRCData) = d.header

"""
    extendedheader(data::MRCData) -> MRCExtendedHeader

Get extended header.
"""
extendedheader(d::MRCData) = d.extendedheader

"""
    updateheader!(data::MRCData; statistics = true)

Update the header stored in `data` from the data and its extended header.
Set `statistics=false` to avoid computing summary statistics from the data.
"""
function updateheader!(d::MRCData; statistics=true)
    h = header(d)

    # update size
    s = size(d)
    for i in eachindex(s)
        setproperty!(h, (:nx, :ny, :nz)[i], s[i])
    end

    # update space group
    if length(s) == 2
        h.nz = 1
        h.mz = 1
        if h.ispg < 0 || h.ispg > 230
            h.ispg = 0
        end
    elseif length(s) == 3 && h.ispg == 0
        h.ispg == 1
    end

    # misc
    h.mode = TYPE_TO_MODE[eltype(d)]
    h.nsymbt = bytesize(extendedheader(d))
    h.nlabl = sum(s -> length(rstrip(s)) > 0, h.label)

    # update statistics
    if statistics
        h.dmin, h.dmax = extrema(d)
        dmean = mean(d)
        h.dmean, h.rms = dmean, stdm(d, dmean; corrected=false)
    end

    return h
end

# Array overloads
@inline Base.parent(d::MRCData) = d.data

@inline Base.getindex(d::MRCData, idx::Int...) = getindex(parent(d), idx...)

@inline Base.setindex!(d::MRCData, val, idx::Int...) = setindex!(parent(d), val, idx...)

Base.size(d::MRCData) = size(parent(d))
Base.size(d::MRCData, i) = size(parent(d), i)

Base.ndims(d::MRCData) = ndims(parent(d))

Base.length(d::MRCData) = length(parent(d))

Base.iterate(d::MRCData, idx...) = iterate(parent(d), idx...)

Base.lastindex(d::MRCData) = lastindex(parent(d))

function Base.similar(d::MRCData)
    return typeof(d)(header(d), extendedheader(d), similar(parent(d)))
end

function Base.IndexStyle(::Type{A}) where {T,N,EH,D,A<:MRCData{T,N,EH,D}}
    return Base.IndexStyle(D)
end

"""
    eachmaprow(d::MRCData)

Return an iterator over map rows.
"""
eachmaprow(d::MRCData) = eachslice(d; dims=header(d).mapr)

"""
    eachmapcol(d::MRCData)

Return an iterator over columns.
"""
eachmapcol(d::MRCData) = eachslice(d; dims=header(d).mapc)

"""
    eachmapsection(d::MRCData)

Return an iterator over sections.
"""
eachmapsection(d::MRCData) = eachslice(d; dims=header(d).maps)

"""
    eachstackunit(d::MRCData)

Return an iterator over elements of stacks.
"""
function eachstackunit(d::MRCData)
    h = header(d)
    mrcdtype = mrcdatatype(h)
    return if mrcdtype ∈ (Volume, Image)
        (d for _ in 1:1)
    elseif mrcdtype == ImageStack
        eachmapsection(d)
    else
        error("Volume stacks not currently supported.")
    end
end


text:
using LinearAlgebra
using TensorToolbox

include("../utils/utils.jl")

function cost_function(cost, X, Xr)
    if cost == "LS"
        return norm(X-Xr)

    elseif cost == "KL"
        dkl = get_dkl(X, Xr)
        #return sum( X .* log.( X ./ Xr ) ) - sum(X) + sum(Xr)
        return dkl
    end
end

function reconst(S, A)
    N = ndims(S)
    Xr = S
    for n=1:N
        Xr = ttm(Xr, A[n], n)
    end
    return Xr
end

function init_based_SVD(X, reqrank)
    D = ndims(X)
    hosvd_X = hosvd(X, reqrank=reqrank)
    S = abs.(hosvd_X.cten)
    f = hosvd_X.fmat
    A = []
    for i = 1:D
        push!(A, abs.(f[i]))
    end

    return S, A
end

function init_random(tensor_size, reqrank)
    A = []
    N = length(reqrank)
    for n = 1:N
        An = rand(tensor_size[n], reqrank[n])
        push!(A, An)
    end
    S = rand( reqrank... )
    return S, A
end

"""
NMF Initializer. See Table 4 in the paper.
"""
function init_for_NMF(X, R; eps=1.0e-12, pre_iter_max_inloop=40)
    (m, l) = size(X)
    A = rand(m, R)
    S = rand(R, l)

    for iter = 1:pre_iter_max_inloop
        A .= max.(X*S', eps)
        A .= A .* ((X*S') ./ (A*S*S'))
        S .= max.(A'*X, eps)
        S .= S .* ( (A'*X) ./ (A'*A*S) )
    end
    return S, A
end

"""
Initializer based on NMF. See Table 5 in the paper.
"""
function init_based_NMF(X, reqrank; eps=1.0e-12, pre_iter_max=40, pre_iter_max_inloop=40)
    N = ndims(X)
    A = []
    for n=1:N
        _, An = init_for_NMF(tenmat(X, n), reqrank[n], pre_iter_max_inloop=pre_iter_max_inloop)
        push!(A, An)
    end

    S = rand( reqrank... )
    for iter = 1:pre_iter_max
        for n=1:N
            SAn = S
            for m in [1:n-1; n+1:N]
                SAn = ttm(SAn, A[m], m)
            end
            SAn = tenmat(SAn, n)
            A[n] = A[n] .* ( ( tenmat(X,n) * SAn' ) ./ (A[n] * SAn * SAn' ))
        end
        S = X
        for m=1:N
            S = ttm(S, A[m]', m)
        end
        S .= max.(S, eps)

        numerator = X
        denominator = S
        for m=1:N
            numerator = ttm(numerator, A[m]', m)
            denominator = ttm(denominator, A[m]'*A[m], m)
        end
        S .= S .* ( numerator ./ denominator )
    end
    return S, A
end

function update_An_LS(A, n, SAn, X)
    An = A[n] .* ( ( tenmat(X,n) * SAn' ) ./ (A[n] * SAn * SAn' ))
   return An
end

function update_An_KL(A, n, SAn, X)
    tensor_size = size(X)
    z = sum(SAn, dims=2)
    An = A[n] .* ((( tenmat(X,n) ./ ( A[n] * SAn ) ) * SAn' ) ./ ( ones(tensor_size[n]) * z' ))
    return An
end

function update_S_LS(S, A, X)
    N = ndims(S)
    numerator = X
    denominator = S
    for m=1:N
        numerator = ttm(numerator, A[m]', m)
        denominator = ttm(denominator, A[m]'*A[m], m)
    end
    S = S .* ( numerator ./ denominator )
    return S
end

function update_S_KL(S, A, X)
    N = ndims(S)
    numerator = X ./ reconst(S, A)
    denominator = ones( size(X)... )
    for m=1:N
        numerator = ttm(numerator, A[m]', m)
        denominator = ttm(denominator, A[m]', m)
    end
    S = S .* ( numerator ./ denominator )
    return S
end

function update(X, S, A, cost, max_iter=100, verbose=true, verbose_interval=20, tol=1.0e-5)
    N = ndims(S)

    cnt_iter = 0
    prev_LSerror = 0.0
    while(cnt_iter < max_iter)
        ############
        # update A #
        ############
        for n=1:N
            # get SAn
            SAn = S
            for m in [1:n-1;n+1:N]
                SAn = ttm(SAn, A[m], m)
            end
            SAn = tenmat(SAn, n)
            if cost == "LS"
                A[n] .= update_An_LS(A, n, SAn, X)
            elseif cost == "KL"
                A[n] .= update_An_KL(A, n, SAn, X)
            end
        end

        ########################
        # update Core tensor S #
        ########################
        if cost == "LS"
            S .= update_S_LS(S, A, X)
        elseif cost == "KL"
            S .= update_S_KL(S, A, X)
        end

        Xr = reconst(S, A)
        if verbose && (cnt_iter % verbose_interval == 0)
            error = cost_function(cost, X, Xr)
            println("$cnt_iter $error")
        end

        LSerror = cost_function("LS", X, Xr)
        if cnt_iter > 1 && abs(LSerror - prev_LSerror) < tol
            break
        end
        prev_LSerror = LSerror
        cnt_iter += 1
    end
    return S, A
end

"""
Non-Negative Tucker Decomposition
proposed by Young-Deok Kim et al.
See [original paper](https://ieeexplore.ieee.org/document/4270403)
[Implemented in also python](
http://tensorly.org/stable/_modules/tensorly/decomposition/_tucker.html#non_negative_tucker)

# Aruguments
- `X` : input non-negative tensor
- `reqrank` : Target Tucker rank, array
- `init_method` : initial values, "NMF" or "random"
- `cost` : cost function, "LS" or "KL"
- `verbose` : true or false
- `pre_iter_max` : iter_max of initialization based on "NMF"
- `pre_iter_max_inloop` : iter_max of initialization of "NMF"

The default values of `tol` and `max_iter` are decided by Tensorly implementation.
"""
function NTD(X, reqrank ;
        cost="LS", init_method="SVD", max_iter=10, verbose=true, verbose_interval=50,
        pre_iter_max=40, pre_iter_max_inloop=40, tol=1.0e-5)

    @assert ndims(X) == length(reqrank)
    tensor_size = size(X)

    #A[n] \in R^( tensor_size[n] \times reqrank[n] )
    #S \in R^(reqrank[1] \times ... \times reqrank[N])
    if init_method == "random"
        S, A = init_random(tensor_size, reqrank)
    elseif init_method == "NMF"
        S, A = init_based_NMF(X, reqrank, pre_iter_max=pre_iter_max, pre_iter_max_inloop=pre_iter_max_inloop)
    elseif init_method == "SVD"
        S, A = init_based_SVD(X, reqrank)
    else
        error("no init method ", init_method)
    end

    S, A = update(X, S, A, cost, max_iter, verbose, verbose_interval, tol)
    Xr = reconst(S, A)
    return S, A, Xr
end

max_stars_repo_name: gkazunii/Legendre-tucker-rank-reduction

max_issues_repo_name: gkazunii/Legendre-tucker-rank-reduction

max_forks_repo_name: gkazunii/Legendre-tucker-rank-reduction


text:
import PathWeightSampling
using Test

system = PathWeightSampling.cooperative_chemotaxis_system(dtimes = 0:0.1:10)
conf = PathWeightSampling.generate_full_configuration(system)
@test length(conf.s_traj.u[1]) == 1
@test length(conf.r_traj.u[1]) == 40
@test length(conf.x_traj.u[1]) == 2

mconf = PathWeightSampling.generate_configuration(system)
@test length(mconf.s_traj.u[1]) == 1
@test length(mconf.x_traj.u[1]) == 2

for u in conf.r_traj.u
    @test sum(u) == sum(conf.r_traj.u[1]) # this is true because the total number of receptors must be a constant
end

for u in conf.x_traj.u
    @test sum(u) == sum(conf.x_traj.u[1]) # this is true because the total number of CheY proteins must be a constant
end

cond_ens = PathWeightSampling.ConditionalEnsemble(system)
for t = 0:0.1:10
    local u0 = PathWeightSampling.sample_initial_condition(cond_ens)
    u_new, weight = PathWeightSampling.propagate(mconf, cond_ens, u0, (0.0, t))
    @test size(u0) == size(u_new)
    @test sum(u0[2:end]) == sum(u_new[2:end]) # this is true because the total number of receptors must be a constant
    @test !isinf(weight)
end

algorithm = PathWeightSampling.SMCEstimate(128)
result = PathWeightSampling.mutual_information(system, algorithm, num_samples = 1)
@test !all(isinf.(result.MutualInformation[1]))


text:
println("# DPM TESTS...")

# run demo code at first
println(" * try demo...")

using DataFrames

# number of samples per cluster
N = 200

# create clusters
Data = DataFrame( x = randn(N), y = randn(N), class = "cluster1" )

# append second cluster
append!(Data, DataFrame( x = randn(N) + 5, y = randn(N) + 5, class = "cluster2" ));

D = 2 # 2 dimensional data
N = 400 # number of data points

# data matrix
X = zeros(N, D)

X[:,1] = convert(Array, Data[:x])
X[:,2] = convert(Array, Data[:y])

# init base distribution parameters
mu0 = mean(X, 1)
kappa0 = 9.0
nu0 = 5.0
Sigma0 = eye(D) * 10

# base distribution and concentration parameter (Gaussian with Normal Inverse Wishart Prior)
H = WishartGaussian(vec(mu0), kappa0, nu0, Sigma0)

# train Dirichlet Process Mixture Model
models = train(DPM(H), Gibbs(), RandomInitialisation(k = 10), X);
train(DPM(H), Gibbs(maxiter = 1000), KMeansInitialisation(), X)

# use $\alpha = 10.0$
train(DPM(H, α = 10.0), Gibbs(), KMeansInitialisation(), X)

# use Random Initialisation
train(DPM(H), Gibbs(), RandomInitialisation(), X)

# use 50 iterations and 10 burnin runs
train(DPM(H), Gibbs(maxiter = 50, burnin = 10),
	KMeansInitialisation(), X)

# access the log likelihood values of all iteration
LLH = map(model -> model.energy, models)

# access the weights of all iterations
W = map(model -> model.weights, models)

# access the distributions of all iterations
D = map(model -> model.distributions, models)

println("# FINISHED")


text:
"""
    struct MultiFieldFEFunction <: CellField
      # private fields
    end
"""
struct MultiFieldFEFunction{T<:MultiFieldCellField} <: FEFunction
  single_fe_functions::Vector{<:SingleFieldFEFunction}
  free_values::AbstractArray
  fe_space::MultiFieldFESpace
  multi_cell_field::T

  function MultiFieldFEFunction(
    free_values::AbstractVector,
    space::MultiFieldFESpace,
    single_fe_functions::Vector{<:SingleFieldFEFunction})

    multi_cell_field = MultiFieldCellField(map(i->i.cell_field,single_fe_functions))
    T = typeof(multi_cell_field)

    new{T}(
      single_fe_functions,
      free_values,
      space,
      multi_cell_field)
  end
end

CellData.get_data(f::MultiFieldFEFunction) = get_data(f.multi_cell_field)
CellData.get_triangulation(f::MultiFieldFEFunction) = get_triangulation(f.multi_cell_field)
CellData.DomainStyle(::Type{MultiFieldFEFunction{T}}) where T = DomainStyle(T)
FESpaces.get_free_dof_values(f::MultiFieldFEFunction) = f.free_values
FESpaces.get_fe_space(f::MultiFieldFEFunction) = f.fe_space

function FESpaces.get_cell_dof_values(f::MultiFieldFEFunction)
  msg = """\n
  This method does not make sense for multi-field
  since each field can be defined on a different triangulation.
  Pass a triangulation in the second argument to get the DOF values
  on top of the corresponding cells.
  """
  trians = map(get_triangulation,f.fe_space.spaces)
  trian = first(trians)
  @check all(map(t->have_compatible_domains(t,trian),trians)) msg
  get_cell_dof_values(f,trian)
end

function FESpaces.get_cell_dof_values(f::MultiFieldFEFunction,trian::Triangulation)
  uhs = f.single_fe_functions
  blockmask = [ _can_be_restricted_to(get_triangulation(uh),trian) for uh in uhs ]
  active_block_ids = findall(blockmask)
  active_block_data = Any[ get_cell_dof_values(uhs[i],trian) for i in active_block_ids ]
  nblocks = length(uhs)
  lazy_map(BlockMap(nblocks,active_block_ids),active_block_data...)
end

function FESpaces.get_cell_dof_values(f::MultiFieldFEFunction,trian::SkeletonTriangulation)
  cell_values_plus = get_cell_dof_values(f,trian.plus)
  cell_values_minus = get_cell_dof_values(f,trian.minus)
  lazy_map(BlockMap(2,[1,2]),cell_values_plus,cell_values_minus)
end

"""
    num_fields(m::MultiFieldFEFunction)
"""
num_fields(m::MultiFieldFEFunction) = length(m.single_fe_functions)
Base.iterate(m::MultiFieldFEFunction) = iterate(m.single_fe_functions)
Base.iterate(m::MultiFieldFEFunction,state) = iterate(m.single_fe_functions,state)
Base.getindex(m::MultiFieldFEFunction,field_id::Integer) = m.single_fe_functions[field_id]


text:
using Test, TORA

@testset "TORA" begin
    # To be written...
end


text:
include("euler/euler.jl")
import StatsBase: countmap
using .BigIntegers: BigInteger
using BenchmarkTools
BenchmarkTools.DEFAULT_PARAMETERS.samples = 100

function compute(n::Int, modulo::Int)::BigInteger
    function get_digits(k::Int, limit::Int)::Vector{Vector{Int}}
        result = Vector{Vector{Int}}()
        if k ≤ 0 return result end
        if k ≤ limit append!(result, [[k]]) end
        for i ∈ 1:limit for j ∈ get_digits(k - i, i) append!(result, [vcat(i, j)]) end end
        return result
    end
    factorials = Dict{Int, BigInteger}()
    all_digits = BigInteger(invmod(9, modulo)) * BigInteger((powermod(10, n, modulo) - 1)) % modulo
    total = BigInteger(0)
    for digit ∈ 1:9
        for other_digits ∈ get_digits(digit, digit)
            len_digits = length(other_digits)
            if n - len_digits - 1 < 0 continue end
            digits = countmap(vcat(digit, other_digits))
            for m ∈ keys(digits)
                partial_digits = copy(digits)
                partial_digits[m] -= 1
                if !haskey(factorials, len_digits) factorials[len_digits] = prod(BigInteger(n - i - 1) for i ∈ 0:len_digits - 1) end
                multinomial = factorials[len_digits]
                for p ∈ values(partial_digits) multinomial ÷= factorial(p) end
                total += all_digits * m * multinomial
            end
        end
        total %= modulo
    end
    return total
end

compute(3, 10_000_000_000_000_000)

compute(7, 10_000_000_000_000_000)

compute(2_020, 10_000_000_000_000_000)

@benchmark compute(2_020, 10_000_000_000_000_000)

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
import Random: shuffle!

struct Context{I<:Integer, J<:AbstractArray{<:Integer,1}}
    n::I
    xs::J

    function Context(n, xs)
        n > 0 || error("context window must be a positive integer - window size $n provided")
        length(xs) > 1 || error("cannot iterate over an x of 1 or less")
        return new{typeof(n), typeof(xs)}(n, xs)
    end
end

function Base.iterate(c::Context{I}, state = 1) where I<:Integer
    if state < 1 || state > length(c.xs) + 1
        error("Invalid iterator state parameter")
    elseif state == length(c.xs) + 1
        return nothing
    else

        if state ≤ c.n
            state == 1 ? lv = I[] : lv = view(c.xs, 1:state-1)
        else
            lv = view(c.xs, state-c.n:state-1)
        end
        if state ≥ length(c.xs) - c.n + 1
            state == length(c.xs) ? uv = I[] : uv = view(c.xs, state+1:length(c.xs))
        else
            uv = view(c.xs, state+1:state+c.n)
        end
        word = c.xs[state]
        v = vcat(lv, uv)
        return ((word, v), state + 1)
    end
end

function Base.length(c::Context)
    return length(c.xs)
end


text:
using BinaryBuilder

sources = [
    "src"
]

# Bash recipe for building across all platforms
script = raw"""
mkdir build && cd build
cmake -DCMAKE_BUILD_TYPE=Release -DCMAKE_TOOLCHAIN_FILE=/opt/$target/$target.toolchain -DCMAKE_CXX_FLAGS="-march=x86-64" -DCMAKE_INSTALL_PREFIX=${prefix} -DCMAKE_FIND_ROOT_PATH=${prefix} -DJulia_PREFIX=${prefix} ../testlib
VERBOSE=ON cmake --build . --config Release --target install
"""

version_number = get(ENV, "TRAVIS_TAG", "")
if version_number == ""
    version_number = "v0.99.0"
end

# This transforms the dependency script to use the binaries generated in the same Travis run. This code should be skipped on a stand-alone project
### begin of part to be skipped ###
@show ws_root = dirname(dirname(dirname(dirname(@__FILE__))))
@show prod_dir = joinpath(ws_root, "products")
@show generated_script_name = joinpath(prod_dir,"build_libcxxwrap-julia-1.0.$(version_number).jl")
local_dep_script = ""
open(generated_script_name, "r") do generated_script
    while !eof(generated_script)
        l = readline(generated_script; keep=true)
        if startswith(l, "bin_prefix = \"https://github.com")
            l = "bin_prefix = \"$prod_dir\"\n"
        end
        global local_dep_script *= l
    end
end
### end of part to be skipped ###

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = Platform[]
_abis(p) = (:gcc7,:gcc8)
_archs(p) = (:x86_64, :i686)
_archs(::Type{Linux}) = (:x86_64,)
for p in (Linux,Windows)
    for a in _archs(p)
        for abi in _abis(p)
            push!(platforms, p(a, compiler_abi=CompilerABI(abi,:cxx11)))
        end
    end
end
push!(platforms, MacOS(:x86_64))

# The products that we will ensure are always built
products = prefix -> [
    LibraryProduct(prefix, "testlib", :testlib),
]

# Dependencies that must be installed before this package can be built
dependencies = [
   "https://github.com/JuliaPackaging/JuliaBuilder/releases/download/v1.0.0-2/build_Julia.v1.0.0.jl",
   BinaryBuilder.InlineBuildDependency(local_dep_script) # Replace this with build script from https://github.com/JuliaInterop/libcxxwrap-julia/releases
]

build_tarballs(ARGS, "testlib", VersionNumber(version_number), sources, script, platforms, products, dependencies)

max_stars_repo_name: dmikushin/libcxxwrap-julia

max_issues_repo_name: mrdaly/ColorCode

max_forks_repo_name: mrdaly/ColorCode


text:
# image utils

"""
Convert a Vector{UInt8} to an image of RGB Vector{N0f8}.
"""
function rgbUint8ToRgb(width::Int64, height::Int64, data::Vector{UInt8})::Array{RGB{N0f8}, 2}
  r=data[1:3:end];g=data[2:3:end];b=data[3:3:end]
  imgForm = rand(RGB{N0f8}, height, width)
  for w in 1:width
      for h in 0:(height-1)
          l = w + h * width
          imgForm[h+1,w] = RGB{N0f8}(r[l]/255.0, g[l]/255.0, b[l]/255.0)
      end
  end
  return imgForm
end

"""
Converts an RGB image to a JPEG image, returns buffer of UInt8 data.
"""
function rgbToJpeg(rgb::Array{RGB{N0f8}, 2})::Vector{UInt8}
  io = IOBuffer()
  save(Stream(format"JPEG",io), rgb)
  bytes = take!(io)
  return bytes
end

"""
Converts an RGB image to a PNG image, returns buffer of UInt8 data.
"""
function rgbToPng(rgb::Array{RGB{N0f8}, 2})::Vector{UInt8}
  io = IOBuffer()
  save(Stream(format"PNG",io), rgb)
  bytes = take!(io)
  return bytes
end

function imshowhackpng(im)
  filename = joinpath("/tmp","tempimgcaesar.png")
  imf = open(filename, "w")
  write(imf, im)
  close(imf)
  run(`eog $(filename)`)
end
function imshowhack(img::Array{T}) where {T <: Colorant}
  filename = joinpath("/tmp","caesarimshowhack.png")
  ImageMagick.save_(filename, img)
  run(`eog $(filename)`)
end

function roi(img, row, col; fov=50)
  r, c = size(img)
  left = round(Int,col-fov)
  right = round(Int,col+fov)
  top = round(Int,row-fov)
  bottom = round(Int,row+fov)
  if c < left
    left = c
  end
  if left < 1
    left=1
  end
  if c < right
    right = c
  end
  if right < 1
    right=1
  end
  if c < top
    top = c
  end
  if top < 1
    top=1
  end
  if c < bottom
    bottom = c
  end
  if bottom < 1
    bottom=1
  end

  img[top:bottom,left:right]
end


text:
using .GZip

# if GZip is loaded in, overwrite openvcf to this method
function openvcf(infile::String)
    if endswith(infile, ".vcf") || endswith(infile, ".bcf")
        return open(infile, "r")
    elseif endswith(infile, ".vcf.gz") || endswith(infile, ".bcf.gz")
        return GZip.open(infile, "r")
    else
        throw(ArgumentError("The filename must end with .vcf/.bcf or .vcf.gz/.bcf.gz"))
    end
end


text:
using Plots, StatsPlots, LaTeXStrings

techcolors = [:lightblue :cyan :yellow :darkgreen :lime :gray :orange :brown :blue]

function plot_generation_dispatch(p_gnt, h_hnt, G, n, T, region_n, technology_g, κ, C_E, κ′, C′_E)
    colors = techcolors
    p = Plots.plot(
        legend=:outertopright,
        size=(800, 600),
        title = "Hourly generation dispatch by technology in $(region_n[n])\nRenewables share = $(round(κ′,digits=3)) ≥ $κ\nCO2 reduction = $(round(C′_E,digits=3)) ≥ $C_E",
        titlefontsize = 10
    )
    for g in G
        plot!(p, T, [p_gnt[g, n, t] for t in T],
              color = colors[g],
              alpha=0.3,
              xlabel=L"t",
              ylabel=L"p_{g,n,t}\,\mathrm{[MWh]}",
              label=technology_g[g])
    end
    plot!(p, T, [h_hnt[1, n, t] for t in T],
          color = colors[9],
          alpha=0.3,
          label="hydro")
    return p
end

function plot_balance_stacked(p_gnt, p̄_gn, h_hnt, h̄_hn, HRmax_n, G, n, T, region_n, technology_g, κ, C_E, κ′, C′_E)
    colors = techcolors
    p = Plots.plot(
        legend=:outertopright,
        size=(800, 600),
        title = "Hourly generation dispatch by technology in $(region_n[n])\nRenewables share = $(round(κ′,digits=3)) ≥ $κ\nCO2 reduction = $(round(C′_E,digits=3)) ≥ $C_E",
        titlefontsize = 10
    )
    for g in G
        plot!(p, T, [p_gnt[g, n, t] for t in T],
              color = colors[g],
              alpha=0.3,
              xlabel=L"t",
              ylabel=L"p_{g,n,t}\,\mathrm{[MWh]}",
              label=technology_g[g])
        plot!(p, T, [p̄_gn[g, n] for t in T],
              color = colors[g],
              label="")
    end
    plot!(p, T, [h_hnt[1, n, t] for t in T],
          color = colors[9],
          alpha=0.3,
          label="hydro")
    plot!(p, T, [h̄_hn[n] + HRmax_n[n] for t in T],
          color = colors[9],
          label="")
    return p
end

function plot_generation_capacities(p̄_gn, h̄_hn, HRmax_n, G, n, region_n, technology_g, κ, C_E, κ′, C′_E)
    StatsPlots.bar([G;9], [[p̄_gn[g, n] for g in G]; [h̄_hn[n] + HRmax_n[n]]],
        xticks=([G;9], [technology_g;"hydro"]),
        ylabel=L"\bar{p}_{g,n}\,\mathrm{[MW]}",
        title = "Generation capacity by technology in $(region_n[n])\nRenewables share = $(round(κ′,digits=3)) ≥ $κ\nCO2 reduction = $(round(C′_E,digits=3)) ≥ $C_E",
        titlefontsize = 10,
        color=permutedims(techcolors),
        alpha=0.7,
        legend=false)
end

function plot_generation_capacities_stacked(p̄_gn, h̄_hn, HRmax_n, N, H, region_n, technology_g, κ, C_E, κ′, C′_E)
    p̄_gn
    H_tot = sum(h̄_hn[h,:] for h in H) .+ HRmax_n
    dispatches = permutedims([p̄_gn; permutedims(H_tot)])
    StatsPlots.groupedbar(dispatches,
        bar_position = :stack,
        bar_width=0.7,
        xticks=(N, region_n),
        ylabel=L"\bar{p}_{g,n}\,\mathrm{[MW]}",
        title = "Generation capacity by region and technology\n Renewables share = $(round(κ′,digits=3)) ≥ $κ\n CO2 reduction = $(round(C′_E,digits=3)) ≥ $C_E",
        titlefontsize = 10,
        color=techcolors,
        labels=permutedims([technology_g;"hydro"]),
        lw = 0,
        alpha=0.7)
end

function plot_transmission_flow(f_lt, f̄_l, l, L, T, region_n, κ, C_E, κ′, C′_E)
    p = Plots.plot(T, [f_lt[l, t] for t in T],
             xlabel=L"t",
             ylabel=L"f_{l,t}\,\mathrm{[MWh]}",
             title = "Hourly transmission flow between $(region_n[L[l][1]]) and $(region_n[L[l][2]])\nRenewables share = $(round(κ′,digits=3)) ≥ $κ\nCO2 reduction = $(round(C′_E,digits=3)) ≥ $C_E",
             titlefontsize = 10,
             alpha=0.7,
             legend=false,
             size=(800, 600))
    plot!(p, T, [f̄_l[l] for t in T], color=:green)
    plot!(p, T, [-f̄_l[l] for t in T], color=:green)
    return p
end

function plot_transmission_bars(f_lt, L, L_ind, T, region_n, κ, C_E, κ′, C′_E)
    lines = Array{AbstractString, 1}(undef, length(L))
    for i in L_ind
        lines[i] = "$(region_n[L[i][1]])-$(region_n[L[i][2]])"
    end
    f_l = sum(f_lt[:,t] for t in T)
    p = StatsPlots.plot(size=(800, 600),
             ylabel=L"\sum_t f_{l,t}\,\mathrm{[MWh]}",
             title = "Transmission flow by line\nRenewables share = $(round(κ′,digits=3)) ≥ $κ\nCO2 reduction = $(round(C′_E,digits=3)) ≥ $C_E",
             titlefontsize = 10,
             xticks=(L_ind, lines),
             xrotation = 90,
             legend=false,
             )
    bar!(L_ind, [f_l[l] for l in L_ind], alpha = 0.7, lw = 0)
    return p
end

function plot_transmission_capacities(f̄_l, L, L_ind, region_n, κ, C_E, κ′, C′_E)
    lines = Array{AbstractString, 1}(undef, length(L))
    for i in L_ind
        lines[i] = "$(region_n[L[i][1]])-$(region_n[L[i][2]])"
    end
    StatsPlots.bar(L_ind, [f̄_l[l] for l in L_ind],
        title = "Transmission capacity by line\nRenewables share = $(round(κ′,digits=3)) ≥ $κ\nCO2 reduction = $(round(C′_E,digits=3)) ≥ $C_E",
        titlefontsize = 10,
        xticks=(L_ind, lines),
        xrotation = 90,
        ylabel=L"\bar{f}_l\,\mathrm{[MW]}",
        legend=false,
        alpha = 0.7,
        size=(800, 600),
        lw = 0)
end

function plot_storage_level(b_snt, b̄_sn, S, n, T, κ, C_E, κ′, C′_E)
    p = Plots.plot(legend=:outertopright,
             title = "Hourly storage levels by storage technology\nRenewables share = $(round(κ′,digits=3)) ≥ $κ\nCO2 reduction = $(round(C′_E,digits=3)) ≥ $C_E",
             titlefontsize = 10,
        )
    for s in S
        plot!(p, T, [b_snt[s, n, t] for t in T],
              xlabel=L"t",
              ylabel=L"b_{s,n,t}\,\mathrm{[MWh]}",
              labels="battery",
              alpha=0.7,
              size=(800, 600))
        plot!(p, T, [b̄_sn[s, n] for t in T], labels="battery capacity")
    end
    return p
end

function plot_storage_capacities(b̄_sn, N, region_n, κ, C_E, κ′, C′_E)
    # @show b̄_sn, N, region_n
    capacities = permutedims(b̄_sn)
    StatsPlots.groupedbar(capacities,
        bar_position = :stack,
        bar_width=0.7,
        ylabel=L"\bar{p}_{s,n}\,\mathrm{[MWh]}",
        xticks=(N, region_n),
        labels="battery",
        title = "Storage capacity by region\nRenewables share = $(round(κ′,digits=3)) ≥ $κ\nCO2 reduction = $(round(C′_E,digits=3)) ≥ $C_E",
        titlefontsize = 10,
        lw = 0,
        alpha=0.7)
end

function plot_loss_of_load(σ_nt, N, T, region_n, κ, C_E, κ′, C′_E)
    p = Plots.plot(
        legend=:outertopright,
        size=(800, 600),
        title = "Loss of load by region\nRenewables share = $(round(κ′,digits=3)) ≥ $κ\nCO2 reduction = $(round(C′_E,digits=3)) ≥ $C_E",
        titlefontsize = 10
    )
    for n in N
        plot!(p, T, [σ_nt[n, t] for t in T],
              xlabel=L"t",
              ylabel=L"\sigma_{n,t}\,\mathrm{[MWh]}",
              label="$(region_n[n])",
              alpha=0.7)
    end
    return p
end

function plot_box(p_gnt, h_hnt, G, n, region_n, technology_g, κ, C_E, κ′, C′_E)
    colors = techcolors
    p = StatsPlots.plot(
        size=(800, 600),
        xticks=([G;9], [technology_g;"hydro"]),
        ylabel=L"p_{g,n,t}\,\mathrm{[MWh]}",
        title = "Hourly generation dispatch by technology in $(region_n[n])\nRenewables share = $(round(κ′,digits=3)) ≥ $κ\nCO2 reduction = $(round(C′_E,digits=3)) ≥ $C_E",
        titlefontsize = 10
    )
    for g in G
        StatsPlots.boxplot!([g], [p_gnt[g,n,:]],
            color = colors[g],
            alpha=0.7,
            label=false
        )
    end
    StatsPlots.boxplot!([9], [h_hnt[1, n,:]],
        color = colors[9],
        alpha=0.7,
        label=false
    )
    return p
end

function plot_box_all(p_gnt, h_hnt, G, H, technology_g, κ, C_E, κ′, C′_E)
    colors = techcolors
    p = StatsPlots.plot(
        size=(800, 600),
        xticks=([G;9], [technology_g;"hydro"]),
        ylabel=L"\sum_n p_{g,n,t}\,\mathrm{[MWh]}",
        title = "Hourly generation dispatch by technology\nRenewables share = $(round(κ′,digits=3)) ≥ $κ\nCO2 reduction = $(round(C′_E,digits=3)) ≥ $C_E",
        titlefontsize = 10
    )
    for g in G
        p_gt = permutedims(sum(p_gnt[g, :, :], dims=1))
        StatsPlots.boxplot!([g], [p_gt],
            color = colors[g],
            alpha=0.7,
            label=false
        )
    end
    h_t = permutedims(sum(sum(h_hnt[h,:,:] for h in H), dims=1))
    StatsPlots.boxplot!([9], [h_t],
        color = colors[9],
        alpha=0.7,
        label=false
    )
    return p
end

function plot_dispatch_bars(max_dem_n, p_gnt, h_hnt, D_nt, N, T, H, region_n, technology_g, κ, C_E, κ′, C′_E)
    p_gn = sum(p_gnt[:,:,t] for t in T)
    h_n = sum(h_hnt[h,:,t] for h in H, t in T)
    D_n = sum(max_dem_n[:].*D_nt[:,t] for t in T)
    dispatches = permutedims([p_gn; permutedims(h_n)])
    groupedbar(dispatches,
        bar_position = :stack,
        bar_width=0.7,
        ylabel=L"\sum_t p_{g,n,t}\,\mathrm{[MWh]}",
        xticks=(N, region_n),
        title = "Generation dispatch by region and technology\n Renewables share = $(round(κ′,digits=3)) ≥ $κ\n CO2 reduction = $(round(C′_E,digits=3)) ≥ $C_E",
        titlefontsize = 10,
        color=techcolors,
        labels=permutedims([technology_g;"hydro"]),
        lw = 0,
        alpha=0.7,
        size = (800,600))
    bar!(N, D_n, bar_width=0.03, fillcolor=repeat([:black], length(N)), label="demand")
end

# Overload functions for signature: parameters::Params, model::Model, ...

"""Plot objective value and individual objective values."""
function plot_objective_values(objectives::Union{Dict{String, Float64}, Dict{String, Any}})
    fs = keys(objectives)
    vs = [objectives[i] for i in fs]
    title = string("Objective value: ", round(sum(vs), digits = 2))
    bar(vs,
        xlabel="Objective function",
        ylabel="EUR",
        title=title,
        legend=false,
        size = (800,600))
end

"""Plot generation dispatch."""
function plot_generation_dispatch(parameters::Params, variables::Union{Dict{String, Float64}, Dict{String, Array{Float64}}}, expressions::Union{Dict{String, Float64}, Dict{String, Any}}, n::Integer)
    plot_generation_dispatch(variables["p_gnt"], variables["h_hnt"], parameters.G, n, parameters.T, parameters.region_n,
                             parameters.technology_g, parameters.κ, parameters.C_E, expressions["κ′"], expressions["C′_E"])
end

"""Plot generation capacities."""
function plot_generation_capacities(parameters::Params, variables::Union{Dict{String, Float64}, Dict{String, Array{Float64}}}, expressions::Union{Dict{String, Float64}, Dict{String, Any}}, n::Integer)
    plot_generation_capacities(variables["p̄_gn"], variables["h̄_hn"], parameters.HRmax_n, parameters.G, n, parameters.region_n,
                               parameters.technology_g, parameters.κ, parameters.C_E, expressions["κ′"], expressions["C′_E"])
end

"""Plot stacked generation capacities as a stacked graph."""
function plot_generation_capacities_stacked(parameters::Params, variables::Union{Dict{String, Float64}, Dict{String, Array{Float64}}}, expressions::Union{Dict{String, Float64}, Dict{String, Any}})
    plot_generation_capacities_stacked(variables["p̄_gn"], variables["h̄_hn"], parameters.HRmax_n, parameters.N, parameters.H,
                       parameters.region_n, parameters.technology_g, parameters.κ, parameters.C_E, expressions["κ′"], expressions["C′_E"])
end

"""Plot transmission flow."""
function plot_transmission_flow(parameters::Params, variables::Union{Dict{String, Float64}, Dict{String, Array{Float64}}}, expressions::Union{Dict{String, Float64}, Dict{String, Any}}, l::Integer)
    plot_transmission_flow(variables["f_lt"], variables["f̄_l"], l, parameters.L, parameters.T, parameters.region_n, parameters.κ, parameters.C_E, expressions["κ′"], expressions["C′_E"])
end

"""Plot transmission bars."""
function plot_transmission_bars(parameters::Params, variables::Union{Dict{String, Float64}, Dict{String, Array{Float64}}}, expressions::Union{Dict{String, Float64}, Dict{String, Any}})
    plot_transmission_bars(variables["f_lt"], parameters.L, parameters.L_ind, parameters.T,
                           parameters.region_n, parameters.κ, parameters.C_E, expressions["κ′"], expressions["C′_E"])
end

"""Plot transmission capacities."""
function plot_transmission_capacities(parameters::Params, variables::Union{Dict{String, Float64}, Dict{String, Array{Float64}}}, expressions::Union{Dict{String, Float64}, Dict{String, Any}})
    plot_transmission_capacities(variables["f̄_l"], parameters.L, parameters.L_ind, parameters.region_n, parameters.κ, parameters.C_E, expressions["κ′"], expressions["C′_E"])
end

"""Plot storage level."""
function plot_storage_level(parameters::Params, variables::Union{Dict{String, Float64}, Dict{String, Array{Float64}}}, expressions::Union{Dict{String, Float64}, Dict{String, Any}}, n::Integer)
    plot_storage_level(variables["b_snt"], variables["b̄_sn"], parameters.S, n, parameters.T, parameters.κ, parameters.C_E, expressions["κ′"], expressions["C′_E"])
end

"""Plot storage capacities."""
function plot_storage_capacities(parameters::Params, variables::Union{Dict{String, Float64}, Dict{String, Array{Float64}}}, expressions::Union{Dict{String, Float64}, Dict{String, Any}})
    plot_storage_capacities(variables["b̄_sn"], parameters.N, parameters.region_n,
                            parameters.κ, parameters.C_E, expressions["κ′"], expressions["C′_E"])
end

"""Plot loss of load."""
function plot_loss_of_load(parameters::Params, variables::Union{Dict{String, Float64}, Dict{String, Array{Float64}}}, expressions::Union{Dict{String, Float64}, Dict{String, Any}})
    plot_loss_of_load(variables["σ_nt"], parameters.N, parameters.T, parameters.region_n, parameters.κ, parameters.C_E, expressions["κ′"], expressions["C′_E"])
end

"""Plot generation box."""
function plot_box(parameters::Params, variables::Union{Dict{String, Float64}, Dict{String, Array{Float64}}}, expressions::Union{Dict{String, Float64}, Dict{String, Any}}, n::Integer)
    plot_box(variables["p_gnt"], variables["h_hnt"], parameters.G, n, parameters.region_n,
             parameters.technology_g, parameters.κ, parameters.C_E, expressions["κ′"], expressions["C′_E"])
end

function plot_box_all(parameters::Params, variables::Union{Dict{String, Float64}, Dict{String, Array{Float64}}}, expressions::Union{Dict{String, Float64}, Dict{String, Any}})
    plot_box_all(variables["p_gnt"], variables["h_hnt"], parameters.G, parameters.H,
              parameters.technology_g, parameters.κ, parameters.C_E, expressions["κ′"], expressions["C′_E"])
end

function plot_dispatch_bars(parameters::Params, variables::Union{Dict{String, Float64}, Dict{String, Array{Float64}}}, expressions::Union{Dict{String, Float64}, Dict{String, Any}})
    plot_dispatch_bars(parameters.max_dem_n, variables["p_gnt"], variables["h_hnt"], parameters.D_nt, parameters.N, parameters.T, parameters.H,
                       parameters.region_n, parameters.technology_g, parameters.κ, parameters.C_E, expressions["κ′"], expressions["C′_E"])
end


text:
module amplitudes_compass
# using masses: mπ, mπ2
using QuadGK  # needed for waves notmalization
using DelimitedFiles
using DalitzPlotAnalysis: λ, change_basis, Z, WignerDϵ, WignerD, Wignerd, ClebschGordon

export swap_kin_parameters
export fρ, ff2, fρ3, fσ, ff0_980, ff0_1500, BlttWskpf
export compass_jmels_basis_psi
export get_wavelist, get_wavenames, get_wavebasis
export get_threshold_mask
#
export read_compass_SDM

mπ = 0.13956755; mπ2 = mπ^2;
mK = 0.493677; mK2 = mK^2;
mK0 = 0.497614; mK02 = mK0^2;

export mπ, mπ2

BlttWskpf = [z->z/(1+z),
                   z->z^2/(9+3z+z^2),
                   z->z*z*z/(z*z*z+6z*z+45z+225),
                   z->z*z*z*z/(z*z*z*z+10z*z*z+135z*z+1575z+11025),
                   z->z*z*z*z*z/(z*z*z*z*z+15z*z*z*z+315z*z*z+6300z*z+99225z+893025),
                   z->z*z*z*z*z*z/(z*z*z*z*z*z+21z*z*z*z*z+630z*z*z*z+18900z*z*z+496125z*z+9823275z+108056025)];

################### SIGMA #################
# AMP Table 1, M solution: f_2^2
# AMP Table 1, M solution: f_1^1 and f_2^1
# Last item is Katchaev modification
_a = [0.1131, 0.0]

# AMP Table 1, M solution: c_11^0
# last item - Katchaev modification
_c = [0.0337, -0.3185, -0.0942, -0.5927, 0.0]
_sP = [-0.0074, 0.9828] .+ 1e-7 # for numerical stability

function _fσ(s::Number)
    (s < 4mπ2) && return 0.0;

    rho00 = sqrt(λ(s, mπ2, mπ2))/s;
    scale = (s / (4 * ((mK+mK0)/2.)^2)) - 1.;

    M00 = 0.0;
    for i = 1:length(_sP)
        M00 += _a[i] / (s - _sP[i]);
    end
    for i = 1:length(_c)
        M00 += scale^(i-1) *_c[i];
    end

  return 1.0/(M00 - 1im*rho00);
end

function _ff0_980(s::Number)
    m = 0.965;
    gPi = 0.165;
    rK = 4.21;
    qsq = λ(s,mπ2,mπ2)/(4*s);
    qsqK = λ(s,mK2,mK2)/(4*s);
    return 1.0/(m*m-s-2.0im*(sqrt(qsq)*gPi+sqrt(convert(Complex{Float64},qsqK))*gPi*rK)/sqrt(s));
end

function _ff0_1500(s::Number)
    m = 1.507;
    G = 0.109;
    return m*G/(m*m-s-1im*m*G);
end

function _fρ(s::Number)
    mρ = 0.7685; Γρ = 0.1507;
    # break up momentum
    p  = sqrt(λ(mπ2,mπ2,s)/(4*s));
    p0 = sqrt(λ(mπ2,mπ2,mρ^2))/(2*mρ);
    # ph.sp.
    ρ  = 1/(8*π)*2*p/sqrt(s)
    ρ0 = 1/(8*π)*2*p0/mρ
    # extra factor due to the spin of ρ
    R = 4.94 # it was 5
    ff = p^2/p0^2*(1.0/R^2+p0^2)/(1.0/R^2+p^2)
    mΓ = mρ*Γρ*p/p0*ff  #changed from ρ/ρ0 to p/p0
    return sqrt(ff)/(mρ^2-s -1.0im*mΓ) #
end

function _ff2(s::Number)
    mπ = 0.13956755;
    mπ2 = mπ^2;
    m = 1.274;  # it was 1.2754;
    G = 0.185;  # it was 0.1852;
    qsq_R = 1.0/4.94^2; # it was 5
    qsq = λ(s,mπ2, mπ2)/(4*s);
    qsq0 = λ(m^2,mπ2, mπ2)/(4*m^2);
    ff = sqrt(qsq/qsq0)*m/sqrt(s)*BlttWskpf[2](qsq/qsq_R)/BlttWskpf[2](qsq0/qsq_R);
    return m*G/(m*m-s-1im*m*G*ff) * sqrt(BlttWskpf[2](qsq/qsq_R));
end

#
function _fρ3(s::Number)
  mπ = 0.13956755; mπ2 = mπ^2;
  m = 1.690;
  G = 0.190;
  qsq_R = 1/4.94^2;
  qsq = λ(s,mπ2, mπ2)/(4*s);
  return sqrt(m*G*sqrt(s))/(m*m-s-1im*m*G) * sqrt(BlttWskpf[3](qsq/qsq_R));
end

for name in ["fρ3", "ff2", "fρ", "ff0_1500", "fσ", "ff0_980"]
    @eval $(Symbol("norm_"*name)) = sqrt(quadgk(
        x->abs2($(Symbol("_"*name))(x))*sqrt((1-4mπ2/x)),4mπ2, Inf)[1]/(16*π^2))
    @eval $(Symbol(name))(σ) = $(Symbol("_"*name))(σ)/$(Symbol("norm_"*name))
end

function ZϵDD(J::Int64,M::Int64,ϵ::Bool,L::Int64,l::Int64, Dϵ, D)
    rng_lm = min(l,J);
    # println("$J,$M,$ϵ")
    sum(ClebschGordon(L,0,l,lm,J,lm)*sqrt((2*L+1)*(2*l+1))*
                Dϵ[J+1,M+1,1+ϵ,4+lm]*D[1+l,4+lm]
        for lm=-min(l,J):min(l,J))
end

isobarsV = [fσ,fρ,ff2,fρ3]
isobarsS = [fσ,ff0_980,ff0_1500]

# constract COMPASS basis
function get_wavelist(path_to_wavelist; path_to_thresholds="", M3pi=3.0)
    (! isfile(path_to_wavelist)) && error("Cannot find $(path_to_wavelist)!")

    wavesInFile = readdlm(path_to_wavelist)

    !(isfile(path_to_thresholds)) && warn("Do not consider thresholds!")
    thresholds_filter = isfile(path_to_thresholds) ?
        get_threshold_mask(path_to_thresholds, M3pi, size(wavesInFile,1)) : fill(true,  size(wavesInFile,1))

    wavesfile = wavesInFile[thresholds_filter,:]

    return wavesfile
end

function get_wavenames(wavelist::Array)
    return vcat(wavelist[:,2]...)::Vector{SubString{String}}
end

function get_wavenames(path_to_wavelist; path_to_thresholds="", M3pi=3.0)
    wavelist = get_wavelist(path_to_wavelist; path_to_thresholds=path_to_thresholds, M3pi=M3pi);
    get_wavenames(wavelist)
end

function get_threshold_mask(path_to_thresholds, M3pi, returnedSize)
    thresholds = fill(0.0 , returnedSize)
    thf = readdlm(path_to_thresholds)
    for v in zip(thf[:,1],thf[:,2])
        thresholds[Int64(v[1])] = v[2]
    end
    thresholds_mask = thresholds .< M3pi
    return thresholds_mask
end

function get_wavebasis(wavelist::Array)
    basis = []
    for i in 1:size(wavelist,1)
        wl = wavelist[i,:]
        # special case
        if wl[2] == "FLAT"
            flat(σ1,cosθ1,ϕ1,cosθ23,ϕ23,s) = 1.0+0.0im
            push!(basis,flat)
            continue
        end
        # special case
        wn, name, J, P, M, ϵ, S, L = wl
        push!(basis, get_compass_jmels_basis_psi(J,M,1*(P==ϵ),L,S))
    end
    return basis
end

function get_compass_jmels_basis_psi(J,M,Pϵ,L,S)
    function wave(σ1,cosθ1,ϕ1,cosθ23,ϕ23,s)
        compass_jmels_basis_psi(QNs=(J,M,Pϵ,L,S),
                                τ1=(σ1,cosθ1,ϕ1,cosθ23,ϕ23),
                                s=s)
    end
    return wave
end

function compass_jmels_basis_psi(;QNs::NTuple{5,Int}=error("give quantum numbers for the wave (J,M,Pϵ,L,S)"),
                                   τ1::NTuple{5,Float64}=error("give kinematic variables (σ1,cosθ1,ϕ1,cosθ23,ϕ23)"),
                                   s::Float64=error("give s=M3π^2"))
    J,M,Pϵ,L,S = QNs
    m1sq, m2sq, m3sq = fill(mπ2,3)
    ###############
    fi = (S ≥ 0) ? isobarsV[S+1] : isobarsS[1-S]
    S = (S ≥ 0 ? S : 0)
    ###############
    τ3 = Array{Float64, 1}(undef, 5)
    τ3 .= change_basis(τ1...,mπ2,mπ2,mπ2,s)
    (abs(τ3[2]) ≈ 1.0) && (τ3[2] = sign(τ3[1])*1.0);
    (abs(τ3[4]) ≈ 1.0) && (τ3[4] = sign(τ3[3])*1.0);
    (abs(τ3[2]) > 1.0 || abs(τ3[4]) > 1.0) && error("Something is wrong! (abs($(τ3[1])) > 1.0 || abs($(τ3[3]) > 1.0)")
    #############
    σ1 = τ1[1]
    σ3 = τ3[1]
    R = 1/0.2024;
    bw1 = (L == 0) ? 1.0 : BlttWskpf[L](λ(s,σ1,m1sq)/(4s)*R^2)
    bw3 = (L == 0) ? 1.0 : BlttWskpf[L](λ(s,σ3,m3sq)/(4s)*R^2)
    #############
    # τ1_rev = collect(τ1); τ1_rev[4] *= -1; τ1_rev[5] += π  # compass `convension` (inconsistent, btw)
    return Z(J,M,Pϵ,L,S,τ1[2],τ1[3],-τ1[4],τ1[5]+π)*fi(σ1)*sqrt(bw1) +
           Z(J,M,Pϵ,L,S,τ3[2],τ3[3], τ3[4],τ3[5])*fi(σ3)*sqrt(bw3)
end

# It is a bit sloppy to repeat the same piece of threshold reading twice
function read_compass_SDM(path_inclusing_sdmXXdot; path_to_wavelist="", path_to_thresholds="", M3pi=3.0)
    (! isfile(path_to_wavelist)) && error("Cannot find $(path_to_wavelist)!")

    wavesInFile = readdlm(path_to_wavelist)

    thresholds = fill(0.0,size(wavesInFile,1))
    if isfile(path_to_thresholds)
        thf = readdlm(path_to_thresholds)
        for v in zip(thf[:,1],thf[:,2])
            thresholds[Int64(v[1])] = v[2]
        end
    else
        warn("Do not consider thresholds!")
    end
    thresholds_filter = thresholds .< M3pi
    SDM = readdlm(path_inclusing_sdmXXdot*"re") + 1im*readdlm(path_inclusing_sdmXXdot*"im")

    return SDM[thresholds_filter,thresholds_filter]
end

function swap_kin_parameters(s,τ1...)
    τ3 = collect(change_basis(τ1...,mπ2,mπ2,mπ2,s))
    τ3[4] *= -1 # angle of the oposite particle
    τ3[5] += π # angle of the oposite particle
    (τ3[5] > π) && (τ3[5] -= 2π)
    vcat([s],τ3)
end

end


text:
function _handle_pmin(pwf_data::Dict, bus_i::Int, dict_dger)
    if haskey(pwf_data, "DGER") && haskey(dict_dger, bus_i)
        bus = dict_dger[bus_i]
        return bus["MINIMUM ACTIVE GENERATION"]
    end
    return 0.0
end

function _handle_pmax(pwf_data::Dict, bus_i::Int, dict_dger)
    if haskey(pwf_data, "DGER") && haskey(dict_dger, bus_i)
        bus = dict_dger[bus_i]
        return bus["MAXIMUM ACTIVE GENERATION"]
    end
    return 99999.0
end

function _create_dict_dger(data::Dict)
    dict_dger = Dict{Int, Any}()

    for (k,v) in data
        sub_data = Dict{String, Any}()
        sub_data["MINIMUM ACTIVE GENERATION"] = v["MINIMUM ACTIVE GENERATION"]
        sub_data["MAXIMUM ACTIVE GENERATION"] = v["MAXIMUM ACTIVE GENERATION"]
        dict_dger[v["NUMBER"]] = sub_data
    end
    return dict_dger
end

function _pwf2pm_generator!(pm_data::Dict, pwf_data::Dict, bus::Dict)
    sub_data = Dict{String,Any}()

    sub_data["gen_bus"] = bus["NUMBER"]
    sub_data["gen_status"] = 1
    sub_data["pg"] = bus["ACTIVE GENERATION"]
    sub_data["qg"] = bus["REACTIVE GENERATION"]
    sub_data["vg"] = pm_data["bus"]["$(bus["NUMBER"])"]["vm"]
    sub_data["mbase"] = _handle_base_mva(pwf_data)

    dict_dger = haskey(pwf_data, "DGER") ? _create_dict_dger(pwf_data["DGER"]) : nothing
    sub_data["pmin"] = _handle_pmin(pwf_data, bus["NUMBER"], dict_dger)
    sub_data["pmax"] = _handle_pmax(pwf_data, bus["NUMBER"], dict_dger)

    sub_data["qmin"] = haskey(bus, "MINIMUM REACTIVE GENERATION") ? bus["MINIMUM REACTIVE GENERATION"] : 0.0
    sub_data["qmax"] = haskey(bus, "MAXIMUM REACTIVE GENERATION") ? bus["MAXIMUM REACTIVE GENERATION"] : 0.0

    # Default Cost functions
    sub_data["model"] = 2
    sub_data["startup"] = 0.0
    sub_data["shutdown"] = 0.0
    sub_data["ncost"] = 2
    sub_data["cost"] = [1.0, 0.0]

    sub_data["source_id"] = ["generator", sub_data["gen_bus"], "1 "]
    sub_data["index"] = length(pm_data["gen"]) + 1

    idx = string(sub_data["index"])
    pm_data["gen"][idx] = sub_data
end

function _pwf2pm_generator!(pm_data::Dict, pwf_data::Dict)

    if !haskey(pwf_data, "DGER")
        @warn("DGER not found, setting pmin as 0.0 MW and pmax as 99999.0 MW")
    end

    pm_data["gen"] = Dict{String, Any}()
    if haskey(pwf_data, "DBAR")
        for (i,bus) in pwf_data["DBAR"]
            if bus["ACTIVE GENERATION"] != 0.0 || bus["REACTIVE GENERATION"] != 0.0 || bus["TYPE"] in [bus_type_raw_to_pwf[bus_type_str_to_num["PV"]], bus_type_raw_to_pwf[bus_type_str_to_num["Vθ"]]]
                _pwf2pm_generator!(pm_data, pwf_data, bus)
            end
        end
    end
end


text:
# This file is a part of HistTools.jl, licensed under the MIT License (MIT).

import Test

Test.@testset "Package HistTools" begin
    include("test.jl")
end # testset