text:
module AlignmentStatistics

# package code goes here
using FastaIO, StatsBase, HypothesisTests, MultipleTesting,
DataFrames, Distributions, GaussDCA

export AAindex1_to_Dict,
read_AAindex1,
export_fasta,
get_sequence_labels,
split_fasta_sequences,
split_sequences,
clean_sequences,
sequence_composition,
sequences_to_numerical_properties,
PDB_AA_to_1_letter_code,
sequence_lengths,
symbols_in_sequences,
majority_consensus,
extract_label_element,
sorted_label_frequencies,
rectangularize_alignment
# binomial_CIs,
# compute_distances_dca_scores_table,
# dca_with_reference_sequence,
# Dirichlet_K3,
# Fisher_test_sequence_sets,
# PDB_to_single_chain_fasta,
# random_rows_from_msa,
# reference_sequence_column_labels,
# sequence_composition_difference,

PDB_AA_to_1_letter_code = Dict(
"ALA" => "A",
"ARG" => "R",
"ASN" => "N",
"ASP" => "D",
"CYS" => "C",
"GLN" => "Q",
"GLU" => "E",
"GLY" => "G",
"HIS" => "H",
"ILE" => "I",
"LEU" => "L",
"LYS" => "K",
"MET" => "M",
"PHE" => "F",
"PRO" => "P",
"SER" => "S",
"THR" => "T",
"TRP" => "W",
"TYR" => "Y",
"VAL" => "V"
)

"""
Input: fasta data read with FastaIO function readfasta

Output: sequences split up into single characters.

"""
function split_fasta_sequences(fasta::Array{Any,1})
    return map(i -> convert(Array{String,1},split(fasta[i][2],"")), 1:length(fasta))
end

"""
Can apply the following operations to sequences read with FastaIO.readfasta:

1. remove last char of each sequence (often "*"),

2. remove sequences that contain any of a set of given substrings, e.g. "X".

3. remove sequences that have a length different from a required length

Input:

- fasta array as read with FastaIO.readfasta

- optionally: should last character be removed of each sequence (default: remove_last_char=false); removing the last character is done first (if requested)

- optionally: remove sequences with any of the given strings (default: remove_seqs_with=["X", "\#", "*"])

- optionally: remove sequences with a length different from a required length (default: required_length = 0, i.e. option switched off)

    - Consider also using the rectangularize_alignment command as a way to produce alignment with homogeneous lengths

Output:

- Cleaned-up fasta (works with a copy of the original data)

"""
function clean_sequences(fasta::Array{Any,1};
                         remove_last_char::Bool=false,
                         remove_seqs_with::Array{String,1}=["X","#","*"],
                         required_length::Int64=0)
    n_seqs=length(fasta)
    new_fasta = copy(fasta)
    clean_seqs = fill(true, n_seqs)
    n_removal_triggers = length(remove_seqs_with)
    if remove_last_char
        for i in 1:n_seqs
            new_fasta[i] = (new_fasta[i][1], new_fasta[i][2][1:(end-1)])
        end
    end
    for i in 1:n_seqs
        for j in 1:n_removal_triggers
            if contains(new_fasta[i][2],remove_seqs_with[j])
                clean_seqs[i] = false
                break
            end
        end
    end
    if required_length != 0
        for i in 1:n_seqs
            for j in 1:n_removal_triggers
                if (clean_seqs[i] == true) &
                   (length(new_fasta[i][2])!=required_length)
                    clean_seqs[i] = false
                    break
                end
            end
        end
    end
    return new_fasta[clean_seqs]
end

# """
# Input:

# - fasta_msa: MSA as read with readfasta of FastaIO

# - n_select: number of rows to be selected

# Output:

# - MSA in fasta format consisting of n_select rows of fasta_msa
# """
# function random_rows_from_msa(fasta_msa::Array{Any,1}, n_select::Int64)
#     n_rows = length(fasta_msa)
#     if n_select >= n_rows
#         error("requested more random rows than available")
#     end
#     fasta_msa[sample(1:n_rows, n_select, replace=false)]
# end

"""
Input: array of strings (=sequences)

Output: sequences split up into single characters.

    Hint: if you have read the sequences from a fasta file with readfasta (package FastaIO), use split_fasta_sequences instead of split_sequences.

"""
function split_sequences(seqs::Any)
    return map(i -> convert(Array{String,1},split(seqs[i],"")), 1:length(seqs))
end

# """
#     Input: fasta data (in general multiple sequences) read with FastaIO function readfasta.

#     Output: vector of sequences lengths.
# """
# function fasta_lengths(fasta::Array{Any,1})
#     return map(i -> length(fasta[i][2]), 1:length(fasta))
# end

"""
Takes an alignment, as read by FastaIO.readall, and appends gap symbols '-' to shorter sequences to make all sequences the same length. Returns this rectangularized sequence array.
"""
function rectangularize_alignment(raw_ali::Array{Any,1})
    seqs = map(x -> convert(Array{String,1},split(x[2],"")), raw_ali[1:end])
    maxlen = maximum(map(x -> length(x), seqs))

    #Fill everything with gaps to make the array rectangular (otherwise we have dangling ends):
    Seqs = fill("-",length(seqs),maxlen);

    for i in 1:length(seqs)
        L = length(seqs[i])
        for j in 1:L
            Seqs[i,j] = seqs[i][j]
        end
    end

    return Seqs
end

"""
Returns sequence labels from sequences read with FastaIO.readall (without the leading '>').
"""
function get_sequence_labels(raw_ali::Array{Any,1})
    return convert(Array{String,1},map(x -> x[1], raw_ali[1:end]))
end

"""
Input: rectangularized sequence array.

Output: array of three components:

(1) consensus sequence where at each position the majority symbol in the respective alignment column is given;

(2) absolute frequency of majority symbols;

(3) relative frequency of majority symbols (absolute / n_rows).
"""
function majority_consensus(seqs::Array{String,2})

    n_rows, n_cols = size(seqs)

    maxelem = fill("-", n_cols)
    maxcount = zeros(Int64, n_cols)

    for i in 1:n_cols
        coltab = countmap(seqs[:,i])
        vals = collect(values(coltab))
        maxval = maximum(vals)
        maxcount[i] = maxval
        maxelem[i] = collect(keys(coltab))[vals .== maxval][1]
    end

    return maxelem, maxcount, (convert(Array{Float64,1},maxcount) ./ n_rows)
end

"""
Input: sequence (ASCIIArray)

Output: dictionary of sequence composition
"""
function sequence_composition(seq::Array{String,1})
    return countmap(seq)
end

"""
Input:

(1) String array of sequence labels (eg from FASTA header)

(2) delimiter character

(3) number of element to be extracted

Output: array of extracted label elements (Strings)
"""
function extract_label_element(labels::Array{String,1}, delim::String, elno::Int64)
    return map(x -> convert(String,x[elno]), map(x -> split(x,delim), labels))
end

"""
# Input: array of labels.

# Output: 2D array with labels in first row and frequencies in second, sorted according to decreasing frequency.
# """
function sorted_label_frequencies(labels::Array{String, 1})
    label_counts = countmap(labels)
    p = sortperm(collect(values(label_counts)),rev=true)
    return permute!(collect(keys(label_counts)), p), permute!(collect(values(label_counts)), p)
end

"""
Input:

(1) symbols: sequence of symbols (String) of the same length as each of the sequences in seqs.

(2) seqs: rectangular set of sequences (2D String array)

Output:

(1) absolute frequency of symbols in the respective columns of seqs.

(2) relative frequency of symbols in respective columns of seqs (absolute / n_rows).

"""
function symbols_in_sequences(symbols::Array{String,1}, seqs::Array{String,2})

    n_rows, n_cols = size(seqs)

    symbol_count = zeros(Int64, n_cols)

    for i in 1:n_cols
        symbol_count[i] = count(x -> x == symbols[i], seqs[:,i])
    end

    return symbols, symbol_count, (convert(Array{Float64,1},symbol_count) ./ n_rows)

end

# """
# Input: two sequence sets, seqsA and seqsB, and a sequence of reference symbols, ref_syms.

# Optionally:
# - set parameter 'FDR' true to correct for multiple testing by control of false discovery rate (method of Benjamini and Hochberg); default value is false.
# - choose 'tail' of Fisher test between :left, :both, and :right

# Output: p-values from Fisher's exact test carried out for each column;
#     we test with 2x2 contingency table:
#     n_ref_sym_i(seqsAi), n_rows_A - n_ref_sym_i(seqsAi)
#     n_ref_sym_i(seqsBi), n_rows_B - n_ref_sym_i(seqsBi)
# """
# function Fisher_test_sequence_sets(
#     seqsA::Array{String,2},
#     seqsB::Array{String,2},
#     ref_syms::Array{String,1};
#     FDR::Bool=false,
#     tail=:both
# )
#     rowsA, colsA = size(seqsA)
#     rowsB, colsB = size(seqsB)
#     len_ref = length(ref_syms)

#     if colsA != colsB || colsA != len_ref || colsB != len_ref
#         error("Fisher_test_sequence_sets: inconsistent lengths of sequences.")
#     end

#     pvals = zeros(len_ref)

#     #go through all columns and carry out test
#     for i in 1:len_ref
#         x = ref_syms[i]
#         nAx = count(y -> y == x, seqsA[:,i])
#         nBx = count(y -> y == x, seqsB[:,i])
#         if nAx + nBx != 0 && (nAx != rowsA || nBx != rowsB)
#             pvals[i] = pvalue(FisherExactTest(nAx, rowsA-nAx, nBx, rowsB-nBx), tail=tail)
#         else
#             pvals[i] = 1.0
#         end
#     end

#     if FDR==true #if user wishes FDR correction for multiple testing
#         return adjust(pvals, BenjaminiHochberg())
#     else #no correction
#         return pvals
#     end

# end

# """
#     Computes column labels of a multiple sequence alignment based on a reference sequence

#     Input:

#     - label of reference sequence (FASTA header without "")

#     - array of all labels (1D ASCII string array)

#     - sequence array (2D ASCII array)

#     Output: Labels for each column:

#     - columns in which the reference sequence has a letter are labelled with the number of this letter in the sequence (without gaps)

#     - columns in which the reference sequence has a gap are labelled with the number of the closest left neighbor letter and the number of gap symbols to the current column
# """

# function reference_sequence_column_labels(ref_label::String,
#                                           labels::Array{String,1},
#                                           seqs::Array{String,2})
#     ref_seq = vec(seqs[labels .== ref_label,:])
#     ali_len = length(ref_seq)
#     res_num = 0
#     gap_num = 1
#     ref_col = fill(" ", ali_len)
#     for i in 1:ali_len
#         if ref_seq[i] == "-"
#             ref_col[i] = join([string(res_num),".",string(gap_num)],"")
#             gap_num += 1
#         else
#             gap_num = 1
#             res_num += 1
#             ref_col[i] = string(res_num)
#         end
#     end
#     return ref_col
# end

"""
     Input:
     - filename for fasta export
     - label vector (1D String array)
     - sequence vector (2D String array)

     Output:
     writes a fasta file with given labels as headers
"""
function export_fasta(
                                filename::String,
                                labels::Array{String,1},
                                seqs::Array{String,2}
                                )
    n_seqs = length(labels)
    to_export = Array{Any}(n_seqs)
    for i in 1:n_seqs
        to_export[i] = (labels[i], join(collect(seqs[i,:]),""))
    end
    writefasta(filename, to_export)
end

"""
     Input:
     - filename for fasta export
     - label vector (1D String array)
     - sequence vector (array of String arrays)

     Output:
     writes a fasta file with given labels as headers
"""
function export_fasta(
                                filename::String,
                                labels::Array{String,1},
                                seqs::Array{String,1}
                                )
    n_seqs = length(labels)
    to_export = Array{Any}(n_seqs)
    for i in 1:n_seqs
        to_export[i] = (labels[i], seqs[i])
    end
    writefasta(filename, to_export)
end

# """
#     Input:
#     - array of successes x (Int64)
#     - array of trials n (Int64)

#     Output: three arrays
#     - fractions of successes
#     - lower boundaries of 95% confidence intervals (CIs)
#     - upper boundaries of 95% confidence intervals

#     If no x = n or x = 0, the CIs have a width of zero
# """
# function fractions_and_binomial_CIs(x::Array{Int64,1}, n::Array{Int64,1})

#     len = length(x)
#     if (len != length(n))
#         error("lengths of x and n do not match")
#     elseif sum(x .> n) > 0
#         error("x > n")
#     end

#     fraction = zeros(len)
#     lowerCI = zeros(len)
#     upperCI = zeros(len)

#     for j in 1:len
#         if ((!isnan(n[j])) & (n[j] != 0))
#             fraction[j] = x[j]/n[j]
#             if x[j] == n[j]
#                 (lowerCI[j],upperCI[j]) = (1.0,1.0)
#             elseif x[j] == 0
#                 (lowerCI[j],upperCI[j]) = (0.0,0.0)
#             else
#                 (lowerCI[j],upperCI[j]) = ci(BinomialTest(x[j],n[j]))
#             end
#         else
#             fraction[j] = NA
#         end
#     end

#     return fraction, lowerCI, upperCI
# end

"""
   Input:

     - translation dictionary (AA -> numerical property)

     - sequences (array of arrays of strings)

   Output:

     - sequences translated to array of array of numerical properties

"""
function sequences_to_numerical_properties(trans::Dict, seqs::Array{Array{String,1},1})
    nseqs = length(seqs)
    props = Array{Any}(nseqs)
    for i in 1:nseqs
        props[i] = map(aa -> trans[aa], seqs[i])
    end
    return convert(Array{Array{Float64,1}},props)
end

"""
  Input:

     - AAindex1 data frame (see function read_AAindex1())

     - index to be used

  Output:

     - dictionary relating AAs with numbers (e.g. "A" => 3.95, ...)
"""
function AAindex1_to_Dict(AAindex1::DataFrame, ix::Int64)
    nAAs = 20
    offset = 2
    AAs = map(x->string(x), names(AAindex1)[3:22])
    trans =
         Dict(AAs[i] => convert(Float64,AAindex1[ix, i+offset]) for i in 1:nAAs)
    return trans
end

"""
   Input: none (reads file AAindex1.csv from data directory of
         package AlignmentStatistics)

   Output: DataFrame of AAindex1

"""
function read_AAindex1()
    readtable(Pkg.dir("AlignmentStatistics","data","AAindex1.csv"))
end

# """
# Input:

#  - Three element Float64 array alpha, e.g. probabilities or counts of an AA to have charge -1, 0, 1.
#    All elements should be positive.

#  - Optional: minimum probability pmin, maximum probability pmax, number of steps np.

# Output:

#  - Array of probability densities computed for a Dirichlet distribution with concentration parameter
#      vector alpha, evaluated at (p_m1 = pmin, ..., pmax; p_0 = pmin, ..., pmax; p_1 = 1.0-p_m1-p_0).
#      This can e.g. correspond to the probabilities for the three charge states -1 (m1), 0 and 1.

# """
# function Dirichlet_K3(alpha::Array{Float64,1}; pmin::Float64=0.01, pmax::Float64=0.99, np::Int64=100)
#     if length(alpha) != 3
#         error("needs alpha of length 3")
#     elseif sum(alpha .<= 0.0) != 0
#         error("all elements of alpha must be > 0")
#     elseif (pmin <= 0.) | (pmax <= 0.)
#         error("pmin and pmax must be positive")
#     elseif pmax <= pmin
#         error("pmax must be > pmin")
#     end

#     p_m1 = linspace(pmin,pmax,np)
#     p_0 = linspace(pmin,pmax,np)
#     pDir = zeros(np,np)
#     for i in 1:np
#         pm1 = p_m1[i]
#         for j in 1:np
#             p0 = p_0[j]
#             if pm1 + p0 >= 1.
#                 break
#             end
#             p1 = 1.0 - pm1 - p0
#             pDir[i,j] = pdf(Dirichlet(alpha),[pm1;p0;p1])
#         end
#     end
#     return collect(p_m1), collect(p_0), pDir
# end

# AAs = ["A";"C";"D";"E";"F";"G";"H";"I";"K";"L";"M";"N";"P";"Q";"R";"S";"T";"V";"W";"Y"]
# DNAs = ["A";"C";"G";"T"]
# RNAs = ["A";"C";"G";"U"]

# """
# Input:
#          - string s (either String or Array{String,1})
#          - alphabet

# Output:
#          DataFrame with columns:
#          - alphabet: elements of the alphabet
#          - abs: absolute frequency of each element of the alphabet in s
#          - rel: relative frequency of each element of the alphabet in s

#          Note that the relative frequency does add up to 1.
#          The absolute frequency is lower than the number of characters in s if s
#              contains characters that are not elements of the alphabet.
# """
# function sequence_composition(s::Any, alphabet::Array{String,1})
#     if typeof(s)==String
#         a = split(s,"")
#     elseif typeof(s) != Array{String,1}
#         error("input has to be String or Array{String,1}")
#     end
#     n = length(alphabet)
#     c = zeros(Int64,n)
#     i = 1
#     result = DataFrame(alphabet = alphabet , abs = zeros(Int64,n), rel = zeros(n))
#     for letter in alphabet
#         c[i] = count(x -> x==letter, a)
#         i+=1
#     end
#     result[:abs] = c
#     result[:rel] = c/sum(c)
#     result
# end

# """
# Input: two DataFrames as returned by function sequence_composition, one DataFrame
#     for each sequence s1 and s2

# Output: \sum_{i \in alphabet} |f_i(s_1)-f_i(s_2)| with relative frequencies f_i

# """
# function sequence_composition_difference(c1::DataFrame, c2::DataFrame)
#     if c1[:alphabet] != c2[:alphabet]
#         error("different alphabets")
#     end
#     sum(abs(c1[:rel] .- c2[:rel]))
# end

"""
Sequence length distribution (as table) in a fasta array, read with FastaIO.readfasta.

Input:

- fasta array

Output:

- table of sequence lengths, generated with countmap

"""
function sequence_lengths(fasta::Array{Any,1})
    countmap(map(x->length(fasta[x][2]), 1:length(fasta)))
end

# """
# Input:

#          - pdb_input: pdb file name

#          - fasta_out_name: fasta output file name

#          - fasta header: header of fasta output file

# Output:

#          - fasta file with sequence extracted from pdb input and interpreted as single chain
# """
# function PDB_to_single_chain_fasta(
#                                    pdb_input::String,
#                                    fasta_out_name::String,
#                                    fasta_header::String
#                                    )
#     pdb = read(pdb_input,PDB)
#     CAs = collectatoms(pdb, calphaselector)
#     seq = (map(i -> PDB_AA_to_1_letter_code[CAs[i].res_name], 1:length(CAs)))'
#     export_fasta(fasta_out_name, [fasta_header], seq)
# end

# """
# Given a PDB of a protein (interpreted as single chain) and a MSA of sequences related to this protein, the function computes a table of distances and Direct Coupling Analysis (DCA) scores for all Calpha pairs.

# Input:

#     - pdb_file: name of PDB file

#     - msa_file: name of MSA file

#     - clean_up_msa=true (default): apply some clean-up on the MSA, e.g. give all columns the same lengths, remove rows (=sequences) with non-AA symbols. To switch this off, use clean_up_msa=false.

#     - remove_last_char=false (default): sometimes fasta blocks in MSA are finished with a "*" symbol. If remove_last_char=true, the last symbol is removed from all fasta blocks (one per block).

#     - dca_type="gDCA_FRN" (default): use package GaussDCA with Frobenius norm as score. Alternative are currently: "gDCA_DI" (direct information), and a file name from which an externally computed DCA output can be read (text array of tuples (i,j,score)) [the latter is currently not available].

#     - contact_threshold=8 (default): contact threshold in Angstrom

# Output:

#     - data frame with columns i, j, dij (=Calpha-Calpha distance), score (=DCA score), contact, true_positive_rate
# """
# function compute_distances_dca_scores_table(
#                 pdb_file::String,
#                 msa_file::String;
#                 clean_up_msa::Bool=true,
#                 remove_last_char::Bool=false, #remove last char of each fasta block
#                 dca_type::String="gDCA_FRN",
#                 contact_threshold::Float64=8.0
#             )

#     #read pdb file and extract a fasta file -> write fasta
#     pdb_fasta_file = pdb_file*".fa"
#     PDB_to_single_chain_fasta(pdb_file, pdb_fasta_file, pdb_file)

#     #if requested: clean up MSA and write file with cleaned-up
#     if clean_up_msa
#         seq = readfasta(msa_file)
#         clean_seq = clean_sequences(seq, remove_last_char=remove_last_char)
#         rect_seq = rectangularize_alignment(clean_seq)
#         msa_file = msa_file*"_clean.fa"
#         export_fasta(msa_file, get_sequence_labels(clean_seq), rect_seq)
#     end

#     #generate MSA from MSA file and pdb-fasta, with pdb-fasta as last entry
#     msa_out = (splitdir(pdb_file)[end])*(splitdir(msa_file)[end])
#     run(pipeline(`mafft --add $pdb_fasta_file $msa_file`, stdout=msa_out))

#     #compute DCA based on current MSA
#     dca_out = 0
#     if dca_type == "gDCA_FRN"
#         dca_out = gDCA(msa_out)
#     elseif dca_type == "gDCA_DI"
#         dca_out = gDCA(msa_out, pseudocount=0.2, score=:DI)
#     else # the DCA has been computed beforehand; dca_type is interpreted as name
#          # of file that contains tuples (i,j,score)
#         error("external DCA currently not available")
#         external = readdlm(dca_type,';')
#         n = length(external)
#         dca_out = Array{Any,1}(n)
#         for ix in 1:n
#             s = split(external[ix],['(',',',')'])
#             dca_out[ix]= (parse(s[2]),parse(s[3]),parse(s[4]))
#         end
#     end

#     #put DCA results in data frame for later join
#     n_dca = length(dca_out)
#     dca_df = DataFrame(
#         i_dca = map(i->dca_out[i][1],1:n_dca),
#         j_dca = map(i->dca_out[i][2],1:n_dca),
#         score = map(i->dca_out[i][3],1:n_dca)
#     )

#     #prepare computation of Calpha-Calpha distances
#     struct = read(pdb_file,PDB)
#     struct_ca = collectatoms(struct,calphaselector)

#     #find mapping of pdb-fasta to MSA columns necessary to map distances to DCA scores
#     struct_msa = readfasta(msa_out)
#     rect_struct_msa = rectangularize_alignment(struct_msa)

#          #the structure-sequence is the last sequence in the MSA
#     seq_struct_ali = collect(rect_struct_msa[end,:])
#     seq_struct = seq_struct_ali[find(x->x!="-", seq_struct_ali)]
#     indices = zeros(Int64,length(seq_struct))
#     j = 1
#     for i in 1:length(seq_struct_ali)
#         if (seq_struct_ali[i]!="-")
#             indices[j]=i
#             j += 1
#         end
#     end
#     #indices contains now the indices of the columns in which the pdb-fasta residues are

#     #generate data frame of Calpha-Calpha distances
#     n_ca = length(struct_ca)
#     n_pairs = round(Int64,n_ca*(n_ca-1)/2)

#     Dij = DataFrame()
#     Dij[:ix] = zeros(Int64,n_pairs)
#     Dij[:jx] = zeros(Int64,n_pairs)
#     Dij[:i_dca] = zeros(Int64,n_pairs)
#     Dij[:j_dca] = zeros(Int64,n_pairs)
#     Dij[:dij]=zeros(n_pairs)
#     k=1
#     for i in 1:(n_ca-1)
#         for j in (i+1):n_ca
#             dr = struct_ca[i].coords-struct_ca[j].coords
#             Dij[k,:dij] = sqrt(dot(dr,dr))
#             Dij[k,:ix] = i
#             Dij[k,:jx] = j
#             Dij[k,:i_dca] = indices[i]
#             Dij[k,:j_dca] = indices[j]
#             k += 1
#         end
#     end

#     #join data frame with distances and DCA scores
#     Dij = join(Dij, dca_df, on=[:i_dca,:j_dca])

#     #sort according to decreasing DCA scores
#     Dij = sort(Dij, cols = [:score], rev=true)

#     #add contact columns
#     Dij[:contact] = (Dij[:dij] .<= contact_threshold)

#     #compute true positive rate up to distance pair i, j
#     Dij[:true_positive_rate] = cumsum(Dij[:contact]) ./ (1:length(Dij[:contact]))

#     #output
#     Dij
# end

# """
# Given an MSA of sequences of a protein and a reference sequence in that MSA given by the fasta header (without initial ">") the function computes DCA scores and provides the corresponding pair indices of the reference sequence in the output data frame.

# Input:

#     - msa_file: name of MSA file

#     - reference_header: fasta header of reference sequence (without initial ">")

#     - clean_up_msa=true (default): apply some clean-up on the MSA, e.g. give all columns the same lengths, remove rows (=sequences) with non-AA symbols. To switch this off, use clean_up_msa=false.

#     - remove_last_char=false (default): sometimes fasta blocks in MSA are finished with a "*" symbol. If remove_last_char=true, the last symbol is removed from all fasta blocks (one per block).

#     - dca_type="gDCA_FRN" (default): use package GaussDCA with Frobenius norm as score. Alternative are currently: "gDCA_DI" (direct information), and a file name from which an externally computed DCA output can be read (text array of tuples (i,j,score)) [the latter is currently not available].

#     - data frame with columns iref, jref, imsa, jmsa, score (=DCA score),
# """
# function dca_with_reference_sequence(
#                 reference_label::String,
#                 msa_file::String;
#                 clean_up_msa::Bool=true,
#                 remove_last_char::Bool=false, #remove last char of each fasta block
#                 dca_type::String="gDCA_FRN"
#             )

#     #if requested: clean up MSA and write file with cleaned-up
#     if clean_up_msa
#         seq = readfasta(msa_file)
#         clean_seq = clean_sequences(seq, remove_last_char=remove_last_char)
#         rect_seq = rectangularize_alignment(clean_seq)
#         msa_file = msa_file*"_clean.fa"
#         export_fasta(msa_file, get_sequence_labels(clean_seq), rect_seq)
#     end

#     #compute DCA based on MSA
#     dca_out = 0
#     if dca_type == "gDCA_FRN"
#         dca_out = gDCA(msa_file)
#     elseif dca_type == "gDCA_DI"
#         dca_out = gDCA(msa_file, pseudocount=0.2, score=:DI)
#     else # the DCA has been computed beforehand; dca_type is interpreted as name
#          # of file that contains tuples (i,j,score)
#         error("external DCA currently not available")
#         external = readdlm(dca_type,';')
#         n = length(external)
#         dca_out = Array{Any,1}(n)
#         for ix in 1:n
#             s = split(external[ix],['(',',',')'])
#             dca_out[ix]= (parse(s[2]),parse(s[3]),parse(s[4]))
#         end
#     end

#     #put DCA results in data frame for later join
#     n_dca = length(dca_out)
#     dca_df = DataFrame(
#         imsa = map(i->dca_out[i][1],1:n_dca),
#         jmsa = map(i->dca_out[i][2],1:n_dca),
#         score = map(i->dca_out[i][3],1:n_dca)
#     )

#     #rectangularize MSA and compute reference
#     msa = readfasta(msa_file)
#     rect_msa = rectangularize_alignment(msa)
#     labels = get_sequence_labels(msa)
#     ref_cols = reference_sequence_column_labels(reference_label, labels, rect_msa)

#     dca_df[:iref] = ref_cols[dca_df[:imsa]]
#     dca_df[:jref] = ref_cols[dca_df[:jmsa]]

#     dca_df

# end

end # module


text:
using NonuniformResampling1D
using Plots
using Statistics

xin = 1:24
yin = [
    0.10221248850845677,
    -0.3017964617372699,
    0.6313812022052889,
    -0.469567525619184,
    -0.742550821690563,
    0.538340882717059,
    0.7645542960785265,
    -0.42847475695892334,
    -0.2030496164626363,
    -0.6961115394075117,
    0.9553667194266671,
    0.6773297980833843,
    0.6595211409431194,
    0.18762849994721864,
    0.9676546988790669,
    -0.6599170737336646,
    0.17851166014365583,
    0.7417238170162541,
    -0.7312025977875327,
    -0.4015347276663681,
    -0.685997771780074,
    -0.7558107010474338,
    0.12236047934097583,
    -0.8627888068515903]
xout = [4.6, 5.7, 8.2, 13.3, 17.3, 18.2, 19.1, 20, 20.8]
red_ind = 3
blue_ind = 6
red_color = :coral
blue_color = :dodgerblue

w = 1.35
win = hann_window(w)
N = 300
win_y = win.(range(0, w, length=N))
required_points_per_slice = 4

red_left_border = xout[red_ind]*(1-w)+w*xout[red_ind-1]
red_right_border = xout[red_ind]*(1-w)+w*xout[red_ind+1]
blue_left_border = xout[blue_ind]*(1-w)+w*xout[blue_ind-1]
blue_right_border = xout[blue_ind]*(1-w)+w*xout[blue_ind+1]

yout = nuresample(xin, yin, xout, win, required_points_per_slice=required_points_per_slice)

function vertical_lines()
    vline!([xout[red_ind]], c=red_color, linewidth=1)
    vline!([red_left_border, red_right_border], c=red_color, line=:dash)
    vline!([xout[blue_ind]], c=blue_color, linewidth=1)
    vline!([blue_left_border, blue_right_border], c=blue_color, line=:dash)
    # unit lines make the image too crowded
    #vline!([xout[red_ind-1], xout[red_ind+1], xout[blue_ind-1], xout[blue_ind+1]], c=:black, line=:dot)
    hline!([0], c=:gray, z_order=:back)
end

# blue upsampling + points
upsample_step = min(xout[blue_ind]-xout[blue_ind-1], xout[blue_ind+1]-xout[blue_ind])/
    required_points_per_slice*w
upsample_x = range(blue_left_border+upsample_step/2, step=upsample_step,
stop=blue_right_border-upsample_step/2)
in_step = Float64(StepRangeLen(xin).step)
upsample_y = [NonuniformResampling1D.interpolate_point(
    StepRangeLen(xin), yin, up_x, in_step, in_step, lanczos_window()) for up_x in upsample_x]

# assume xout around blue is symmetrical
blue_weights = @. win(abs(upsample_x-xout[blue_ind])/(xout[blue_ind]-xout[blue_ind-1]))
blue_weighted = blue_weights .* upsample_y

plt_in = begin
    plot(xin, yin, m=:circle, c=:black, legend=false, xaxis=false)
    # red window
    red_x_left = range(stop=xout[red_ind], length=N, step=(xout[red_ind]-xout[red_ind-1])/N*w)
    plot!(red_x_left, reverse(win_y), c=red_color)
    red_x_right = range(start=xout[red_ind], length=N, step=(xout[red_ind+1]-xout[red_ind])/N*w)
    plot!(red_x_right, win_y, c=red_color)
    #blue window
    blue_x_left = range(stop=xout[blue_ind], length=N, step=(xout[blue_ind]-xout[blue_ind-1])/N*w)
    plot!(blue_x_left, reverse(win_y), c=blue_color)
    blue_x_right = range(start=xout[blue_ind], length=N, step=(xout[blue_ind+1]-xout[blue_ind])/N*w)
    plot!(blue_x_right, win_y, c=blue_color)
    plot!(upsample_x, upsample_y, c=:deepskyblue, m=:x, markersize=2)

    annotate!([
        (4, -1.1, ("adaptive, asymmetric\nwindow width", 9, red_color)),
        (21.5, .8, ("ad-hoc\nupsampling", 9, blue_color))
    ])

    vertical_lines()
    title!("Input", titlefontsize=13)
end

plt_windowed = begin
    red_start_ind = NonuniformResampling1D.find_first_above_or_equal(red_left_border, StepRangeLen(xin))
    red_middle_ind = NonuniformResampling1D.find_last_below(xout[red_ind], StepRangeLen(xin))
    red_stop_ind = NonuniformResampling1D.find_last_below_or_equal(red_right_border, StepRangeLen(xin))

    red_left_x = xin[red_start_ind:red_middle_ind]
    red_right_x = xin[red_middle_ind+1:red_stop_ind]

    red_left_weights = @. win(abs(red_left_x-xout[red_ind])/(xout[red_ind]-xout[red_ind-1]))
    red_right_weights = @. win(abs(red_right_x-xout[red_ind])/(xout[red_ind+1]-xout[red_ind]))

    red_y = yin[red_start_ind:red_stop_ind]
    red_weighted = [red_left_weights; red_right_weights] .* red_y
    plot(red_start_ind:red_stop_ind, red_weighted,
        m=:circle, c=red_color, legend=false, markerstrokewidth=0)

    plot!(upsample_x, blue_weighted, m=:circle, c=blue_color, markerstrokewidth=0)

    vertical_lines()
    title!("Example Slices Windowed", titlefontsize=13)
end

plt_out = begin
    plot(xout, yout, m=:circle, c=:black, legend=false)
    scatter!([xout[red_ind]], [yout[red_ind]], c=red_color, markerstrokewidth=0)
    scatter!([xout[blue_ind]], [yout[blue_ind]], c=blue_color, markerstrokewidth=0)

    quiver!([xout[red_ind], xout[red_ind]], [-.5, -1.0], quiver=(-[(xout[red_ind]-xout[red_ind-1]), (xout[red_ind]-red_left_border)], [0,0]), c=:black)
    annotate!([
        (mean(xout[red_ind-1:red_ind]), -.7, ("Unit", 9)),
        (mean([xout[red_ind], red_left_border]), -1.2, ("Slice", 9))
    ])

    vertical_lines()
    title!("Output", titlefontsize=13)
end

plot(plt_in, plt_windowed, plt_out, layout=(3,1), link=:x, xticks=false, xaxis=false, yaxis=false)
ylims!(-1.1, 1.1)
savefig("explanation_diagram.svg")

ext: jl


text:
@testset "Not closed error"  begin
    s = "<!--"
    @test_throws FP.FranklinParserException md_blockifier(s)
end

@testset "Comment - no nesting" begin
    s = """
        <!--ABC-->
        """
    blocks = s |> md_blockifier
    @test FP.content(blocks[1]) == "ABC"
    @test isempty(blocks[1].inner_tokens)
    s = """
        <!--A<!--B-->
        """
    blocks = s |> md_blockifier
    @test blocks[1].inner_tokens[1].name == :COMMENT_OPEN
    @test FP.content(blocks[1]) == "A<!--B"
    s = """
        <!--A<!--B-->C-->
        """
    blocks = s |> md_blockifier
    @test FP.content(blocks[1]) == "A<!--B"
    @test blocks[2].name == :P_BREAK
end

@testset "Other - no nesting" begin
    blocks = "~~~ABC~~~" |> md_blockifier
    @test FP.content(blocks[1]) == "ABC"
    blocks = """
        +++
        ABC
        +++
        """ |> md_blockifier
    # code
    @test strip(FP.content(blocks[1])) == "ABC"
    blocks = "```julia ABC``` ````julia ABC```` `````julia ABC`````" |> md_blockifier
    @test FP.content(blocks[1]) // "julia ABC"
    @test FP.content(blocks[2]) // "julia ABC"
    @test FP.content(blocks[3]) // "julia ABC"
    blocks = "`A` ``A`` ``` A```" |> md_blockifier
    @test FP.content(blocks[1]) == "A"
    @test FP.content(blocks[2]) == "A"
    @test FP.content(blocks[3]) // "A"
    blocks = "```! ABC```" |> md_blockifier
    @test FP.content(blocks[1]) // "! ABC"
    # headers
    blocks = """
        # abc
        ## def
        ### ghi
        #### klm
        ##### nop
        ###### qrs
        """ |> md_blockifier
    @test strip(FP.content(blocks[1])) == "abc"
    @test strip(FP.content(blocks[2])) == "def"
    @test strip(FP.content(blocks[3])) == "ghi"
    @test strip(FP.content(blocks[4])) == "klm"
    @test strip(FP.content(blocks[5])) == "nop"
    @test strip(FP.content(blocks[6])) == "qrs"
    # footnotes
    blocks = """
        abc [^1] and [^ab] ef
        [1]:
        """
end

@testset "Braces - nesting" begin
    blocks = "{ABC}" |> md_blockifier
    @test FP.content(blocks[1]) == "ABC"
    blocks = "{ABC{DEF}H{IJ}K{L{M}O}P}" |> md_blockifier
    @test FP.content(blocks[1]) == "ABC{DEF}H{IJ}K{L{M}O}P"
end

@testset "Other - nesting" begin
    blocks = """
        @@abc DEF {GHI} KLM @@
        """ |> md_blockifier
    @test strip(FP.content(blocks[1])) == "DEF {GHI} KLM"
end

@testset "@def" begin
    blocks = """
        @def a = 5
        """ |> md_blockifier
    @test strip(FP.content(blocks[1])) == "a = 5"
    # XXX this is now two separate things and will error eventually
    blocks = """
        @def a = [1,
            2]
        """ |> md_blockifier
    @test isapproxstr(FP.content(blocks[1]), "a=[1,")
end

@testset "Double brace" begin
    blocks = """
        {{abc {g} def}}
        """ |> md_blockifier
    @test FP.content(blocks[1]) == "abc {g} def"
end

@testset "Ambiguous double braces" begin
    blocks = raw"""
        \a{\b{c}} \a{{b}\c{d}}
        """ |> md_blockifier
    @test blocks[1].ss == raw"\a"
    @test blocks[2].ss == raw"{\b{c}}"
    @test blocks[3].ss == raw"\a"
    @test blocks[4].ss == raw"{{b}\c{d}}"
end

@testset "Mix" begin
    s = """
        <!--~~~ABC~~~-->
        ~~~<!--ABC-->~~~
        """
    blocks = s |> md_blockifier
    @test FP.content(blocks[1]) == "~~~ABC~~~"
    @test FP.content(blocks[2]) == "<!--ABC-->"

    b = """
        +++
        a = "~~~567"
        b = "-->"
        +++
        """ |> md_blockifier
    @test strip(FP.content(b[1])) == "a = \"~~~567\"\nb = \"-->\""

    b = """
        ```foo
        `A` ``B``
        ```
        """ |> md_blockifier
    @test FP.content(b[1]) // "foo\n`A` ``B``"

    b = """
        @@abc
            @@def
                ```julia
                hello
                ```
            @@
            {ABC}
        @@
        """ |> md_blockifier
    @test isapproxstr(FP.content(b[1]), "@@def ```julia hello```@@ {ABC}")
end

@testset "maths" begin
    b = raw"""
        A $B$ C
        """ |> md_blockifier
    @test b[1].name == :MATH_INLINE
    @test FP.content(b[1]) == "B"
    b = raw"""
        A $$B$$ C
        """ |> md_blockifier
    @test b[1].name == :MATH_DISPL_A
    @test FP.content(b[1]) == "B"
    b = raw"""
        A \[B\] C
        """ |> md_blockifier
    @test b[1].name == :MATH_DISPL_B
    @test FP.content(b[1]) == "B"
    # b = raw"""
    #     A _$>_B_$<_ C
    #     """ |> md_blockifier
    # @test b[1].name == :MATH_I
    b = raw"""
        A $$B$$ C $D$
        """ |> md_blockifier
    @test b[1].name == :MATH_DISPL_A
    @test b[2].name == :MATH_INLINE
end

# -----------------------------
# Added to take over commonmark

@testset "emphasis" begin
    # no nest
    b = """
        a *b* c _d_ e _ f g * h
        """ |> md_blockifier
    @test b[1].name == :EMPH_EM
    @test FP.content(b[1]) == "b"
    @test b[2].name == :EMPH_EM
    @test FP.content(b[2]) == "d"
    @test b[3].name == :P_BREAK

    # nest (we only recover the outer block of course)
    b = """
        **abc_def_ghi**
        """ |> md_blockifier
    @test b[1].name == :EMPH_STRONG
    @test FP.content(b[1]) == "abc_def_ghi"
    @test b[2].name == :P_BREAK
end

@testset "cov" begin
    t = "a *bc* d _ef_ g" |> toks
    b = FP.find_blocks(t)
    @test b[1] // "*bc*"
    @test b[2] // "_ef_"
end

@testset "env" begin
    s = raw"""
        \begin{abc}
        foo
        \end{abc}
        """
    t = s |> toks
    b = FP.find_blocks(t)
    @test b[1].ss // s

    s = raw"""
        \begin{abc}
        \begin{def}
        1
        \end{def}
        2
        \end{abc}
        """
    t = s |> toks
    b = FP.find_blocks(t)
    @test b[1].ss // s

    s = raw"""
        \begin{abc}
        foo
        \end{abc}
        def
        """
    t = s |> toks
    b = FP.find_blocks(t)

    s = raw"""
        \newenvironment{foo}{bar:}{:baz}
        \begin{foo}
        abc
        \end{foo}
        ABC
        """
    t = s |> toks
    b = FP.find_blocks(t)
    @test b[5].ss // "\\begin{foo}\nabc\n\\end{foo}"
    @test b[5].inner_tokens[end].name == :CU_BRACKET_CLOSE
end

@testset "env issue" begin
    s = raw"""
        \begin{abc}
        foo
        """
    t = s |> toks
    @test_throws FP.FranklinParserException FP.find_blocks(t)
end


text:
module Streams

export Stream, closebody, isaborted,
       header, hasheader,
       setstatus, setheader

import ..HTTP
using ..Sockets
using ..IOExtras
using ..Messages
import ..Messages: header, hasheader, setheader,
                   writeheaders, writestartline
import ..ConnectionPool: getrawstream
import ..@require, ..precondition_error
import ..@ensure, ..postcondition_error
import ..@debug, ..DEBUG_LEVEL

mutable struct Stream{M <: Message, S <: IO} <: IO
    message::M
    stream::S
    writechunked::Bool
    readchunked::Bool
    warn_not_to_read_one_byte_at_a_time::Bool
    ntoread::Int
    nwritten::Int
end

"""
    Stream(::Request, ::IO)

Creates a `HTTP.Stream` that wraps an existing `IO` stream.

 - `startwrite(::Stream)` sends the `Request` headers to the `IO` stream.
 - `write(::Stream, body)` sends the `body` (or a chunk of the body).
 - `closewrite(::Stream)` sends the final `0` chunk (if needed) and calls
   `closewrite` on the `IO` stream. When the `IO` stream is a
   [`HTTP.ConnectionPool.Transaction`](@ref), calling `closewrite` releases
   the [`HTTP.ConnectionPool.Connection`](@ref) back into the pool for use by the
   next pipelined request.

 - `startread(::Stream)` calls `startread` on the `IO` stream then
    reads and parses the `Response` headers.  When the `IO` stream is a
   [`HTTP.ConnectionPool.Transaction`](@ref), calling `startread` waits for other
   pipelined responses to be read from the [`HTTP.ConnectionPool.Connection`](@ref).
 - `eof(::Stream)` and `readavailable(::Stream)` parse the body from the `IO`
    stream.
 - `closeread(::Stream)` reads the trailers and calls `closeread` on the `IO`
    stream.  When the `IO` stream is a [`HTTP.ConnectionPool.Transaction`](@ref),
    calling `closeread` releases the readlock and allows the next pipelined
    response to be read by another `Stream` that is waiting in `startread`.
    If a complete response has not been received, `closeread` throws `EOFError`.
"""
Stream(r::M, io::S) where {M, S} = Stream{M,S}(r, io, false, false, true, 0, 0)

header(http::Stream, a...) = header(http.message, a...)
setstatus(http::Stream, status) = (http.message.response.status = status)
setheader(http::Stream, a...) = setheader(http.message.response, a...)
getrawstream(http::Stream) = getrawstream(http.stream)

Sockets.getsockname(http::Stream) = Sockets.getsockname(getrawstream(http))
function Sockets.getpeername(http::Stream)
    # TODO: MbedTLS only forwards getsockname(::SSLContext)
    # so we use IOExtras.tcpsocket to reach into the MbedTLS internals
    # for now to keep compatibility with older MbedTLS versions.
    # return Sockets.getpeername(getrawstream(http))
    return Sockets.getpeername(IOExtras.tcpsocket(getrawstream(http)))
end

IOExtras.isopen(http::Stream) = isopen(http.stream)

# Writing HTTP Messages

messagetowrite(http::Stream{Response}) = http.message.request
messagetowrite(http::Stream{Request}) = http.message.response

IOExtras.iswritable(http::Stream) = iswritable(http.stream)

function IOExtras.startwrite(http::Stream)
    if !iswritable(http.stream)
        startwrite(http.stream)
    end
    m = messagetowrite(http)
    if !hasheader(m, "Content-Length") &&
       !hasheader(m, "Transfer-Encoding") &&
       !hasheader(m, "Upgrade") &&
       (m isa Request || (m.request.version >= v"1.1" && bodylength(m) > 0))

        http.writechunked = true
        setheader(m, "Transfer-Encoding" => "chunked")
    else
        http.writechunked = ischunked(m)
    end
    buf = IOBuffer()
    writeheaders(buf, m)
    n = write(http.stream, take!(buf))
    http.nwritten = 0 # should not include headers
    return n
end

function Base.unsafe_write(http::Stream, p::Ptr{UInt8}, n::UInt)
    if n == 0
        return 0
    end
    if !iswritable(http) && isopen(http.stream)
        startwrite(http)
    end
    nw = if !http.writechunked
        unsafe_write(http.stream, p, n)
    else
        write(http.stream, string(n, base=16), "\r\n") +
        unsafe_write(http.stream, p, n) +
        write(http.stream, "\r\n")
    end
    http.nwritten += nw
    return nw
end

"""
    closebody(::Stream)

Write the final `0` chunk if needed.
"""
function closebody(http::Stream)
    if http.writechunked
        http.writechunked = false
        try
            write(http.stream, "0\r\n\r\n")
        catch end
    end
end

function IOExtras.closewrite(http::Stream{Response})
    if !iswritable(http)
        return
    end
    closebody(http)
    closewrite(http.stream)
end

function IOExtras.closewrite(http::Stream{Request})

    if iswritable(http)
        closebody(http)
        closewrite(http.stream)
    end

    if hasheader(http.message, "Connection", "close") ||
       hasheader(http.message, "Connection", "upgrade") ||
       http.message.version < v"1.1" &&
      !hasheader(http.message, "Connection", "keep-alive")

        @debug 1 "✋  \"Connection: close\": $(http.stream)"
        close(http.stream)
    end
end

# Reading HTTP Messages

IOExtras.isreadable(http::Stream) = isreadable(http.stream)

Base.bytesavailable(http::Stream) = min(ntoread(http),
                                        bytesavailable(http.stream))

function IOExtras.startread(http::Stream)

    if !isreadable(http.stream)
        startread(http.stream)
    end

    readheaders(http.stream, http.message)
    handle_continue(http)

    http.readchunked = ischunked(http.message)
    http.ntoread = bodylength(http.message)

    return http.message
end

"""
100 Continue
https://tools.ietf.org/html/rfc7230#section-5.6
https://tools.ietf.org/html/rfc7231#section-6.2.1
"""
function handle_continue(http::Stream{Response})
    if http.message.status == 100
        @debug 1 "✅  Continue:   $(http.stream)"
        readheaders(http.stream, http.message)
    end
end

function handle_continue(http::Stream{Request})
    if hasheader(http.message, "Expect", "100-continue")
        if !iswritable(http.stream)
            startwrite(http.stream)
        end
        @debug 1 "✅  Continue:   $(http.stream)"
        writeheaders(http.stream, Response(100))
    end
end

function Base.eof(http::Stream)
    if !headerscomplete(http.message)
        startread(http)
    end
    if http.ntoread == 0
        return true
    end
    if eof(http.stream)
        return true
    end
    return false
end

@inline function ntoread(http::Stream)

    if !headerscomplete(http.message)
        startread(http)
    end

    # Find length of next chunk
    if http.ntoread == unknown_length && http.readchunked
        http.ntoread = readchunksize(http.stream, http.message)
    end

    return http.ntoread
end

@inline function update_ntoread(http::Stream, n)

    if http.ntoread != unknown_length
        http.ntoread -= n
    end

    if http.readchunked
        if http.ntoread == 0
            http.ntoread = unknown_length
        end
    end

    @ensure http.ntoread >= 0
end

function Base.readavailable(http::Stream, n::Int=typemax(Int))

    ntr = ntoread(http)
    if ntr == 0
        return UInt8[]
    end

    bytes = read(http.stream, min(n, ntr))
    update_ntoread(http, length(bytes))
    return bytes
end

Base.read(http::Stream, n::Integer) = readavailable(http, Int(n))

function Base.read(http::Stream, ::Type{UInt8})

    if http.warn_not_to_read_one_byte_at_a_time
        @warn "Reading one byte at a time from HTTP.Stream is inefficient.\n" *
              "Use: io = BufferedInputStream(http::HTTP.Stream) instead.\n" *
              "See: https://github.com/BioJulia/BufferedStreams.jl" stacktrace()
        http.warn_not_to_read_one_byte_at_a_time = false
    end

    if ntoread(http) == 0
        throw(EOFError())
    end
    update_ntoread(http, 1)

    return read(http.stream, UInt8)
end

function http_unsafe_read(http::Stream, p::Ptr{UInt8}, n::UInt)::Int

    ntr = UInt(ntoread(http))
    if ntr == 0
        return 0
    end

    unsafe_read(http.stream, p, min(n, ntr + (http.readchunked ? 2 : 0)))
                                             # If there is spare space in `p`
                                             # read two extra bytes
    n = min(n, ntr)                          # (`\r\n` at end ofchunk).
    update_ntoread(http, n)
    return n
end

function Base.readbytes!(http::Stream, buf::AbstractVector{UInt8},
                                       n=length(buf))
    @require n <= length(buf)
    return http_unsafe_read(http, pointer(buf), UInt(n))
end

function Base.unsafe_read(http::Stream, p::Ptr{UInt8}, n::UInt)
    nread = 0
    while nread < n
        if eof(http)
            throw(EOFError())
        end
        nread += http_unsafe_read(http, p + nread, n - nread)
    end
    nothing
end

Base.readbytes!(http::Stream, buf::Base.BufferStream, n=bytesavailable(http)) = Base.readbytes!(http, buf.buffer, n)

function Base.readbytes!(http::Stream, buf::IOBuffer, n=bytesavailable(http))
    Base.ensureroom(buf, n)
    unsafe_read(http, pointer(buf.data, buf.size + 1), n)
    buf.size += n
end

# used in julia v1.0
function Base.readbytes!(http::Stream, buf::IOStream, n=bytesavailable(http))
    nread = 0
    while nread < n
        nread += write(buf, readavailable(http, n - nread))
    end
    nread
end

function Base.read(http::Stream)
    buf = PipeBuffer()
    if ntoread(http) == unknown_length
        while !eof(http)
            readbytes!(http, buf)
        end
    else
        while !eof(http)
            readbytes!(http, buf, ntoread(http))
        end
    end
    return take!(buf)
end

"""
    isaborted(::Stream{Response})

Has the server signaled that it does not wish to receive the message body?

"If [the response] indicates the server does not wish to receive the
 message body and is closing the connection, the client SHOULD
 immediately cease transmitting the body and close the connection."
[RFC7230, 6.5](https://tools.ietf.org/html/rfc7230#section-6.5)
"""
function isaborted(http::Stream{Response})

    if iswritable(http.stream) &&
       iserror(http.message) &&
       hasheader(http.message, "Connection", "close")
        @debug 1 "✋  Abort on $(sprint(writestartline, http.message)): " *
                 "$(http.stream)"
        @debug 2 "✋  $(http.message)"
        return true
    end
    return false
end

incomplete(http::Stream) =
    http.ntoread > 0 && (http.readchunked || http.ntoread != unknown_length)

function IOExtras.closeread(http::Stream{Response})

    if hasheader(http.message, "Connection", "close")
        # Close conncetion if server sent "Connection: close"...
        @debug 1 "✋  \"Connection: close\": $(http.stream)"
        close(http.stream)
    else

        # Discard body bytes that were not read...
        while !eof(http)
            readavailable(http)
        end

        if incomplete(http)
            # Error if Message is not complete...
            close(http.stream)
            throw(EOFError())
        elseif isreadable(http.stream)
            closeread(http.stream)
        end
    end

    return http.message
end

function IOExtras.closeread(http::Stream{Request})
    if incomplete(http)
        # Error if Message is not complete...
        close(http.stream)
        throw(EOFError())
    end
    if isreadable(http)
        closeread(http.stream)
    end
end

end #module Streams


text:
# Neighborhood Strength Driven Label Propagation Algorithm
function nsdlpa{V,T<:Real}(g::AbstractGraph{V}, λ::T=1.0, CN::Vector{Int}=triangle(g))
    @graph_requires g edge_map incidence_list vertex_map

    # calculate number of common neighbors
    N = num_vertices(g)
    membership = [1:N]
    label_counters = zeros(T, N)
    dominant_labels = Int[]
    nonzero_labels = Int[]

    running = true
    STEP = 0
    while running
        STEP += 1
        running = false

        # Shuffle the node ordering vector
        X = shuffle(collect(vertices(g)))

        # In the prescribed order, loop over the vertices and reassign labels
        for v in X
            v_idx = vertex_index(v, g)

            # Clear dominant_labels and nonzero_labels
            dominant_labels = Int[]
            nonzero_labels = Int[]

            # recount
            max_count = zero(T)

            for w in out_edges(v, g)
                k = membership[vertex_index(target(w, g), g)]

                # skip if it has no label yet
                k != 0 || continue

                was_zero = label_counters[k] == zero(T)
                label_counters[k] += 1 + λ*CN[edge_index(w, g)]
                if was_zero && label_counters[k] != zero(T)
                    # counter just became nonzero
                    push!(nonzero_labels, k)
                end
                if max_count < label_counters[k]
                    max_count = label_counters[k]
                    resize!(dominant_labels, 1)
                    dominant_labels[1] = k
                elseif max_count == label_counters[k]
                    push!(dominant_labels, k)
                end
            end

            if !isempty(dominant_labels)

                # Select randomly from the dominant labels
                k = dominant_labels[rand(1:length(dominant_labels))]

                # Check if the current label of the node is also dominant
                if label_counters[membership[v_idx]] != max_count
                    # Nope, we need at least one more iteration
                    running = true
                end

                # Update label of the current node
                membership[v_idx] = k
            end
            # Clear the nonzero elements in label_counters
            for i in nonzero_labels
                label_counters[i] = zero(T)
            end
        end
    end
    permute_labels!(membership), STEP
end

function nsdlpa1{V,T<:Real}(g::AbstractGraph{V}, λ::T=1.0, CN::Vector{Int}=triangle(g), DS::Vector{Int}=quadrangle(g))
    @graph_requires g edge_map incidence_list vertex_map

    # calculate number of common neighbors
    N = num_vertices(g)
    membership = [1:N]
    label_counters = zeros(T, N)
    proximity_sum = zeros(T, N)

    running = true
    STEP = 0
    while running
        STEP += 1
        running = false

        # Shuffle the node ordering vector
        X = shuffle(collect(vertices(g)))

        # In the prescribed order, loop over the vertices and reassign labels
        for v in X
            v_idx = vertex_index(v, g)

            # Clear dominant_labels and nonzero_labels
            dominant_labels = Int[]
            nonzero_labels = Int[]
            similarity = Float64[]

            # recount
            max_count = zero(T)

            for w in out_edges(v, g)
                k = membership[vertex_index(target(w, g), g)]

                # skip if it has no label yet
                k != 0 || continue

                was_zero = label_counters[k] == zero(T)
                label_counters[k] += 1 + λ*CN[edge_index(w, g)]
                proximity_sum[k] += DS[edge_index(w, g)]
                if was_zero && label_counters[k] != zero(T)
                    # counter just became nonzero
                    push!(nonzero_labels, k)
                end
                if max_count < label_counters[k]
                    max_count = label_counters[k]
                    resize!(dominant_labels, 1)
                    resize!(similarity, 1)
                    dominant_labels[1] = k
                    similarity[1] = proximity_sum[k]
                elseif max_count == label_counters[k]
                    push!(dominant_labels, k)
                    push!(similarity, proximity_sum[k])
                end
            end

            if !isempty(dominant_labels)

                max_similarity = maximum(similarity)
                dominant_labels = dominant_labels[similarity.<=max_similarity]

                # Select randomly from the dominant labels
                k = dominant_labels[rand(1:length(dominant_labels))]

                # Check if the current label of the node is also dominant
                if label_counters[membership[v_idx]] != max_count
                    # Nope, we need at least one more iteration
                    running = true
                end

                # Update label of the current node
                membership[v_idx] = k
            end
            # Clear the nonzero elements in label_counters
            for i in nonzero_labels
                label_counters[i] = zero(T)
                proximity_sum[i] = zero(T)
            end
        end
    end
    permute_labels!(membership), STEP
end

max_stars_repo_name: afternone/LPA

max_issues_repo_name: afternone/LPA

max_forks_repo_name: afternone/LPA


text:
using
  Documenter,
  FourierFlows

makedocs(
   modules = [FourierFlows],
   doctest = false, clean = true,
 checkdocs = :all,
    format = :html,
   authors = "Gregory L. Wagner and Navid C. Constantinou",
  sitename = "FourierFlows.jl",
     pages = Any[
              "Home" => "index.md",
              "Code Basics" => "basics.md",
              "Forcing" => "forcing.md",
              "Modules" => Any[
                "modules/kuramotosivashinsky.md",
                "modules/traceradvdiff.md",
              ],
              "DocStrings" => Any[
              "man/types.md",
              "man/functions.md"]
             ]
)

deploydocs(
       repo = "github.com/FourierFlows/FourierFlows.jl.git",
     target = "build",
      julia = "0.7",
     osname = "linux",
       deps = nothing,
       make = nothing
 )


text:
module Denoising
export
    # denoising
    RelErrorShrink,
    SureShrink,
    denoiseall,
    relerrorthreshold

using
    Wavelets,
    LinearAlgebra,
    Statistics,
    Plots

using
    ..DWT,
    ..SWT,
    ..ACWT,
    ..Utils

# ========== Threshold Determination Methods ==========
"""
    RelErrorShrink(th, t) <: DNFT

Relative Error Shrink method used in their paper "Efficient Approximation and Denoising of
Graph Signals using the Multiscale Basis Dictionary" for IEEE Transactions on Signal and
Information Processing over Networks, Vol 0, No. 0, 2016.

# Attributes
- `th::Wavelets.Threshold.THType`: (Default: `Wavelets.HardTH()`) Threshold type.
- `t::AbstractFloat`: (Default: 1.0) Threshold size.

# Examples
```julia
using Wavelets, WaveletsExt

RelErrorShrink()                    # Using default th and t values
RelErrorShrink(SoftTH())            # Using default t value
RelErrorShrink(HardTH(), 0.5)       # Using user input th and t values
```

**See also:** [`SureShrink`](@ref), [`VisuShrink`](@ref)
"""
struct RelErrorShrink <: DNFT
    # Attributes
    th::Wavelets.Threshold.THType
    t::AbstractFloat
    # Constructor
    RelErrorShrink(th = HardTH(), t = 1.0) = new(th, t)
end

"""
    SureShrink(th, t) <: DNFT

Stein's Unbiased Risk Estimate (SURE) Shrink

# Attributes
- `th::Wavelets.Threshold.THType`: (Default: `Wavelets.HardTH()`) Threshold type.
- `t::AbstractFloat`: (Default: 1.0) Threshold size.

**See also:** [`RelErrorShrink`](@ref), [`VisuShrink`](@ref)
"""
struct SureShrink <: DNFT
    # Attributes
    th::Wavelets.Threshold.THType
    t::AbstractFloat
    # Constructor
    SureShrink(th, t) = new(th, t)
end

"""
    SureShrink(xw[, redundant, tree, th])

Struct constructor for `SureShrink` based on the signal coefficients `xw`.

# Arguments
- `xw::AbstractArray{T} where T<:Number`: Decomposed signal.
- `redundant::Bool`: (Default: `false`) Whether the transform type of `xw` is a redundant
  transform. Autocorrelation and stationary wavelet transforms are examples of redundant
  transforms.
- `tree::Union{BitVector, Nothing}`: (Default: `nothing`) The basis tree for decomposing
  `xw`. Must be provided if `xw` is decomposed using `wpt`, `swpd`, or `acwpd`.
- `th::Wavelets.Threshold.THType`: (Default: `HardTH()`) Threshold type.

# Returns
`::SureShrink`: SUREShrink object.

# Examples
```julia
SureShrink(xw)                  # `xw` is output of dwt, wpt
SureShrink(xw, true)            # `xw` is output of sdwt, acdwt, swpt, acwpt
SureShrink(xw, true, tree)      # `xw` is output of swpd, acwpd
```

**See also:** [`SureShrink`](@ref), [`surethreshold`](@ref)
"""
function SureShrink(xw::AbstractArray{T},
                    redundant::Bool = false,
                    tree::Union{BitVector, Nothing} = nothing,
                    th::Threshold.THType = HardTH()) where T<:Number
    t = surethreshold(xw, redundant, tree)
    return SureShrink(th, t)
end

"""
    VisuShrink(n, th)

Extension to the `VisuShrink` struct constructor from `Wavelets.jl`.

# Arguments
- `n::Integer`: Signal length.
- `th::Wavelets.Threshold.THType`: Threshold type.

# Returns
`::Wavelets.Threshold.VisuShrink`: VisuShrink object.

# Examples
```julia
using Wavelets, WaveletsExt

VisuShrink(128, SoftTH())
```
"""
function Wavelets.Threshold.VisuShrink(n::Integer, th::Wavelets.Threshold.THType)
    return VisuShrink(th, sqrt(2*log(n)))
end

"""
    surethreshold(coef, redundant[, tree])

Determination of the `t` value used for `SureShrink`. `t` is defined as the threshold value
when the standard deviation of the noisy signal is 1.

# Arguments
- `coef::AbstractArray{T} where T<:Number`: Coefficients of decomposed signal.
- `redundant::Bool`: Whether the transform type of `xw` is a redundant transform.
  Autocorrelation and stationary wavelet transforms are examples of redundant transforms.
- `tree::Union{BitVector, Nothing}`: (Default: `nothing`) The basis tree for decomposing
  `xw`. Must be provided if `xw` is decomposed using `wpt`, `swpd`, or `acwpd`.

# Returns
`::AbstractFloat`: `t` value used for `SureShrink`.

**See also:** [`SureShrink`](@ref)
"""
function surethreshold(coef::AbstractArray{T},
                       redundant::Bool,
                       tree::Union{BitVector,Nothing} = nothing) where T<:Number
    # extract necessary coefficients
    if !redundant                              # dwt or wpt
        y = coef
    elseif redundant && isa(tree, Nothing)     # sdwt, acdwt, swpt, or acwpt
        y = reshape(coef, :)
    else                                       # swpd or acwpd
        leaves = getleaf(tree,:binary)
        y = reshape(coef[:, leaves], :)
    end
    a = sort(abs.(y)).^2
    b = cumsum(a)
    n = length(y)
    c = collect(reverse(0:(n-1)))
    s = b + c .* a
    risk = (n .- (2*(1:n)) + s)/n
    i = argmin(risk)
    return sqrt(a[i])
end

# ========== Noise Estimation ==========
"""
    noisest(x, redundant[, tree])

Extension to the `noisest` function from `Wavelets.jl`. Estimates the standard deviation of
a signal's noise assuming that the noise is distributed normally. This function is generally
used in combination with `VisuShrink` and `SureShrink` in the `denoise`/`denoiseall`
functions.

# Arguments
- `x::AbstractArray{T}`: Decomposed signal.
- `redundant::Bool`: Whether the transform type of `xw` is a redundant transform.
  Autocorrelation and stationary wavelet transforms are examples of redundant transforms.
- `tree::Union{BitVector, Nothing}`: (Default: `nothing`) The basis tree for decomposing
  `xw`. Must be provided if `xw` is decomposed using `wpt`, `swpd`, or `acwpd`.

# Returns
`::AbstractFloat`: Estimated standard deviation of the noise of the signal.

# Examples
```julia
using Wavelets, WaveletsExt

x = randn(128)
wt = wavelet(WT.haar)

# noise estimate for dwt transformation
y = dwt(x, wt)
noise = noisest(y, false)

# noise estimate for wpt transformation
tree = maketree(x, :full)
y = wpt(x, wt, tree)
noise = noisest(y, false, tree)

# noise estimate for sdwt transformation
y = sdwt(x, wt)
noise = noisest(y, true)

# noise estimate for swpd transformation
y = swpd(x, wt)
noise = noisest(y, true, tree)
```

**See also:** [`relerrorthreshold`](@ref), [`VisuShrink`](@ref), [`SureShrink`](@ref)
"""
function Wavelets.Threshold.noisest(x::AbstractArray{T},
                                    redundant::Bool,
                                    tree::Union{BitVector,Nothing} = nothing) where T<:Number
    # Sanity check
    @assert isdyadic(size(x,1))

    # Get detail coefficients
    if !redundant && isa(tree, Nothing)        # regular dwt
        n₀ = size(x,1) ÷ 2
        dr = x[(n₀+1):end]
    elseif !redundant && isa(tree, BitVector)  # regular wpt
        dr = x[finestdetailrange(x, tree)]
    elseif redundant && isa(tree, Nothing)     # redundant dwt
        dr = x[:,end]
    else                                       # redundant wpd
        dr = x[finestdetailrange(x, tree, true)...]
    end
    return Wavelets.Threshold.mad!(dr)/0.6745
end

"""
    relerrorthreshold(coef, [redundant, tree, elbows; makeplot])

Takes in a set of expansion coefficients, 'plot' the threshold vs relative error curve and
select the best threshold value based on the elbow method. If one wants to see the resulting
plot from this computation, simply set `makeplot=true`.

# Arguments
- `coef::AbstractArray{T} where T<:Number`: Decomposed signal.
- `redundant::Bool`: (Default: `false`) Whether the transform type of `xw` is a redundant
  transform. Autocorrelation and stationary wavelet transforms are examples of redundant
  transforms.
- `tree::Union{BitVector, Nothing}`: (Default: `nothing`) The basis tree for decomposing
  `xw`. Must be provided if `xw` is decomposed using `swpd` or `acwpd`.
- `elbows::Integer`: (Default: 2) Number of elbows used to determine the best threshold
  value.

# Keyword Arguments
- `makeplot::Bool`: (Default: `false`) Whether to return the plot that was used to determine
  the best threshold value.

# Returns
- `::AbstractFloat`: Best threshold value.
- `::GR.Plot`: Plot that was used to determine the best threshold value. Only returned if
  `makeplot = true`.

# Examples
```julia
x = randn(128)
wt = wavelet(WT.haar)

# noise estimate for dwt transformation
y = dwt(x, wt)
noise = relerrorthreshold(y, false)

# noise estimate for wpt transformation
tree = maketree(x, :full)
y = wpt(x, wt, tree)
noise = relerrorthreshold(y, false, tree)

# noise estimate for sdwt transformation
y = sdwt(x, wt)
noise = relerrorthreshold(y, true)

# noise estimate for swpd transformation
y = swpd(x, wt)
noise = relerrorthreshold(y, true, tree)
```

**See also:** [`noisest`](@ref), [`RelErrorShrink`](@ref)
"""
function relerrorthreshold(coef::AbstractArray{T},
                           redundant::Bool = false,
                           tree::Union{BitVector,Nothing} = nothing,
                           elbows::Integer = 2;
                           makeplot::Bool = false) where T<:Number
    # Sanity check
    @assert elbows >= 1
    # extract necessary coefficients
    if !redundant                              # dwt or wpt
        c = coef
    elseif redundant && isa(tree, Nothing)     # sdwt
        c = reshape(coef, :)
    else                                       # swpd
        leaves = getleaf(tree,:binary)
        c = reshape(coef[:, leaves], :)
    end
    # the magnitudes of the coefficients
    x = sort(abs.(c), rev = true)
    # compute the relative error curve
    r = orth2relerror(c)
    # shift the data points
    push!(x, 0)
    pushfirst!(r, r[1])
    # reorder the data points
    xmax = maximum(x)
    ymax = maximum(r)
    x = x[end:-1:1]/xmax
    y = r[end:-1:1]/ymax
    # compute elbow method
    ix = Vector{Integer}(undef, elbows)
    A = Vector{Vector{T}}(undef, elbows)
    v = Vector{Vector{T}}(undef, elbows)
    ix[1], A[1], v[1] = findelbow(x,y)          # First elbow point
    for i in 2:elbows                           # Second elbow point and beyond
        @inbounds (ix[i], A[i], v[i]) = findelbow(x[1:ix[i-1]], y[1:ix[i-1]])
    end
    # plot relative error curve
    if makeplot
        p = relerrorplot(x*xmax, y*ymax, ix, A, v)
        return x[ix[end]]*xmax, p
    else
        return x[ix[end]]*xmax
    end
end

"""
    orth2relerror(orth)

Given a vector 'orth' of orthonormal expansion coefficients, return a vector of relative
approximation errors when retaining the 1,2,...,N largest coefficients in magnitude.

# Arguments
- `orth::AbstractVector{T} where T<:Number`: Vector of coefficients.

# Returns
`::Vector{T}`: Relative errors.

**See also:** [`RelErrorShrink`](@ref), [`relerrorthreshold`](@ref), [`findelbow`](@ref)
"""
function orth2relerror(orth::AbstractVector{T}) where T<:Number
    # sort the coefficients
    orth = sort(orth.^2, rev = true)
    # compute the relative errors
    return ((abs.(sum(orth) .- cumsum(orth))).^0.5) / sum(orth).^0.5
end

"""
    findelbow(x, y)

Given the x and y coordinates of a curve, return the elbow.

# Arguments
- `x::AbstractVector{T} where T<:Number`: x-coordinates.
- `y::AbstractVector{T} where T<:Number`: y-coordinates.

# Returns
- `::Integer`: Index of the elbow point.
- `::Vector{T}`: Length of adjacent sides.
- `::Vector{T}`: The y-coordinates going in the direction of (x₁, y₁) to (xₙ, yₙ)

**See also:** [`RelErrorShrink`](@ref), [`relerrorthreshold`](@ref), [`orth2relerror`](@ref)
"""
function findelbow(x::AbstractVector{T}, y::AbstractVector{T}) where T<:Number
    # a unit vector pointing from (x1,y1) to (xN,yN)
    v = [x[end] - x[1], y[end] - y[1]]
    v = v/norm(v,2)
    # subtract (x1,y1) from the coordinates
    xy = [x.-x[1] y.-y[1]]
    # the hypothenuse
    H = reshape((sum(xy.^2, dims = 2)).^0.5, :)
    # the adjacent side
    A = xy * v
    # the opposite side
    O = abs.(H.^2 - A.^2).^0.5
    # return the largest distance
    return (findmax(O)[2], A, v)
end

"""
    relerrorplot(x, y, ix, A, v)

Relative error plot used for threshold determination using the elbow rule.

# Arguments
- `x::Vector{T} where T<:Number`: x-coordinates.
- `y::Vector{T} where T<:Number`: y-coordinates.
- `ix::Vector{<:Integer}`: Indices for elbow points.
- `A::Vector{S} where {S<:Vector{T}, T<:Number}`: Length of adjacent sides.
- `v::Vector{S} where {S<:Vector{T}, T<:Number}`: The y-coordinates going in the direction
  of (x[1], y[1]) to (x[ix], y[ix])

# Returns
`::Plot`: Relative error plot.

**See also:** [`relerrorthreshold`](@ref), [`findelbow`](@ref)
"""
function relerrorplot(x::Vector{T},
                      y::Vector{T},
                      ix::Vector{<:Integer},
                      A::Vector{S},
                      v::Vector{S}) where {T<:Number, S<:Vector{T}}
    @assert length(ix) == length(A) == length(v)
    elbows = length(ix)
    # rescale x and y values
    xmax = maximum(x)
    ymax = maximum(y)
    # relative error line
    p = plot(x, y, lw = 2, color = :blue, legend = false)
    plot!(p, xlims = (0, 1.004*xmax), ylims = (0, 1.004*ymax))
    for i in 1:elbows
        col = i + 1
        # diagonal line
        endpoint = i > 1 ? ix[i-1] : length(x)
        plot!([x[1], 1.004*x[endpoint]], [y[1], 1.004*y[endpoint]], lw = 2,
            color = col)
        # perpendicular line
        dropto = [x[1], y[1]] + A[i][ix[i]]*(v[i].*[xmax, ymax])
        plot!(p, [x[ix[i]], dropto[1]], [y[ix[i]], dropto[2]], lw = 2,
            color = col)
        # highlight point
        scatter!(p, [x[ix[i]]], [y[ix[i]]], color = col)
    end
    # add plot labels
    plot!(p, xlabel = "Threshold", ylabel = "Relative Error")
    return p
end

# TODO: add additional threshold determination methods

# ========== Denoising ==========
"""
    denoise(x, inputtype, wt[; L=maxtransformlevels(size(x,1)),
        tree=maketree(size(x,1), L, :dwt), dnt=VisuShrink(size(x,1)),
        estnoise=noisest, smooth=:regular])

Extension of the `denoise` function from `Wavelets.jl`. Denoise a signal of
input type `inputtype`.

# Arguments:
- `x::AbstractArray{<:Number}`: input signals/coefficients.
- `inputtype::Symbol`: input type of `x`. Current accepted types of inputs are
    - `:sig`: original signals; `x` should be a 2-D array with each column
        representing a signal.
    - `:dwt`: `dwt`-transformed signal coefficients; `x` should be a 1-D array
        with each column representing the coefficients of a signal.
    - `:wpt`: `wpt`-transformed signal coefficients; `x` should be a 1-D array
        with each column representing the coefficients of a signal.
    - `:sdwt`: `sdwt`-transformed signal coefficients; `x` should be a 2-D array
        with each column representing the coefficients of a node.
    - `:swpd`: `swpd`-transformed signal coefficients; `x` should be a 2-D array
        with each column representing the coefficients of a node.
    - `:acwt`: `acwt`-transformed signal coefficients from
        AutocorrelationShell.jl; `x` should be a 2-D array with each column
        representing the coefficients of a node.
    - `:acwpt`: `acwpt`-transformed signal coefficients from
        AutocorrelationShell.jl; `x` should be a 2-D array with each column
        representing the coefficients of a node.
- `wt::Union{DiscreteWavelet, Nothing}`: the discrete wavelet to be used for
    decomposition (for input type `:sig`) and reconstruction. `nothing` can
    be supplied if no reconstruction is necessary.
- `L::Integer`: the number of decomposition levels. Necessary for input types
    `:sig`, `:dwt`, and `:sdwt`. Default value is set to be
    `maxtransformlevels(size(x,1))`.
- `tree::BitVector`: the decomposition tree of the signals. Necessary for input
    types `:wpt` and `:swpd`. Default value is set to be
    `maketree(size(x,1), L, :dwt)`.
- `dnt::DNFT`: denoise type. Default type is set to be `VisuShrink(size(x,1))`.
- `estnoise::Union{Function, Vector{<:Number}}`: noise estimation. Input can be
    provided as a function to estimate noise in signal, or a vector of estimated
    noise. Default is set to be the `noisest` function.
- `smooth::Symbol`: the smoothing method used. `:regular` smoothing thresholds
    all given coefficients, whereas `:undersmooth` smoothing does not threshold
    the lowest frequency subspace node of the wavelet transform. Default is set
    to be `:regular`.

**See also:** [`denoiseall`](@ref), [`noisest`](@ref),
    [`relerrorthreshold`](@ref)
"""
function Wavelets.Threshold.denoise(x::AbstractArray{T},
        inputtype::Symbol,
        wt::Union{DiscreteWavelet,Nothing};
        L::Integer=maxtransformlevels(size(x,1)),
        tree::BitVector=maketree(size(x,1), L, :dwt),
        dnt::S=VisuShrink(size(x,1)),
        estnoise::Union{Function,Number}=noisest,   # can be precomputed
        smooth::Symbol=:regular) where {T<:Number, S<:DNFT}

    @assert smooth ∈ [:undersmooth, :regular]
    @assert inputtype ∈ [:sig, :dwt, :wpt, :sdwt, :swpd, :acwt, :acwpt]

    # wavelet transform if inputtype == :sig
    if inputtype == :sig
        wt === nothing && error("inputtype=:sig not supported with wt=nothing")
        x = dwt(x, wt, L)
        inputtype = :dwt
    end

    if inputtype == :dwt                    # regular dwt
        # noise estimation
        σ = isa(estnoise, Function) ? estnoise(x, false, nothing) : estnoise
        # thresholding
        if smooth == :regular
            x̃ = threshold(x, dnt.th, σ*dnt.t)
        else    # :undersmooth
            n₀ = nodelength(size(x,1), L)
            x̃ = [x[1:n₀]; threshold(x[(n₀+1):end], dnt.th, σ*dnt.t)]
        end
        # reconstruction
        y = wt === nothing ? x̃ : idwt(x̃, wt, L)

    elseif inputtype == :wpt                # regular wpt
        # noise estimation
        σ = isa(estnoise, Function) ? estnoise(x, false, tree) : estnoise
        # thresholding
        if smooth == :regular
            x̃ = threshold(x, dnt.th, σ*dnt.t)
        else    # :undersmooth
            crng = coarsestscalingrange(x, tree, false)
            rng = setdiff(1:size(x,1), crng)
            x̃ = [x[crng]; threshold(x[rng], dnt.th, σ*dnt.t)]
        end
        # reconstruction
        y = wt === nothing ? x̃ : iwpt(x̃, wt, tree)

    elseif inputtype == :sdwt               # stationary dwt
        @assert ndims(x) > 1
        # noise estimation
        σ = isa(estnoise, Function) ? estnoise(x, true, nothing) : estnoise
        # thresholding
        if smooth == :regular
            x̃ = copy(x)
            threshold!(x̃, dnt.th, σ*dnt.t)
        else    # :undersmooth
            temp = x[:,2:end]
            x̃ = [x[:,1] threshold!(temp, dnt.th, σ*dnt.t)]
        end
        # reconstruction
        y = wt === nothing ? x̃ : isdwt(x̃, wt)

    elseif inputtype == :swpd               # stationary wpd
        @assert ndims(x) > 1
        # noise estimation
        σ = isa(estnoise, Function) ? estnoise(x, true, tree) : estnoise
        # thresholding
        if smooth == :regular
            leaves = findall(getleaf(tree,:binary))
            x̃ = copy(x)
            @inbounds x̃[:, leaves] = threshold!(x[:,leaves], dnt.th, σ*dnt.t)
        else    # :undersmooth
            x̃ = copy(x)
            leaves = findall(getleaf(tree,:binary))
            _, coarsestnode = coarsestscalingrange(x, tree, true)
            rng = setdiff(leaves, coarsestnode)
            @inbounds x̃[:,rng] = threshold!(x[:,rng], dnt.th, σ*dnt.t)
        end
        # reconstruction
        y = wt === nothing ? x̃ : iswpd(x̃, wt, tree)
    elseif inputtype == :acwt               # autocorrelation dwt
        @assert ndims(x) > 1
        # noise estimation
        σ = isa(estnoise, Function) ? estnoise(x, true, nothing) : estnoise
        # thresholding
        if smooth == :regular
            x̃ = copy(x)
            threshold!(x̃, dnt.th, σ*dnt.t)
        else    # :undersmooth
            temp = x[:,2:end]
            x̃ = [x[:,1] threshold!(temp, dnt.th, σ*dnt.t)]
        end
        # reconstruction
        y = iacdwt(x̃)

    else                                    # autocorrelation wpt
        @assert ndims(x) > 1
        # noise estimation
        σ = isa(estnoise, Function) ? estnoise(x, true, tree) : estnoise
        # thresholding
        if smooth == :regular
            leaves = findall(getleaf(tree,:binary))
            x̃ = copy(x)
            @inbounds x̃[:, leaves] = threshold!(x[:,leaves], dnt.th, σ*dnt.t)
        else    # :undersmooth
            x̃ = copy(x)
            leaves = findall(getleaf(tree,:binary))
            _, coarsestnode = coarsestscalingrange(x, tree, true)
            rng = setdiff(leaves, coarsestnode)
            @inbounds x̃[:,rng] = threshold!(x[:,rng], dnt.th, σ*dnt.t)
        end
        # reconstruction
        y = iacwpd(x̃, tree)
    end
    return y
end

"""
    denoiseall(x, inputtype, wt[; L=maxtransformlevels(size(x,1)),
        tree=maketree(size(x,1), L, :dwt), dnt=VisuShrink(size(x,1)),
        estnoise=noisest, bestTH=nothing, smooth=:regular])

Denoise multiple signals of input type `inputtype`.

# Arguments:
- `x::AbstractArray{<:Number}`: input signals/coefficients.
- `inputtype::Symbol`: input type of `x`. Current accepted types of inputs are
    - `:sig`: original signals; `x` should be a 2-D array with each column
        representing a signal.
    - `:dwt`: `dwt`-transformed signal coefficients; `x` should be a 2-D array
        with each column representing the coefficients of a signal.
    - `:wpt`: `wpt`-transformed signal coefficients; `x` should be a 2-D array
        with each column representing the coefficients of a signal.
    - `:sdwt`: `sdwt`-transformed signal coefficients; `x` should be a 3-D array
        with each 2-D slice representing the coefficients of a signal.
    - `:swpd`: `swpd`-transformed signal coefficients; `x` should be a 3-D array
        with each 2-D slice representing the coefficients of a signal.
    - `:acwt`: `acwt`-transformed signal coefficients from
        AutocorrelationShell.jl; `x` should be a 3-D array with each 2-D slice
        representing the coefficients of a signal.
    - `:acwpt`: `acwpt`-transformed signal coefficients from
        AutocorrelationShell.jl; `x` should be a 3-D array with each 2-D slice
        representing the coefficients of a signal.
- `wt::Union{DiscreteWavelet, Nothing}`: the discrete wavelet to be used for
    decomposition (for input type `:sig`) and reconstruction. `nothing` can
    be supplied if no reconstruction is necessary.
- `L::Integer`: the number of decomposition levels. Necessary for input types
    `:sig`, `:dwt`, and `:sdwt`. Default value is set to be
    `maxtransformlevels(size(x,1))`.
- `tree::BitVector`: the decomposition tree of the signals. Necessary for input
    types `:wpt` and `:swpd`. Default value is set to be
    `maketree(size(x,1), L, :dwt)`.
- `dnt::DNFT`: denoise type. Default type is set to be `VisuShrink(size(x,1))`.
- `estnoise::Union{Function, Vector{<:Number}}`: noise estimation. Input can be
    provided as a function to estimate noise in signal, or a vector of estimated
    noise. Default is set to be the `noisest` function.
- `bestTH::Union{Function, Nothing}`: method to determine the best threshold
    value for a group of signals. If `nothing` is given, then each signal will
    be denoised by its respective best threshold value determined from the
    parameters `dnt` and `estnoise`; otherwise some function can be passed
    to determine the best threshold value from a vector of threshold values, eg:
    `mean` and `median`. Default is set to be `nothing`.
- `smooth::Symbol`: the smoothing method used. `:regular` smoothing thresholds
    all given coefficients, whereas `:undersmooth` smoothing does not threshold
    the lowest frequency subspace node of the wavelet transform. Default is set
    to be `:regular`.

**See alse:** [`denoise`](@ref), [`noisest`](@ref), [`relerrorthreshold`](@ref)
"""
function denoiseall(x::AbstractArray{T1},
        inputtype::Symbol,
        wt::Union{DiscreteWavelet,Nothing};
        L::Integer=maxtransformlevels(size(x,1)),
        tree::BitVector=maketree(size(x,1), L, :dwt),
        dnt::S=VisuShrink(size(x,1)),
        estnoise::Union{Function,Vector{T2}}=noisest,   # can be precomputed
        bestTH::Union{Function,Nothing}=nothing,
        smooth::Symbol=:regular) where {T1<:Number, T2<:Number, S<:DNFT}

    @assert inputtype ∈ [:sig, :dwt, :wpt, :sdwt, :swpd, :acwt, :acwpt]
    @assert smooth ∈ [:regular, :undersmooth]
    @assert ndims(x) > 1
    n = size(x, 1)              # signal length
    N = size(x)[end]            # number of signals
    # wavelet transform if inputtype == :sig
    if inputtype == :sig
        wt === nothing && error("inputtype=:sig not supported with wt=nothing")
        x = hcat([dwt(x[:,i], wt, L) for i in 1:N]...)
        inputtype = :dwt
    end

    # denoise each signal
    y = Array{T1}(undef, (n,N))
    if bestTH === nothing           # using individual threshold values
        for i in axes(x, ndims(x))
            xᵢ = ndims(x) == 2 ? x[:,i] : x[:,:,i]
            # noise estimation
            σ = isa(estnoise, Function) ? estnoise : estnoise[i]
            # denoising
            @inbounds y[:,i] = denoise(xᵢ, inputtype, wt, L=L, tree=tree,
                dnt=dnt, estnoise=σ, smooth=smooth)
        end
    else                            # using summary threshold value
        # noise estimation
        if isa(estnoise, Function)
            σ = Vector{AbstractFloat}(undef, N)
            for i in axes(x, ndims(x))
                xᵢ = ndims(x) == 2 ? x[:,i] : x[:,:,i]
                if inputtype == :dwt
                    @inbounds σ[i] = estnoise(xᵢ, false, nothing)
                elseif inputtype == :wpt
                    @inbounds σ[i] = estnoise(xᵢ, false, tree)
                elseif inputtype == :sdwt
                    @inbounds σ[i] = estnoise(xᵢ, true, nothing)
                else
                    @inbounds σ[i] = estnoise(xᵢ, true, tree)
                end
            end
            σ = bestTH(σ)
        else
            σ = bestTH(estnoise)
        end
        # denoising
        for i in axes(x, ndims(x))
            xᵢ = ndims(x) == 2 ? x[:,i] : x[:,:,i]
            @inbounds y[:,i] = denoise(xᵢ, inputtype, wt, L=L, tree=tree,
                dnt=dnt, estnoise=σ, smooth=smooth)
        end
    end
    return y
end

end # end module


text:
include("double_integrator.jl")

## Discrete Double Integrator model
struct DiscreteModel{NM} <: RobotDynamics.DiscreteDynamics
    cfg::ForwardDiff.JacobianConfig{
        Nothing,
        Float64,
        NM,
        Tuple{
            Vector{ForwardDiff.Dual{Nothing,Float64,NM}},
            Vector{ForwardDiff.Dual{Nothing,Float64,NM}},
        },
    }
    cache::FiniteDiff.JacobianCache{
        Vector{Float64},
        Vector{Float64},
        Vector{Float64},
        UnitRange{Int64},
        Nothing,
        Val{:forward}(),
        Float64,
    }
    function DiscreteModel()
        n,m = 4,2
        y = zeros(n)
        z = zeros(n + m)
        cfg = ForwardDiff.JacobianConfig(nothing, y, z)
        cache = FiniteDiff.JacobianCache(z, y)
        new{length(cfg.seeds)}(cfg, cache)
    end
end
RD.state_dim(::DiscreteModel) = 4
RD.control_dim(::DiscreteModel) = 2

function RD.discrete_dynamics(model::DiscreteModel, x, u, t, dt)
    SA[
        x[1] + x[3] * dt,
        x[2] + x[4] * dt,
        x[3] + u[1] * dt,
        x[4] + u[2] * dt,
    ]
end

function RD.discrete_dynamics!(model::DiscreteModel, xn, x, u, t, dt)
    xn[1] = x[1] + x[3] * dt
    xn[2] = x[2] + x[4] * dt
    xn[3] = x[3] + u[1] * dt
    xn[4] = x[4] + u[2] * dt
    return nothing
end

function RD.jacobian!(model::DiscreteModel, J, xn, x, u, t, dt)
    J .= 0
    J[1,1] = 1
    J[1,3] = dt
    J[2,2] = 1
    J[2,4] = dt
    J[3,3] = 1
    J[3,5] = dt
    J[4,4] = 1
    J[4,6] = dt
    return nothing
end

function RD.dynamics_error_jacobian!(model::DiscreteModel, J2, J1, y2, y1, z2::RD.AbstractKnotPoint, z1::RD.AbstractKnotPoint)
    n = RD.state_dim(model)
    RD.jacobian!(model, J1, y2, z1)
    J2 .= 0
    for i = 1:n
        J2[i,i] = -1.0
    end
    return nothing
end

# Explicit Dynamics derivatives
#   These are automatically defined by @autodiff
function RobotDynamics.jacobian!(::StaticReturn, ::ForwardAD, model::DiscreteModel, J, y, z)
    f(_z) = RD.evaluate(model, RD.getstate(z, _z), RD.getcontrol(z, _z), RD.getparams(z))
    J .= ForwardDiff.jacobian(f, getdata(z))
end

function RobotDynamics.jacobian!(::InPlace, ::ForwardAD, model::DiscreteModel, J, y, z)
    f!(_y, _z) = RD.evaluate!(model, _y, RD.getstate(z, _z), RD.getcontrol(z, _z), RD.getparams(z))
    J .= ForwardDiff.jacobian!(J, f!, y, getdata(z), model.cfg)
end

function RobotDynamics.jacobian!(::StaticReturn, ::FiniteDifference, model::DiscreteModel, J, y, z)
    f!(_y,_z) = _y .= RD.evaluate(model, RD.getstate(z, _z), RD.getcontrol(z, _z), RD.getparams(z))
    FiniteDiff.finite_difference_jacobian!(J, f!, RD.getdata(z), model.cache)
end

function RobotDynamics.jacobian!(::InPlace, ::FiniteDifference, model::DiscreteModel, J, y, z)
    f!(_y,_z) = RD.evaluate!(model, _y, RD.getstate(z, _z), RD.getcontrol(z, _z), RD.getparams(z))
    FiniteDiff.finite_difference_jacobian!(J, f!, RD.getdata(z), model.cache)
end

# Dynamics error (implicit dynamics) derivatives
#   These are automatically defined by @autodiff for DiscreteDynamics models
function RD.dynamics_error_jacobian!(::StaticReturn, ::ForwardAD, model::DiscreteModel, J2, J1, y2, y1, z2::KnotPoint, z1::KnotPoint)
    f1(_z) = RD.dynamics_error(model, z2, RD.StaticKnotPoint(z1, _z))
    J1 .= ForwardDiff.jacobian(f1, RD.getdata(z1))

    f2(_z) = RD.dynamics_error(model, RD.StaticKnotPoint(z2, _z), z1)
    J2 .= ForwardDiff.jacobian(f2, RD.getdata(z2))
    return nothing
end

function RD.dynamics_error_jacobian!(::InPlace, ::ForwardAD, model::DiscreteModel, J2, J1, y2, y1, z2::KnotPoint, z1::KnotPoint)
    f1!(_y, _z) = RD.dynamics_error!(model, _y, y1, z2, RD.StaticKnotPoint(z1, _z))
    ForwardDiff.jacobian!(J1, f1!, y2, RD.getdata(z1), model.cfg)

    f2!(_y, _z) = RD.dynamics_error!(model, _y, y1, RD.StaticKnotPoint(z2, _z), z1)
    ForwardDiff.jacobian!(J1, f1!, y2, RD.getdata(z2), model.cfg)
    return nothing
end

function RD.dynamics_error_jacobian!(::StaticReturn, ::FiniteDifference, model::DiscreteModel, J2, J1, y2, y1, z2::KnotPoint, z1::KnotPoint)
    f1!(_y, _z) = _y .= RD.dynamics_error(model, z2, RD.StaticKnotPoint(z1, _z))
    FiniteDiff.finite_difference_jacobian!(J1, f1!, RD.getdata(z1), model.cache)

    f2!(_y, _z) = _y .= RD.dynamics_error(model, RD.StaticKnotPoint(z2, _z), z1)
    FiniteDiff.finite_difference_jacobian!(J2, f2!, RD.getdata(z2), model.cache)
    return nothing
end

function RD.dynamics_error_jacobian!(::InPlace, ::FiniteDifference, model::DiscreteModel, J2, J1, y2, y1, z2::KnotPoint, z1::KnotPoint)
    f1!(_y, _z) = RD.dynamics_error!(model, _y, y1, z2, RD.StaticKnotPoint(z1, _z))
    FiniteDiff.finite_difference_jacobian!(J1, f1!, RD.getdata(z1), model.cache)

    f2!(_y, _z) = RD.dynamics_error!(model, _y, y1, RD.StaticKnotPoint(z2, _z), z1)
    FiniteDiff.finite_difference_jacobian!(J2, f2!, RD.getdata(z2), model.cache)
    return nothing
end

model = DiscreteModel()
test_model(model)

# Test dynamics error
z1 = KnotPoint(randn(model)..., 0, 0.1)
x2 = RD.evaluate(model, z1)
z2 = KnotPoint(x2, RD.control(z1)*0, 0.1, z1.dt)
RD.control(z1)
@test RD.dynamics_error(model, z2, z1) == zeros(4)

# Test dynamics error jacobian
n,m = size(model)
J1 = zeros(n, n+m); J2 = copy(J1)
y1 = zeros(n);      y2 = copy(y1)
J10 = copy(J1);     J20 = copy(J2)
RD.dynamics_error_jacobian!(model, J20, J10, y2, y1, z2, z1)
RD.dynamics_error_jacobian!(StaticReturn(), ForwardAD(), model, J2, J1, y2, y1, z2, z1)
@test J10 ≈ J1
@test J20 ≈ J2
RD.dynamics_error_jacobian!(StaticReturn(), RD.UserDefined(), model, J2, J1, y2, y1, z2, z1)
@test J10 ≈ J1
@test J20 ≈ J2

test_error_allocs(model)

## Discrete model w/ autogen
@autodiff struct DiscreteModelAuto <: RobotDynamics.DiscreteDynamics end
RD.state_dim(::DiscreteModelAuto) = 4
RD.control_dim(::DiscreteModelAuto) = 2

function RD.discrete_dynamics(model::DiscreteModelAuto, x, u, t, dt)
    SA[
        x[1] + x[3] * dt,
        x[2] + x[4] * dt,
        x[3] + u[1] * dt,
        x[4] + u[2] * dt,
    ]
end

function RD.discrete_dynamics!(model::DiscreteModelAuto, xn, x, u, t, dt)
    xn[1] = x[1] + x[3] * dt
    xn[2] = x[2] + x[4] * dt
    xn[3] = x[3] + u[1] * dt
    xn[4] = x[4] + u[2] * dt
    return nothing
end

function RD.jacobian!(model::DiscreteModelAuto, J, xn, x, u, t, dt)
    J .= 0
    J[1,1] = 1
    J[1,3] = dt
    J[2,2] = 1
    J[2,4] = dt
    J[3,3] = 1
    J[3,5] = dt
    J[4,4] = 1
    J[4,6] = dt
    return nothing
end

model = DiscreteModelAuto()
test_model(model)
test_error_allocs(model)


text:
lines = readlines("input.txt")

function main(lines)
    init = "shiny gold" => 1

    part1 = let tree = Dict{String, Vector{String}}()
        for line in lines
            outer, innards = split(line, " bags contain ")
            for inner in filter(!isempty, split(innards, r" ?bags?[.,] ?"))
                bagtype = replace(inner, r"\d+ ?" => "")
                if haskey(tree, bagtype)
                    push!(tree[bagtype], outer)
                else
                    tree[bagtype] = [outer]
                end
            end
        end

        stack = String[first(init)]
        works = Set{String}()
        while !isempty(stack)
            top = pop!(stack)
            if top ∉ works
                if haskey(tree, top)
                    append!(stack, tree[top])
                end
                push!(works, top)
            end
        end
        length(works) - 1
    end

    part2 = let tree = Dict{String, Vector{Pair{String, Int}}}()
        for line in lines
            outer, innards = split(line, " bags contain ")
            for inner in filter(!isempty, split(innards, r" ?bags?[.,] ?"))
                m = match(r"(\d+)\s([\w\s]+)", inner)
                if isnothing(m)
                    tree[outer] = []
                else
                    n = parse(Int, m[1])
                    bagtype = m[2]
                    if haskey(tree, outer)
                        push!(tree[outer], bagtype => n)
                    else
                        tree[outer] = [bagtype => n]
                    end
                end
            end
        end

        stack = Pair{String, Int}[init]
        numbags = 0
        while !isempty(stack)
            bag, count = pop!(stack)
            numbags += count
            if haskey(tree, bag)
                for (b, c) in tree[bag]
                    push!(stack, b => count * c)
                end
            end
        end

        numbags - last(init)
    end

    @info "Solutions" part1 part2
end

main(lines)

max_stars_repo_name: dglmoore/advent-of-code

max_issues_repo_name: dglmoore/advent-of-code

max_forks_repo_name: dglmoore/advent-of-code


text:
using PayOff
using Test

@testset "PayOff.jl" begin
    # Write your tests here.
end


text:
@testset "Model Methods: Fast Base (Quadruped)" begin
	s = ContactImplicitMPC.get_simulation("quadruped", "flat_2D_lc", "flat")
	model = s.model
	env = s.env

	# Setup variables
	T = Float64
	nq = model.nq
	q0s = rand(nq)
	q̇0s = rand(nq)

	# Testing fast methods
	@test norm(ContactImplicitMPC.M_fast(model, q0s) - ContactImplicitMPC.M_func(model, q0s), Inf) < 1.0e-8
	@test norm(ContactImplicitMPC.B_fast(model, q0s) - ContactImplicitMPC.B_func(model, q0s), Inf) < 1.0e-8
	@test norm(ContactImplicitMPC.A_fast(model, q0s) - ContactImplicitMPC.A_func(model, q0s), Inf) < 1.0e-8
	@test norm(ContactImplicitMPC.C_fast(model, q0s, q̇0s) - ContactImplicitMPC.C_func(model, q0s, q̇0s), Inf) < 1.0e-8

	nq = model.nq
	nu = model.nu
	nw = model.nw
	nc = model.nc
	nb = nc * friction_dim(env)

	hs = 0.1 * ones(1)
	q0s = rand(nq)
	q1s = rand(nq)
	u1s = rand(nu)
	u1s = rand(nu)
	w1s = rand(nw)
	γ1s = rand(nc)
	b1s = rand(nb)
	λ1s = rand(nc * dim(env))
	q2s = rand(nq)

	# Testing dynamics methods
	@test norm(s.con.d(hs, q0s, q1s, u1s, w1s, λ1s, q2s) -
		ContactImplicitMPC.dynamics(model, env, hs, q0s, q1s, u1s, w1s, λ1s, q2s), Inf) < 1.0e-8
end