text:
using DelimitedFiles
using Distributions
using IterTools
using NonUniformRandomVariateGeneration
using Printf
using Random
using SpecialFunctions
using Statistics
using StatsBase

"""
This is a test file. It runs exactly the same as pmdi but returns some internals.
Should not be used.
"""
function __pmdi(dataFiles, dataTypes, N::Int64, particles::Int64,
    ρ::Float64, iter::Int64;
     featureSelect::Union{String, Nothing} = nothing,
     dataNames = nothing)

    K       = length(dataFiles) # No. of datasets
    n_obs   = size(dataFiles[1], 1)

    # Set names if not specified
    if dataNames == nothing
        dataNames = ["K$i" for i in 1:K]
    end

    @assert length(dataTypes) == K "Number of datatypes not equal to number of datasets"
    @assert length(dataNames) == K "Number of data names not equal to number of datasets"
    @assert all(x->x==n_obs, [size(dataFiles[k])[1] for k = 1:K]) "Datasets don't have same number of observations. Each row must correspond to the same underlying observational unit across datasets."
    @assert (ρ < 1) && (ρ > 0) "ρ must be between 0 and 1"
    @assert (N <= n_obs) & (N > 1) "Number of clusters must be greater than 1 and not greater than the number of observations, suggest using floor(log(n)) = $(floor(Int64, log(n_obs)))"
    @assert particles > 1 "Conditional particle filter requires 2 or more particles"

    # Initialise the hyperparameters
    M = ones(Float64, K) .* 2 # Mass parameter
    γc = rand(Gamma(1.0 / N, 1), N, K) .+ eps(Float64) # Component weights
    Φ = K > 1 ? rand(Gamma(1, 0.2), Int64(K * (K - 1) * 0.5)) : zeros(1) # Dataset concordance measure

    # Initialise allocations randomly according to γc
    s = Matrix{Int64}(undef, n_obs, K)
    for k = 1:K
        s[:, k] = sampleCategorical(n_obs, γc[:, k])
    end

    # Get a matrix of all combinations of allocations
    c_combn = Matrix{Int64}(undef, N ^ K, K)
    for k in 1:K
        c_combn[:, K - k + 1] = div.(0:(N ^ K - 1), N ^ (K - k)) .% N .+ 1
    end

    # The corresponding gammas
    Γc = Matrix{Float64}(undef, N ^ K, K)
    log_γ = log.(γc)
    for k = 1:K
        Γc[:, k] = view(log.(γc[:, k]), c_combn[:, k])
    end

    # Which Φ value is activated by each of the above combinations
    Φ_index = K > 1 ? Matrix{Bool}(undef, N ^ K, Int64(K * (K - 1) / 2)) : fill(1, (N, 1))
    if K > 1
        i = 1
        for k1 in 1:(K - 1)
            for k2 in (k1 + 1):K
                Φ_index[:, i] = (c_combn[:, k1] .== c_combn[:, k2])
                i += 1
            end
        end
    end

    # Normalising constant
    Z = update_Z(Φ, Φ_index, Γc)
    v = update_v(n_obs, Z)

    # Particle weights
    logweight = zeros(Float64, particles)
    # Mutation weights
    logprob = Matrix{Float64}(undef, N * particles + 1, K)
    fprob = Vector{Float64}(undef, N)
    cluster_update = zeros(Bool, N * particles + 1)
    n_operations = 0

    # Feature selection index
    featureFlag = [rand(Bool, size(dataFiles[k], 2)) for k in 1:K]
    if featureSelect == nothing
        for k in 1:K
            featureFlag[k] .= true
        end
    else
        featureNames = ["$(dataNames[k])_d$d" for k in 1:K for d in 1:size(dataFiles[k], 2)]
        writedlm(featureSelect, reshape(featureNames, 1, length(featureNames)), ',')
        featureFile = open(featureSelect, "a")
        writedlm(featureFile, [featureFlag...;]', ',')
    end

    # Feature select probabilities
    featureProb = [zeros(Float64, size(dataFiles[k], 2)) for k in 1:K]
    featureNull = [zeros(Float64, size(dataFiles[k], 2)) for k in 1:K]
    nullCluster = [dataTypes[k](dataFiles[k]) for k in 1:K]
    for k = 1:K
        for i = 1:n_obs
            cluster_add!(nullCluster[k], dataFiles[k][i, :], ones(Bool, size(dataFiles[k], 2)))
        end
        featureNull[k] = - calc_logmarginal(nullCluster[k])
    end

    # particle matches the cluster labels to the cluster IDs
    particle = ones(Int, N, particles, K)
    particle_id = ones(Int, particles, K)
    fprob_dict = Matrix{Float64}(undef, N + 1, particles)
    fprob_done = Vector{Bool}(undef, particles)
    new_id = zeros(Int, N, particles, K)
    logprob = Matrix{Float64}(undef, N * particles + 1, K)
    logprob_particle = view(logprob, particle)

    # A vector containing all of the clusters
    clusters = [Vector{dataTypes[k]}(undef, N * particles + 1) for k in 1:K]
    # Keep track of how many copies of a cluster exists
    clusters_counts = zeros(Int, N * particles + 1, K)
    clusters_counts[1, :] .= particles * N

    sstar_id = Matrix{Int64}(undef, particles, K)
    sstar = zeros(Int64, particles, n_obs, K)

    order_obs = collect(1:n_obs)
    n1 = floor(Int64, ρ * n_obs)

    @inbounds for it in 1:iter
        fill!(clusters_counts, 0)
        clusters_counts[1, :] .= particles * N
        fill!(new_id, 0)
        fill!(particle_id, 1)
        for i in eachindex(particle)
            particle[i] = 1
        end
        shuffle!(order_obs)

        # Update hyperparameters
        if K > 1
            update_Φ!(Φ, v, s, Φ_index, γc, K, Γc)
        end
        update_γ!(γc, Φ, v, M, s, Φ_index, c_combn, Γc, N, K)
        log_γ = log.(γc)

        Π = γc ./ sum(γc, dims = 1)
        Z = update_Z(Φ, Φ_index, Γc)
        v = update_v(n_obs, Z)
        update_M!(M, γc, K, N)

        for k = 1:K
            clusters[k][1] = dataTypes[k](dataFiles[k])
            clust_ids = Dict{Int, Int}()
            id = 2
            us = unique(s[order_obs[1:(n1 - 1)], k])
            for u in us
                clusters[k][id] = dataTypes[k](dataFiles[k])
                clusters_counts[id, k] = particles
                clusters_counts[1, k] -= particles
                clust_ids[u] = id
                particle[u, :, k] .= id
                id += 1
            end
            for i in order_obs[1:(n1 - 1)]
                # id = findall((in)(s[i, k]), us)[1] + 1
                id = clust_ids[s[i, k]]
                sstar[:, i, k] .= s[i, k]
                cluster_add!(clusters[k][id], dataFiles[k][i, :], featureFlag[k])
            end
        end

        for i in order_obs[n1:n_obs]
            for k in 1:K
                fill!(cluster_update, false)
                fill!(fprob_done, false)
                particle_k = view(particle, :, :, k)
                logprob_particle = view(logprob, particle_k, k)
                obs = view(dataFiles[k], i, :)
                sstar_id_k = view(sstar_id, :, k)
                Π_k = view(Π, :, k)
                for id in 1:maximum(particle_k)
                    logprob[id, k] = calc_logprob(obs, clusters[k][id], featureFlag[k])
                    n_operations += 1
                end
                # Draw the new allocations
                curr_id = 0
                for p in 1:particles
                    id = particle_id[p, k]
                    if fprob_done[id]
                        for n in 1:N
                            @inbounds fprob[n] = fprob_dict[n, id]
                        end
                        logweight[p] += fprob_dict[end, id]
                    else
                        for n in 1:N
                             fprob[n] = logprob_particle[n, p]
                         end
                         max_logprob = maximum(fprob)
                         for n in 1:N
                             fprob[n] -= max_logprob
                             fprob[n] = exp(fprob[n])
                             fprob[n] *= Π_k[n]
                        end
                        fprob = cumsum!(fprob, fprob)
                        logprob_inc = log(fprob[N]) + max_logprob
                        fprob_dict[end, id] =  logprob_inc
                        logweight[p] += logprob_inc
                        for n in 1:N
                            fprob[n] = fprob[n] / fprob[N]
                        end
                        fprob_dict[1:N, particle_id[p, k]] = fprob
                        fprob_done[particle_id[p, k]] = true
                    end
                    # Set reference trajectory
                    if p != 1
                        new_s = 1
                        u = rand()
                        for i in 1:(N - 1)
                            if fprob[new_s] > u
                                break
                            else
                                new_s += 1
                            end
                        end
                    else
                        new_s = s[i, k]
                    end
                    sstar_id_k[p] = particle[new_s, p, k]
                    sstar[p, i, k] = new_s
                    if new_id[new_s, id, k] == 0
                        curr_id += 1
                        new_id[new_s, id, k] = curr_id
                        particle_id[p, k] = curr_id
                    else
                        particle_id[p, k] = new_id[new_s, id, k]
                    end
                end
                # Add observation to new cluster
                max_k = maximum(particle_k)
                for p in sstar_id_k
                    if !cluster_update[p]
                        cluster_update[p] = true
                        ncopies = count(x -> x == p, sstar_id_k)
                        if ncopies == clusters_counts[p, k]
                            # If the origin cluster still exists somewhere
                            # Need to create a new cluster
                            # with obs added to it
                            id = p
                        else
                            id = max_k + 1
                            clusters_counts[p, k] -= ncopies
                            clusters_counts[id, k] = ncopies
                            # clusters[k][id] = deepcopy(clusters[k][p])
                            # clusters[k][id] = copy_particle(clusters[k][p], dataFiles[k])
                            clusters[k][id] = deepcopy(clusters[k][p])
                            max_k += 1
                        end
                        cluster_add!(clusters[k][id], obs, featureFlag[k])
                        if id !== p
                            for part in 1:particles
                                s_id = sstar[part, i, k]
                                if particle[s_id, part, k] == p
                                    particle[s_id, part, k] = id
                                end
                            end
                        end
                    end
                end
            end
            if K > 1
                Φ_upweight!(logweight, sstar[:, i, :], K, Φ, particles)
            end

            # Resampling
            if calc_ESS(logweight) <=  0.5 * particles
                partstar = draw_partstar(logweight, particles)
                logweight .= 1.0
                for k in 1:K
                    particle[:, :, k] = particle[:, partstar, k]
                    particle_id[:, k] = particle_id[partstar, k]
                    sstar[:, :, k] = sstar[partstar, :, k]
                    # Reset clusters_counts
                    clusters_counts[:, k] .= 0
                    # Renumber the clusters to ensure elements don't grow too large
                    particle_k = view(particle, :, :, k)
                    for (i, id) in enumerate(sort(unique(particle[:, :, k])))
                        if id !== i
                            for j in eachindex(particle_k)
                                if particle_k[j] == id
                                    particle_k[j] = i
                                end
                            end
                            clusters[k][i] = deepcopy(clusters[k][id])
                        end
                        clusters_counts[i, k] = count(x -> x == i, particle_k)
                    end
                end
            end
        end

        # Select a single particle
        p_weights = similar(logweight)
        max_logweight = maximum(logweight)
        for i in eachindex(p_weights)
            p_weights[i] = exp(logweight[i] - max_logweight)
        end
        p_star = sample(1:particles, Weights(p_weights))

        # Feature selection
        ## Create a null particle with every obs in one cluster
        if featureSelect != nothing
            ## Compare this to the marginal likelihood for each cluster
            for k = 1:K
                featureProb[k] = featureNull[k] .+ 0
                occupiedClusters = unique(sstar[p_star, :, k])
                for clust in occupiedClusters
                    clust_members = findindices(sstar[p_star, :, k], clust)
                    clust_params = dataTypes[k](dataFiles[k])
                    for obs in clust_members
                        cluster_add!(clust_params, dataFiles[k][obs, :], ones(Bool, size(dataFiles[k], 2)))
                    end
                    featureProb[k] += calc_logmarginal(clust_params)
                end
                featureFlag[k] = (1 .- 1 ./ (exp.(featureProb[k] .+ 1))) .> rand(length(featureProb[k]))
                # featurePosterior[k] += featureFlag[k] ./ (iter + 1)
            end
        end

        logweight .= 1.0
        s[:] = sstar[p_star, :, :]
        # Match up labels across datasets
        align_labels!(s, Φ, γc, N, K)

    end
    if featureSelect != nothing
        close(featureFile)
    end
    return particle, clusters, clusters_counts, n_operations
end


text:
# Read and write files

export
    checkfilepath,
    readdssp,
    readpdb,
    fixstring,
    spaceatomname,
    writepdb,
    readensemble,
    writeensemble,
    writeensemblescores,
    writeprojections,
    writepcview,
    writepcviews,
    writeintarray,
    writefloatarray,
    writestringarray,
    readpocketpoints,
    readligsite,
    readpdblines,
    writeclusterpoints

"Checks if an output filepath is valid, and throws an error if not."
function checkfilepath(out_filepath::AbstractString)
    directory = splitdir(out_filepath)[1]
    @assert directory == "" || isdir(directory) "Not a valid file location: \"$out_filepath\""
end

"""
Reads a DSSP file. If atoms are provided, creates additional '-' entries for residues not present in the DSSP file.
Returns a dictionary where key is residue number and chain, and value is secondary structure.
"""
function readdssp(dssp_filepath::AbstractString, atoms::Array{Atom,1})
    dssps = readdssp(dssp_filepath)
    res_pres = collect(keys(dssps))
    for atom in atoms
        res = "$(atom.res_n)$(atom.chain_id)"
        if !(res in res_pres)
            dssps[res] = '-'
            push!(res_pres, res)
        end
    end
    return dssps
end

function readdssp(dssp_filepath::AbstractString)
    expanded_path = expanduser(dssp_filepath)
    @assert isfile(expanded_path) "Not a valid filepath: \"$expanded_path\""
    dssps = Dict{String, Char}()
    reading = false
    open(expanded_path, "r") do dssp_file
        for line in eachline(dssp_file)
            if reading && length(line) >= 17 && line[14] != '!'
                res = "$(strip(line[7:10]))$(line[12])"
                sec_str = line[17]
                dssps[res] = sec_str
            elseif !reading && length(line) >= 12 && line[1:12] == "  #  RESIDUE"
                reading = true
            end
        end
    end
    if length(dssps) == 0
        println("Nothing read from DSSP file")
    else
        println("Read ", length(dssps), " residues from DSSP file")
    end
    return dssps
end

"""
Parse a PDB file.
Returns a list of atoms.
hetatm determines whether HETATM records are read in.
"""
function readpdb(in_filepath::AbstractString; hetatm::Bool=false)
    expanded_path = expanduser(in_filepath)
    @assert isfile(expanded_path) "Not a valid filepath: \"$expanded_path\""
    atoms = Atom[]
    open(expanded_path, "r") do in_file
        if hetatm
            line_starts = ["ATOM  ", "HETATM"]
        else
            line_starts = ["ATOM  "]
        end
        for line in eachline(in_file)
            if length(line) >= 54 && line[1:6] in line_starts
                # Atom name string is columns 13-16 with leading/trailing spaces removed
                atom_name = strip(line[13:16])
                # Residue name string is columns 18-20
                res_name = line[18:20]
                # Chain character is column 22
                chain_id = line[22]
                # Residue number integer is columns 23-26 converted to an integer
                res_n = parse(Int, line[23:26])
                # x, y and z coordinates are columns 31-38, 39-46 and 47-54 respectively converted to floats
                x = parse(Float64, line[31:38])
                y = parse(Float64, line[39:46])
                z = parse(Float64, line[47:54])
                # Element is determined from atom name, not from columns 77-78
                element = inferelement(atom_name)
                # Remove disorder
                if (line[17] == ' ' || line[17] == 'A')
                    push!(atoms, Atom(atom_name, res_name, chain_id, res_n, [x, y, z], element))
                end
            end
        end
    end
    println("Read ", length(atoms), " atoms from PDB file")
    return atoms
end

"""
Converts a string, character or number to a fixed length string.
Be careful - could e.g. cut 0.5 to "0".
Also does not round correctly.
"""
function fixstring(val_in, new_length::Integer; fill_char::Char=' ')
    string_out = string(val_in)
    string_out = string_out[1:min(length(string_out), new_length)]
    while length(string_out) < new_length
        string_out = "$fill_char$string_out"
    end
    return string_out
end

"""
Space an atom name correctly in the four-character limit.
The element is in position 2.
Returns a string of length 4.
"""
function spaceatomname(atom_name::AbstractString, element::AbstractString)
    chars = length(atom_name)
    @assert chars <= 4 "Atom name is greater than four characters"
    if element != "-"
        fi = findfirst(x -> x == element[1], atom_name)
        cent_ind = fi == nothing ? 0 : fi
    else
        cent_ind = 1
    end
    @assert cent_ind <= 2 "Atom name is too long to align correctly"
    @assert chars - cent_ind <= 3 "Atom name is too long to align correctly"
    if cent_ind == 1 && chars < 4
        out_string = " "
    else
        out_string = ""
    end
    out_string = "$out_string$atom_name"
    while length(out_string) < 4
        out_string = "$out_string "
    end
    return out_string
end

"""
Write a list of atoms as a PDB file.
The atoms are not reordered so it is assumed they are in a sensible order, e.g. the order they were read in.
Only ATOM records are written; no TER labels.
Could write a brief header for generated models.
"""
function writepdb(out_filepath::AbstractString, atoms::Array{Atom,1})
    expanded_path = expanduser(out_filepath)
    checkfilepath(expanded_path)
    open(expanded_path, "w") do out_file
        println(out_file, "REMARK   Generated by ExProSE on ", Libc.strftime(time()))
        for (atom_counter, atom) in enumerate(atoms)
            line = String[
                "ATOM  ",
                fixstring(atom_counter, 5),
                " ",
                spaceatomname(atom.atom_name, atom.element),
                " ",
                fixstring(atom.res_name, 3),
                " ",
                string(atom.chain_id),
                fixstring(atom.res_n, 4),
                "      ",
                fixstring(atom.coords[1], 6),
                "  ",
                fixstring(atom.coords[2], 6),
                "  ",
                fixstring(atom.coords[3], 6),
                "  1.00  0.00          ",
                fixstring(atom.element, 2),
                "  "
            ]
            @assert length(join(line)) == 80 "PDB line not correct length"
            println(out_file, line...)
        end
    end
end

"Read an ensemble of structures and return a ModelledEnsemble."
function readensemble(out_prefix::AbstractString, n_to_read::Integer)
    @assert n_to_read > 0 "n_to_read must be positive"
    atoms = readpdb("$out_prefix.1")
    strucs = ModelledStructure[]
    for i in 1:n_to_read
        curr_filepath = "$out_prefix.$i"
        @assert isfile(curr_filepath) "Not a valid filepath: \"$curr_filepath\""
        curr_atoms = readpdb(curr_filepath)
        # Placeholder score of -1.0 given to structures read in
        struc = ModelledStructure(-1.0, atomcoords(curr_atoms))
        push!(strucs, struc)
    end
    println("Read in protein ensemble")
    return ModelledEnsemble(atoms, strucs)
end

"Write an ensemble out as separate PDB files."
function writeensemble(out_prefix::AbstractString, ensemble::ModelledEnsemble)
    checkfilepath(out_prefix)
    atoms = ensemble.atoms
    strucs = ensemble.strucs
    @assert length(atoms) > 0 "No atoms associated with ensemble"
    @assert length(strucs) > 0 "The ensemble does not contain any structures"
    for (i, struc) in enumerate(strucs)
        out_filepath = "$out_prefix.$i"
        struc_atoms = atomcoords(atoms, struc.coords)
        writepdb(out_filepath, struc_atoms)
    end
    println("Wrote ensemble to PDB file(s) with prefix \"$out_prefix\"")
end

"Write out the scores associated with the structures in an ensemble."
function writeensemblescores(out_filepath::AbstractString, ensemble::ModelledEnsemble)
    expanded_path = expanduser(out_filepath)
    checkfilepath(expanded_path)
    strucs = ensemble.strucs
    @assert length(strucs) > 0 "The ensemble does not contain any structures"
    open(expanded_path, "w") do out_file
        for struc in strucs
            println(out_file, round(struc.score, digits=1))
        end
    end
    println("Wrote ensemble scores to file \"$expanded_path\"")
end

"""
Write out the projections of structures onto the eigenvectors of a PCA.
pcs_to_write is a list of the PCs to write, e.g. [1, 2].
All PCs written by default.
"""
function writeprojections(out_filepath::AbstractString, pcs::Array{Float64}, pcs_to_write::Array{Int,1}=collect(1:size(pcs, 1)))
    expanded_path = expanduser(out_filepath)
    checkfilepath(expanded_path)
    @assert size(pcs, 1) >= maximum(pcs_to_write) "Chose to write out a PC higher than the number of PCs present"
    open(expanded_path, "w") do out_file
        for i in 1:size(pcs, 2)
            vals_out = getindex(pcs[:,i], pcs_to_write)
            string_out = join([string(val) for val in vals_out], "\t")
            println(out_file, string_out)
        end
    end
    println("Projections written to file \"$expanded_path\"")
end

"Write a PyMol script to view the atom motions of a principal component."
function writepcview(out_filepath::AbstractString,
                    pca::PCA,
                    atoms::Array{Atom,1},
                    pc::Integer;
                    weighting::Real=50.0,
                    line_colour::AbstractString="magenta")
    expanded_path = expanduser(out_filepath)
    checkfilepath(expanded_path)
    @assert pc > 0 "Selected principal component must be positive"
    @assert pc <= size(pca.evecs, 2) "Selected principal component is too high"
    inds_to_use = calphaindices(atoms)
    @assert size(pca.evecs, 1) == 3 * length(inds_to_use) "PCA and atom list have different sizes"
    disps = pca.evecs[:,pc]
    open(expanded_path, "w") do out_file
        println(out_file, "set dash_gap, 0")
        println(out_file, "set dash_color, ", line_colour)
        for (i, j) in enumerate(inds_to_use)
            pc_coords = atoms[j].coords + disps[3i-2:3i] * weighting
            println(out_file, "pseudoatom pt", 2i-1, ", pos=[", round(atoms[j].coords[1], digits=2), ", ", round(atoms[j].coords[2], digits=2), ", ", round(atoms[j].coords[3], digits=2), "]")
            println(out_file, "pseudoatom pt", 2i, ", pos=[", round(pc_coords[1], digits=2), ", ", round(pc_coords[2], digits=2), ", ", round(pc_coords[3], digits=2), "]")
            println(out_file, "distance dist", i, ", /pt", 2i-1, ", /pt", 2i)
        end
        println(out_file, "hide everything, pt*")
        println(out_file, "hide labels")
    end
end

"Wrapper function for writepcview that produces visualisation scripts for the top `n_to_write` PCs."
function writepcviews(out_prefix::AbstractString,
                    pca::PCA,
                    atoms::Array{Atom,1},
                    n_to_write::Integer)
    expanded_prefix = expanduser(out_prefix)
    checkfilepath(expanded_prefix)
    if n_to_write <= size(pca.evecs, 2)
        n_to_use = n_to_write
    else
        n_to_use = size(pca.evecs, 2)
        println("n_to_write is larger than the number of eigenvectors - proceeding with all possible eigenvectors")
    end
    for i in 1:n_to_use
        out_filepath = "$expanded_prefix$i.pml"
        writepcview(out_filepath, pca, atoms, i)
    end
    println("PyMol script to view principal components written to file(s) \"$(expanded_prefix)x.pml\"")
end

"Write an array of integers to an output file."
function writeintarray(out_filepath::AbstractString, ints::Array{Int,1})
    expanded_path = expanduser(out_filepath)
    checkfilepath(expanded_path)
    open(expanded_path, "w") do out_file
        for i in ints
            println(out_file, i)
        end
    end
end

"Write an array of floats to an output file."
function writefloatarray(out_filepath::AbstractString, vals::Array{Float64,1}; dec_places::Integer=3)
    expanded_path = expanduser(out_filepath)
    checkfilepath(expanded_path)
    open(expanded_path, "w") do out_file
        for val in vals
            println(out_file, round(val, digits=dec_places))
        end
    end
end

"Write an array of strings to an output file."
function writestringarray(out_filepath::AbstractString, vals::Array{<:AbstractString,1})
    expanded_path = expanduser(out_filepath)
    checkfilepath(expanded_path)
    open(expanded_path, "w") do out_file
        for val in vals
            println(out_file, val)
        end
    end
end

"""
Read a custom pocket point file in PDB format where the residue number is the pocket number.
This file can be generated using a script from the LIGSITEcs pocket_r.pdb and pocket_all.pdb files.
Returns a dictionary where key is pocket number and value is point coords.
"""
function readpocketpoints(pdb_filepath::AbstractString)
    expanded_path = expanduser(pdb_filepath)
    @assert isfile(expanded_path) "Not a valid filepath: \"$expanded_path\""
    pock_points = Dict{Int, Array{Float64}}()
    counter = 0
    open(expanded_path, "r") do pdb_file
        for line in eachline(pdb_file)
            if startswith(line, "ATOM  ") || startswith(line, "HETATM")
                counter += 1
                x = parse(Float64, line[31:38])
                y = parse(Float64, line[39:46])
                z = parse(Float64, line[47:54])
                pocket_n = parse(Int, line[23:26])
                if haskey(pock_points, pocket_n)
                    pock_points[pocket_n] = hcat(pock_points[pocket_n], [x, y, z])
                else
                    pock_points[pocket_n] = [x, y, z]
                end
            end
        end
    end
    println("Read ", counter, " pocket points from pocket points PDB file")
    return pock_points
end

"""
Read a LIGSITEcs PDB output file.
Returns the pocket centre coordinate list and list of pocket volumes.
"""
function readligsite(pdb_filepath::AbstractString)
    expanded_path = expanduser(pdb_filepath)
    @assert isfile(expanded_path) "Not a valid filepath: \"$expanded_path\""
    xs = Float64[]
    ys = Float64[]
    zs = Float64[]
    vols = Int[]
    counter = 0
    open(expanded_path, "r") do pdb_file
        for line in eachline(pdb_file)
            if startswith(line, "ATOM  ") || startswith(line, "HETATM")
                counter += 1
                x = parse(Float64, line[31:38])
                y = parse(Float64, line[39:46])
                z = parse(Float64, line[47:54])
                vol = parse(Int, line[23:26])
                push!(xs, x)
                push!(ys, y)
                push!(zs, z)
                push!(vols, vol)
            end
        end
    end
    coords = transpose(hcat(xs, ys, zs))
    println("Read ", counter, " pocket centres from LIGSITEcs output PDB file")
    return coords, vols
end

"Read the ATOM lines from a PDB file into an array."
function readpdblines(pdb_filepath::AbstractString)
    expanded_path = expanduser(pdb_filepath)
    @assert isfile(expanded_path) "Not a valid filepath: \"$expanded_path\""
    point_lines = String[]
    open(expanded_path, "r") do in_file
        for line in eachline(in_file)
            if startswith(line, "ATOM  ")
                push!(point_lines, chomp(line))
            end
        end
    end
    return point_lines
end

"Write output file with assignment in residue number column"
function writeclusterpoints(out_filepath::AbstractString,
                    point_lines::Array{<:AbstractString,1},
                    assignments::Array{Int,1})
    expanded_path = expanduser(out_filepath)
    checkfilepath(expanded_path)
    open(expanded_path, "w") do out_file
        println(out_file, "REMARK    Pocket number in residue number column")
        println(out_file, "REMARK    Pocket numbers found using cluster-ligsite")
        for (i, line) in enumerate(point_lines)
            print(out_file, line[1:22], lpad(string(assignments[i]), 4), line[27:end])
        end
    end
    println("Wrote pocket number assignments to file \"$expanded_path\"")
end

ext: jl


text:
const CONTEXTS = (Training(), Validation(), Inference())

"""
    checktask_core(task, sample, model; device = identity)
    checktask_core(task; device = identity)

Check if `task` conforms to the [core interface](docs/interfaces/core.md).
`sample` and `model` are used for testing. If you have implemented the testing
interface and don't supply these as arguments, `mocksample(task)` and
`mockmodel(task)` will be used.
"""
function checktask_core(
    task;
    model = mockmodel(task),
    sample = mocksample(task),
    devicefn = identity,
)
    Test.@testset "Core interface" begin
        Test.@testset "`encode`" begin
            for context in CONTEXTS
                Test.@test_nowarn encodesample(task, context, sample)
            end
        end
        Test.@testset "Model compatibility" begin
            for context in CONTEXTS
                x, _ = encodesample(task, context, sample)
                Test.@test_nowarn ŷ = _predictx(task, model, x, devicefn)
            end
        end
        Test.@testset "`decodeypred" begin
            for context in CONTEXTS
                x, _ = encodesample(task, context, sample)
                ŷ = _predictx(task, model, x, devicefn)
                Test.@test_nowarn decodeypred(task, context, ŷ)
            end
        end
    end
end

function _predictx(method, model, x, device = identity)
    if shouldbatch(method)
        x = DataLoaders.collate([x])
    end
    ŷs = device(model)(device(x))
    if shouldbatch(method)
        ŷ = ŷs[((:) for _ in 1:ndims(ŷs)-1)..., 1]
    end
    return ŷ
end


text:
include("cli_wrappers.jl")

main()


text:
#!/usr/bin/env julia
include("eur_dmon.jl")
push!(PGFPlotsX.CUSTOM_PREAMBLE, raw"\usepackage{fontspec}")
push!(PGFPlotsX.CUSTOM_PREAMBLE, raw"\setmainfont{Asana-Math}")
nord = ["DK", "FI", "NO", "SE"]
months = [12, 1, 2, 3, 4]

for month in months
    swe_month = swe_months[month]
    monthstr = "M$(lpad("$(month)", 2, '0'))"
    pgfsave("../../images/nord_mort_6018_$(monthstr).svg", plot_geos_month(nord, monthstr))
    print(*("![Dödstal nordiska länder $(swe_month) 1960–2018.]",
        "(../../images/nord_mort_6018_$(monthstr).svg)",
        "{#fig:nord_mort_6018_$(monthstr) width=100%}\n\n"))
end

max_stars_repo_name: klpn/static-dust

max_issues_repo_name: klpn/static-dust

max_forks_repo_name: klpn/static-dust


text:
#TODO make test for this
"""
    Example of 'intake_info' = Dict("gt" => Dict("ub" => 100.0, "c" => 10.0))
"""
function apply_bound!(metnet::MetNet, xi, intake_info::Dict = Dict();
        emptyfirst = false, ignore_miss = false)
    xi < 0 && error("xi must be positive")

    emptyfirst && empty!(metnet.intake_info)
    for (exch, exch_info) in intake_info
        ignore_miss && !(exch in metnet.rxns) && continue
        metnet.intake_info[exch] = exch_info
    end

    for (exch, exch_info) in metnet.intake_info

        !haskey(exch_info, "c") && error("'c' not defined for exchange '$(rxns(metnet, exch))'");

        c = exch_info["c"]
        c < 0.0 && error("negative 'c' ($c) in exchange '$(rxns(metnet, exch))'")

        if haskey(exch_info, "ub")
            ub = exch_info["ub"]
            ub!(metnet, exch, stst_bound(c, xi, ub))
        elseif haskey(exch_info, "lb")
            lb = exch_info["lb"]
            lb!(metnet, exch, stst_bound(c, xi, lb))
        else
            error("neither 'lb' or 'ub' defined for exchange '$(rxns(metnet, exch))'")
        end
    end
    return metnet
end

max_stars_repo_name: josePereiro/Chemostat

max_issues_repo_name: josePereiro/Chemostat

max_forks_repo_name: josePereiro/Chemostat


text:
export problem

"""
`problems(terms...)`

Constructs a problem.

# Example

```julia

julia> x = Variable(4)
Variable(Float64, (4,))

julia> A, b = randn(10,4), randn(10);

julia> p = problem(ls(A*x-b), norm(x) <= 1)

```

"""
function problem(terms::Vararg)
	cf = ()
	for i = 1:length(terms)
		cf = (cf...,terms[i]...)
	end
	return cf
end


text:
using FullyThreadedTree
using BenchmarkTools

function sumvol(tree::Tree)
    val = 0.0
    for cell ∈ cells(tree, filter = active)
        val += volume(cell)
    end
    val
end
function sumvol(vec::Vector{Tree{N}}) where N
    val = 0.0
    for cell ∈ filter(active, vec)
        val += volume(cell)
    end
    val
end
function sumvol_faces(tree::Tree)
    val = 0.0
    for face ∈ all_faces(tree, filter = active)
        val += cell_volume(face)
    end
    val
end
function sumvol_faces(vec::Vector{Face{N}}) where N
    val = 0.0
    for face ∈ filter(active, vec)
        val += cell_volume(face)
    end
    val
end

function compare_performance(dim, nr_steps)
    tree = Tree(zeros(dim))
    child = tree
    for step ∈ 1:nr_steps
        refine!(child)
        child = child.children[mod(step, 1<<dim) + 1]
    end

    println("Collecting cells")
    display(@benchmark collect_cells($tree))

    println("\nCollecting faces")
    display(@benchmark collect_faces($tree))

    cells = collect_cells(tree)
    faces = collect_faces(tree)

    println("\nCells using a vector of Trees")
    display(@benchmark sumvol($cells))

    println("\nFaces using a vector of Faces")
    display(@benchmark sumvol_faces($faces))

    println("\nCells using the Tree data structure")
    display(@benchmark sumvol($tree))

    println("\nFaces using the Tree data structure")
    display(@benchmark sumvol_faces($tree))

    return tree
end


text:
using ForwardDiff: Dual

# import Flatten: ulflatten, reconstruct
# flatten(x::Unitful.Quantity) = (x.val,)
# reconstruct(::T, data, n) where T <: Unitful.Quantity = (unit(T) * data[n],), n + 1
# retype(::T, data, n) where T <: Unitful.Quantity = (unit(T) * data[n],), n + 1
# update!(::T, data, n) where T <: Unitful.Quantity = (unit(T) * data[n],), n + 1

dualize_records(o, du, u) = begin
    vars1 = dualize_vars(o.records[1].vars[1], zero(u[1]))
    vars2 = dualize_vars(o.records[2].vars[1], zero(u[1]))
    Records(o.params, vars=(vars1, vars2), val=zero(du[1]))
end

dualize_vars(vars, val) = begin
    flatvars = ulflatten(Vector, vars)
    dualvars = ForwardDiff.seed!(fill(val, length(flatvars)), flatvars)
    reconstruct(vars, dualvars)
end

(o::Organism)(du::AbstractVector{<:Dual}, u::AbstractVector{<:Dual}, p::AbstractVector{<:Real}, t::Number) = begin
    println("Dual/Real")
    u1 = u .* mol
    o1 = reconstruct(o, p)
    du1 = du .* (mol/hr)
    rec = dualize_records(o, du1, u)
    o1(du1, u1, t * (hr), define_organs(o1.params, rec, o1, t))
    du .= ustrip.(du1)
end
(o::Organism)(du::AbstractVector{<:Dual}, u::AbstractVector{<:Dual}, p::AbstractVector{<:Dual}, t::Number) = begin
    println("Dual/Dual")
    u1 = u .* mol
    o1 = reconstruct(o, p)
    du1 = [zero(p[1]) .* (mol/hr) for d in du]
    rec = dualize_records(o, du1, u)
    define_organs(o1.params, rec, o1, t)
    check_params.(o1.organs)
    o1(du1, u1, t * (hr), )
    du .= ustrip.(du1)
end
(o::Organism)(du::AbstractVector{<:Real}, u::AbstractVector{<:Real}, p::AbstractVector{<:Dual}, t::Number) = begin
    println("Real/Dual")
    u1 = u .* mol
    o1 = reconstruct(o, p)
    du1 = [zero(p[1]) .* (mol/hr) for d in du]
    rec = dualize_records(o, du1, p)
    o1(du1, u1, t * (hr), define_organs(o1.params, rec, o1, t))
    ustrip.(du1)
end

function deb_function(sol)
   tot_loss = 0.0
   if any((s.retcode != :Success for s in sol))
     tot_loss = Inf
   else
     # calculation for the loss here
   end
   tot_loss
end


text:
module Helper

params = Dict(
    :default =>
    Dict(:x0 => 0.0, :steps => 1000000, :dt => 0.1, :sigma => 1.0, :seed => 1),
)

function init(name)
    param = copy(params[name])
    x = zeros(Float64, param[:steps])
    x[:, 1] = param[:x0]
    param[:x] = x
    param
end

end

max_stars_repo_name: tkf/comparatist

max_issues_repo_name: tkf/comparatist

max_forks_repo_name: tkf/comparatist


text:
using Gtk, Graphics, Cairo, Colors
using svgPath

include("samples.jl")
path = parsePathData(Computer)
#path = ScalePath(path, 2)

c = @GtkCanvas()
win = GtkWindow(c, "Canvas")
@guarded draw(c) do widget
    ctx = getgc(c)
    h = height(c)
    w = width(c)
    rectangle(ctx, 50,50, 50, 50)
    set_source_rgb(ctx, 1, 1, 1)
    fill(ctx)

    set_line_width(ctx, 1)
    set_source_rgb(ctx, 0, 0, 0)
    # stroke(ctx)
    drawPath(ctx, path)
    #println(get_current_point(ctx))
    #println(copy_path_flat(ctx))

    stroke(ctx)
end
show(c)


text:
"""
        saw(decisionMat, weights, fns)

Apply SAW (Simple Additive Weighting) method for a given matrix and weights.
This method also known as WSM (Weighted Sum Model)

# Arguments:
 - `decisionMat::DataFrame`: n × m matrix of objective values for n candidate (or strategy) and m criteria
 - `weights::Array{Float64, 1}`: m-vector of weights that sum up to 1.0. If the sum of weights is not 1.0, it is automatically normalized.
 - `fns::Array{Function, 1}`: m-vector of functions to be applied on the columns.

# Description
saw() applies the SAW method to rank n strategies subject to m criteria which are supposed to be
either maximized or minimized.

# Output
- `::SawResult`: SawResult object that holds multiple outputs including scores, rankings, and best index.

# Examples
```julia-repl
julia> decmat = [4.0  7  3  2  2  2  2;
                 4.0  4  6  4  4  3  7;
                 7.0  6  4  2  5  5  3;
                 3.0  2  5  3  3  2  5;
                 4.0  2  2  5  5  3  6];
julia> df = makeDecisionMatrix(decmat)
5×7 DataFrame
 Row │ Crt1     Crt2     Crt3     Crt4     Crt5     Crt6     Crt7
     │ Float64  Float64  Float64  Float64  Float64  Float64  Float64
─────┼───────────────────────────────────────────────────────────────
   1 │     4.0      7.0      3.0      2.0      2.0      2.0      2.0
   2 │     4.0      4.0      6.0      4.0      4.0      3.0      7.0
   3 │     7.0      6.0      4.0      2.0      5.0      5.0      3.0
   4 │     3.0      2.0      5.0      3.0      3.0      2.0      5.0
   5 │     4.0      2.0      2.0      5.0      5.0      3.0      6.0

julia> fns = convert(Array{Function,1}, [maximum for i in 1:7])
7-element Array{Function,1}:
 maximum (generic function with 16 methods)
 maximum (generic function with 16 methods)
 maximum (generic function with 16 methods)
 maximum (generic function with 16 methods)
 maximum (generic function with 16 methods)
 maximum (generic function with 16 methods)
 maximum (generic function with 16 methods)

julia> weights = [0.283, 0.162, 0.162, 0.07, 0.085, 0.162, 0.076];

julia> result = saw(df, weights, fns);

julia> result.scores
5-element Array{Float64,1}:
 0.5532285714285714
 0.7134857142857142
 0.8374285714285714
 0.5146571428571429
 0.5793428571428572

julia> result.bestIndex
3
```

# References
Afshari, Alireza, Majid Mojahed, and Rosnah Mohd Yusuff. "Simple additive weighting approach to
personnel selection problem." International Journal of Innovation, Management and Technology
1.5 (2010): 511.
"""
function saw(decisionMat::DataFrame, weights::Array{Float64,1}, fns::Array{Function,1})::SawResult

    n, p = size(decisionMat)

    normalizedDecisionMat = similar(decisionMat)

        w = unitize(weights)

    colminmax = zeros(Float64, p)
    @inbounds for i in 1:p
        colminmax[i] = decisionMat[:, i] |> fns[i]
        if fns[i] == maximum
            normalizedDecisionMat[:, i] = decisionMat[:, i] ./ colminmax[i]
        elseif fns[i] == minimum
            normalizedDecisionMat[:, i] = colminmax[i] ./ decisionMat[:, i]
        else
            @error fns[i]
            error("Function not found")
        end
    end

    scores = w * normalizedDecisionMat |> rowsums

    rankings = scores |> sortperm |> reverse

    bestIndex = rankings |> first

    result = SawResult(
        decisionMat,
        normalizedDecisionMat,
        w,
        scores,
        rankings,
        bestIndex
    )

    return result
end

"""
        saw(setting)

Apply SAW (Simple Additive Weighting) method for a given matrix and weights.
This method also known as WSM (Weighted Sum Model)

# Arguments:
 - `setting::MCDMSetting`: MCDMSetting object.

# Description
saw() applies the SAW method to rank n strategies subject to m criteria which are supposed to be
either maximized or minimized.

# Output
- `::SawResult`: SawResult object that holds multiple outputs including scores, rankings, and best index.
"""
function saw(setting::MCDMSetting)::SawResult
    saw(
        setting.df,
        setting.weights,
        setting.fns
    )
end

max_stars_repo_name: drvinceknight/JMcDM

max_issues_repo_name: drvinceknight/JMcDM

max_forks_repo_name: drvinceknight/JMcDM


text:
module AppServer

using Sockets
using GitHub
using HTTP
using Base64
using Pkg
using Logging

using ..Config
using ..Utils
using ..Runner
using ..ContinuousBenchmarks: set_benchmark_config_file

const event_queue = Channel{Any}(1024)
const httpsock = Ref{Sockets.TCPServer}()

const report_repo = Config.get_config("github.report_repo")
const report_branch = Config.get_config("github.report_branch", "master")
const bot_user = Config.get_config("github.user")

function bm_from_comment(data)
    data.payload["action"] != "created" && return

    bot_auth = GitHub.authenticate(Config.get_config("github.token"))

    comment_url = data.payload["comment"]["html_url"]
    comment_body = data.payload["comment"]["body"]
    user = data.payload["comment"]["user"]["login"]
    issue_no = data.payload["issue"]["number"]
    issue_url = data.payload["issue"]["html_url"]
    issue_repo = data.repository

    occursin("@$bot_user", comment_body) || return
    branches = benchmark_branches(comment_body)
    if isempty(branches)
        params = Dict("body" => "@$user Yes Sir?")
        create_comment(issue_repo, issue_no, :issue; params=params, auth=bot_auth)
        return
    end
    try
        gitbranches(Config.get_config("target.project_dir", "."), branches)
    catch ex # not all of the branches exist
        params = Dict("body" => "@$user $ex")
        create_comment(issue_repo, issue_no, :issue; params=params, auth=bot_auth)
        return
    end

    name = bm_name(branches)

    # create the branch
    base_br = GitHub.reference(report_repo, "heads/$report_branch"; auth=bot_auth)
    params = Dict("ref" => "refs/heads/" * name, "sha" => base_br.object["sha"])
    GitHub.create_reference(report_repo; params=params, auth=bot_auth)

    # commit bm info file on the new branch
    content = base64encode(Utils.bm_file_content(user, issue_url, comment_url, branches))
    params = Dict("branch" => name,
                  "message" => name,
                  "content" => content)
    commit = GitHub.create_file(report_repo, "jobs/$(name).toml"; params=params, auth=bot_auth)

    # create pull request
    params = Dict("title" => name,
                  "head" => name,
                  "base" => report_branch,
                  "body" => Utils.bm_pr_content(comment_url))
    pr = GitHub.create_pull_request(report_repo; params=params, auth=bot_auth)

    params = Dict("body" => Utils.bm_reply0_content(user, pr.html_url.uri))
    create_comment(issue_repo, issue_no, :issue; params=params, auth=bot_auth)
end

function handle_events(e)
    @debug("Handle Event: ", e)
    e.kind == "issue_comment" && return bm_from_comment(e)
    e.kind == "push" && Config.get_config("server.bm_runner") == "self" &&
        return Runner.trigger(e)
end

function event_handler(event::WebhookEvent)
    global event_queue
    @debug("Received event for $(event.repository.full_name)")
    push!(event_queue, event)
    @debug("event pushed")
    return HTTP.Messages.Response(200)
end

function comment_handler(event::WebhookEvent, phrase::RegexMatch)
    global event_queue
    @debug("Received event for $(event.repository.full_name), phrase: $phrase")
    push!(event_queue, event)
    @debug("event pushed")
    return HTTP.Messages.Response(200)
end

function recover(name, keep_running, do_action, handle_exception;
                 backoff=0, backoffmax=120, backoffincrement=1)
    while keep_running()
        try
            do_action()
            backoff = 0
        catch ex
            exception_action = handle_exception(ex)
            if exception_action == :exit
                @warn("Stopping", name)
                return
            else # exception_action == :continue
                bt = get_backtrace(ex)
                @error("Recovering from unknown exception", name, ex, bt, backoff)
                sleep(backoff)
                backoff = min(backoffmax, backoff+backoffincrement)
            end
        end
    end
    @warn("Stopped", name)
end

function request_processor()
    do_action() = begin
        event = take!(event_queue)
        tsk = @task handle_events(event)
        schedule(tsk)
    end
    handle_exception(ex) =
        (isa(ex, InvalidStateException) && (ex.state == :closed)) ? :exit : :continue
    keep_running() = isopen(event_queue)
    recover("request_processor", keep_running, do_action, handle_exception)
end

function github_webhook(
    http_ip=Config.get_config("server.http_ip"),
    http_port=Config.get_config("server.http_port")
)
    app_auth = GitHub.JWTAuth(
        Config.get_config("github.app_id"), Config.get_config("github.priv_pem"))
    # trigger = Regex("@$bot_user")
    # listener = GitHub.CommentListener(
    #    comment_handler,
    #    trigger;
    #    check_collab=false,
    #    auth=app_auth,
    #    secret=Config.get_config("github.secret"))

    listener = GitHub.EventListener(
        event_handler,
        auth=app_auth,
        secret=Config.get_config("github.secret"),
        events=["push", "issue_comment"])
    httpsock[] = Sockets.listen(IPv4(http_ip), http_port)

    do_action() = GitHub.run(listener, httpsock[], IPv4(http_ip), http_port, verbose=true)
    handle_exception(ex) = (isa(ex, Base.IOError) && (ex.code == -103)) ? :exit : :continue
    keep_running() = isopen(httpsock[])

    @info("GitHub Webhook starting...", http_ip, http_port)
    recover("github_webhook", keep_running, do_action, handle_exception)
end

function status_monitor()
    stop_file = config["server"]["stop_file"]
    while isopen(event_queue)
        sleep(5)
        flush(stdout); flush(stderr);
        # stop server if stop is requested
        if isfile(stop_file)
            @warn("Server stop requested.")
            flush(stdout); flush(stderr)
            # stop accepting new requests
            close(httpsock[])
            # wait for queued requests to be processed and close queue
            while isready(event_queue)
                yield()
            end
            close(event_queue)
            rm(stop_file; force=true)
        end
    end
end

function main()
    if isempty(get(ENV, "ENV", "")) && isempty(get(ENV, "TBB_CONFIG", ""))
        println("Usage: \n\t" *
                "export TBB_CONFIG=/path/to/config/file # or export ENV=dev|prod|...\n\t" *
                "julia -e 'using ContinuousBenchmarks; ContinuousBenchmarks.AppServer.main()' ")
        return
    end

    global_logger(SimpleLogger(stdout, Config.log_level()))

    project_dir = Config.get_config("target.project_dir")
    benchmark_config_file = Config.get_config("target.benchmark_config_file")
    set_benchmark_config_file(joinpath(project_dir, benchmark_config_file))

    @info("Starting server...")
    handler_task = @async request_processor()
    monitor_task = @async status_monitor()
    github_webhook()
    wait(handler_task)
    wait(monitor_task)
    @warn("Server stopped.")
end

end


text:
module UnsafeAtomicsTests

include("test_orderings.jl")
include("test_core.jl")

end  # module UnsafeAtomicsTests


text:
# methods
abstract type AbstractESSMethod end

"""
    ESSMethod <: AbstractESSMethod

The `ESSMethod` uses a standard algorithm for estimating the
effective sample size of MCMC chains.

It is is based on the discussion by [^Vehtari2019] and uses the
biased estimator of the autocovariance, as discussed by [^Geyer1992].
In contrast to Geyer, the divisor `n - 1` is used in the estimation of
the autocovariance to obtain the unbiased estimator of the variance for lag 0.

[^Geyer1992]: Geyer, C. J. (1992). Practical Markov Chain Monte Carlo. Statistical Science, 473-483. <https://projecteuclid.org/euclid.ss/1177011137>.

[^Vehtari2019]: Vehtari, A., Gelman, A., Simpson, D., Carpenter, B., & Bürkner, P. C. (2021). Rank-normalization, folding, and localization: An improved ``\\widehat {R}`` for assessing convergence of MCMC. Bayesian Analysis. <https://arxiv.org/pdf/1903.08008.pdf>.
"""
struct ESSMethod <: AbstractESSMethod end

"""
    FFTESSMethod <: AbstractESSMethod

The `FFTESSMethod` uses a standard algorithm for estimating
the effective sample size of MCMC chains.

It is is based on the discussion by [^Vehtari2019] and uses the
biased estimator of the autocovariance, as discussed by [^Geyer1992].
In contrast to Geyer, the divisor `n - 1` is used in the estimation of
the autocovariance to obtain the unbiased estimator of the variance for lag 0.

In contrast to [`ESSMethod`](@ref), this method uses fast Fourier transforms
(FFTs) for estimating the autocorrelation.
"""
struct FFTESSMethod <: AbstractESSMethod end

"""
    BDAESSMethod <: AbstractESSMethod

The `BDAESSMethod` uses a standard algorithm for estimating the effective sample size of
MCMC chains.

It is is based on the discussion by [^Vehtari2019] and uses the
variogram estimator of the autocorrelation function discussed in [^Gelman2013].

[^Gelman2013]: Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., & Rubin, D. B. (2013). Bayesian data analysis. CRC press.
"""
struct BDAESSMethod <: AbstractESSMethod end

# caches
struct ESSCache{T,S}
    samples::Matrix{T}
    chain_var::Vector{S}
end

struct FFTESSCache{T,S,C,P,I}
    samples::Matrix{T}
    chain_var::Vector{S}
    samples_cache::C
    plan::P
    invplan::I
end

mutable struct BDAESSCache{T,S,M}
    samples::Matrix{T}
    chain_var::Vector{S}
    mean_chain_var::M
end

function build_cache(::ESSMethod, samples::Matrix, var::Vector)
    # check arguments
    niter, nchains = size(samples)
    length(var) == nchains || throw(DimensionMismatch())

    return ESSCache(samples, var)
end

function build_cache(::FFTESSMethod, samples::Matrix, var::Vector)
    # check arguments
    niter, nchains = size(samples)
    length(var) == nchains || throw(DimensionMismatch())

    # create cache for FFT
    T = complex(eltype(samples))
    n = nextprod([2, 3], 2 * niter - 1)
    samples_cache = Matrix{T}(undef, n, nchains)

    # create plans of FFTs
    fft_plan = plan_fft!(samples_cache, 1)
    ifft_plan = plan_ifft!(samples_cache, 1)

    return FFTESSCache(samples, var, samples_cache, fft_plan, ifft_plan)
end

function build_cache(::BDAESSMethod, samples::Matrix, var::Vector)
    # check arguments
    nchains = size(samples, 2)
    length(var) == nchains || throw(DimensionMismatch())

    return BDAESSCache(samples, var, mean(var))
end

update!(cache::ESSCache) = nothing

function update!(cache::FFTESSCache)
    # copy samples and add zero padding
    samples = cache.samples
    samples_cache = cache.samples_cache
    niter, nchains = size(samples)
    n = size(samples_cache, 1)
    T = eltype(samples_cache)
    @inbounds for j in 1:nchains
        for i in 1:niter
            samples_cache[i, j] = samples[i, j]
        end
        for i in (niter + 1):n
            samples_cache[i, j] = zero(T)
        end
    end

    # compute unnormalized autocovariance
    cache.plan * samples_cache
    @. samples_cache = abs2(samples_cache)
    cache.invplan * samples_cache

    nothing
end

function update!(cache::BDAESSCache)
    # recompute mean of within-chain variances
    cache.mean_chain_var = mean(cache.chain_var)

    return
end

function mean_autocov(k::Int, cache::ESSCache)
    # check arguments
    samples = cache.samples
    niter, nchains = size(samples)
    0 ≤ k < niter || throw(ArgumentError("only lags ≥ 0 and < $niter are supported"))

    # compute mean of unnormalized autocovariance estimates
    firstrange = 1:(niter - k)
    lastrange = (k + 1):niter
    s = mean(1:nchains) do i
        if eltype(samples) isa BlasReal
            # call into BLAS if possible
            x = dot(samples, firstrange, samples, lastrange)
            firstrange = firstrange .+ niter
            lastrange = lastrange .+ niter
            return x
        else
            # otherwise use views
            return dot(view(samples, firstrange, i), view(samples, lastrange, i))
        end
    end

    # normalize autocovariance estimators by `niter - 1` instead
    # of `niter - k` to obtain
    # - unbiased estimators of the variance for lag 0
    # - biased but more stable estimators for all other lags as discussed by
    #   Geyer (1992)
    return s / (niter - 1)
end

function mean_autocov(k::Int, cache::FFTESSCache)
    # check arguments
    niter, nchains = size(cache.samples)
    0 ≤ k < niter || throw(ArgumentError("only lags ≥ 0 and < $niter are supported"))

    # compute mean autocovariance
    # we use biased but more stable estimators as discussed by Geyer (1992)
    samples_cache = cache.samples_cache
    chain_var = cache.chain_var
    return mean(1:nchains) do i
        real(samples_cache[k + 1, i]) / real(samples_cache[1, i]) * chain_var[i]
    end
end

function mean_autocov(k::Int, cache::BDAESSCache)
    # check arguments
    samples = cache.samples
    niter, nchains = size(samples)
    0 ≤ k < niter || throw(ArgumentError("only lags ≥ 0 and < $niter are supported"))

    # compute mean autocovariance
    n = niter - k
    idxs = 1:n
    s = mean(1:nchains) do j
        return sum(idxs) do i
            abs2(samples[i, j] - samples[k + i, j])
        end
    end

    return cache.mean_chain_var - s / (2 * n)
end

"""
    ess(chains::Chains; kwargs...)

Estimate the effective sample size and the potential scale reduction.
"""
function ess(
    chains::Chains;
    sections = _default_sections(chains),
    kwargs...
)
    # subset the chain
    _chains = Chains(chains, _clean_sections(chains, sections))

    # estimate the effective sample size and rhat
    ess, rhat = ess_rhat(_chains.value.data; kwargs...)

    # convert to namedtuple
    nt = merge((parameters = names(_chains),), (ess = ess, rhat = rhat))

    return ChainDataFrame("ESS", nt)
end

function ess_rhat(
    chains::AbstractArray{<:Union{Missing,Real},3};
    method::AbstractESSMethod = ESSMethod(),
    maxlag::Int = 250
)
    # compute size of matrices (each chain is split!)
    niter = size(chains, 1) ÷ 2
    nparams = size(chains, 2)
    nchains = 2 * size(chains, 3)
    ntotal = niter * nchains

    # do not compute estimates if there is only one sample or lag
    maxlag = min(maxlag, niter - 1)
    maxlag > 0 || return fill(missing, nparams), fill(missing, nparams)

    # define caches for mean and variance
    U = typeof(zero(eltype(chains)) / 1)
    T = promote_type(eltype(chains), typeof(zero(eltype(chains)) / 1))
    chain_mean = Array{T}(undef, 1, nchains)
    chain_var = Array{T}(undef, nchains)
    samples = Array{T}(undef, niter, nchains)

    # compute correction factor
    correctionfactor = (niter - 1) / niter

    # define cache for the computation of the autocorrelation
    esscache = build_cache(method, samples, chain_var)

    # define output arrays
    ess = Vector{T}(undef, nparams)
    rhat = Vector{T}(undef, nparams)

    # for each parameter
    for (i, chains_slice) in enumerate(eachslice(chains; dims = 2))
        # check that no values are missing
        if any(x -> x === missing, chains_slice)
            rhat[i] = missing
            ess[i] = missing
            continue
        end

        # split chains
        copyto_split!(samples, chains_slice)

        # calculate mean of chains
        mean!(chain_mean, samples)

        # calculate within-chain variance
        @inbounds for j in 1:nchains
            chain_var[j] = var(view(samples, :, j); mean = chain_mean[j], corrected = true)
        end
        W = mean(chain_var)

        # compute variance estimator var₊, which accounts for between-chain variance as well
        var₊ = correctionfactor * W + var(chain_mean; corrected = true)
        inv_var₊ = inv(var₊)

        # estimate the potential scale reduction
        rhat[i] = sqrt(var₊ / W)

        # center the data around 0
        samples .-= chain_mean

        # update cache
        update!(esscache)

        # compute the first two autocorrelation estimates
        # by combining autocorrelation (or rather autocovariance) estimates of each chain
        ρ_odd = 1 - inv_var₊ * (W - mean_autocov(1, esscache))
        ρ_even = one(ρ_odd) # estimate at lag 0 is known

        # sum correlation estimates
        pₜ = ρ_even + ρ_odd
        sum_pₜ = pₜ

        k = 2
        while k < maxlag
            # compute subsequent autocorrelation of all chains
            # by combining estimates of each chain
            ρ_even = 1 - inv_var₊ * (W - mean_autocov(k, esscache))
            ρ_odd = 1 - inv_var₊ * (W - mean_autocov(k + 1, esscache))

            # stop summation if p becomes non-positive
            Δ = ρ_even + ρ_odd
            Δ > zero(Δ) || break

            # generate a monotone sequence
            pₜ = min(Δ, pₜ)

            # update sum
            sum_pₜ += pₜ

            # update indices
            k += 2
        end

        # estimate the effective sample size
        τ = 2 * sum_pₜ - 1
        ess[i] = ntotal / τ
    end

    return ess, rhat
end

"""
	copyto_split!(out::AbstractMatrix, x::AbstractMatrix)

Copy the elements of matrix `x` to matrix `out`, in which each column of `x` is split.

If the number of rows in `x` is odd, the sample at index `(size(x, 1) + 1) / 2` is dropped.
"""
function copyto_split!(out::AbstractMatrix, x::AbstractMatrix)
    # check dimensions
    nrows_out, ncols_out = size(out)
    nrows_x, ncols_x = size(x)
    ncols_out == 2 * ncols_x ||
        throw(DimensionMismatch("the output matrix must have twice as many columns as the input matrix"))
    nrows_out == nrows_x ÷ 2 ||
        throw(DimensionMismatch("the output matrix must have half as many rows as as the input matrix"))

    jout = 0
    offset = iseven(nrows_x) ? nrows_out : nrows_out + 1
    @inbounds for j in 1:ncols_x
        jout += 1
        for i in 1:nrows_out
            out[i, jout] = x[i, j]
        end

        jout += 1
        ix = offset
        for i in 1:nrows_out
            ix += 1
            out[i, jout] = x[ix, j]
        end
    end

    out
end


text:
"""
    LogNumber{S,L}

Represent a real or complex number in terms of its sign and logarithm.

This is a hack around the limitation that standard floating-point
types link the bitcount of the exponent to that of the mantissa, which
makes it impossible to represent very large or small values without
simultaneously increasing the precision of the representation. The
drawback of this hack is that values x close to 1 have a larger relative
accuracy than values far away from 1 since

    1 - exp(l*(1±ε)) / exp(l) == 1 - exp(±l*ε) ≈ ∓l*ε.

It remains to be seen whether this has any practical relevance.

A better workaround would be to create a floating-point type with
arbitrary exponent and mantissa.
"""
struct LogNumber{S,L} <: Number
    sign::S
    logabs::L
end

LogNumber(x::Number) where {S,L} = LogNumber(sign(x),logabs(x))
LogNumber{S,L}(x::Number) where {S,L} = LogNumber{S,L}(sign(x),logabs(x))
Base.AbstractFloat(x::LogNumber) = x.sign*exp(x.logabs)

Base.sign(x::LogNumber) = x.sign

"""
    logabs(x) -> log(abs(x))
"""
logabs(x::Number) = log(abs2(x))/2
logabs(x::LogNumber) = x.logabs

Base.promote_rule(::Type{LogNumber{S,L}}, ::Type{T}) where {S,L,T<:Number} =
    LogNumber{promote_type(S,typeof(sign(one(T)))), promote_type(L,typeof(logabs(one(T))))}

Base.:*(x::LogNumber,y::LogNumber) = LogNumber(x.sign * y.sign, x.logabs + y.logabs)
Base.:/(x::LogNumber,y::LogNumber) = LogNumber(x.sign * conj(y.sign), x.logabs - y.logabs)


text:
# # [Overdetermined Beam](@ref overdetermined)
#
# This example shows how to predict the behavior of a beam which is clamped at one end and
# simply supported at the other end when subjected to a linear distributed load.
#
# ![](../assets/overdetermined-drawing.svg)
#
#-
#md # !!! tip
#md #     This example is also available as a Jupyter notebook:
#md #     [`overdetermined.ipynb`](@__NBVIEWER_ROOT_URL__/examples/overdetermined.ipynb).
#-

using GXBeam, LinearAlgebra

nelem = 16

## create points
L = 1
x = range(0, L, length=nelem+1)
y = zero(x)
z = zero(x)
points = [[x[i],y[i],z[i]] for i = 1:length(x)]

## index of endpoints for each beam element
start = 1:nelem
stop = 2:nelem+1

## create compliance matrix for each beam element
EI = 1e7
compliance = fill(Diagonal([0, 0, 0, 0, 1/EI, 0]), nelem)

## create assembly
assembly = Assembly(points, start, stop, compliance=compliance)

## set prescribed conditions
prescribed_conditions = Dict(
    ## simply supported left endpoint
    1 => PrescribedConditions(uz=0),
    ## clamped right endpoint
    nelem+1 => PrescribedConditions(ux=0, uy=0, uz=0, theta_x=0, theta_y=0,
        theta_z=0)
)

## create distributed load
qmax = 1000
distributed_loads = Dict()
for i = 1:nelem
    distributed_loads[i] = DistributedLoads(assembly, i; s1=x[i],
        s2=x[i+1], fz = (s) -> qmax*s)
end

## perform a static analysis
system, converged = static_analysis(assembly;
    prescribed_conditions = prescribed_conditions,
    distributed_loads = distributed_loads,
    linear = true)

## post-process the results
state = AssemblyState(system, assembly;
    prescribed_conditions = prescribed_conditions)
#!jl nothing #hide

# An analytical solution to this over-determined problem can be found using the method of
# superposition using the analytical solutions to a cantilever with a linearly distributed
# load and a cantilever with an end load.

## construct analytical solution
x_a = range(0.0, L, length=1000)
w_a = @. qmax*(1-x_a)^2/(120*EI)*(4 - 8*(1-x_a) + 5*(1-x_a)^2 - (1-x_a)^3)
slope_a = @. -qmax*(1-x_a)/(120*EI)*(8 - 24*(1-x_a) + 20*(1-x_a)^2 - 5*(1-x_a)^3)
M_a = @. qmax/120*(8 - 48*(1-x_a) + 60*(1-x_a)^2 - 20*(1-x_a)^3)

theta_a = -atan.(slope_a)

## adjust coordinate system of analytical solution
M_a = -M_a
#!jl nothing #hide

# Plotting the results reveals that the analytical and computational solutions show
# excellent agreement.

using Plots
#md using Suppressor #hide
pyplot()
#!jl nothing #hide

#-

#md @suppress_err begin #hide
#md local x #hide

plot(
    xlim = (0.0, 1.0),
    xticks = 0.0:0.2:1.0,
    xlabel = "x (m)",
    ylabel = "Linear Deflection (m)",
    grid = false,
    overwrite_figure=false
    )

x = [assembly.points[ipoint][1] + state.points[ipoint].u[1] for ipoint =
    1:length(assembly.points)]
deflection = [state.points[ipoint].u[3] for ipoint = 1:length(assembly.points)]
plot!(x_a, w_a, label="Analytical")
scatter!(x, deflection, label="GXBeam")
#!nb plot!(show=true)
#md savefig("../assets/overdetermined-linear-deflection.svg") #hide
#md closeall() #hide
#md end #hide
#md nothing #hide

#md # ![](../assets/overdetermined-linear-deflection.svg)

#-

#md @suppress_err begin #hide
#md local x #hide

plot(
    xlim = (0.0, 1.0),
    xticks = 0.0:0.2:1.0,
    xlabel = "x (m)",
    ylabel = "Angular Deflection (rad)",
    grid = false,
    overwrite_figure=false
    )

x = [assembly.points[ipoint][1] + state.points[ipoint].u[1]
    for ipoint = 1:length(assembly.points)]
theta = [4*atan.(state.points[ipoint].theta[2]/4)
    for ipoint = 1:length(assembly.points)]
plot!(x_a, theta_a, label="Analytical")
scatter!(x, theta, label="GXBeam")
#!nb plot!(show=true)
#md savefig("../assets/overdetermined-angular-deflection.svg") #hide
#md closeall() #hide
#md end #hide
#md nothing #hide

#md # ![](../assets/overdetermined-angular-deflection.svg)

#-

#md @suppress_err begin #hide
#md local x #hide

plot(
    xlim = (0.0, 1.0),
    xticks = 0.0:0.2:1.0,
    xlabel = "x (m)",
    ylabel = "Bending Moment (\$Nm\$)",
    grid = false,
    overwrite_figure=false
    )

x = [assembly.elements[ielem].x[1] + state.elements[ielem].u[1] for ielem =
    1:length(assembly.elements)]
M = [state.elements[ielem].M[2] for ielem = 1:length(assembly.elements)]
plot!(x_a, M_a, label="Analytical")
scatter!(x, M, label="GXBeam")
#!nb plot!(show=true)
#md savefig("../assets/overdetermined-bending-moment.svg") #hide
#md closeall() #hide
#md end #hide
#md nothing #hide

#md # ![](../assets/overdetermined-bending-moment.svg)

# Note that we could have easily performed a nonlinear analysis for this problem by
# setting `linear=false`.


text:
export MostPopular

"""

    MostPopular(data::DataAccessor)

Recommend most popular items.
"""
struct MostPopular <: Recommender
    data::DataAccessor
    scores::AbstractVector

    function MostPopular(data::DataAccessor)
        n_item = size(data.R, 2)
        new(data, vector(n_item))
    end
end

isbuilt(recommender::MostPopular) = isfilled(recommender.scores)

function build!(recommender::MostPopular)
    n_item = size(recommender.data.R, 2)

    for i in 1:n_item
        recommender.scores[i] = count(!iszero, recommender.data.R[:, i])
    end
end

function ranking(recommender::MostPopular, u::Int, i::Int)
    check_build_status(recommender)
    recommender.scores[i]
end


text:
using ExponentialUtilities, YaoArrayRegister

export TimeEvolution, time_evolve

"""
    TimeEvolution{N, TT, GT} <: PrimitiveBlock{N}

TimeEvolution, where GT is block type. input matrix should be hermitian.

!!!note:
    `TimeEvolution` contructor check hermicity of the input block by default, but sometimes it can be slow. Turn off the check manually by specifying optional parameter `check_hermicity = false`.
"""
mutable struct TimeEvolution{N, T, Tt, Hamilton <: AbstractBlock{N}} <: PrimitiveBlock{N}
    H::BlockMap{Complex{T}, Hamilton}
    dt::Tt
    tol::T

    function TimeEvolution(
        H::BlockMap{Complex{T}, TH},
        dt::Tt, tol::T; check_hermicity::Bool=true) where {N, Tt, T, TH <: AbstractBlock{N}}
        (check_hermicity && !ishermitian(H)) && error("Time evolution Hamiltonian has to be a Hermitian")
        return new{N, T, Tt, TH}(H, dt, tol)
    end
end

"""
    TimeEvolution(H, dt[; tol::Real=1e-7])

Create a [`TimeEvolution`](@ref) block with Hamiltonian `H` and time step `dt`. The
`TimeEvolution` block will use Krylove based `expv` to calculate time propagation.

Optional keywords are tolerance `tol` (default is `1e-7`)
`TimeEvolution` block can also be used for
[imaginary time evolution](http://large.stanford.edu/courses/2008/ph372/behroozi2/) if dt is complex.
"""
TimeEvolution(H::AbstractBlock, dt; tol::Real=1e-7, check_hermicity=true) =
    TimeEvolution(BlockMap(H), dt, tol, check_hermicity=check_hermicity)

TimeEvolution(M::BlockMap, dt; tol::Real, check_hermicity=true) =
    TimeEvolution(M, dt, tol, check_hermicity=check_hermicity)

time_evolve(M::BlockMap, dt; kwargs...) = TimeEvolution(M, dt; kwargs...)
time_evolve(M::AbstractBlock, dt; kwargs...) = TimeEvolution(M, dt; kwargs...)
time_evolve(dt; kwargs...) = @λ(M->time_evolve(M, dt; kwargs...))

function mat(::Type{T}, te::TimeEvolution{N}) where {T, N}
    A = Matrix{T}(te.H.block)
    return exp(-im*T(te.dt) * A)
end

function apply!(reg::ArrayReg, te::TimeEvolution)
    st = state(reg)
    dt = real(te.dt) == 0 ? imag(te.dt) : -im*te.dt
    @inbounds for j in 1:size(st, 2)
        v = view(st, :, j)
        Ks = arnoldi(te.H, v; tol=te.tol, ishermitian=true)
        expv!(v, dt, Ks)
    end
    return reg
end

cache_key(te::TimeEvolution) = (te.dt, cache_key(te.H))

# parametric interface
niparams(::Type{<:TimeEvolution}) = 1
getiparams(x::TimeEvolution) = x.dt
setiparams!(r::TimeEvolution, param::Real) = (r.dt = param; r)

function Base.:(==)(lhs::TimeEvolution, rhs::TimeEvolution)
    return lhs.H == rhs.H && lhs.dt == rhs.dt
end

function Base.adjoint(te::TimeEvolution)
    return TimeEvolution(te.H, -adjoint(te.dt); tol=te.tol)
end
Base.copy(te::TimeEvolution) = TimeEvolution(te.H, te.dt, tol=te.tol)

function YaoBase.isunitary(te::TimeEvolution)
    iszero(imag(te.dt)) || return false
    return true
end


text:
importall Bukdu

type VueController <: ApplicationController
end

index(::VueController) = "Hello Bukdu"

Router() do
    get("/", VueController, index)
end

Endpoint() do
    plug(Plug.Static, at= "/", from= "public")
    plug(Plug.Logger)
    plug(Router)
end

Bukdu.start(8080)

# wait()

# Bukdu.stop()


text:
@testset "ContingencyTables" begin

    @testset "Creation and Getters" begin

        for alphabet in (UngappedAlphabet(),
                         GappedAlphabet(),
                         ReducedAlphabet("(AILMV)(NQST)(RHK)(DE)(FWY)CGP"))
            for N in 1:3
                table = ContingencyTable(Float64, Val{N}, alphabet) # zeros in MIToS 1.0

                @test size(table) == (Int[length(alphabet) for i in 1:N]...,)
                @test length(table) == length(alphabet)^N
                @test length(getmarginals(table)) == length(alphabet)* N
                @test size(getmarginals(table)) == (length(alphabet), N)
                @test sum(gettable(table)) == 0.0
                @test sum(table) == 0.0 # == gettotal(table)
                @test sum(table.temporal) == 0.0
                @test sum(getmarginals(table)) == 0.0
                @test gettotal(table) == 0.0
                @test collect(table) == gettablearray(table) # Iteration interface in MIToS 1.0
                @test isa(gettablearray(table), Array{Float64,N})
                @test isa(getmarginalsarray(table), Array{Float64,2})
            end
        end

        @testset "Similar" begin
            for alphabet in (UngappedAlphabet(),
                            GappedAlphabet(),
                            ReducedAlphabet("(AILMV)(NQST)(RHK)(DE)(FWY)CGP"))
                for N in 1:3
                    table = ContingencyTable(Float64, Val{N}, alphabet) # zeros in MIToS 1.0

                    @test table == similar(table)
                    @test typeof(table) == typeof(similar(table))
                    @test table == similar(table, BigFloat)
                    @test typeof(table) != typeof(similar(table, BigFloat))
                    @test BigFloat == eltype(similar(table, BigFloat))
                end
            end
        end
    end

    @testset "Indexing" begin

        for alphabet in (UngappedAlphabet(),
                         GappedAlphabet(),
                         ReducedAlphabet("(AILMV)(NQST)(RHK)(DE)(FWY)CGP"))

            table = ContingencyTable(Float64, Val{2}, alphabet)

            if isa(getalphabet(table), ReducedAlphabet)
                @test table[1,3] == 0.0
                table[1,3] = 100.0
                i = 2 * length(getalphabet(table)) + 1 # using one index
                @test table[i] == 100.0
                table[i] = 10.0
                @test table[Residue('A'), Residue('R')] == 10.0 # using two indices
                table[Residue('A'), Residue('R')] = 20.0
                @test table["AILMV", "RHK"] == 20.0
                table["AILMV", "RHK"] = 30.0
                @test table[1,3] == 30.0
            else
                @test table[1,2] == 0.0
                table[1,2] = 100.0
                i = length(getalphabet(table)) + 1 # using one index
                @test table[i] == 100.0
                table[i] = 10.0
                @test table[Residue('A'), Residue('R')] == 10.0 # using two indices
                table[Residue('A'), Residue('R')] = 20.0
                @test table["A", "R"] == 20.0
                table["A", "R"] = 30.0
                @test table[1,2] == 30.0
             end
        end
    end

    @testset "Update" begin

        for alphabet in (UngappedAlphabet(),
                         GappedAlphabet(),
                         ReducedAlphabet("(AILMV)(NQST)(RHK)(DE)(FWY)CGP"))
            for N in 1:3

                table = ContingencyTable(Float64, Val{N}, alphabet)

                fill!(table.temporal, 1.0)
                @test sum(table.temporal) == 22.0^N
                @test sum(table) == 0.0
                @test sum(getmarginals(table)) == 0.0
                @test gettotal(table) == 0.0

                Information._update!(table)
                @test sum(table.temporal) == 22.0^N
                if isa(getalphabet(table),ReducedAlphabet)
                    @test table[1] == 5.0^N
                    @test sum(table) == 20.0^N
                    @test sum(getmarginals(table)) == N*(20.0^N)
                    @test gettotal(table) == 20.0^N
                    if N == 1
                        @test vec(getarray(getmarginals(table))) ==
                            vec(getarray(gettable(table)))
                    elseif N == 2
                        @test getmarginals(table)[1] ==
                            sum(5.0*n for n in [5,4,3,2,3,1,1,1])
                    elseif N == 2
                        @test getmarginals(table)[1] ==
                            sum(5.0*n*20.0 for n in [5,4,3,2,3,1,1,1])
                    end
                else
                    @test table[1] == 1.0
                    @test getmarginals(table)[1] == length(alphabet)^(N-1)
                    @test sum(table) == length(alphabet)^N
                    @test sum(getmarginals(table)) == float(N)*(length(alphabet)^N)
                    @test gettotal(table) == float(length(alphabet))^N
                end
            end
        end
    end

    @testset "Fill" begin

        for alphabet in (UngappedAlphabet(),
                         GappedAlphabet(),
                         ReducedAlphabet("(AILMV)(NQST)(RHK)(DE)(FWY)CGP"))
            for N in 1:3

                table = ContingencyTable(Float64, Val{N}, alphabet)

                fill!(table, AdditiveSmoothing(1.0))
                @test table[1] == 1.0
                @test getmarginals(table)[1] == length(alphabet)^(N-1)
                @test sum(table) == length(alphabet)^N
                @test sum(getmarginals(table)) == float(N)*(length(alphabet)^N)
                @test gettotal(table) == float(length(alphabet))^N
            end
        end
    end

    @testset "Pseudocount" begin

        @test AdditiveSmoothing(1.0) ==  one(AdditiveSmoothing{Float64})
        @test AdditiveSmoothing(0.0) == zero(AdditiveSmoothing{Float64})

        for alphabet in (UngappedAlphabet(),
                         GappedAlphabet(),
                         ReducedAlphabet("(AILMV)(NQST)(RHK)(DE)(FWY)CGP"))

            for N in 1:3

                table = ContingencyTable(Float64, Val{N}, alphabet)

                apply_pseudocount!(table, zero(AdditiveSmoothing{Float64}))

                @test sum(table.temporal) == 0.0
                @test sum(table) == 0.0
                @test sum(getmarginals(table)) == 0.0
                @test gettotal(table) == 0.0

                apply_pseudocount!(table, AdditiveSmoothing(1.0))

                @test table[1] == 1.0
                @test getmarginals(table)[1] == length(alphabet)^(N-1)
                @test sum(table) == length(alphabet)^N
                @test sum(getmarginals(table)) == float(N)*(length(alphabet)^N)
                @test gettotal(table) == float(length(alphabet))^N

                table = ContingencyTable(Float64, Val{N}, alphabet)
                apply_pseudocount!(table, 1.0)

                @test table[1] == 1.0
                @test getmarginals(table)[1] == length(alphabet)^(N-1)
                @test sum(table) == length(alphabet)^N
                @test sum(getmarginals(table)) == float(N)*(length(alphabet)^N)
                @test gettotal(table) == float(length(alphabet))^N
            end
        end
    end

end


text:
VarianceFunction(::Type{Binomial{T}}, μ̂::Vector{<:Real}) where {T<:Real} = μ̂ .* (1 .- μ̂)

VarianceFunction(::Type{Normal{T}}, μ̂::Vector{<:Real}) where {T<:Real} = 1

VarianceFunction(::Type{Poisson{T}}, μ̂::Vector{<:Real}) where {T<:Real} = μ̂

pearson(d::UnivariateDistribution, Y, μ̂) = (Y .- μ̂) ./ sqrt(VarianceFunction(typeof(d), μ̂))

function phi(e)
    ϕ, m, n = 0, length(e), maximum(length.(e))
    for i in 1:m
        for j in 1:length(e[i])
            ϕ += e[i][j]^2
        end
    end
    ϕ / (n * m)
end

# D' = X'A where A is diagonal matrix with Var(Y_ij|X_ij), i.e., variance function
updateD(d, X, μ̂) = X .* VarianceFunction(typeof(d), μ̂)

updateV(corstr, μ̂) = sqrt(diagm(μ̂)) * corstr() * sqrt(diagm(μ̂))


text:
# Copyright (c) 2020-2021, Bank of Canada
# All rights reserved.

@testset "TSeries" begin
    # test constructors
    s = TSeries(20Q1, collect(10.0 .+ (1:12)))
    @test typeof(s) === TSeries{Quarterly,Float64,Array{Float64,1}}
    @test size(s) == (12,)
    @test axes(s) == (20Q1:22Q4,)
    @test length(s) == 12

    t = TSeries(Int, 5)   # from type and number
    @test typeof(t) === TSeries{Unit,Int,Vector{Int}} && t.firstdate == 1U && length(t.values) == 5
    t = TSeries(UInt8, 4 .+ (1:5)) # from type and Int range
    @test typeof(t) === TSeries{Unit,UInt8,Vector{UInt8}} && t.firstdate == 5U && length(t.values) == 5
    t = TSeries(Float32, 1:5, undef) # from type, range and undef
    @test typeof(t) === TSeries{Unit,Float32,Vector{Float32}} && t.firstdate == 1U && length(t.values) == 5

    # constructing with similar()
    t = similar(ones(Float64, 5), (2Q1:4Q4))
    @test typeof(t) === TSeries{Quarterly,Float64,Vector{Float64}} && t.firstdate == 2Q1 && length(t.values) == 12

    # indexing
    @test s[1] == 11.0
    @test s[12] == 22.0
    @test s[1:3] == [11.0, 12.0, 13.0]
    @test s[1:2:12] == collect(10.0 .+ (1:2:12))
    @test s[s .< 13] == [11.0, 12.0]
    #
    @test s[20Q1] == 11.0
    @test s[begin] == s.values[1]
    @test s[end] == s.values[end]
    @test s[begin:begin + 3] isa typeof(s)
    @test s[begin:begin + 3].values == s.values[begin:begin + 3]
    @test (@. 13 < s < 16 ) isa TSeries{frequencyof(s),Bool}
    #
    @test_throws ArgumentError s[2:end]  # can't mix Int indexing with begin/end
    @test_throws ArgumentError s[begin:4]
    #
    @test_throws ArgumentError s[1U]  # wrong frequency
    @test_throws ArgumentError s[1Y:3Y]  # wrong frequency
    @test_throws ArgumentError s[2Y] = 5  # wrong frequency

    q = copy(s)
    s[19Q1] = 5  # outside range
    @test s.values ≈ [5, NaN, NaN, NaN, q.values...] nans=true
    s[end + 3] = 3  # outside range
    @test s.values ≈ [5, NaN, NaN, NaN, q.values..., NaN, NaN, 3] nans=true

    @test_throws ArgumentError s[20Y:21Y] = [2, 3]  # wrong frequency

    i = TSeries(20Y, ones(Int32, 5))
    i[17Y] = -1
    @test i.values == [-1, typenan(Int32), typenan(Int32), ones(5)...]

    @test_throws ArgumentError resize!(i, 17U:24U)  # wrong frequency
    @test_throws ArgumentError copyto!(i, 17U:24U, i)  # wrong frequency
    @test_throws ArgumentError copyto!(i, s)  # wrong frequency

    # rangeof with drop
    let myts = TSeries(20Q1:21Q4,1)
        @test rangeof(myts,drop= 2) == 20Q3:21Q4
        @test rangeof(myts,drop=-2) == 20Q1:21Q2
    end

end

@testset "Bcast" begin
    t = TSeries(5U:10U, rand(6))

    # we can broadcast with a singleton
    r = t .+ 5
    @test typeof(r) == typeof(t) && eachindex(r) == eachindex(t) && all(r.values .== t.values .+ 5)

    # we can broadcast with another TSeries of identical range
    r = t .+ TSeries(5U, collect(1:6))
    @test typeof(r) == typeof(t) && eachindex(r) == eachindex(t) && all(r.values .== t.values .+ (1:6))

    # we can broadcast with another TSeries of different range
    r = t .+ TSeries(4U, collect(1:6))
    @test typeof(r) == typeof(t) && eachindex(r) == 5U:9U && all(r.values .== t.values[1:end - 1] .+ (2:6))

    # we can broadcast with a Vector
    r = t .+ collect(1:6)
    @test typeof(r) == typeof(t) && eachindex(r) == eachindex(t) && all(r.values .== t.values .+ (1:6))

    # broadcast with a vector of the wrong dimension throws a DimensionMismatch
    @test_throws DimensionMismatch t .+ collect(1:4)
    @test_throws DimensionMismatch t .+ collect(1:8)

    # broadcast with a TSeries of another frequency throws a mixed_freq_error
    @test_throws ArgumentError t .+ TSeries(2020Q1 .+ (0:5), 0.75)

    # we can .= correctly
    t .= 1.0
    @test all(t.values .== 1.0)

    t .+= collect(1:6)
    @test eachindex(t) == (5U:10U) && all(t.values .== 1.0 .+ (1:6))

    s = similar(t, 3U:12U)
    fill!(s, 0.0)
    s .= t .- 2
    @test eachindex(s) === (3U:12U) && all(s[5U:10U].values .== (1:6) .- 1.0)
    @test all(s[3U:4U].values .== 0.0) && all(s[11U:end] .== 0.0)

    # dot-assign when the rhs is a vector
    @test_throws DimensionMismatch s .= ones(length(s) + 1)
    @test_throws DimensionMismatch s .= ones(length(s) - 1)
    @test (s .= ones(length(s)); rangeof(s) === 3U:12U && all(s.values .== 1.0))

    # dot-assign into a range
    t[begin + 2:end - 2] .= 2
    @test t.values == [2, 3, 2, 2, 6, 7]

    t[end + 2:end + 4] .= 8
    @test t.values ≈ [2, 3, 2, 2, 6, 7, NaN, 8, 8, 8] nans = true

    # setindex with BitArray and broadcasting
    t[t .< 7] .= 0
    @test t.values ≈ [0, 0, 0, 0, 0, 7, NaN, 8, 8, 8] nans = true

    t[2:4] .= 1
    @test t.values ≈ [0, 1, 1, 1, 0, 7, NaN, 8, 8, 8] nans = true

end

ts_u = TSeries(5)
ts_v = TSeries(3:5)
ts_m = TSeries(mm(2018, 1), collect(1.0:12.0))
ts_q = TSeries(qq(2018, 1):qq(2020, 4), collect(1:12))
ts_y = TSeries(yy(2018), collect(1:12))

@testset "TSeries 1" begin

    @test ts_u.firstdate == 1U
    # @test ts_u.values == 1:5

    @test ts_m.firstdate == mm(2018, 1)
    @test ts_m.values == collect(1.0:12.0)

    @test ts_q.firstdate == qq(2018, 1)
    @test ts_q.values == collect(1.0:12.0)

    @test ts_y.firstdate == yy(2018)
    @test ts_y.values == collect(1.0:12.0)

    # Make sure if lengths are different we get an error
    @test_throws ArgumentError TSeries(1U:5U, 1:6)

    let t = TSeries(2U, rand(5))
        @test firstdate(t) == 2U
        @test lastdate(t) == 6U
        @test length(t.values) == 5
        @test length(t) == 5
    end

    let t = TSeries(1991Q1:1992Q4), s = TSeries(2222Y:2225Y, undef), r = TSeries(1006M3:1009M5, 0.3), e = TSeries(2000Y:1995Y, 7)
        @test isempty(e)
        @test length(r) == 10 + 12 + 12 + 5
        @test length(t) == length(t.values) == 8
        @test firstindex(t) == firstdate(t) == 1991Q1
        @test lastindex(t) == lastdate(t) == 1992Q4
    end

end

@testset "Int indexing" begin
    let t = TSeries(4U:8U, rand(5))
        @test t.firstdate == 4U && lastdate(t) == 8U
        # test access
        @test t[1] isa Number
        @test t[1] == t.values[1]
        @test t[2:4] isa Vector{Float64} # TSeries{frequencyof(t),Float64,Vector{Float64}}
        @test t[2:4] == t.values[2:4]
        @test t[[1,3,4]] isa Vector{Float64}
        @test t[[1,3,4]] == t.values[[1,3,4]]
        # test assignment
        @test begin
            t[2:4] .= 2.5
            t.values[2:4] == fill(2.5, 3)
        end
        @test 5 == (t[3] = 5)
        @test t.values == [first(t), 2.5, 5.0, 2.5, last(t)]
    end
end

@testset "Views" begin
    let t = TSeries(2010M1, rand(20))
        @test axes(t) == (2010M1 - 1 .+ (1:20),)
        @test Base.axes1(t) == 2010M1 - 1 .+ (1:20)
        z = similar(t)
        @test z isa typeof(t)
        @test z.firstdate == t.firstdate
        @test z != t
        z = copy(t)
        @test z == t
        z[begin .+ (2:4)] .+= 0.2
        z[begin .+ (3:4)] = [3,4]
        @test z != t
        z = view(t, 2:5)
        c = view(t, 2010M2:2010M5)
@test c == z
        @test z == t[begin .+ (1:4)]
        z[[1,3]] += [0.5, 0.5]
        z[[2,4]] = [1,1.5]
        @test c == z
        @test z == t[begin .+ (1:4)]
    end
end

@testset "show" begin
    for (nrow, fd) = zip([3, 4, 5, 6, 7, 8, 22, 23, 24, 25, 26, 30],
        Iterators.cycle((qq(2010, 1), mm(2010, 1), yy(2010), 1U)))
        let io = IOBuffer()
            t = TSeries(fd, rand(24))
            show(IOContext(io, :displaysize => (nrow, 80)), MIME"text/plain"(), t)
            @test length(readlines(seek(io, 0))) == max(2, min(length(t) + 1, nrow - 3))
        end
    end
    for fd = (2020Q1, 2020M1, 2020Y, 2020U)
        let io = IOBuffer()
            t = TSeries(fd, Float64[])
            show(IOContext(io, :displaysize => (10, 80)), MIME"text/plain"(), t)
            @test startswith(readlines(seek(io, 0))[1], "Empty")
        end
    end
end

@testset "math" begin
    tq = TSeries(2020Q1, rand(12))
    tm = TSeries(2020M1, copy(tq.values))
    tu = TSeries(11U, copy(tq.values))

    tmp = tq .* 5
    @test 5tq isa TSeries{Quarterly}
    @test tq * 5 isa TSeries{Quarterly}
    @test 5tq == tmp
    @test tq * 5 == tmp
    @test (5tm).values == (5tq).values && 5tm ≠ 5tq

    @test (tq + tq) == (tq .+ tq)

    # approx works when ranges don't match
    sq = resize!(copy(tq), 2020Q3:2022Q2)
    @test sq ≠ tq && sq ≈ tq

    # adding two TSeries of same frequency and range works
    @test tq + 5tq ≈ 6tq
    # adding TSeries and
    @test_throws MethodError 5 + tq
    @test_throws MethodError tq + 5
    @test 5 .+ tq == tq .+ 5  # broadcasting works
    @test_throws ArgumentError tq + 5tm   # different frequencies not allowed

end

@testset "Monthly" begin
    @test ts_m[mm(2018, 1):mm(2018, 12)] == ts_m
    @test ts_m[mm(2018, 1):mm(2018, 12)].firstdate == mm(2018, 1)

    # access outside of ts boundaries
    @test_throws BoundsError ts_m[mm(2017, 1):mm(2019, 12)] == ts_m
    # @test ts_m[mm(2017, 1):mm(2019, 12)].firstdate == ts_m.firstdate

    # partially out of boundary
    @test_throws BoundsError ts_m[mm(2017, 1):mm(2018, 6)] == ts_m[mm(2018, 1):mm(2018, 6)]
    @test_throws BoundsError ts_m[mm(2017, 1):mm(2018, 6)] == ts_m[mm(2018, 1):mm(2018, 6)]

    @test_throws BoundsError ts_m[mm(2017, 1):mm(2018, 6)] == ts_m[1:6]

    @test_throws BoundsError ts_m[mm(2018, 6):mm(2019, 12)] == ts_m[mm(2018, 6):mm(2018, 12)]
    @test_throws BoundsError ts_m[mm(2018, 6):mm(2019, 12)] == ts_m[6:12]

    # fully out of boundary
    @test_throws BoundsError ts_m[mm(2017, 1)] === nothing
    @test_throws BoundsError ts_m[mm(2017, 1):mm(2017, 3)] === nothing
end

@testset "Quarterly" begin
    @test ts_q[qq(2018, 1):qq(2020, 4)] == ts_q

    # access outside of ts boundaries
    @test_throws BoundsError ts_q[qq(2017, 1):qq(2021, 4)] == ts_q

    # partially out of boundary
    @test_throws BoundsError ts_q[qq(2017, 1):qq(2018, 4)] == ts_q[qq(2018, 1):qq(2018, 4)]
    @test_throws BoundsError ts_q[qq(2017, 1):qq(2018, 4)] == ts_q[1:4]

    @test_throws BoundsError ts_q[qq(2018, 4):qq(2021, 4)] == ts_q[qq(2018, 4):qq(2020, 4)]
    @test_throws BoundsError ts_q[qq(2018, 4):qq(2021, 4)] == ts_q[4:12]

    # fully out of boundary
    @test_throws BoundsError ts_q[qq(2017, 1)] === nothing
    @test_throws BoundsError ts_q[qq(2017, 1):qq(2017, 3)] === nothing
end

@testset "Yearly" begin
    @test ts_y[yy(2018):yy(2029)] == ts_y

    # access outside of ts boundaries
    @test_throws BoundsError ts_y[yy(2017):yy(2017) + 100] == ts_y

    # partially out of boundary
    @test_throws BoundsError ts_y[yy(2017):yy(2018)] == ts_y[yy(2018):yy(2018)]
    @test_throws BoundsError ts_y[yy(2017):yy(2021)] == ts_y[1:4]

    @test_throws BoundsError ts_y[yy(2018):yy(2100)] == ts_y[yy(2018):yy(2029)]
    @test_throws BoundsError ts_y[yy(2021):yy(2100)] == ts_y[4:12]

    # fully out of boundary
    @test_throws BoundsError ts_y[yy(2017)] === nothing
    @test_throws BoundsError ts_y[yy(2010):yy(2017)] === nothing
end

# ts_m = TSeries(mm(2018, 1), collect(1.0:12.0))
@testset "Setting" begin
    begin
        ts_m = TSeries(mm(2018, 1):mm(2018, 12), collect(1.0:12.0))
        @test_throws BoundsError ts_m[mm(2019, 2):mm(2019, 4)] = 1;
        ts_m[mm(2019, 2):mm(2019, 4)] .= 1;
        @test ts_m[mm(2019, 2):mm(2019, 4)].values == [1, 1, 1]
        @test ts_m.firstdate == mm(2018, 1)
        @test isequal(ts_m.values, vcat(collect(1.0:12.0), [NaN], [1, 1, 1]))
    end

    begin
        ts_m = TSeries(mm(2018, 1):mm(2018, 12), collect(1.0:12.0))
        @test_throws BoundsError ts_m[mm(2017, 10):mm(2017, 11)] = 1;
        ts_m[mm(2017, 10):mm(2017, 11)] .= 1;
        @test ts_m[mm(2017, 10):mm(2017, 11)].values == [1, 1]
        @test ts_m.firstdate == mm(2017, 10)
        @test isequal(ts_m.values, vcat([1, 1], [NaN], collect(1.0:12.0)))
    end

    begin
        ts_m = TSeries(mm(2018, 1):mm(2018, 12), collect(1.0:12.0))
        ts_m[mm(2019, 2):mm(2019, 4)] = [9, 10, 11];
        @test ts_m[mm(2019, 2):mm(2019, 4)].values == [9, 10, 11]
        @test ts_m.firstdate == mm(2018, 1)
@test isequal(ts_m.values, vcat(collect(1.0:12.0), [NaN], [9, 10, 11]))
    end

    begin
        ts_m = TSeries(mm(2018, 1):mm(2018, 12), collect(1.0:12.0))
        ts_m[mm(2017, 10):mm(2017, 11)] = [9, 10];
        @test ts_m[mm(2017, 10):mm(2017, 11)].values == [9, 10]
        @test ts_m.firstdate == mm(2017, 10)
        @test isequal(ts_m.values, vcat([9, 10], [NaN], collect(1.0:12.0)))
    end
end

@testset "Addition" begin
    x = TSeries(1U, [7, 7, 7])
    y = TSeries(3U, [2, 4, 5])
    @test x + y == TSeries(3U, [9])

    x = TSeries(1U, [7, 7, 7])
    y = TSeries(2U, [2, 4, 5])
    @test x + y == TSeries(2U, [9, 11])
end

@testset "Iris" begin
    # IRIS based assignment of values from other TSeries
    x = TSeries(qq(2020, 1), zeros(3));
    y = TSeries(qq(2020, 1), ones(3));
    x[qq(2020, 1):qq(2020, 2)] = y;
    @test x == TSeries(qq(2020, 1), [1, 1, 0])

    let t = TSeries(5U:10U, ones)
        s = t[6U:8U]
        @test rangeof(s) === 6U:8U
        s .= 0.7
        t[6U:8U] = s
        @test t.values == [1, 0.7, 0.7, 0.7, 1, 1]

        s .= 0.8
        t[6U:7U] = s
        @test t.values == [1, 0.8, 0.8, 0.7, 1, 1]

        t[6U:7U] = [2,3]
        @test t.values == [1, 2, 3, 0.7, 1, 1]

        @test_throws ArgumentError t[6Y:7Y] = s  # mixed frequency in indexing range
        @test_throws ArgumentError t[6U:7U] = TSeries(6Y, s.values)  # mixed frequency of src and dest
    end

        # IRIS related: shift
    x = TSeries(qq(2020, 1), zeros(3));
    @test shift(x, 1) == TSeries(qq(2019, 4), zeros(3))

    shift!(x, 1)
    @test x == TSeries(qq(2019, 4), zeros(3))

    # # IRIS related: nanrm!
    # x = TSeries(qq(2020, 1), [NaN, 123, NaN]);
    # nanrm!(x)
    # @test x == TSeries(qq(2020, 2), [123])

    # TODO
    # - pct
    # - apct

end

@testset "TS.math" begin
    let x = TSeries(2000Y:2010Y, ones)
        # lags and leads
        y = cumsum(x)
        @test rangeof(y) === rangeof(x)
        z = y
        @test z === y
        y1 = lag(y)
        @test rangeof(y1) === 1 .+ (rangeof(y))  && y1.values == y.values
        lag!(y)
        @test y1 == y && y1 !== y
        @test z === y
        y2 = lead(y, 2)
        @test rangeof(y2) === -2 .+ (rangeof(y)) && y2.values == y.values
        lead!(y, 3)
        @test lead(y2) == y && lead(y2) !== y
        @test z === y
        # opeartions
        @test x + x + 3x == 5x
        @test_throws ArgumentError x + x.values
        @test_throws ArgumentError x.values + x
    end
end

@testset "axes.range" begin
    @test axes(1U:5U) == axes(1:5)
    @test Base.axes1(2020Y:2030Y) == Base.OneTo(11)
end

@testset "overlay" begin
    # from FAME manual
    A, B = (TSeries(87Y, [1, 2, NaN, 4]), TSeries(87Y, [NaN, 6, 7, 8]))
    @test overlay(A, B) == TSeries(87Y, [1,2,7,4])
    @test overlay(B, A) == TSeries(87Y, [1,6,7,8])
    @test overlay(86Y:92Y, A, B) ≈ TSeries(86Y, [NaN, 1, 2, 7, 4, NaN, NaN]) nans = true
    @test (C = overlay(A,B); overlay(C,A).values == C.values)
end

@testset "fconvert" begin
    t = TSeries(5U, collect(1:10))
    @test fconvert(Unit, t) === t
    @test_throws ErrorException fconvert(Quarterly, t)

    q = TSeries(5Q1, 1.0collect(1:10))
    @test_throws ErrorException  fconvert(Unit, q)
    mq = fconvert(Monthly, q)
    @test typeof(mq) === TSeries{Monthly, Float64, Vector{Float64}}
    yq = fconvert(Yearly, q)
    @test typeof(yq) === TSeries{Yearly, Float64, Vector{Float64}}

    for i = 1:11
        @test rangeof(fconvert(Yearly, TSeries(1M1 .+ (i:50)))) == 2Y:4Y
        @test rangeof(fconvert(Yearly, TSeries(1M1 .+ (0:47+i)))) == 1Y:4Y
    end
    for i = 1:3
        @test rangeof(fconvert(Yearly, TSeries(1Q1 .+ (i:50)))) == 2Y:12Y
        @test rangeof(fconvert(Yearly, TSeries(1Q1 .+ (0:47+i)))) == 1Y:12Y
    end
    for i = 1:11
        @test rangeof(fconvert(Quarterly, TSeries(1M1 .+ (i:50)))) == 1Q2+div(i-1,3):5Q1
        # @test rangeof(fconvert(Quarterly, TSeries(1M1 .+ (0:47+i)))) == 1Y:4Y
    end

end

@testset "strip" begin
    let rng_x = 2000Y:2010Y, x = TSeries(rng_x, ones)
        x[2011Y:2015Y] .= NaN
        x[1995Y:1999Y] .= NaN
        @test (rangeof(strip(x)) == rng_x)
        @test (TimeSeriesEcon.strip!(x); rangeof(x) == rng_x)
    end
end

@testset "recursive" begin
    ts = TSeries(1U, zeros(0))
    ts[1U] = ts[2U] = 1.0
    @rec 3U:10U ts[t] = ts[t-1]+ts[t-2]
            @test ts.values == [1.0,1,2,3,5,8,13,21,34,55]
    t = zeros(10,7)
    r = rand(1, 7)
    t[1, :] = r
    @rec s=2:10 t[s,:] = t[s-1,:] .* s
    @test t ≈ factorial.(1:10) * r
    #
    s = ones(15);
    @rec i=3:15 s[i] = s[i-1] + s[i-2]
    @test s == Float64[1,1,2,3,5,8,13,21,34,55,89,144,233,377,610]
    s = TSeries(2020Q1:2021Q4)
    s[begin] = 0
    @test_throws UndefVarError @rec 2020Q1:2021Q2 s[i+1] = 1.0 + s[i]
    @rec i=firstdate(s):2023Q2 s[i+1] = 1.0 + s[i]
    @test rangeof(s) == 2020Q1:2023Q3
    @test values(s) == Float64[0:14...]
    resize!(s, 2020Q1:2020Q2)
    @rec firstdate(s)+1:2022Q3 s[t+1] = 2.0 * s[t] + s[t-1]
    @test rangeof(s) == 2020Q1:2022Q4
    @test values(s) == Float64[0,1,2,5,12,29,70,169,408,985,2378,5741]
end


text:
module ValueHistoriesLogger

# Imports
using Base.CoreLogging:
    global_logger, LogLevel, Info

import Base.CoreLogging:
    AbstractLogger, handle_message, shouldlog, min_enabled_level,
	catch_exceptions

using ValueHistories

# Exports
export MVLogger, data

# Core Code
"""
	MVLogger <: AbstractLogger

A `Logger` that will save any variable logged to a `MVHistory` object from
`ValueHistories` package.

Any message such as
`@info \"text\" key1=val1 key2=val2`
will result into the keys `text/key1` and `text/key2` with the corresponding
data being logged to the history.

The special key `log_step_increment=1` is always included. If redefined then
this represents by how much the step counter should be incremented.
"""
mutable struct MVLogger <: AbstractLogger
    hist::MVHistory
    min_level::LogLevel
    global_step::Int
end

"""
	MVLogger(level=Info, step=0)

Creates a Multi Value Logger `MVLogger` that will log any message above level
`LogLevel.Info` and current `step` set to 0.
"""
MVLogger(level=Info; step=0)  = MVLogger(MVHistory(), level, step)

"""
	data(logger)

Returns the `MVHistory` object to which the logger logged data.
"""
data(logger::MVLogger) = logger.hist

"""
	set_step(logger, iter)

Sets the step of `logger` to `iter`.
"""
set_step(lg::MVLogger, iter::Int) = lg.global_step = iter

"""
	increment_step(logger, incr)

Increments the step of `logger` by `incr` and returns the new step value.
"""
increment_step(lg::MVLogger, iter::Int) = lg.global_step += iter

"""
	step(logger)

Returns the current step.
"""
step(lg::MVLogger) = lg.global_step

shouldlog(lg::MVLogger, level, _module, group, id) = true

min_enabled_level(lg::MVLogger) = lg.min_level

catch_exceptions(lg::MVLogger) = false

function handle_message(lg::MVLogger, level, message, _module, group, id,
                        filepath, line; maxlog=nothing, kwargs...)
    i_step = 1
	#=
    for (key, val) in pairs(kwargs)
        if key == :log_step_increment
            i_step = val
            break
        end
    end
    iter = increment_step(lg, i_step)

	!isempty(message) ? message = message*"/" : message = message

    for (key, val) in pairs(kwargs)
        key == :log_step_increment && continue

        push!(lg.hist, Symbol(message, key), iter, deepcopy(val))
    end
	=#
	!isempty(message) ? message = message*"/" : message = message
	# new
	if !isempty(kwargs)
		data = Vector{Pair{String,Any}}()

		# ∀ (k-v) pairs, decompose values into objects that can be serialized
		for (key,val) in pairs(kwargs)
			# special key describing step increment
			if key == :log_step_increment
				i_step = val
				continue
			end

			preprocess(message*"$key", val, data)
		end
		iter = increment_step(lg, i_step)

		# Serialize every object
		for (name,val) in data
			push!(lg.hist, Symbol(name), iter, deepcopy(val))
		end
	end
end

function preprocess(name, val::T, data) where T
    if isstructtype(T)
        fn = fieldnames(T)
        for f=fn
            prop = getfield(val, f)
            preprocess(name*"/$f", prop, data)
        end
	else
		push!(data, name=>val)
    end
	data
end

function preprocess(name, val::Union{AbstractArray, Complex}, data)
	push!(data, name=>val)
end

ValueHistories.getindex(h::MVHistory, s::String) = h[Symbol(s)]

end # module