text:
"""
    nw_align(s1::String, s2::String; edge_reduction = 0.99)

Returns aligned strings using the Needleman-Wunch Algorithm (quadratic),
with end gaps penalized slightly less. edge_reduction is a multiplier (usually
less than one) on gaps on end of strings.
"""
function nw_align(s1::String, s2::String; edge_reduction = 0.99, mismatch_cost = -1.0)
    # edge_reduction is a multiplicative score that gets multiplied to
    # the penalties along the edges, to prefer terminal gaps.
    ins_cost = -1.0  # -> : horizontal
    del_cost = -1.0  # V : vertical
    match_cost = 1.0

    s1arr = collect(s1)  # vertical
    s2arr = collect(s2)  # horizontal
    arr = zeros(length(s1arr)+1, length(s2arr)+1)
    traceArr = zeros(Int, length(s1arr), length(s2arr))

    # this will need to be generalized when we want to allow overhang
    arr[:, 1] = edge_reduction*del_cost*(0:length(s1arr))
    arr[1,:] = edge_reduction*ins_cost*(0:length(s2arr))

    for i in 2:length(s1arr)+1
        for j in 2:length(s2arr)+1
            if s1arr[i-1] == s2arr[j-1]
                diag = arr[i-1, j-1] + match_cost
            else
                diag = arr[i-1, j-1] + mismatch_cost
            end

            # to handle the lower edge penalties.
            delMult = (i == length(s1arr)+1) ? edge_reduction : 1
            insMult = (j == length(s2arr)+1) ? edge_reduction : 1

            ins = arr[i-1, j]+(ins_cost*insMult)
            del = arr[i, j-1]+(del_cost*delMult)
            scores = [diag, del, ins]
            best = indmax(scores)
            arr[i, j] = scores[best]
            traceArr[i-1, j-1] = best
        end
    end
    alignedScore = arr[end, end]

    # the trace endpoing will need to be generalized if we want to
    # allow overhang.
    trI, trJ = length(s1arr), length(s2arr)

    # First compute the trace running backwards. Initialized to the
    # maximum size. With unit penalties, is it possible to tell how
    # big this should be in advance?
    backtrace = Array{Int}(length(s1arr) + length(s2arr))
    # If you hit any boundary, you have to run all the way to the side.
    btInd = 1
    while (trI > 0) && (trJ > 0)
        backtrace[btInd] = traceArr[trI, trJ]
        if backtrace[btInd] == 1
            trI += -1
            trJ += -1
        elseif backtrace[btInd] == 2
            trJ += -1
        elseif backtrace[btInd] == 3
            trI += -1
        end
        btInd += 1
    end
    # If you hit the boundaries not at the top left corner.
    while trI > 0
        backtrace[btInd] = 3
        btInd += 1
        trI += -1
    end
    while trJ > 0
        backtrace[btInd] = 2
        btInd += 1
        trJ += -1
    end

    backtrace = backtrace[1:btInd-1]
    ali1arr = Array{Char}(length(backtrace))
    ali2arr = Array{Char}(length(backtrace))
    ind1 = 1
    ind2 = 1
    for i in 1:length(backtrace)
        tr = backtrace[length(backtrace)-(i-1)]
        if tr == 1
            ali1arr[i] = s1arr[ind1]
            ali2arr[i] = s2arr[ind2]
            ind1 += 1
            ind2 += 1
        elseif tr == 2
            ali2arr[i] = s2arr[ind2]
            ali1arr[i] = '-'
            ind2 += 1
        elseif tr == 3
            ali1arr[i] = s1arr[ind1]
            ali2arr[i] = '-'
            ind1 += 1
        end
    end
    return join(ali1arr), join(ali2arr)
end

#--------Banded alignment---------

"""
    nw_align(s1::String, s2::String, banded::Float64)

Wrapper for `nw_align` and `banded_nw_align`. A larger `banded` value makes alignment slower
but more accurate.
"""
function nw_align(s1::String, s2::String, banded::Float64)
    if banded <= 0
        return nw_align(s1, s2)
    else
        return banded_nw_align(s1, s2, band_coeff = banded)
    end
end

# methods for converting from rectangular coordinates to diagonal band coordinates and back

"""Sets value in band at `i`, `j` to `val`, where `i` and `j` are in square matrix coords, `dim_diff` = ncols - nrows"""
function add_to_band!(band, val, i::Int, j::Int, bandwidth::Int, dim_diff::Int)
    # smallest coord is slice, largest is +/- offset from middle of slice + offset for diff in string lengths
    ii = min(i, j)
    jj = j - i + bandwidth + 1 + max(0, dim_diff)
    band[ii, jj] = val
end

"""Returns value from band where `i` and `j` are in square matrix coords, `dim_diff` = ncols - nrows"""
function get_band_val(band, i::Int, j::Int, bandwidth::Int, dim_diff::Int)
    # smallest coord is slice, largest is +/- offset from middle of slice + offset for diff in string lengths
    ii = min(i, j)
    jj = j - i + bandwidth + 1 + max(0, dim_diff)
    return band[ii, jj]
end

"""Checks if given square matrix coords are within width of band (ignore length)"""
function in_band(i, j, bandwidth, dim_diff)
    ii = min(i, j)
    jj = j - i + bandwidth + 1 + max(0, dim_diff)
    return !( ii < 1 || jj < 1 || jj > 2*bandwidth + 1 + abs(dim_diff) )
end

"""
    banded_nw_align(s1::String, s2::String; edge_reduction = 0.99, band_coeff = 1)

Like nw_align, but sub quadratic by only computing values within a band around the center diagonal.
One 'band' of radius 3 = (4,1), (3,1), (2,1), (1,1), (1,2), (1,3), (1,4), aka upside-down L shape.
band_coeff = 1 is sufficient to get same alignments as nw_align for 10% diverged sequences ~97% of the time;
increase this value for more conservative alignment with longer computation time.
Radius of band = `bandwidth` = `band_coeff` * sqrt(avg seq length)
"""
function banded_nw_align(s1::String, s2::String; edge_reduction = 0.99, band_coeff = 1)
    # calculate band width
    avg_len = (length(s1) + length(s2)) / 2
    bandwidth = Int(ceil(band_coeff * sqrt(avg_len)))
    # edge_reduction is a multiplicative score that gets multiplied to
    # the penalties along the edges, to prefer terminal gaps.
    ins_cost = -1.0  # -> : horizontal
    del_cost = -1.0  # V : vertical
    mismatch_cost = -1.0
    match_cost = 1.0

    s1arr = collect(s1)  # vertical
    s2arr = collect(s2)  # horizontal

    # dimension difference added to band length on one side
    # to give padding for large sequence length differences
    dim_diff = length(s1arr) - length(s2arr)
    maxlen = max(length(s1), length(s2))
    arr = zeros(maxlen + 1, 2*bandwidth + 1 + abs(dim_diff))
    traceArr = zeros(Int, size(arr))
    for i in 1:bandwidth+1
        add_to_band!(arr, edge_reduction*del_cost*(i-1), i, 1, bandwidth, dim_diff)
        add_to_band!(arr, edge_reduction*ins_cost*(i-1), 1, i, bandwidth, dim_diff)
    end

    for i in 2:length(s1arr)+1
        # keep j within 2*bandwidth + 1 band of i, but keep bandwidth sized buffer on edges
        lo = 0
        # to keep a big enough padding around the edges when there's a big length mismatch
        if i <= bandwidth + 1 + max(0, dim_diff)
            lo = 2
        else
            lo = i - bandwidth - max(0, dim_diff)
        end
        for j in lo:length(s2arr)+1
            # keep j within band
            if !(in_band(i, j, bandwidth, dim_diff))
                break
            end

            if in_band(i-1, j-1, bandwidth, dim_diff) && s1arr[i-1] == s2arr[j-1]
                diag = get_band_val(arr, i-1, j-1, bandwidth, dim_diff) + match_cost
            else
                diag = get_band_val(arr, i-1, j-1, bandwidth, dim_diff) + mismatch_cost
            end

            # to handle the lower edge penalties.
            delMult = (i == length(s1arr)+1) ? edge_reduction : 1
            insMult = (j == length(s2arr)+1) ? edge_reduction : 1

            ins = in_band(i-1, j, bandwidth, dim_diff) ?
                    (get_band_val(arr, i-1, j, bandwidth, dim_diff)+(ins_cost*insMult)) : NaN
            del = in_band(i, j-1, bandwidth, dim_diff) ?
                    (get_band_val(arr, i, j-1, bandwidth, dim_diff)+(del_cost*delMult)) : NaN

            scores = [diag, del, ins]
            best = indmax(scores)
            add_to_band!(arr, scores[best], i, j, bandwidth, dim_diff)
            add_to_band!(traceArr, best, i-1, j-1, bandwidth, dim_diff)
        end
    end
    alignedScore = get_band_val(arr, length(s1arr)+1, length(s2arr)+1, bandwidth, dim_diff)

    trI, trJ = length(s1arr), length(s2arr)
    backtrace = Array{Int}(length(s1arr) + length(s2arr))
    btInd = 1
    while (trI > 0) && (trJ > 0)
        backtrace[btInd] = get_band_val(traceArr, trI, trJ, bandwidth, dim_diff)
        if backtrace[btInd] == 1
            trI += -1
            trJ += -1
        elseif backtrace[btInd] == 2
            trJ += -1
        elseif backtrace[btInd] == 3
            trI += -1
        else
            error("Bad trace value: $(backtrace[btInd]): ($trI, $trJ)")
        end
        btInd += 1
    end
    # If you hit the boundaries not at the top left corner.
    while trI > 0
        backtrace[btInd] = 3
        btInd += 1
        trI += -1
    end
    while trJ > 0
        backtrace[btInd] = 2
        btInd += 1
        trJ += -1
    end

    backtrace = backtrace[1:btInd-1]
    ali1arr = Array{Char}(length(backtrace))
    ali2arr = Array{Char}(length(backtrace))
    ind1 = 1
    ind2 = 1
    for i in 1:length(backtrace)
        tr = backtrace[length(backtrace)-(i-1)]
        if tr == 1
            ali1arr[i] = s1arr[ind1]
            ali2arr[i] = s2arr[ind2]
            ind1 += 1
            ind2 += 1
        elseif tr == 2
            ali2arr[i] = s2arr[ind2]
            ali1arr[i] = '-'
            ind2 += 1
        elseif tr == 3
            ali1arr[i] = s1arr[ind1]
            ali2arr[i] = '-'
            ind1 += 1
        end
    end
    return join(ali1arr), join(ali2arr)
end

#----Triplet/Codon Alignment----

"""
    triplet_nw_align(s1::String, s2::String; edge_reduction = 0.99, boundary_mult = 2)

Returns alignment of two sequences where `s1` is a reference with reading frame to be preserved and `s2` is a query sequence.
`boundary_mult` adjusts penalties for gaps preserving the reading frame of `s1`.
This usually works best on range 0 to 3, higher values for more strongly enforced gaps aligned on
reference frame (divisible-by-3 indices)
"""
function triplet_nw_align(s1::String, s2::String; edge_reduction = 0.99, boundary_mult = 2)
    # edge_reduction is a multiplicative score that gets multiplied to
    # the penalties along the edges, to prefer terminal gaps.
    ins_cost = -1.0  # -> : horizontal
    del_cost = -1.0  # V : vertical

    triple_ins_cost = -0.5
    triple_del_cost = -0.5
    triple_ins_score(qv) = triple_ins_cost
    triple_del_score(qv) = triple_del_cost

    mismatch_cost = -0.0
    match_cost = 1.0

    s1arr = collect(s1)  # vertical
    s2arr = collect(s2)  # horizontal
    arr = zeros(length(s1arr)+1, length(s2arr)+1)
    traceArr = zeros(Int, length(s1arr), length(s2arr))

    # this will need to be generalized when we want to allow overhang
    arr[:, 1] = edge_reduction*del_cost*(0:length(s1arr))
    arr[1, :] = edge_reduction*ins_cost*(0:length(s2arr))

    for i in 2:length(s1arr)+1
        for j in 2:length(s2arr)+1
            if s1arr[i-1] == s2arr[j-1]
                diag = arr[i-1, j-1] + match_cost
            else
                diag = arr[i-1, j-1] + mismatch_cost
            end

            # to handle the lower edge penalties.
            delMult = (i == length(s1arr)+1) ? edge_reduction : 1
            insMult = (j == length(s2arr)+1) ? edge_reduction : 1
            del_bndry_mult = ((i-1) % 3 == 0) ? -boundary_mult : 1
            ins_bndry_mult = ((j-1) % 3 == 0) ? -boundary_mult : 1

            ins = arr[i-1, j]+(ins_cost*insMult)
            del = arr[i, j-1]+(del_cost*delMult)
            # TODO: factor in qv scores
            s1qv = NaN
            s2qv = NaN
            triple_ins = i > 3 ? (arr[i-3, j]+(triple_ins_score(s1qv) * ins_bndry_mult)) : NaN
            triple_del = j > 3 ? (arr[i, j-3]+(triple_del_score(s2qv) * del_bndry_mult)) : NaN
            scores = [diag, del, ins, triple_del, triple_ins]
            best = indmax(scores)
            arr[i, j] = scores[best]
            traceArr[i-1, j-1] = best
        end
    end
    alignedScore = arr[end, end]

    trI, trJ = length(s1arr), length(s2arr)
    # First compute the trace running backwards
    backtrace = Array{Int}(length(s1arr) + length(s2arr))
    btInd = 1
    len = 0
    while (trI > 0) && (trJ > 0)
        backtrace[btInd] = traceArr[trI, trJ]
        if backtrace[btInd] == 1
            trI += -1
            trJ += -1
            len += 1
        elseif backtrace[btInd] == 2
            trJ += -1
            len += 1
        elseif backtrace[btInd] == 3
            trI += -1
            len += 1
        elseif backtrace[btInd] == 4
            trJ += -3
            len += 3
        elseif backtrace[btInd] == 5
            trI += -3
            len += 3
        else
            error("uh oh: $(traceArr[trI, trJ])")
        end
        btInd += 1
    end
    # If you hit the boundaries not at the top left corner.
    while trI > 0
        backtrace[btInd] = 3
        btInd += 1
        trI += -1
        len += 1
    end
    while trJ > 0
        backtrace[btInd] = 2
        btInd += 1
        trJ += -1
        len += 1
    end

    backtrace = backtrace[1:btInd-1]
    ali1arr = Array{Char}(len+3)
    ali2arr = Array{Char}(len+3)
    ind1 = 1
    ind2 = 1
    i = 1
    for tr in reverse(backtrace)
        if tr == 1
            ali1arr[i] = s1arr[ind1]
            ali2arr[i] = s2arr[ind2]
            ind1 += 1
            ind2 += 1
            i += 1
        elseif tr == 2
            ali2arr[i] = s2arr[ind2]
            ali1arr[i] = '-'
            ind2 += 1
            i += 1
        elseif tr == 3
            ali1arr[i] = s1arr[ind1]
            ali2arr[i] = '-'
            ind1 += 1
            i += 1
        elseif tr == 4
            ali2arr[i:i+2] = s2arr[ind2:ind2+2]
            ali1arr[i:i+2] = ['-', '-', '-']
            ind2 += 3
            i += 3
        elseif tr == 5
            ali1arr[i:i+2] = s1arr[ind1:ind1+2]
            ali2arr[i:i+2] = ['-', '-', '-']
            ind1 += 3
            i += 3
        else error("uh oh: tr = $tr")
        end
    end
    # to fix: weird boundary checking -- sometimes get off by 1 in length
    return join(ali1arr[1:i-1]), join(ali2arr[1:i-1])
end

codon_nw_align = triplet_nw_align

#-------Local Alignment--------

"""
    local_align(ref::String, query::String; mismatch_score = -1,
                match_score = 1, gap_penalty = -1,
                rightaligned=true, refend = false)

Aligns a query sequence locally to a reference. If true, `rightaligned` keeps the
right ends of each sequence in final alignment- otherwise they are trimmed;
`refend` keeps the beginning/left end of `ref`.
If you want to keep both ends of both strings, use nw_align.
For best alignments use the default score values.
"""
function local_align(ref::String, query::String; mismatch_score = -1,
                     match_score = 1, gap_penalty = -1,
                     rightaligned=true, refend = false)
    s1 = ref
    s2 = query
    # Populate Scores
    scoremat = zeros(length(s1)+1, length(s2)+1)
    scorematdirs = zeros(length(s1), length(s2))
    hiscore = 0
    hiscore_ij = (0, 0)
    for i in 2:size(scoremat)[1]
        for j in 2:size(scoremat)[2]
            # indexes indications: 1 = left and up, 2 = up, 3 = left, 4 = terminate (0)
            score_ij = [
                        scoremat[i-1, j-1] + (s1[i-1] == s2[j-1] ? match_score : mismatch_score),
                        scoremat[i-1, j] + gap_penalty,
                        scoremat[i, j-1] + gap_penalty,
                        0
                    ]
            # record index of max in matrix, track current value and location of high score
            ind = indmax(score_ij)
            scoremat[i, j] = score_ij[ind]
            scorematdirs[i-1, j-1] = ind
            if score_ij[ind] > hiscore
                hiscore = score_ij[ind]
                hiscore_ij = (i-1, j-1)
            end
        end
    end

    # Traceback Path
    i, j = size(scorematdirs)
    if !rightaligned
        i, j = hiscore_ij
    end
    backtrace = Array{UInt8}(length(s1) + length(s2))
    step = 0
    while i > 0 && j > 0
        direction = scorematdirs[i, j]
        if direction == 4
            break
        end
        step += 1
        backtrace[step] = direction
        if direction == 1
            # both chars
            i -= 1
            j -= 1
        elseif direction == 2
            # s1 char, s2 gap
            i -= 1
        else  # direction == 3
            # s1 gap, s2 char
            j -= 1
        end
    end
    # run along edges to end
    while refend && i > 0
        step += 1
        backtrace[step] = 2
        i -= 1
    end

    # Construct Strings
    ali1arr = Array{Char}(step)
    ali2arr = Array{Char}(step)
    ind1 = 0
    ind2 = 0
    for k in 1:step
        dir = backtrace[k]
        if dir == 1
            ali1arr[k] = s1[end-ind1]
            ali2arr[k] = s2[end-ind2]
            ind1 += 1
            ind2 += 1
        elseif dir == 2
            ali1arr[k] = s1[end-ind1]
            ali2arr[k] = '-'
            ind1 += 1
        elseif dir == 3
            ali2arr[k] = s2[end-ind2]
            ali1arr[k] = '-'
            ind2 += 1
        else
            error("error in local_align: alignment score matrix has invalid values")
        end
    end
    return join(reverse(ali1arr)), join(reverse(ali2arr))
end

#--------Kmer_seeded_align internals---------

"""
Sets value of key in dictionary to given index, unless the key already exists,
in which case value is set to -1.
"""
function unique_key(dicto::Dict{String, Int}, keyo::String, indo::Int)
    if haskey(dicto, keyo)
        dicto[keyo] = -1
    else
        dicto[keyo] = indo
    end
end

# Obsolete
function clean_unique_key(dicto::Dict{String, Int})
    for i in keys(dicto)
        if dicto[i] == -1
            delete!(dicto, i)
        end
    end
end

"""
Returns a list of indices of matches of unique words between `s1` and `s2`,
sorted by index of `s1`.
"""
function sorted_matches(s1, s2, wordlength, skip, aligncodons)
    word_dict1 = Dict{String, Int}()
    word_dict2 = Dict{String, Int}()
    # we can space one of these out, but not both
    for i in 1:skip:(length(s1)-(wordlength-1))  # bounds checked
        if aligncodons
            unique_key(word_dict1, translate_to_aa(s1[i:i+(wordlength-1)]), i)
        else
            unique_key(word_dict1, s1[i:i+(wordlength-1)], i)  # bounds checked
        end
    end
    for i in 1:(length(s2)-(wordlength-1))  # bounds checked
        if aligncodons
            unique_key(word_dict2, translate_to_aa(s2[i:i+(wordlength-1)]), i)
        else
            unique_key(word_dict2, s2[i:i+(wordlength-1)], i)  # bounds checked
        end
    end

    intersection = Dict{String, UInt8}()
    for (word, ind) in word_dict1
        if ind > 0 && haskey(word_dict2, word) && word_dict2[word] > 0
            intersection[word] = 0
        end
    end
    common = collect(keys(intersection))
    matches = zeros(Int, length(common), 2)
    for i in 1:length(common)
        matches[i, 1] = word_dict1[common[i]]
        matches[i, 2] = word_dict2[common[i]]
    end
    return sortrows(matches, by=x->(x[1]))
end

"""
Returns a list of indices of matches of unique words between `s1` and `s2`,
sorted by index of `s1`. Words match by amino acid encoding, in any reading frame.
"""
function sorted_aa_matches(str1, str2, wordlength)
    word_dict1 = Dict{String, Int}()
    word_dict2 = Dict{String, Int}()
    # get amino acid sequence in each reference frame
    s1s = generate_aa_seqs(str1)
    s2s = generate_aa_seqs(str2)
    for (offset, s1) in enumerate(s1s)
        for i in 1:(length(s1) - (div(wordlength, 3) - 1))
            # offset acounts for the shift in index due to being in a different reference frame.
            unique_key(word_dict1, s1[i:i+(div(wordlength, 3) - 1)], i*3 + offset - 3)
        end
    end
    for (offset, s2) in enumerate(s2s)
        for i in 1:(length(s2) - (div(wordlength, 3) - 1))
            unique_key(word_dict2, s2[i:i+(div(wordlength, 3) - 1)], i*3 + offset - 3)
        end
    end
    intersection = Dict{String, UInt8}()
    for (word, ind) in word_dict1
        if ind > 0 && haskey(word_dict2, word) && word_dict2[word] > 0
            intersection[word] = 0
        end
    end
    common = collect(keys(intersection))
    matches = zeros(Int, length(common), 2)
    for i in 1:length(common)
        matches[i, 1] = word_dict1[common[i]]
        matches[i, 2] = word_dict2[common[i]]
    end
    return sortrows(matches, by=x->(x[1]))
end

"""
Returns true if indices of word matches in second sequence (second column of
`matches`) is strictly increasing, else returns false.
"""
function matches_are_inconsistent(matches)
    # `matches` must already be sorted
    if size(matches)[1] <= 1
        return false
    end
    return minimum(matches[2:length(matches[:, 2]), 2] -
                   matches[1:(length(matches[:, 2]) - 1), 2]) < 1
end

"""
Some edge cases commonly arise where a kmer match starts before
the previous kmer match ends, but the two sequences still
mismatch.

To handle these, we make sure that either the word match in both sequences is spaced
`skip` after the previous word match, or both matches have spacing greater than
`wordlength`.
"""
function clean_matches(matches, wordlength, skip)
    # `matches` must already by sorted
    cleanedmatches = Vector{Int}[]
    push!(cleanedmatches, matches[1,:])
    currentvec = matches[1,:]
    for i in 2:length(matches[:, 2])
        # The condition below involves a kludge. It should really be
        # considering the case that skip == diff, rather than skip <
        # diff. Fix later.
        if !(skip < matches[i, 1] - currentvec[1] < wordlength ||
             (matches[i, 2] - currentvec[2] != skip &&
              matches[i, 2] - currentvec[2] < wordlength))
            push!(cleanedmatches, matches[i,:])
            currentvec = matches[i,:]
        end
    end
    return hcat(cleanedmatches...)'
end

"""
Returns ranges of overlapping word matches.
"""
function merge_overlapping(matches, wordlength, skip)
    range_inds = Vector{Int}[]
    push!(range_inds,[1, 1])
    current_range_ind_ind = 1
    for i in 2:length(matches[:, 2])
        # TODO: THESE BOUNDS (especially +wordlength) NEED TO BE CHECKED.
        # if ((matches[i, 1] < matches[i-1, 1] + (wordlength - 1) &&
        #      (matches[i, 2] < matches[i-1, 2] + (wordlength - 1))))
        if (matches[i, 1] - matches[i-1, 1] == skip) &&
            (matches[i, 2] - matches[i-1, 2] == skip)
            range_inds[current_range_ind_ind][2] = i
        else
            push!(range_inds,[i, i])
            current_range_ind_ind += 1
        end
    end
    return range_inds
end

"""
Returns actual word matches from indices of the matches.
"""
function get_matches(s1, s2, clean, range_inds, wordlength)
    s1matches = [clean[range_inds[i], 1] + [0, wordlength - 1]
                 for i in 1:length(range_inds)]
    s2matches = [clean[range_inds[i], 2] + [0, wordlength - 1]
                 for i in 1:length(range_inds)]
    match1 = [s1[i[1]:i[2]] for i in s1matches]
    match2 = [s2[i[1]:i[2]] for i in s2matches]
    return match1, match2
end

"""
Returns actual word mismatches (intervals between word matches) from indices of matches.
"""
function get_mismatches(s1, s2, matches, range_inds, wordlength)
    s1mismatches = [matches[(range_inds[i-1][2]):range_inds[i][1], 1] .+ [wordlength, -1]
                    for i in 2:length(range_inds)]
    s2mismatches = [matches[(range_inds[i-1][2]):range_inds[i][1], 2] .+ [wordlength,-1]
                    for i in 2:length(range_inds)]
    s1mismatches = vcat([[1, matches[1, 1]-1]],
                        s1mismatches,
                        [[matches[range_inds[length(range_inds)][2], 1]+wordlength, length(s1)]])
    s2mismatches = vcat([[1, matches[1, 2]-1]],
                        s2mismatches,
                        [[matches[range_inds[length(range_inds)][2], 2]+wordlength, length(s2)]])

    mismatch1 = [s1[i[1]:i[2]] for i in s1mismatches]
    mismatch2 = [s2[i[1]:i[2]] for i in s2mismatches]
    return mismatch1, mismatch2
end

"""
Find longest increasing subsequence of second column of given array.
Used to resolve bad orders of word matches while preserving as many matches as possible.
Because the first column (match locations in the first sequence) is sorted, corresponding
matches in the second column must be sorted, so we get the maximum number of such matches.
"""
function longest_incr_subseq(arr::Array{Int, 2})
    len = size(arr)[1]
    prevs = ones(Int, len)
    curr = ones(Int, len)
    maxlen = 0
    for i in 1:len
        # binary search for longest seq so far that's less than current value.
        lo = 1
        hi = maxlen
        while lo <= hi
            mid = Int(ceil((lo+hi)/2))
            if arr[curr[mid], 2] < arr[i, 2]
                lo = mid + 1
            else
                hi = mid - 1
            end
        end
        # update sequences seen and corresponding previous indices
        if lo > 1
            prevs[i] = curr[lo-1]
        else
            prevs[i] = 1
        end
        curr[lo] = i
        if lo > maxlen
            maxlen = lo
        end
    end
    # build subsequence from indices
    subseq = zeros(Int, maxlen, 2)
    k = curr[maxlen]
    for i in maxlen:-1:1
        subseq[i, :] = arr[k, :]
        k = prevs[k]
    end
    return subseq
end

#--------end of Kmer_seeded_align internals---------

"""
    kmer_seeded_align(s1::String, s2::String;
                      wordlength = 30,
                      skip = 10,
                      aligncodons = false,
                      banded = 1.0,
                      debug::Bool = false)

Returns aligned strings, where alignment is first done with larger word matches and
then (possibly banded) Needleman-Wunsch on intermediate intervals.
`skip` gives a necessary gap between searched-for words in `s1`.
For best results, use the default `wordlength` and `skip` values.
See `nw_align` for explanation of `banded`.
"""
function kmer_seeded_align(s1::String, s2::String;
                           wordlength = 30,
                           skip = 10,
                           aligncodons = false,
                           banded = 1.0,
                           debug::Bool = false)

    # ToDo (i think we tried this and it had insignificant speed up):
    # 1: Make this recurse. So instead fo calling nw_align on the set
    # of mismatches strings, it calls itself, but with a lower word
    # length. This will help for really noisy sequences.

    if (length(s1) < wordlength || length(s2) < wordlength || s1 == "" || s2 == "")
        return nw_align(s1, s2)
    end

    sorted = sorted_matches(s1, s2, wordlength, skip, aligncodons)
    if size(sorted)[1] == 0
        return nw_align(s1, s2, banded)
    end
    if matches_are_inconsistent(sorted)
        # try to extract an in-order subsequence to resolve inconsistency
        sorted = longest_incr_subseq(sorted)
        if matches_are_inconsistent(sorted)
            println("Notice: Word matching produced inconsistent ordering.",
                    " Consider a larger word size. Returning full DP alignment.",
                    "\nNote: this message shouldn't print with code changes. If it does,
                    check kner_seeded_align helper methods.")
            return nw_align(s1, s2, banded)
        end
    end

    clean = clean_matches(sorted, wordlength, skip)
    range_inds = merge_overlapping(clean, wordlength, skip)
    # get matches
    match1, match2 = get_matches(s1, s2, clean, range_inds, wordlength)
    # get mismatches
    mismatch1, mismatch2 = get_mismatches(s1, s2, clean, range_inds, wordlength)

    alignedStrings = collect(nw_align(mismatch1[1], mismatch2[1]))
    for i in 1:(length(match1))
        alignedStrings[1] = alignedStrings[1] * match1[i]
        alignedStrings[2] = alignedStrings[2] * match2[i]
        # al1, al2 = align(mismatch1[i+1], mismatch2[i+1],-1.0,-1.0, 20)
        al1, al2 = nw_align(mismatch1[i+1], mismatch2[i+1])
        alignedStrings[1] = alignedStrings[1] * al1
        alignedStrings[2] = alignedStrings[2] * al2
    end

    if debug
        if (degap(alignedStrings[1]) != degap(s1) || degap(alignedStrings[2]) != degap(s2))
            error("Aligned strings do not match original strings")
        end
    end
    return alignedStrings
end

"""
    triplet_kmer_seeded_align(s1::String, s2::String;
                              wordlength = 30,
                              skip = 9,
                              boundary_mult = 2,
                              alignedcodons = true,
                              debug::Bool=false)

Returns aligned strings, where alignment is first done with word matches and
then Needleman-Wunsch on intermediate intervals, prefering to preserve the
reading frame of the first arg `s1`.
`skip` gives a necessary gap between searched-for words in `s1`.
For best results, use the default `wordlength` and `skip` values.
See `triplet_nw_align` for explanation of `boundary_mult`.
"""
function triplet_kmer_seeded_align(s1::String, s2::String;
                           wordlength = 30,
                           skip = 9,
                           boundary_mult = 2,
                           alignedcodons = true,
                           debug::Bool=false)
    if !(wordlength % 3 == skip % 3 == 0)
        error("wordlength and skip must be divisible by 3")
    end

    if s1 == "" || s2 == ""
        return triplet_nw_align(s1, s2, boundary_mult=boundary_mult)
    end
    sorted = sorted_matches(s1, s2, wordlength, skip, alignedcodons)
    if size(sorted)[1] == 0
        return triplet_nw_align(s1, s2, boundary_mult=boundary_mult)
    end
    if matches_are_inconsistent(sorted)
        println("Notice: Word matching produced inconsistent ordering.",
                " Consider a larger word size. Returning full DP alignment.")
        return triplet_nw_align(s1, s2, boundary_mult=boundary_mult)
    end

    clean = clean_matches(sorted, wordlength, skip)
    range_inds = merge_overlapping(clean, wordlength, skip)
    # get matches
    match1, match2 = get_matches(s1, s2, clean, range_inds, wordlength)
    # get mismatches
    mismatch1, mismatch2 = get_mismatches(s1, s2, clean, range_inds, wordlength)

    alignedStrings = collect(triplet_nw_align(mismatch1[1], mismatch2[1], boundary_mult=boundary_mult))
    for i in 1:(length(match1))
        alignedStrings[1] = alignedStrings[1] * match1[i]
        alignedStrings[2] = alignedStrings[2] * match2[i]
        # al1, al2 = align(mismatch1[i+1], mismatch2[i+1],-1.0,-1.0, 20)
        al1, al2 = triplet_nw_align(mismatch1[i+1], mismatch2[i+1], boundary_mult=boundary_mult)
        alignedStrings[1] = alignedStrings[1] * al1
        alignedStrings[2] = alignedStrings[2] * al2
    end

    debug = true
    if debug && (degap(alignedStrings[1]) != degap(s1) || degap(alignedStrings[2]) != degap(s2))
        error("Aligned strings do not match original strings")
    end
    return alignedStrings
end

"""
    function local_kmer_seeded_align(s1::String, s2::String;
                                     wordlength = 30,
                                     skip = 10,
                                     trimpadding = 100,
                                     debug::Bool=false)

Returns locally aligned strings, where alignment is first done with word matches and
then Needleman-Wunsch on intermediate intervals.

`s1` is a reference to align to, and `s2` is a query to extract a local match from.
`s2` may be trimmed or expanded with gaps.
Before locally aligning ends of sequences, the ends of `s2` are trimmed to length
`trimpadding` for faster alignment. Increasing this will possible increase alignment accuracy
but effect runtime.
`skip` gives a necessary gap between searched-for words in `s1`.
For best results, use the default `wordlength` and `skip` values.
"""
function local_kmer_seeded_align(s1::String, s2::String;
                                   wordlength = 30,
                                   skip = 10,
                                   trimpadding = 100,
                                   debug::Bool=false)
    if s1 == "" || s2 == ""
        return local_align(s1, s2, rightaligned=false)
    end

    sorted = sorted_matches(s1, s2, wordlength, skip, false)
    if size(sorted)[1] == 0
        return local_align(s1, s2, rightaligned=false)
    end
    if matches_are_inconsistent(sorted)
        sorted = longest_incr_subseq(sorted)
        if matches_are_inconsistent(sorted)
            println("Notice: Word matching produced inconsistent ordering.",
                    " Consider a larger word size. Returning full DP alignment.")
            return local_align(s1, s2, rightaligned=false)
        end
    end

    clean = clean_matches(sorted, wordlength, skip)
    range_inds = merge_overlapping(clean, wordlength, skip)
    # get mismatches
    mismatch1, mismatch2 = get_mismatches(s1, s2, clean, range_inds, wordlength)
    # partition into front, middle, end mismatches, where middle is between beginning and end matching words
    frontmismatch1, mismatch1, endmismatch1 = mismatch1[1], mismatch1[2:end-1], mismatch1[end]
    frontmismatch2, mismatch2, endmismatch2 = mismatch2[1], mismatch2[2:end-1], mismatch2[end]
    # get matches
    match1, match2 = get_matches(s1, s2, clean, range_inds, wordlength)

    alignedStrings = ["", ""]
    for i in 1:(length(match1)-1)
        # append a match
        alignedStrings[1] = alignedStrings[1] * match1[i]
        alignedStrings[2] = alignedStrings[2] * match2[i]
        # append a mismatch
        al1, al2 = nw_align(mismatch1[i], mismatch2[i])
        alignedStrings[1] = alignedStrings[1] * al1
        alignedStrings[2] = alignedStrings[2] * al2
    end
    # append final match
    alignedStrings[1] = alignedStrings[1] * match1[end]
    alignedStrings[2] = alignedStrings[2] * match2[end]

    # trim query ends to trimpadding + length of ends of reference before local alignment to go fast
    if trimpadding > 0
        frontmismatch2 = frontmismatch2[max(1, end-length(frontmismatch1)-trimpadding):end]
        endmismatch2 = endmismatch2[1:end-length(endmismatch1)-trimpadding]
    end
    # do semiglobal alignment with beginning and ends and make them closer to the middle matches
    # reverse end mismatches to make local_align's default right alignment be left alignment
    fali1, fali2 = local_align(frontmismatch1, frontmismatch2, refend=true)
    eali1, eali2 = local_align(reverse(endmismatch1), reverse(endmismatch2), refend=true)
    eali1, eali2 = reverse(eali1), reverse(eali2)
    alignedStrings[1] = fali1*alignedStrings[1]*eali1
    alignedStrings[2] = fali2*alignedStrings[2]*eali2

    if debug
        if (degap(alignedStrings[1]) != degap(s1) || degap(alignedStrings[2]) != degap(s2))
            error("Aligned strings do not match original strings")
        end
    end
    return alignedStrings
end

loc_kmer_seeded_align = local_kmer_seeded_align

"""
    kmer_seeded_edit_dist(s1::String , s2::String;
                          wordlength = 30,
                          skip = 5,
                          aa_matches = false)

Computes levenshtein edit distance with speedups from only computing the dp scoring matrix between word matches.
If aa_matches = true, will attempt to find amino acid matches in any reference frame,
and add the nucleotide Hamming distance of these matches to Levenshtein distances of mismatches.
`skip` gives a necessary gap between searched-for words in `s1`.
For best results, use the default `wordlength` and `skip` values.
"""
function kmer_seeded_edit_dist(s1::String , s2::String;
                               wordlength = 30,
                               skip = 5,
                               aa_matches = false)
    if (length(s1) < wordlength || length(s2) < wordlength || s1 == "" || s2 == "")
        return levenshtein(s1, s2)
    end
    if aa_matches
        if wordlength % 3 != 0
            error("If aa_matches == true, wordlength must be divisible by 3")
        end
        sorted = sorted_aa_matches(s1, s2, wordlength)
        skip = 1
    else
        sorted = sorted_matches(s1, s2, wordlength, skip, false)
    end
    if size(sorted)[1] == 0
        return levenshtein(s1, s2)
    end
    if matches_are_inconsistent(sorted)
        # try to extract an in-order subsequence to resolve inconsistency
        sorted = longest_incr_subseq(sorted)
        if matches_are_inconsistent(sorted)
            println("Notice: Word matching produced inconsistent ordering.",
                    " Consider a larger word size. Returning full edit distance.")
            return levenshtein(s1, s2)
        end
    end
    trim = 0
    if aa_matches
        # trim word matches to avoid double counting some errors on ends of matches
        trim = max(div(wordlength, 5), 1)
        for i in 1:size(sorted)[1]
            sorted[i, 1] += trim
            sorted[i, 2] += trim
        end
    end
    clean = clean_matches(sorted, wordlength, skip)
    range_inds = merge_overlapping(clean, wordlength, skip)
    # 2*trim to accomodate trimmed word matches on each side of the word match
    mismatch1, mismatch2 = get_mismatches(s1, s2, clean, range_inds, wordlength - 2*trim)
    matched_diffs = 0
    if aa_matches
        match1, match2 = get_matches(s1, s2, clean, range_inds, wordlength - 2*trim)
        for i in 1:length(match1)
            for j in 1:length(match1[i])
                if match1[i][j] != match2[i][j]
                    matched_diffs += 1
                end
            end
        end
    end
    return matched_diffs + sum([levenshtein(mismatch1[i], mismatch2[i]) for i in 1:length(mismatch1)])
end

resolve_alignments(alignments::Array{String, 1}; mode = 1) = resolve_alignments(alignments[1], alignments[2], mode = mode)

"""
    resolve_alignments(ref::String, query::String; mode = 1)

Called on aligned strings. Resolves `query` with respect to `ref`.
`mode` = 1 for resolving single indels, `mode` = 2 for resolving single indels and codon insertions in query.
"""
function resolve_alignments(ref::String, query::String; mode = 1)
    if length(ref) != length(query)
        error("Arguments must be aligned")
    end
    for i in 1:length(ref)
        if ref[i] == query[i]
            continue
        end
        # 1%3 deletions, with gap of length 4 or something
#        if ((i%3 == 1 && i < length(query) && query[i+1] != '-') ||
#            (i%3 == 0 && query[i-1] != '-')) && single_mod_three_gap(query, i)
#            query = query[1:i-1] * "N" * query[i+1:end]
        if query[i] == '-' && !triple_gap(query, i)
            query = query[1:i-1] * "N" * query[i+1:end]
        # single insertions in query
        elseif single_gap(ref, i)
            query = query[1:i-1] * "-" * query[i+1:end]
        # single deletions in query
#        elseif single_gap(query, i)
#            query = query[1:i-1] * "N" * query[i+1:end]
        # copy codon of reference
        elseif mode == 2 && i <= length(ref)-3 && ref[i] == ref[i+1] == ref[i+2] == '-'
            query = query[1:i-1] * "---" * query[i+3:end]
        end
    end
    return degap(ref), degap(query)
end

"""
    align_reading_frames(clusters; k = 6, thresh = 0.03, verbose = false)

Takes `clusters` = [consensus_sequences, cluster_sizes], chooses references
out of consensuses that do not have stop codons in the middle, and makes all
consensus sequence reading frames agree. Returns resolved consensus seqs
(`goods`) along with filtered out consensus seqs that are >`thresh` divergent from
nearest reference (`bads`).
`k` = kmer size for computing kmer vectors of sequences.
"""
function align_reading_frames(clusters; k = 6, thresh = 0.03, verbose = false)
    consensuses = clusters[1]
    refs = String[]
    refkmers = []
    # find references
    for sq in consensuses
        tr = translate_to_aa(sq[1:end - length(sq)%3])
        if !('*' in tr[1:end-1])
            push!(refs, sq)
            push!(refkmers, kmer_count(sq, k))
        end
    end
    if length(refs) == 0
        return clusters, [[], []]
    end
    if verbose
        println("Using $(length(refs)) references to align reading frames")
    end
    goods = ([], [])
    bads = ([], [])
    # find closest ref and align
    for (i, sq) in enumerate(consensuses)
        sqkmer = kmer_count(sq, k)
        dists = zeros(length(refs))
        for (j, rfk) in enumerate(refkmers)
            dists[j] = corrected_kmer_dist(sqkmer, rfk, k=k)
        end
        minref_ind = indmin(dists)
        if dists[minref_ind] > thresh
            push!(bads[1], sq)
            push!(bads[2], clusters[2][i])
        else
            ali = triplet_kmer_seeded_align(refs[minref_ind], sq)
            #ali = triplet_kmer_seeded_align(refs[1], sq)
            _, newsq = resolve_alignments(ali...)
            # new reference frame may be really bad
            if (kmer_seeded_edit_dist(newsq, refs[minref_ind], aa_matches=true) > thresh*length(newsq))
                push!(bads[1], sq)
                push!(bads[2], clusters[2][i])
            else
                push!(goods[1], newsq)
                push!(goods[2], clusters[2][i])
            end
        end
    end
    if verbose
        println("Reading frame alignment: $(length(goods[1])) good sequences, $(length(bads[1])) bad sequences")
    end
    return goods, bads
end

align_reference_frames = align_reading_frames

"""
    local_edit_dist(s1::String, s2::String)

Returns the edit distance between two sequences after local alignment
"""
function local_edit_dist(s1::String, s2::String)
    str1, str2 = s1, s2
    if length(s1) > length(s2)
        str1, str2 = s2, s1
    end
    a, b = loc_kmer_seeded_align(str1, str2)
    dst = 0
    for i in 1:length(a)
        if a[i] != b[i] && a[i] != 'N' && b[i] != 'N'
            dst += 1
        end
    end
    return dst
end


text:
using NestedSamplers: default_update_interval

@testset "default helpers" begin
    @test default_update_interval(Proposals.Rejection(), 3) == 1.5
    @test default_update_interval(Proposals.RWalk(), 10) == 3.75
    @test default_update_interval(Proposals.RWalk(walks=10), 10) == 1.5
    @test default_update_interval(Proposals.RStagger(), 10) == 3.75
    @test default_update_interval(Proposals.RStagger(walks=10), 10) == 1.5
    @test default_update_interval(Proposals.Slice(), 30) == 135
    @test default_update_interval(Proposals.Slice(slices=10), 25) == 225
    @test default_update_interval(Proposals.RSlice(), 30) == 10
    @test default_update_interval(Proposals.RSlice(slices=10), 25) == 20
end

spl = Nested(3, 100)
@test spl.proposal isa Proposals.Rejection
@test spl.bounds == Bounds.MultiEllipsoid
@test spl.update_interval == 150
@test spl.enlarge == 1.25
@test spl.min_ncall == 200
@test spl.dlv ≈ log(101/100)

spl = Nested(10, 1000)
@test spl.proposal isa Proposals.RWalk
@test spl.update_interval == 3750

spl = Nested(30, 1500)
@test spl.proposal isa Proposals.Slice
@test spl.update_interval == 202500

ext: jl


text:
######### CmndStan program example  ###########

using CmdStan, Test, Statistics

ProjDir = dirname(@__FILE__)
cd(ProjDir) do

  bernoulli = "
  data {
    int<lower=0> N;
    int<lower=0,upper=1> y[N];
  }
  parameters {
    real<lower=0,upper=1> theta;
  }
  model {
    theta ~ beta(1,1);
      y ~ bernoulli(theta);
  }
  "

  bernoullidata = Dict("N" => 10, "y" => [0, 1, 0, 1, 0, 0, 0, 0, 0, 1])
  global stanmodel, rc, chns, cnames
  stanmodel = Stanmodel(CmdStan.Variational(), name="bernoulli",
    model=bernoulli)

  rc, chns, cnames = stan(stanmodel, bernoullidata, ProjDir, CmdStanDir=CMDSTAN_HOME)

  if rc == 0
    println()
    sdf = read_summary(stanmodel)
    @test 0.0 <= sdf[:theta, :mean][1] <= 2.0
  end
end # cd


text:
include("drop.jl")

include("hard_drop.jl")
include("soft_drop.jl")


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct MultiStartCost <: OperationalCost
        variable::VariableCost
        no_load::Float64
        fixed::Float64
        startup::NamedTuple{(:hot, :warm, :cold), NTuple{3, Float64}}
        shutdn::Float64
        forecasts::InfrastructureSystems.Forecasts
        internal::InfrastructureSystemsInternal
    end

Data Structure Operational Cost Data which includes fixed, variable cost, multiple start up cost and stop costs.

# Arguments
- `variable::VariableCost`: variable cost
- `no_load::Float64`: no load cost
- `fixed::Float64`: fixed cost
- `startup::NamedTuple{(:hot, :warm, :cold), NTuple{3, Float64}}`: startup cost
- `shutdn::Float64`: shutdown cost, validation range: (0, nothing), action if invalid: warn
- `forecasts::InfrastructureSystems.Forecasts`: internal forecast storage
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct MultiStartCost <: OperationalCost
    "variable cost"
    variable::VariableCost
    "no load cost"
    no_load::Float64
    "fixed cost"
    fixed::Float64
    "startup cost"
    startup::NamedTuple{(:hot, :warm, :cold), NTuple{3, Float64}}
    "shutdown cost"
    shutdn::Float64
    "internal forecast storage"
    forecasts::InfrastructureSystems.Forecasts
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function MultiStartCost(variable, no_load, fixed, startup, shutdn, forecasts=InfrastructureSystems.Forecasts(), )
    MultiStartCost(variable, no_load, fixed, startup, shutdn, forecasts, InfrastructureSystemsInternal(), )
end

function MultiStartCost(; variable, no_load, fixed, startup, shutdn, forecasts=InfrastructureSystems.Forecasts(), )
    MultiStartCost(variable, no_load, fixed, startup, shutdn, forecasts, )
end

# Constructor for demo purposes; non-functional.
function MultiStartCost(::Nothing)
    MultiStartCost(;
        variable=VariableCost((0.0, 0.0)),
        no_load=0.0,
        fixed=0.0,
        startup=(hot = START_COST, warm = START_COST,cold = START_COST),
        shutdn=0.0,
        forecasts=InfrastructureSystems.Forecasts(),
    )
end

"""Get MultiStartCost variable."""
get_variable(value::MultiStartCost) = value.variable
"""Get MultiStartCost no_load."""
get_no_load(value::MultiStartCost) = value.no_load
"""Get MultiStartCost fixed."""
get_fixed(value::MultiStartCost) = value.fixed
"""Get MultiStartCost startup."""
get_startup(value::MultiStartCost) = value.startup
"""Get MultiStartCost shutdn."""
get_shutdn(value::MultiStartCost) = value.shutdn

InfrastructureSystems.get_forecasts(value::MultiStartCost) = value.forecasts
"""Get MultiStartCost internal."""
get_internal(value::MultiStartCost) = value.internal

"""Set MultiStartCost variable."""
set_variable!(value::MultiStartCost, val::VariableCost) = value.variable = val
"""Set MultiStartCost no_load."""
set_no_load!(value::MultiStartCost, val::Float64) = value.no_load = val
"""Set MultiStartCost fixed."""
set_fixed!(value::MultiStartCost, val::Float64) = value.fixed = val
"""Set MultiStartCost startup."""
set_startup!(value::MultiStartCost, val::NamedTuple{(:hot, :warm, :cold), NTuple{3, Float64}}) = value.startup = val
"""Set MultiStartCost shutdn."""
set_shutdn!(value::MultiStartCost, val::Float64) = value.shutdn = val

InfrastructureSystems.set_forecasts!(value::MultiStartCost, val::InfrastructureSystems.Forecasts) = value.forecasts = val
"""Set MultiStartCost internal."""
set_internal!(value::MultiStartCost, val::InfrastructureSystemsInternal) = value.internal = val


text:
# Autogenerated wrapper script for OpenJpeg_jll for powerpc64le-linux-gnu
export libopenjp2, opj_compress, opj_decompress, opj_dump

using LittleCMS_jll
using libpng_jll
using Libtiff_jll
JLLWrappers.@generate_wrapper_header("OpenJpeg")
JLLWrappers.@declare_library_product(libopenjp2, "libopenjp2.so.7")
JLLWrappers.@declare_executable_product(opj_compress)
JLLWrappers.@declare_executable_product(opj_decompress)
JLLWrappers.@declare_executable_product(opj_dump)
function __init__()
    JLLWrappers.@generate_init_header(LittleCMS_jll, libpng_jll, Libtiff_jll)
    JLLWrappers.@init_library_product(
        libopenjp2,
        "lib/libopenjp2.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        opj_compress,
        "bin/opj_compress",
    )

    JLLWrappers.@init_executable_product(
        opj_decompress,
        "bin/opj_decompress",
    )

    JLLWrappers.@init_executable_product(
        opj_dump,
        "bin/opj_dump",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# This file contains functions to create new positions.

"""
function make_new_positions

Return a new vector of SPADL.
"""
function make_new_positions(vector_wyscout_data::Vector{WyscoutData})::Vector{WyscoutData}

    for i in eachindex(vector_wyscout_data)
        if length(vector_wyscout_data[i].event.positions) == 2
            vector_wyscout_data[i].event_fixed.start_x = clamp(vector_wyscout_data[i].event.positions[1]["x"], 0, 105)
            vector_wyscout_data[i].event_fixed.start_y = clamp(vector_wyscout_data[i].event.positions[1]["y"], 0, 68)
            vector_wyscout_data[i].event_fixed.end_x = clamp(vector_wyscout_data[i].event.positions[2]["x"], 0, 105)
            vector_wyscout_data[i].event_fixed.end_y = clamp(vector_wyscout_data[i].event.positions[2]["y"], 0, 68)
        elseif length(event_df[i].positions) == 1
            vector_wyscout_data[i].event_fixed.start_x = clamp(vector_wyscout_data[i].event.positions[1]["x"], 0, 105)
            vector_wyscout_data[i].event_fixed.start_y = clamp(vector_wyscout_data[i].event.positions[1]["y"], 0, 68)
            vector_wyscout_data[i].event_fixed.end_x = clamp(vector_wyscout_data[i].event.start_x, 0, 105)
            vector_wyscout_data[i].event_fixed.end_y = clamp(vector_wyscout_data[i].event.start_y, 0, 68)
        end
    end

    return vector_wyscout_data
end


text:
using MAT, DelimitedFiles, SparseArrays

ncol(io) = mapreduce(!isempty, +, eachcol(read(io, "xdata"))) # only works well if the empty column is last

function savepixels(pixelfolder, resfile)
    matopen(resfile) do io
        n = ncol(io)
        ids = Vector{UUID}(undef, n)
        for i in 1:n
            x, y, t = getcoordinates(io, i)
            @assert !isempty(x) "no coordinates in column $i of $n columns in $resfile"
            @assert length(x) == length(y) == length(t) "x, y, and t must have the same length in $resfile"
            id = uuid1()
            writedlm(joinpath(pixelfolder, "$id.csv"), zip(x, y, t))
            ids[i] = id
        end
        ids
    end
end

function getcoordinates(io, i)
    _x = read(io, "xdata")[:, i]
    x = nonzeros(_x)
    y = nonzeros(read(io, "ydata")[:,i])
    fr = read(io, "status")["FrameRate"]
    t = length(x) == 1 ? findfirst(!iszero, _x)/fr : findall(!iszero, _x)./fr
    return x, y, t
end


text:
using LifeTable, DataFrames

function zfrac(a, b)
	if b==0
		return 0
	else
		return a / b
	end
end

function amplt(inframe, sex, rate = "inc")
	age = inframe[1]
	pop = inframe[2]
	acd = inframe[3]
	cd = inframe[4]
	cc = inframe[5]
	if rate == "inc"
		ncol = cc
		dcol = acd .- cd .+ cc
	elseif rate == "mort"
		ncol = cd
		dcol = acd
	end
	df = DataFrame(age = age, pop = pop, dcol = dcol)
	cprop = zfrac.(ncol, dcol)
	lt = periodlifetable(df, sex)
	Dict(:sex=>sex, :cl_frame=>causelife(lt, cprop))
end

function ratechange(inframe, changefac_cd, changefac_cc, changefac_nocd)
	ncd = inframe[4] .* changefac_cd
	ncc = inframe[5] .* changefac_cc
	nacd = (inframe[3].-inframe[4]).*changefac_nocd .+ ncd
	DataFrame(age = inframe[1], pop = inframe[2],
		acd = nacd, cd = ncd, cc = ncc)
end

ptoage(cl_frame, plt_frame) = [cl_frame[:f][1] .* (1-plt_frame[:l]); cl_frame[:f][1]]

function changes(ifdict1, ifdict2)
	if1 = ifdict1[:cl_frame]
	if2 = ifdict2[:cl_frame]
	lts = map(i->periodlifetable(i, ifdict1[:sex], true, "rate"), [if1, if2])
	ptoages1 = ptoage(if1, lts[1])[2:end]
	ptoages2 = ptoage(if2, lts[2])[2:end]
	DataFrame(age = [if1[:age][2:end]; Inf],
		ptoages1 = ptoages1, ptoages2 = ptoages2,
		relchange = ptoages1./ptoages2, abschange = ptoages1.-ptoages2)
end

max_stars_repo_name: klpn/static-dust

max_issues_repo_name: klpn/static-dust

max_forks_repo_name: klpn/static-dust


text:
@testset "Error Handling{$ST}" begin
  msg = GridapDistributedPETScWrappers.PetscErrorMessage(73)
  @test msg == "Object is in wrong state"
  @test_throws GridapDistributedPETScWrappers.PetscError GridapDistributedPETScWrappers.chk(76)
end


text:
# All code in here are directly derived from the C code

using MetaBalls:MetaBalls

ISO_VALUE = 0.1

function edge_lookup(voxel)
    cubeIndex = 0

    # println(ISO_MIN, " ", voxel[1].value, " ", ISO_MAX)

    # if (ISO_MIN < voxel[1].value) && (voxel[1].value < ISO_MAX) cubeIndex |= 1 end
    # if (ISO_MIN < voxel[2].value) && (voxel[2].value < ISO_MAX) cubeIndex |= 2 end
    # if (ISO_MIN < voxel[3].value) && (voxel[3].value < ISO_MAX) cubeIndex |= 4 end
    # if (ISO_MIN < voxel[4].value) && (voxel[4].value < ISO_MAX) cubeIndex |= 8 end
    # if (ISO_MIN < voxel[5].value) && (voxel[5].value < ISO_MAX) cubeIndex |= 16 end
    # if (ISO_MIN < voxel[6].value) && (voxel[6].value < ISO_MAX) cubeIndex |= 32 end
    # if (ISO_MIN < voxel[7].value) && (voxel[7].value < ISO_MAX) cubeIndex |= 64 end
    # if (ISO_MIN < voxel[8].value) && (voxel[8].value < ISO_MAX) cubeIndex |= 128 end

    if (voxel[1].iso_value < ISO_VALUE) cubeIndex |= 1 end
    if (voxel[2].iso_value < ISO_VALUE) cubeIndex |= 2 end
    if (voxel[3].iso_value < ISO_VALUE) cubeIndex |= 4 end
    if (voxel[4].iso_value < ISO_VALUE) cubeIndex |= 8 end
    if (voxel[5].iso_value < ISO_VALUE) cubeIndex |= 16 end
    if (voxel[6].iso_value < ISO_VALUE) cubeIndex |= 32 end
    if (voxel[7].iso_value < ISO_VALUE) cubeIndex |= 64 end
    if (voxel[8].iso_value < ISO_VALUE) cubeIndex |= 128 end

    return cubeIndex
end


text:
using RemoteREPL
using Test
using Sockets

ENV["JULIA_DEBUG"] = "RemoteREPL"

@testset "Protocol header handshake" begin
    # Happy case
    io = IOBuffer()
    RemoteREPL.send_header(io)
    seek(io, 0)
    @test RemoteREPL.verify_header(io)

    # Broken magic number
    io = IOBuffer()
    write(io, "BrokenMagic", RemoteREPL.protocol_version)
    seek(io, 0)
    @test_throws ErrorException RemoteREPL.verify_header(io)
    # Version mismatch
    io = IOBuffer()
    write(io, RemoteREPL.protocol_version, typemax(UInt32))
    seek(io, 0)
    @test_throws ErrorException RemoteREPL.verify_header(io)

    # ser_version=10 in julia 1.5
    for (local_ver, remote_ver) in [(10,13), (13,10)]
        io = IOBuffer()
        seek(io, 0)
        RemoteREPL.send_header(io, remote_ver)
        @test_throws ErrorException RemoteREPL.verify_header(io, local_ver)
    end
end

# Connect to a non-default loopback address to test SSH integration
test_interface = ip"127.111.111.111"

# Detect whether an ssh server is running on the default port 22
use_ssh = if "use_ssh=true" in ARGS
    true
elseif "use_ssh=false" in ARGS
    false
else
    # Autodetct
    try
        socket = Sockets.connect(test_interface, 22)
        # https://tools.ietf.org/html/rfc4253#section-4.2
        id_string = String(readavailable(socket))
        startswith(id_string, "SSH-")
    catch
        false
    end
end

if !use_ssh
    test_interface = Sockets.localhost
end
@info use_ssh ? "Running tests with SSH tunnel" : "Testing without SSH tunnel - localhost only"

# Use non-default port to avoid clashes with concurrent interactive use or testing.
test_port = RemoteREPL.find_free_port(Sockets.localhost)
server_proc = run(`$(Base.julia_cmd()) -e "using Sockets; using RemoteREPL; serve_repl($test_port)"`, wait=false)

try

@testset "RemoteREPL.jl" begin
    local socket = nothing
    max_tries = 4
    for i=1:max_tries
        try
            socket = RemoteREPL.setup_connection(test_interface, test_port,
                                                 tunnel=use_ssh ? :ssh : :none,
                                                 ssh_opts=`-o StrictHostKeyChecking=no`)
            break
        catch exc
            if i == max_tries
                rethrown()
            end
            # Server not yet started - continue waiting
            sleep(2)
        end
    end
    @assert isopen(socket)

    # Some basic tests of the transport and server side and partial client side.
    #
    # More full testing of the client code would requires some tricky mocking
    # of the REPL environment.
    runcommand(cmdstr) = sprint(io->RemoteREPL.run_remote_repl_command(socket, io, cmdstr))

    @test runcommand("asdf = 42") == "42\n"
    @test runcommand("Main.asdf") == "42\n"
    @test !isdefined(Main, :asdf) # asdf not defined locally

    # Output Limiting
    @test '⋮' in runcommand("ones(1000)")

    # Error formatting
    @test occursin(r"DivideError.*Stacktrace"s, runcommand("1÷0"))

    # Semicolon suppresses output
    @test runcommand("asdf;") == ""

    # Execute a single command on a separate connection
    @test RemoteREPL.remote_eval(test_interface, test_port, "asdf") == "42\n"
end

finally
    kill(server_proc)
end


text:
using Test
using KeboolaConnectionTransformation
using JSON
import KeboolaConnectionComponent

function createconfig(data)
    path = mktempdir();
    ENV["KBC_DATADIR"] = path
    open(path * "/config.json", "w") do file
        JSON.print(file, data, 4)
    end
    path
end

function getconfigdata()
    Dict(
        "storage" => Dict(
            "input" => Dict(
                "tables" => [
                    Dict(
                        "source" => "in.c-main.data",
                        "destination" => "sample.csv"
                    )
                ]
            ),
            "output" => Dict(
                "tables" => [
                    Dict(
                        "source" => "sample.csv",
                        "destination" => "out.c-main.data"
                    )
                ]
            )
        ),
        "parameters" => Dict(
        )
    )
end

@testset "basic run" begin
    data = getconfigdata()
    data["parameters"] = Dict(
        "packages" => ["CSV", "LightXML"],
        "script" => ["println(\"Hello\")", "println(\"World\")\nusing LightXML"]
    )
    path = createconfig(data)
    ENV["KBC_DATADIR"] = path
    @test KeboolaConnectionTransformation.run() == 0

    open(joinpath(path, "script.jl"), "r") do file
        script = read(file, String)
        @test script == "println(\"Hello\")\nprintln(\"World\")\nusing LightXML"
    end
    @test pwd() == path

    delete!(ENV, "KBC_DATADIR")
end

@testset "no packages" begin
    data = getconfigdata()
    data["parameters"] = Dict(
        "packages" => [],
        "script" => ["println(\"Hello\")\nprintln(\"World\")"]
    )
    path = createconfig(data)
    ENV["KBC_DATADIR"] = path
    @test KeboolaConnectionTransformation.run() == 0

    open(joinpath(path, "script.jl"), "r") do file
        script = read(file, String)
        @test script == "println(\"Hello\")\nprintln(\"World\")"
    end
    @test pwd() == path

    delete!(ENV, "KBC_DATADIR")
end

@testset "invalid package" begin
    data = getconfigdata()
    data["parameters"] = Dict(
        "packages" => ["CSV", "non existent package"],
        "script" => ["println(\"Hello\")", "println(\"World\")\nusing LightXML"]
    )
    path = createconfig(data)
    ENV["KBC_DATADIR"] = path
    try
        KeboolaConnectionTransformation.run()
    catch e
        @test e isa ArgumentError
        @test e.msg == "Failed to install package 'non existent package' error: `non existent package` is not a valid package name"
    end
    delete!(ENV, "KBC_DATADIR")
end

@testset "invalid script" begin
    data = getconfigdata()
    data["parameters"] = Dict(
        "packages" => [],
        "script" => ["println(\"Hello\")", "this is not a valid code"]
    )
    path = createconfig(data)
    ENV["KBC_DATADIR"] = path
    try
        KeboolaConnectionTransformation.run()
    catch e
        @test e isa ArgumentError
        @test occursin("The transformation script failed: LoadError: syntax: extra token \"is\" after end of expression\nin expression starting at", e.msg)
        @test occursin("script.jl:2", e.msg)
    end
    delete!(ENV, "KBC_DATADIR")
end

max_stars_repo_name: keboola/julia-transformation

max_issues_repo_name: keboola/julia-transformation

max_forks_repo_name: keboola/julia-transformation


text:
# Part 2: Thermomechanical coupling
const GPU_USE = true
using ParallelStencil
using ParallelStencil.FiniteDifferences2D
@static if GPU_USE
    @init_parallel_stencil(CUDA, Float64, 2)
else
    @init_parallel_stencil(Threads, Float64, 2)
end
using Plots, Printf, ImplicitGlobalGrid
import MPI

# macros to avoid array allocation
macro  ∇V(ix,iy)  esc(:(               (Vx[$ix,$iy]-Vx[$ix-1,$iy])/dx + (Vy[$ix,$iy]-Vy[$ix,$iy-1])/dy  )) end
macro ϵxx(ix,iy)  esc(:(         0.5*( (Vx[$ix,$iy]-Vx[$ix-1,$iy])/dx - (Vy[$ix,$iy]-Vy[$ix,$iy-1])/dy) )) end
macro τxx(ix,iy)  esc(:( ηc[$ix,$iy]*( (Vx[$ix,$iy]-Vx[$ix-1,$iy])/dx - (Vy[$ix,$iy]-Vy[$ix,$iy-1])/dy + 2*λb*( (Vx[$ix,$iy]-Vx[$ix-1,$iy])/dx + (Vy[$ix,$iy]-Vy[$ix,$iy-1])/dy )) - P[$ix,$iy] )) end
macro ϵyy(ix,iy)  esc(:(         0.5*(-(Vx[$ix,$iy]-Vx[$ix-1,$iy])/dx + (Vy[$ix,$iy]-Vy[$ix,$iy-1])/dy) )) end
macro τyy(ix,iy)  esc(:( ηc[$ix,$iy]*(-(Vx[$ix,$iy]-Vx[$ix-1,$iy])/dx + (Vy[$ix,$iy]-Vy[$ix,$iy-1])/dy + 2*λb*( (Vx[$ix,$iy]-Vx[$ix-1,$iy])/dx + (Vy[$ix,$iy]-Vy[$ix,$iy-1])/dy )) - P[$ix,$iy] )) end
macro ϵxy(ix,iy)  esc(:(         0.5*(( Vx[$ix,$iy+1]-Vx[$ix,$iy])/dy + (Vy[$ix+1,$iy]-Vy[$ix,$iy])/dx) )) end
macro τxy(ix,iy)  esc(:( ηv[$ix,$iy]*( (Vx[$ix,$iy+1]-Vx[$ix,$iy])/dy + (Vy[$ix+1,$iy]-Vy[$ix,$iy])/dx) )) end
macro  qx(ix,iy)  esc(:( -(T[$ix+1,$iy+1] - T[$ix,$iy+1])/dx )) end
macro  qy(ix,iy)  esc(:( -(T[$ix+1,$iy+1] - T[$ix+1,$iy])/dy )) end
macro   S(ix,iy)  esc(:( 4* ηc[$ix,$iy]* ϵii2[$ix,$iy]^2 )) end
macro  η0(ix,iy)  esc(:( ϵii2[$ix,$iy]^(1/n-1)*exp( -T[$ix,$iy]*(1 /(1+T[$ix,$iy]/T0)) ) )) end

macro  dτP(ix,iy)  esc(:( θp*4.1/max(nx-2,ny-2)*ηc[$ix,$iy]*(1.0+λb) )) end
macro dτVx(ix,iy)  esc(:( θv/4.1*min(dx,dy)^2 /((ηc[$ix,$iy]+ηc[$ix+1,$iy])/2)/(1.0+λb) )) end
macro dτVy(ix,iy)  esc(:( θv/4.1*min(dx,dy)^2 /((ηc[$ix,$iy]+ηc[$ix,$iy+1])/2)/(1.0+λb) )) end
macro  dτT(ix,iy)  esc(:( θT/4.1*min(dx,dy)^2 )) end
"""
  compute the second order invariant of the deviatoric strain rate tensor
"""
@parallel_indices (ix,iy) function compute_2invar!(ϵii2, Vx, Vy, dx, dy, nx, ny)
  if ( ix<=(nx-2) && iy<=(ny-2) )
    ϵii2[ix+1,iy+1] = ( 0.5*(@ϵxx(ix+1,iy+1)^2 + @ϵyy(ix+1,iy+1)^2) + (0.25*(@ϵxy(ix,iy) + @ϵxy(ix+1,iy) + @ϵxy(ix,iy+1) + @ϵxy(ix+1,iy+1)))^2 )^0.5
  end
  return
end
"""
  compute nonlinear viscosity defined at cell centers
"""
@parallel_indices (ix,iy) function compute_viscosity_c!(ηc, ϵii2, n, T, T0, θη, nx, ny)
  if ( ix<=(nx-2) && iy<=(ny-2) )
    ηc[ix+1,iy+1] = exp( θη*log(ηc[ix+1,iy+1]) + (1-θη)*log(@η0(ix+1,iy+1)) )
  end
  return
end
"""
  compute nonlinear viscosity at cell vertices by averaging the values of 4 neighboring cell centers
"""
@parallel_indices (ix,iy) function compute_viscosity_v!(ηv, ηc, size_ηv1, size_ηv2)
  if ( ix<=size_ηv1 && iy<=size_ηv2 )
    ηv[ix,iy] = 0.25*( ηc[ix,iy] + ηc[ix+1,iy] + ηc[ix,iy+1] + ηc[ix+1,iy+1] )
  end
  return
end
"""
  compute residuals of continuity and temperature equations
"""
@parallel_indices (ix,iy) function compute_residuals_P_T!(dPdτ, dTdτ, P, Vx, Vy, T, Tt, ηc, ϵii2, n, T0, dt, dx, dy, nx, ny)
  if ( ix<=(nx-2) && iy<=(ny-2) )
    dPdτ[ix+1,iy+1] = - @∇V(ix+1,iy+1)
    dTdτ[ix+1,iy+1] = (Tt[ix+1,iy+1]-T[ix+1,iy+1])/dt - (@qx(ix+1,iy)-@qx(ix,iy))/dx - (@qy(ix,iy+1)-@qy(ix,iy))/dy + @S(ix+1,iy+1)
  end
  return
end
"""
  compute residuals of x y momentum equations, including damped residuals
"""
@parallel_indices (ix,iy) function compute_residuals_Vx_Vy!(dVxdτ, dVydτ, dVxdτ0, dVydτ0, P, Vx, Vy, ηv, ηc, λb, dampx, dampy, dx, dy, size_Vx1, size_Vx2, size_Vy1, size_Vy2)
  if ( ix<=(size_Vx1-2) && iy<=(size_Vx2-2) )
    dVxdτ0[ix+1,iy+1] = dVxdτ[ix+1,iy+1] + dampx * dVxdτ0[ix+1,iy+1]
     dVxdτ[ix+1,iy+1] = (@τxy(ix+1,iy+1) - @τxy(ix+1,iy))/dy + (@τxx(ix+2,iy+1) - @τxx(ix+1,iy+1))/dx end
  if ( ix<=(size_Vy1-2) && iy<=(size_Vy2-2) )
    dVydτ0[ix+1,iy+1] = dVydτ[ix+1,iy+1] + dampy * dVydτ0[ix+1,iy+1]
     dVydτ[ix+1,iy+1] = (@τxy(ix+1,iy+1) - @τxy(ix,iy+1))/dx + (@τyy(ix+1,iy+2) - @τyy(ix+1,iy+1))/dy end
  return
end
"""
  update pressure and temperature in pseudo time
"""
@parallel_indices (ix,iy) function stepping_pseudo_P_T!(P, T, dPdτ, dTdτ, ηc, θp, θT, λb, dx, dy, nx, ny)
  if ( ix<=(nx-2) && iy<=(ny-2) )
    P[ix+1,iy+1] = P[ix+1,iy+1] + @dτP(ix+1,iy+1)*dPdτ[ix+1,iy+1]
    T[ix+1,iy+1] = T[ix+1,iy+1] + @dτT(ix+1,iy+1)*dTdτ[ix+1,iy+1]
  end
  return
end
"""
  update velocities in pseudo time
"""
@parallel_indices (ix,iy) function stepping_pseudo_Vx_Vy!(Vx, Vy, dVxdτ, dVydτ, dVxdτ0, dVydτ0, ηc, λb, θv, dampx, dampy, dx, dy, size_Vx1, size_Vx2, size_Vy1, size_Vy2)
  if ( ix<=(size_Vx1-2) && iy<=(size_Vx2-2) ) Vx[ix+1,iy+1] = Vx[ix+1,iy+1] + @dτVx(ix+1,iy+1)*( dVxdτ[ix+1,iy+1]+ dampx*dVxdτ0[ix+1,iy+1] )  end
  if ( ix<=(size_Vy1-2) && iy<=(size_Vy2-2) ) Vy[ix+1,iy+1] = Vy[ix+1,iy+1] + @dτVy(ix+1,iy+1)*( dVydτ[ix+1,iy+1]+ dampy*dVydτ0[ix+1,iy+1] ) end
  return
end
"""
  compute the norm of a matrix
"""
@views norm_g(A) = (sum2_l = sum(A.^2); sqrt(MPI.Allreduce(sum2_l, MPI.SUM, MPI.COMM_WORLD)))

@views function TM_2D(; do_visu = false, do_print = false)
  # Physics
  Lx, Ly = 0.86038, 0.86038     # global domain size
  n      = 3                    # power law rheology exponent
  Vbc    = 66.4437              # boundary velocity
  T0     = 49.3269 / n          # initial temperature
  r      = 0.0737               # perturbation radius
  Tamp   = 0.1*T0               # perturbation temperature

  # Numerics
  nx, ny = 32*2^2, 32*2^2       # local grid size
  tol    = 1e-5                 # non-linear tolerance
  ξ      = 10.0                 # physical steps reduction for temperature
  λb     = 1.0                  # numerical bulk viscosity
  θη     = 0.9                  # relaxation factor for viscosity
  ν      = 4.0                  # damping factor for velocity residuals
  nt     = 5                    # number of time steps
  iter_max = 1e6                # max PT iterations
  n_check  = 1e3                # residuals checking step
  n_vis    = 1                  # visualization step
  θp,θv,θT = 0.5, 0.5, 0.5      # PT steps reduction for pressure, velocities and temperature

  # Derived numerics
  me, dims = init_global_grid(nx, ny, 1)  # Initialize 2D MPI
  @static if GPU_USE select_device() end  # select one GPU per MPI local rank (if >1 GPU per node)
  dx, dy = Lx/(nx_g()-2), Ly/(ny_g()-2)   # local grid step
  dt     = ξ*min(dx, dy)^2/4.1            # physical time step
  dampx  = 1 - ν/(nx-2)                   # damping for x-momentum residuals
  dampy  = 1 - ν/(ny-2)                   # damping for y-momentum residuals

  # Array initialisation
  P      = @zeros(nx  , ny  )    # pressure
  T      = @zeros(nx  , ny  )    # temperature
  Tt     = @zeros(nx  , ny  )    # temperature at current physical time step
  ηc     =  @ones(nx  , ny  )    # effective viscosity at cell centroids
  ηv     =  @ones(nx-1, ny-1)    # effective viscosity at cell vertices
  Vx     = @zeros(nx-1, ny  )    # x-velocity
  Vy     = @zeros(nx  , ny-1)    # y-velocity
  ϵii2   = @zeros(nx  , ny  )    # 2nd order invariant
  dPdτ   = copy(P)               # continuity residuals
  dTdτ   = copy(T)               # temperature residuals
  dVxdτ  = copy(Vx)              # x-momentum residuals
  dVydτ  = copy(Vy)              # y-momentum residuals
  dVxdτ0 = copy(dVxdτ)           # damped x-momentum residuals
  dVydτ0 = copy(dVydτ)           # damped y-momentum residuals
  size_Vx1, size_Vx2 = size(Vx,1), size(Vx,2)
  size_Vy1, size_Vy2 = size(Vy,1), size(Vy,2)
  size_ηv1, size_ηv2 = size(ηv,1), size(ηv,2)

  # Initial conditions
  Vx[:,2:end-1] .= Data.Array([ Vbc*(x_g(ix,dx,Vx)-dx/2)/Lx for ix=1:size(Vx,1), iy=2:size(Vx,2)-1])
  Vy[2:end-1,:] .= Data.Array([-Vbc*(y_g(iy,dy,Vy)-dy/2)/Ly for ix=2:size(Vy,1)-1, iy=1:size(Vy,2)])
  Vx[:,[1 end]] .= Vx[:,[2 end-1]]
  Vy[[1 end],:] .= Vy[[2 end-1],:]
  T[2:end-1,2:end-1] .= Data.Array([(x_g(ix,dx,T)-dx/2)^2 + (y_g(iy,dy,T)-dy/2)^2 < r^2 for ix=2:nx-1, iy=2:ny-1] .* Tamp)
  T[:, [1 end]] .= T[:,[2 end-1]]
  T[[1 end], :] .= T[[2 end-1],:]

  # Preparation of visualisation
  if do_visu
    if (me==0) ENV["GKSwstype"]="nul"; if isdir("TM2D_out")==false mkdir("TM2D_out") end; loadpath = "./TM2D_out/"; anim = Animation(loadpath,String[]); println("Animation directory: $(anim.dir)") end
    nx_v, ny_v = (nx-2)*dims[1], (ny-2)*dims[2]
    if (nx_v*ny_v*sizeof(Data.Number) > 0.8*Sys.free_memory()) error("Not enough memory for visualization.") end
    T_v   = zeros(nx_v, ny_v) # global array for visu
    T_inn = zeros(nx-2, ny-2) # no halo local array for visu
    xc, yc = LinRange(dx/2, Lx-dx/2, nx_v), LinRange(dy/2, Ly-dy/2, ny_v) # inner cell centers
  end

  # Physical time stepping
  for it = 1:nt
    Tt .= T     # temperature updates from previous timestep
    t_tic = 0.0; niter = 0
    for iter = 1:iter_max # pseudo transient stepping
      # start timer after warming up
      if (iter==11) t_tic = Base.time(); niter = 0 end

      # compute 2nd order invariant of the deviatoric strain rate tensor
      @parallel compute_2invar!(ϵii2, Vx, Vy, dx, dy, nx, ny)

      # compute viscosity at cell centers and vertices
      @parallel compute_viscosity_c!(ηc, ϵii2, n, T, T0, θη, nx, ny)
      update_halo!(ηc)
      @parallel compute_viscosity_v!(ηv, ηc, size_ηv1, size_ηv2)

      # compute residuals for inner cells
      @parallel compute_residuals_Vx_Vy!(dVxdτ, dVydτ, dVxdτ0, dVydτ0, P, Vx, Vy, ηv, ηc, λb, dampx, dampy, dx, dy, size_Vx1, size_Vx2, size_Vy1, size_Vy2)
      @parallel compute_residuals_P_T!(dPdτ, dTdτ, P, Vx, Vy, T, Tt, ηc, ϵii2, n, T0, dt, dx, dy, nx, ny)

      # update velocities, pressure and temperature
      @parallel stepping_pseudo_Vx_Vy!(Vx, Vy, dVxdτ, dVydτ, dVxdτ0, dVydτ0, ηc, λb, θv, dampx, dampy, dx, dy, size_Vx1, size_Vx2, size_Vy1, size_Vy2) # inner cells
      @parallel stepping_pseudo_P_T!(P, T, dPdτ, dTdτ, ηc, θp, θT, λb, dx, dy, nx, ny)
       T[:, [1 end]] .=  T[:,[2 end-1]] # zero heat flux at top and bottom boundaries
       T[[1 end], :] .=  T[[2 end-1],:] # zero heat flux at left and right boundaries
      Vx[:, [1 end]] .= Vx[:,[2 end-1]] # zero shear stress at top and bottom boundaries
      Vy[[1 end], :] .= Vy[[2 end-1],:] # zero shear stress at left and right boundaries
      update_halo!(T, Vx, Vy)
      niter += 1

      # check residuals
      if iter % n_check == 0
        R_P  = norm_g( Array( dPdτ[2:end,2:end]) )/( (nx_g()-2)*(ny_g()-2) )
        R_Vx = norm_g( Array(dVxdτ[2:end,2:end]) )/( (nx_g()-2)*(ny_g()-2) )
        R_Vy = norm_g( Array(dVydτ[2:end,2:end]) )/( (nx_g()-2)*(ny_g()-2) )
        R_T  = norm_g( Array( dTdτ[2:end,2:end]) )/( (nx_g()-2)*(ny_g()-2) )
        if (me==0 && do_print) println(">>time = $(round(it*dt, sigdigits=3)), iter = $(iter), R_P = $(round(R_P,sigdigits=3)), R_Vx = $(round(R_Vx,sigdigits=3)), R_Vy = $(round(R_Vy,sigdigits=3)), R_T = $(round(R_T,sigdigits=3))") end
        if max(R_P, R_Vx, R_Vy, R_T) < tol break; end
      end
    end # end pseudo transient stepping

    # performance evaluation
    t_toc = Base.time() - t_tic
    A_eff = 10/1e9*nx_g()*ny_g()*sizeof(Float64)  # Effective main memory access per iteration [GB]
    t_it  = t_toc/niter                           # Execution time per iteration [s]
    T_eff = A_eff/t_it                            # Effective memory throughput [GB/s]
    if (me==0) println("-----time = $(round(it*dt, sigdigits=3)) s, T_eff = $(round(T_eff, sigdigits=3)) GB/s-----") end

    # visualization
    if do_visu && (it%n_vis == 0) && me==0
      T_inn .= Array(T[2:end-1,2:end-1]); gather!(T_inn, T_v)
      opts = (aspect_ratio=1, xlims=(xc[1], xc[end]), ylims=(yc[1], yc[end]), xlabel="Lx", ylabel="Ly", title="time = $(round(it*dt, sigdigits=3))")
      display(heatmap(xc, yc, Array(T_v)'; opts...)); frame(anim)
    end

  end # end physical stepping

  # Save animation
  if (do_visu && me==0) gif(anim, "TM_2D.gif", fps = 5)  end

  finalize_global_grid()
  return
end

TM_2D(; do_visu = true, do_print = true)


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule fpzip_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("fpzip")
JLLWrappers.@generate_main_file("fpzip", UUID("dabd510a-92a7-54df-aa64-c306dd35554a"))
end  # module fpzip_jll


text:
function trival((twos, threes, fives))
    BigInt(2)^twos * BigInt(3)^threes * BigInt(5)^fives
end

mutable struct Hammings
    ham532 :: Vector{Tuple{Float64,Tuple{Int,Int,Int}}}
    ham53 :: Vector{Tuple{Float64,Tuple{Int,Int,Int}}}
    ndx532 :: Int
    ndx53 :: Int
    next2 :: Tuple{Float64,Tuple{Int,Int,Int}}
    next3 :: Tuple{Float64,Tuple{Int,Int,Int}}
    next5 :: Tuple{Float64,Tuple{Int,Int,Int}}
    next53 :: Tuple{Float64,Tuple{Int,Int,Int}}
    Hammings() = new(
        Vector{Tuple{Float64,Tuple{Int,Int,Int}}}(),
        Vector{Tuple{Float64,Tuple{Int,Int,Int}}}(),
        1, 1,
        (1.0, (1, 0, 0)), (log(2, 3), (0, 1, 0)),
        (log(2, 5), (0, 0, 1)), (0.0, (0, 0, 0))
    )
end
Base.eltype(::Type{Hammings}) = Tuple{Int,Int,Int}
function Base.iterate(HM::Hammings, st = HM) # :: Union{Nothing,Tuple{Tuple{Float64,Tuple{Int,Int,Int}},Hammings}}
    log2of2, log2of3, log2of5 = 1.0, log(2,3), log(2,5)
    if st.next2[1] < st.next53[1]
        push!(st.ham532, st.next2); st.ndx532 += 1
        last, (twos, threes, fives) = st.ham532[st.ndx532]
        st.next2 = (log2of2 + last, (twos + 1, threes, fives))
    else
        push!(st.ham532, st.next53)
        if st.next3[1] < st.next5[1]
            st.next53 = st.next3; push!(st.ham53, st.next3)
            last, (_, threes, fives) = st.ham53[st.ndx53]; st.ndx53 += 1
            st.next3 = (log2of3 + last, (0, threes + 1, fives))
        else
            st.next53 = st.next5; push!(st.ham53, st.next5)
            last, (_, _, fives) = st.next5
            st.next5 = (log2of5 + last, (0, 0, fives + 1))
        end
    end
    len53 = length(st.ham53)
    if st.ndx53 > (len53 >>> 1)
        nlen53 = len53 - st.ndx53 + 1
        copyto!(st.ham53, 1, st.ham53, st.ndx53, nlen53)
        resize!(st.ham53, nlen53); st.ndx53 = 1
    end
    len532 = length(st.ham532)
    if st.ndx532 > (len532 >>> 1)
        nlen532 = len532 - st.ndx532 + 1
        copyto!(st.ham532, 1, st.ham532, st.ndx532, nlen532)
        resize!(st.ham532, nlen532); st.ndx532 = 1
    end
    _, tri = st.ham532[end]
    tri, st
#    convert(Union{Nothing,Tuple{Tuple{Float64,Tuple{Int,Int,Int}},Hammings}},(st.ham532[end], st))
#   (length(st.ham532), length(st.ham53)), st
end

foreach(x -> print(trival(x)," "), (Iterators.take(Hammings(), 20))); println()
let count = 1691; for t in Hammings() count <= 1 && (println(trival(t)); break); count -= 1 end end
let count = 1000000; for t in Hammings() count <= 1 && (println(trival(t)); break); count -= 1 end end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
function MassMatrix{T<:AbstractFloat}(vp::Array{T,2},nz::Int,nx::Int,ext::Int)

# This function builds the diagonal mass matrix/operator, where the elements
# of the main diagonal contain the slowness-squared (1 over velocity-squared)
# for each point in the modeling region AND the surrounding boundary region.
#
# INPUTS:     vp        - Acoustic velocity model
#             nz        - Number of grid points in z-direction INCLUDING the absorbing boundary region
#             nx        - Number of grid points in x-direction INCLUDING the absorbing boundary region
#             ext       - Thickness (number of grid points) of absorbing boundary region
#
# OUTPUTS:    M         - Diagonal matrix containing the slowness-squared for each point within the model and boundary region

    M = zeros(T,nz,nx)
    M[ext+1:end-ext,ext+1:end-ext] = 1./(vp.^2)
    M[1:ext,1:ext] = 1./(vp[1,1].^2)
    M[1:ext,nx-ext+1:nx] = 1./(vp[1,end].^2)
    M[nz-ext+1:nz,1:ext] = 1./(vp[end,1].^2)
    M[nz-ext+1:nz,nx-ext+1:nx] = 1./(vp[end,end].^2)
    M[ext+1:nz-ext,1:ext] = 1./(kron(vp[:,1],ones(1,ext)).^2)
    M[ext+1:nz-ext,nx-ext+1:nx] = 1./(kron(vp[:,end],ones(1,ext)).^2)
    M[1:ext,ext+1:nx-ext] = 1./(kron(vp[1,:]',ones(ext,1)).^2)
    M[nz-ext+1:nz,ext+1:nx-ext] = 1./(kron(vp[end,:]',ones(ext,1)).^2)

    return spdiagm((M[:]),(0))

end