text:
function test_lbfgs()
  ϵ = eps(Float64)
  rtol = sqrt(ϵ)

  # test limited-memory BFGS
  @testset ExtendedTestSet "LBFGS" begin
    n = 10
    mem = 5
    B = LBFGSOperator(n, mem = mem, scaling = false)
    H = InverseLBFGSOperator(n, mem = mem, scaling = false)

    for t = 1:2 # Run again after reset!
      @test norm(diag(B) - diag(Matrix(B))) <= rtol

      @test B.data.insert == 1
      @test H.data.insert == 1
      @test norm(Matrix(B) - Matrix(1.0I, n, n)) <= ϵ
      @test norm(Matrix(H) - Matrix(1.0I, n, n)) <= ϵ

      # Test that nonpositive curvature can't be added.
      s = simple_vector(Float64, n)
      z = zeros(n)
      push!(B, s, -s)
      @test B.data.insert == 1
      push!(B, s, z)
      @test B.data.insert == 1
      push!(H, s, -s)
      @test H.data.insert == 1
      push!(H, s, z)
      @test H.data.insert == 1

      # Insert a few {s,y} pairs.
      insert = 0
      for i = 1:(mem + 2)
        s = ones(n) * i
        y = [i; ones(n - 1)]
        if dot(s, y) > 1.0e-20
          insert += 1
          push!(B, s, y)
          push!(H, s, y)
        end
      end

      @test B.data.insert == mod(insert, B.data.mem) + 1
      @test H.data.insert == mod(insert, H.data.mem) + 1

      @test check_positive_definite(B)
      @test check_positive_definite(H)

      @test check_hermitian(B)
      @test check_hermitian(H)

      @test norm(diag(B) - diag(Matrix(B))) <= rtol

      @test norm(Matrix(H * B) - Matrix(1.0I, n, n)) <= rtol

      # Testing reset! function
      v = simple_vector(Float64, n)
      @test norm(B * v - v) > rtol
      @test norm(H * v - v) > rtol
      reset!(B)
      reset!(H)
      @test B.data.scaling_factor == 1.0
      @test H.data.scaling_factor == 1.0
      @test norm(B * v - v) < rtol
      @test norm(H * v - v) < rtol
    end

    # test against full BFGS without scaling
    mem = n
    LB = LBFGSOperator(n, mem = mem, scaling = false)
    B = Matrix(1.0I, n, n)

    function bfgs!(B, s, y, damped = false)
      # dense BFGS update
      ys = dot(y, s)
      Bs = B * s
      tol = damped ? (0.2 * dot(s, Bs)) : 1.0e-20
      if ys > tol
        B = B - Bs * Bs' / dot(s, Bs) + y * y' / ys
      end
      return B
    end

    @test norm(Matrix(LB) - B) < rtol * norm(B)
    @test norm(diag(LB) - diag(B)) < rtol * norm(diag(B))

    for k = 1:mem
      s = simple_vector(Float64, n)
      y = simple_vector(Float64, n)
      B = bfgs!(B, s, y)
      LB = push!(LB, s, y)
      @test norm(Matrix(LB) - B) < rtol * norm(B)
      @test norm(diag(LB) - diag(B)) < rtol * norm(diag(B))
    end

    # test damped L-BFGS
    B = LBFGSOperator(n, mem = mem, damped = true, scaling = false, σ₂ = 0.8, σ₃ = Inf)
    H = InverseLBFGSOperator(n, mem = mem, damped = true, scaling = false, σ₂ = 0.8, σ₃ = Inf)

    insert_B = insert_H = 0
    for i = 1:(mem + 2)
      s = simple_vector(Float64, n)
      y = simple_vector(Float64, n)
      ys = dot(y, s)
      g = simple_vector(Float64, n)
      d = -(H * g) # added brackets here
      α = i / mem
      s = α * d
      if ys > 0.2 * dot(s, B * s)
        insert_B += 1
        insert_H += 1
        push!(B, s, y)
        push!(H, s, y, α, g)
      end
    end

    @test B.data.insert == mod(insert_B, B.data.mem) + 1
    @test H.data.insert == mod(insert_H, H.data.mem) + 1

    @test check_positive_definite(B)
    @test check_positive_definite(H)

    @test check_hermitian(B)
    @test check_hermitian(H)

    @test norm(diag(B) - diag(Matrix(B))) <= rtol

    @test norm(Matrix(H * B) - Matrix(1.0I, n, n)) <= rtol

    # test against full BFGS without scaling
    mem = n
    LB = LBFGSOperator(n, mem = mem, damped = true, scaling = false)
    B = Matrix(1.0I, n, n)

    @test norm(Matrix(LB) - B) < rtol * norm(B)
    @test norm(diag(LB) - diag(B)) < rtol * norm(diag(B))

    for k = 1:mem
      s = simple_vector(Float64, n)
      y = simple_vector(Float64, n)
      B = bfgs!(B, s, y, true)
      LB = push!(LB, s, y)
      @test norm(Matrix(LB) - B) < rtol * norm(B)
      @test norm(diag(LB) - diag(B)) < rtol * norm(diag(B))
    end
  end

  @testset ExtendedTestSet "Different precision" begin
    n = 10
    mem = 5
    for T in (Float16, Float32, Float64, BigFloat)
      B = LBFGSOperator(T, n, mem = mem)
      H = InverseLBFGSOperator(T, n, mem = mem)
      s = ones(T, n)
      y = ones(T, n)
      push!(B, s, y)
      push!(H, s, y)
      @test eltype(B) == T
      @test eltype(H) == T
      v = [-(-one(T))^i for i = 1:n]
      @test eltype(B * v) == T
      @test eltype(H * v) == T
    end
  end

  @testset "LBFGS allocations" begin
    n = 100
    mem = 20
    B = LBFGSOperator(n, mem = mem)
    H = InverseLBFGSOperator(n, mem = mem)
    for _ = 1:2:n
      s = rand(n)
      y = rand(n)
      push!(B, s, y)
      push!(H, s, y)
    end
    x = rand(n)
    res = similar(x)
    mul!(res, B, x)  # warmup
    nallocs = @allocated mul!(res, B, x)
    @test nallocs == 0
    mul!(res, H, x)  # warmup
    nallocs = @allocated mul!(res, H, x)
    @test nallocs == 0
    nallocs = @allocated diag!(B, x)
    @test nallocs == 0
  end

  @testset "LBFGS eigenvalues" begin
    n = 100
    mem = 20
    B = LBFGSOperator(n, mem = mem)
    H = InverseLBFGSOperator(n, mem = mem)
    for _ = 1:2:n
      s = rand(n)
      y = rand(n)
      push!(B, s, y)
      push!(H, s, y)
    end
    λ = eigs(B, nev = n)[1]  # call Arpack
    @test minimum(λ) > 0
    λ = eigs(H)[1]
    @test minimum(λ) > 0
  end
end

test_lbfgs()


text:
# # [How to port NiLang to ChainRules](@id port_chainrules)
#
# In [How to port NiLang to Zygote](@ref port_zygote) we showed the way to insert Nilang-based
# gradient as Zygote's pullback/adjoint. Given that [ChainRules](https://github.com/JuliaDiff/ChainRules.jl)
# is now the core of many AD packages including Zygote, extending `ChainRules.rrule` with Nilang
# does the same job, except that it affects all ChainRules-based AD packages and not just Zygote.
#
# We'll use the same example as [How to port NiLang to Zygote](@ref port_zygote), so you might need
# to restart your Julia to get a fresh environment.

using NiLang, NiLang.AD, Zygote, ChainRules

# Let's start from the Julia native implementation of `norm2` function.
function norm2(x::AbstractArray{T}) where T
    out = zero(T)
    for i=1:length(x)
        @inbounds out += x[i]^2
    end
    return out
end

# Zygote is able to generate correct dual function, i.e., gradients, but much slower than the primal
# function `norm2`
using BenchmarkTools
x = randn(1000);
original_grad = norm2'(x)
@benchmark norm2'($x) seconds=1

# The primal function is
@benchmark norm2($x) seconds=1

# Then we have the reversible implementation
@i function r_norm2(out::T, x::AbstractArray{T}) where T
    for i=1:length(x)
        @inbounds out += x[i]^2
    end
end

# The gradient generated by NiLang is much faster, which is comparable to the forward program
@benchmark (~r_norm2)(GVar($(norm2(x)), 1.0), $(GVar(x))) seconds=1

# By defining our custom `rrule` using Nilang's gradient implementation, `Zygote` automaticallly
# gets boosted because it internally uses the available ChainRules ruleset.
# Here we need to create a new symbol here because otherwise Zygote will still use the
# previously generated slow implementation.
norm2_faster(x) = norm2(x)
function ChainRules.rrule(::typeof(norm2_faster), x::AbstractArray{T}) where T
    out = norm2_faster(x)
    function pullback(ȳ)
        ChainRules.NoTangent(), grad((~r_norm2)(GVar(out, ȳ), GVar(x))[2])
    end
    out, pullback
end
@assert norm2_faster'(x) ≈ original_grad

# See, much faster
@benchmark norm2_faster'(x) seconds=1

ext: jl


text:
module KoalaElasticNet

# new:
export ElasticNetRegressor

# needed in this module:
import Koala: Regressor, BaseType, keys_ordered_by_values, SupervisedMachine, params
import Koala: Transformer
import KoalaTransforms
import Lasso
import DataFrames: AbstractDataFrame, DataFrame
import UnicodePlots

# to be extended (but not explicitly rexported):
import Koala: setup, fit, predict
import Koala: default_transformer_X, default_transformer_y, transform, inverse_transform

# development only:
# using ADBUtilities

## Helpers

"""
`getcol(A, j)` gets the jth column of a sparse matrix `A`, returned as
a sparse vector. If `j=0` then the last column is returned. If `j`
exceeds the number of columns of `A`, then a zero vector is returned.

"""
function getcol(A::SparseMatrixCSC{Float64,Int}, j::Int)
    if j == 0
        j = A.n
    end
    j > A.n ? (return SparseVector(A.m,Int[],Float64[])) : nothing
    I = Int[]     # indices for non-zero entries of col
    V = Float64[] # values of non-zero entries of col
    for k in A.colptr[j]:(A.colptr[j+1] - 1)
        push!(I, A.rowval[k])
        push!(V, A.nzval[k])
    end
    return SparseVector(A.m,I,V)
end

# function to determine the largest lambda value conceivable for given
# alpha and training data:
function lambda_max(X, y, alpha)
    Xy = X'*(y .- mean(y)) # dot products of feature columns with centred response
    λmax = abs(Xy[1])
    for i = 2:length(Xy)
        x = abs(Xy[i])
        if x > λmax
            λmax = x
        end
    end
    return λmax/(alpha*length(y))
end

## Model type definitions

struct LinearPredictor <: BaseType
    intercept::Float64
    coefs::SparseVector{Float64,Int64}
end

# Following returns a `DataFrame` with three columns:
#
# column name | description
# :-----------|:-------------------------------------------------
# `:index`    | index of a feature used to train `predictor`
# `:feature`  | corresponding feature label provided by `features`
# `:coef`     | coefficient for that feature in the predictor
#
# The rows are ordered by the absolute value of the coefficients. If
# `rgs` is unfitted, an error is returned.
function coef_info(predictor::LinearPredictor, features)
    coef_given_index = Dict{Int, Float64}()
    abs_coef_given_index = Dict{Int, Float64}()
    v = predictor.coefs # SparseVector
    for k in eachindex(v.nzval)
        coef_given_index[v.nzind[k]] = v.nzval[k]
        abs_coef_given_index[v.nzind[k]] = abs(v.nzval[k])
    end
    df = DataFrame()
    df[:index] = reverse(keys_ordered_by_values(abs_coef_given_index))
    df[:feature] = map(df[:index]) do index
        features[index]
    end
    df[:coef] = map(df[:index]) do index
        coef_given_index[index]
    end
    return df
end

mutable struct ElasticNetRegressor <: Regressor{LinearPredictor}

    lambda::Float64
    n_lambdas::Int
    alpha::Float64
    max_n_coefs::Int
    criterion::Symbol
    lambda_min_ratio::Float64

    # for cv-optimizing lambda:
    n_folds::Int

    function ElasticNetRegressor(lambda, n_lambdas::Int, alpha,
                                 max_n_coefs::Int, criterion::Symbol,
                                 lambda_min_ratio, n_folds::Int)
        if alpha <= 0.0 || alpha > 1.0
            if alpha == 0.0
                throw(Base.error("alpha=0 dissallowed."*
                                 " Consider ridge regression instead."))
            end
            throw(DomainError)
        end
        return new(lambda, n_lambdas, alpha, max_n_coefs, criterion,
                  lambda_min_ratio, n_folds)
    end

end

# lazy keywork constructor
ElasticNetRegressor(;lambda=0.0, n_lambdas::Int=100,
                    alpha=1.0,
                    max_n_coefs::Int=0, criterion::Symbol=:coef,
                    lambda_min_ratio=0.0, n_folds::Int=9) =
                        ElasticNetRegressor(lambda, n_lambdas,
                                            alpha,
                                            max_n_coefs, criterion,
                                            lambda_min_ratio, n_folds)

# `showall` method for `ElasticNetRegressor` machines:
function Base.showall(stream::IO,
                      mach::SupervisedMachine{LinearPredictor, ElasticNetRegressor})
    show(stream, mach)
    println(stream)
    if isdefined(mach,:report) && :feature_importance_curve in keys(mach.report)
        features, importance = mach.report[:feature_importance_curve]
        plt = UnicodePlots.barplot(features, importance,
              title="Feature importance (coefs of linear predictor)")
    end
    dict = params(mach)
    report_items = sort(collect(keys(dict[:report])))
    dict[:report] = "Dict with keys: $report_items"
    dict[:Xt] = string(typeof(mach.Xt), " of shape ", size(mach.Xt))
    dict[:yt] = string(typeof(mach.yt), " of shape ", size(mach.yt))
    delete!(dict, :cache)
    showall(stream, dict)
    println(stream, "\nModel detail:")
    showall(stream, mach.model)
    if isdefined(mach,:report) && :feature_importance_curve in keys(mach.report)
        show(stream, plt)
    end
end

default_transformer_X(model::ElasticNetRegressor) =
    KoalaTransforms.DataFrameToArrayTransformer()
default_transformer_y(model::ElasticNetRegressor) =
    KoalaTransforms.RegressionTargetTransformer()

# Note: For readability we now use `X` and `y` in place of `Xt` and `yt`

struct Cache <: BaseType
    X::Matrix{Float64}
    y::Vector{Float64}
    features::Vector{Symbol}
end

setup(model::ElasticNetRegressor, X, y, scheme_X, parallel, verbosity) =
    Cache(X, y, scheme_X.spawned_features)

# The core elastic net algorithm does not fit to model parameters but
# fits a *sequence* (aka path) of predictors corresponding to a
# specified sequence of `lambda` values (all other model parameters
# fixed). The reason is this: for high `lambda` values, the optimal
# predictor is a trivial (constant) linear predictor, and it is
# efficient to compute predictors for slightly lower values as
# perturbations to this model, and so on as `lambda` is slowly
# decreased.

# The following core function returns the optimal predictor (really a
# *path* of predictors) for a specified path of lambda values,
# `lambdas`. If `lambdas` is empty, it uses an automatically generated
# path instead.

function fit_path(model::ElasticNetRegressor, cache, rows, lambdas)

    X = cache.X[rows,:]
    y = cache.y[rows]

    if isempty(lambdas)
        if model.max_n_coefs == 0
            if model.lambda_min_ratio == 0.0
                predictors =
                    Lasso.fit(Lasso.LassoPath, X, y;
                              nλ=model.n_lambdas, α=model.alpha,
                              standardize=false, # in Koala any standardization is external to fit
                              criterion=model.criterion, verbose=false)
            else
                predictors =
                    Lasso.fit(Lasso.LassoPath, X, y;
                              nλ=model.n_lambdas, α=model.alpha,
                              standardize=false, # in Koala any standardization is external to fit
                              criterion=model.criterion, verbose=false,
                              λminratio=model.lambda_min_ratio)
            end
        else
            if model.lambda_min_ratio == 0.0
                predictors =
                    Lasso.fit(Lasso.LassoPath, X, y;
                              nλ=model.n_lambdas, α=model.alpha,
                              standardize=false, # in Koala any standardization is external to fit
                              maxncoef=model.max_n_coefs, criterion=model.criterion,
                              verbose=false)
            else
                predictors =
                    Lasso.fit(Lasso.LassoPath, X, y;
                              nλ=model.n_lambdas, α=model.alpha,
                              standardize=false, # in Koala any standardization is external to fit
                              maxncoef=model.max_n_coefs, criterion=model.criterion,
                              verbose=false, λminratio=model.lambda_min_ratio)
            end
        end
    else
        reverse!(sort!(lambdas))
        if model.max_n_coefs == 0
            if model.lambda_min_ratio == 0.0
                predictors =
                    Lasso.fit(Lasso.LassoPath, X, y;
                              λ =lambdas, α=model.alpha,
                              standardize=false, # in Koala any standardization is external to fit
                              criterion=model.criterion, verbose=false)
            else
                predictors =
                    Lasso.fit(Lasso.LassoPath, X, y;
                              λ =lambdas, α=model.alpha,
                              standardize=false, # in Koala any standardization is external to fit
                              criterion=model.criterion, verbose=false,
                              λminratio=model.lambda_min_ratio)
            end
        else
            if model.lambda_min_ratio == 0.0
                predictors =
                    Lasso.fit(Lasso.LassoPath, X, y;
                              λ =lambdas, α=model.alpha,
                              standardize=false, # in Koala any standardization is external to fit
                              maxncoef=model.max_n_coefs,
                              criterion=model.criterion, verbose=false)
            else
                predictors =
                    Lasso.fit(Lasso.LassoPath, X, y;
                              λ =lambdas, α=model.alpha,
                              standardize=false, # in Koala any standardization is external to fit
                              maxncoef=model.max_n_coefs,
                              criterion=model.criterion,
                              verbose=false, λminratio=model.lambda_min_ratio)
            end
        end
    end

    return predictors::Lasso.LassoPath

end

# for generating predictions on a pattern, for each lambda in the path
# used to generate `predictors`:
predict(predictors::Lasso.LassoPath, pattern) =
    vec(Lasso.predict(predictors, pattern'))

function fit(model::ElasticNetRegressor, cache, add, parallel, verbosity)

    report = Dict{Symbol, Any}()

    rows = eachindex(cache.y)

    if model.lambda == 0

        n_folds = model.n_folds

        best_lambdas = Array{Float64}(n_folds)
        errors = Array{Float64}(n_folds)
        n_samples = length(rows)

        k = floor(Int, n_samples/n_folds)
        first = 1       # first test_bag index
        second = k
        verbosity < 1 ||
            println("Optimizing regularization parameter using "*
                    "$(n_folds)-fold cross-validation. ")
        for n in 1:n_folds
            test_rows = rows[first:second]
            train_rows = vcat(rows[1:first-1], rows[second+1:end])
#            @dbg length(test_rows) length(train_rows)
            first = first + k
            second = second + k

            # get predictors along an automatically generated lambda-path:
            predictors = fit_path(model, cache, train_rows, Float64[])
            lambdas = predictors.λ
#            @dbg lambdas

            # get RMS errors on the current test set for each lambda in path:
            ss_deviations = zeros(length(lambdas))
            for i in test_rows
                ss_deviations = ss_deviations +
                    map(predict(predictors, cache.X[i,:])) do yhat
                        (yhat - cache.y[i])^2
                    end
            end

            # get lambda with minimum error:
            L = indmin(ss_deviations)
            if L == length(lambdas) && verbosity > 0
                Base.warn("Optimal lambda not reached on descent.\n"*
                          "Consider lowering lambda lambda_min_ratio.")
            end
            best_lambdas[n] = lambdas[L]
            verbosity < 1 || println("fold: $n  optimal lambda: $(lambdas[L])")
            errors[n] = sqrt(ss_deviations[L]/length(test_rows))
        end

        loglambdas = log.(best_lambdas)

        report[:loglambdaopt] = mean(loglambdas)
        report[:loglambdaopt_stde] = std(loglambdas)
        report[:lambdaopt] = exp.(report[:loglambdaopt])
        report[:cv_rmse] = mean(errors)
        report[:cv_rmse_stde] = std(errors)

        verbosity < 1 ||
            println("\nOptimal lambda = $(report[:lambdaopt]). See machine "*
                    "report for details.")

        # prepare for final train:
        lambda = report[:lambdaopt]

    else

        lambda = model.lambda

    end

    # calculate lambda-path `lambdas` for final train
    λmax = lambda_max(cache.X, cache.y, model.alpha)
    λmax >= model.lambda ||
        throw(Base.error("Something wrong here"))
    lambdas = exp.(linspace(log(λmax), log(lambda), model.n_lambdas))

    # make sure last element of lambdas is *exactly* lambda:
    lambdas[end]=lambda

    # final train on all the data:
    predictors = fit_path(model, cache, rows, lambdas)

    # check path finishes at optimal lambda:
    if predictors.λ[end] != lambda && verbosity < 1
        verbosity < 1 || Base.warn("Early stopping of path before optimal "*
                                   "or prescribed lambda reached.")
    end

    # assemble predictor:
    intercept =  predictors.b0[end]
    coefs = getcol(predictors.coefs, 0) # 0 gets last column
    predictor = LinearPredictor(intercept, coefs)

    # report on the relative strength of each feature in the predictor:
    cinfo = coef_info(predictor, cache.features) # a DataFrame object
    u = Symbol[]
    v = Float64[]
    for i in 1:size(cinfo, 1)
        feature, coef = (cinfo[i, :feature], cinfo[i, :coef])
        coef = floor(1000*coef)/1000
        if coef < 0
            label = string(feature, " (-)")
        else
            label = string(feature, " (+)")
        end
        push!(u, label)
        push!(v, abs(coef))
    end
    report[:feature_importance_curve] = u, v

    return predictor, report, cache

end

function predict(predictor::LinearPredictor, pattern::Vector{Float64})

    ret = predictor.intercept
    for i in eachindex(predictor.coefs.nzval)
        ret = ret + predictor.coefs.nzval[i]*pattern[predictor.coefs.nzind[i]]
    end

    return ret

end

predict(model::ElasticNetRegressor, predictor::LinearPredictor,
        X::Matrix{Float64}, parallel, verbosity) =
            Float64[predict(predictor, X[i,:]) for i in 1:size(X,1)]

end # of module


text:
"""
    filtering(;input::AbstractString="", featurelist::AbstractString="", samplelist::AbstractString="", thr1::Number=0.0, thr2::Number=0.0, direct1::AbstractString="+", direct2::AbstractString="+", outdir::AbstractString=".")

This function filters the genes by some standards such as mean or variance of the genes.

Input Arguments
---------
- `input` : A Julia Binary file generated by `csv2bin` function.
- `featurelist` : A row-wise summary data such as. The CSV files are generated by `csv2bin` function.
- `thr` : The threshold to reject low-signal feature.
- `outdir` : The directory specified the directory you want to save the result.

Output Files
---------
- `filtered.zst` : Filtered binary file.
"""
function filtering(;input::AbstractString="", featurelist::AbstractString="", samplelist::AbstractString="", thr1::Number=0.0, thr2::Number=0.0, direct1::AbstractString="+", direct2::AbstractString="+", outdir::AbstractString=".")
    # Feature selection
    if featurelist == "" && samplelist == ""
        error("At least featurelist or samplelist must be specified")
    end
    if featurelist != ""
        featurevec = readcsv(featurelist) # e.g., HVG_pval.csv
    else
        featurevec = Nothing
    end
    if samplelist != ""
        samplevec = readcsv(samplelist) # e.g., Sample_NoCounts.csv
    else
        samplevec = Nothing
    end

    # Setting
    if thr1 isa String
        thr1 = parse(Float64, thr1)
    end
    if thr2 isa String
        thr2 = parse(Float64, thr2)
    end
    N, M = nm(input)
    tmpN = zeros(UInt32, 1)
    tmpM = zeros(UInt32, 1)
    x = zeros(UInt32, M)
    outfile1 = joinpath(outdir, "filtered.zst")
    outfile2 = joinpath(outdir, "filteredFeature.csv")
    outfile3 = joinpath(outdir, "filteredSample.csv")

    # Search
    nr, nc, ix1, ix2 = rcfilter(input, featurevec, samplevec, thr1, thr2, direct1, direct2)
    nr = UInt32(nr)
    nc = UInt32(nc)

    # save
    writecsv(outfile2, ix1)
    writecsv(outfile3, ix2)
    open(outfile1, "w") do file1
        stream1 = ZstdCompressorStream(file1)
        write(stream1, nr) # by rcfilter
        write(stream1, nc) # by rcfilter
        open(input , "r") do file2
            stream2 = ZstdDecompressorStream(file2)
            read!(stream2, tmpN)
            read!(stream2, tmpM)
            progress = Progress(N)
            for n = 1:N
                next!(progress)
                read!(stream2, x)
                    if ix1[n] == 1
                        write(stream1, x[ix2])
                    end
            end
            close(stream2)
        end
        close(stream1)
    end
    print("\n")
end

function rcfilter(input, featurevec, samplevec, thr1, thr2, direct1, direct2)
    ncol = 0
    nrow = 0
    N, M = nm(input)
    idx1 = zeros(UInt32, N)
    idx2 = Int64[]
    tmpN = zeros(UInt32, 1)
    tmpM = zeros(UInt32, 1)
    # Row-wise filter
    if featurevec != Nothing
        open(input, "r") do file
            stream = ZstdDecompressorStream(file)
            read!(stream, tmpN)
            read!(stream, tmpM)
            for n = 1:N
                if direct1 == "+"
                    if featurevec[n] >= thr1
                        nrow += 1
                        idx1[n] = 1
                    end
                elseif direct1 == "-"
                    if featurevec[n] <= thr1
                        nrow += 1
                        idx1[n] = 1
                    end
                end
            end
            close(stream)
        end
    else
        nrow = N
        for n = 1:N
            idx1[n] = 1
        end
    end

    # Column-wise filter
    if samplevec != Nothing
        for m=1:M
            if direct2 == "+"
                if(samplevec[m] >= thr2)
                    ncol += 1
                    push!(idx2, m)
                end
            elseif direct2 == "-"
                if(samplevec[m] <= thr2)
                    ncol += 1
                    push!(idx2, m)
                end
            end
        end
    else
        ncol = M
        for m = 1:M
            push!(idx2, m)
        end
    end
    nrow, ncol, idx1, idx2
end


text:
using Documenter
using GriddingMachine
using Literate
using PkgUtility

# define default docs pages
pages = Pair{Any,Any}[
    "Home" => "index.md",
    "API"  => "API.md"  ,
]

# add preview pages
gen_preview = false;
gen_dir     = joinpath(@__DIR__, "src/generated");
rm(gen_dir, force=true, recursive=true);
mkpath(gen_dir);

if gen_preview
    filename    = joinpath(@__DIR__, "src/preview.jl");
    script      = Literate.script(filename, gen_dir);
    code        = strip(read(script, String));
    mdpost(str) = replace(str, "@__CODE__" => code);
    Literate.markdown(filename, gen_dir, postprocess=mdpost);
    push!(pages, "Data Preview" => "generated/preview.md");
end

@info "Pages to generate:"
pretty_display!(pages);

# format the docs
mathengine = MathJax(Dict(
    :TeX => Dict(
        :equationNumbers => Dict(:autoNumber => "AMS"),
        :Macros => Dict(),
    ),
))

format = Documenter.HTML(
    prettyurls = get(ENV, "CI", nothing) == "true",
    mathengine = mathengine,
    collapselevel = 1,
    assets = ["assets/favicon.ico"]
)

# build the docs
makedocs(
    sitename = "GriddingMachine",
    format = format,
    clean = false,
    modules = [GriddingMachine],
    pages = pages,
)

# function to replace strings (copied from Yujie-W/PAGES)
# TODO move to PkgUtility.jl later
function replace_html(file_name::String)
    list_str   = readlines(file_name);
    file_write = open(file_name, "w");
    for i in eachindex(list_str)
        list_str[i] = replace(list_str[i], "
<details> <summary>" => "");
        list_str[i] = replace(list_str[i], "</summary>
" => "");
        list_str[i] = replace(list_str[i], "
</details>
" => "");
        write(file_write, list_str[i] * "\n");
    end

    return nothing
end

# Replace the strings
file_name = joinpath(@__DIR__, "build/API.html");
if isfile(file_name)
    replace_html(file_name);
else
    file_name = joinpath(@__DIR__, "build/API/index.html");
    if isfile(file_name)
        replace_html(file_name);
    else
        @warn "No file found to work on, please check file names";
    end
end

# deploy the docs to Github gh-pages
deploydocs(
    repo = "github.com/CliMA/GriddingMachine.jl.git",
    target = "build",
    devbranch = "main",
    push_preview = true,
)


text:
using LaserFileSplitter
using Test

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

max_stars_repo_name: zoemcc/LaserFileSplitter

max_issues_repo_name: zoemcc/LaserFileSplitter

max_forks_repo_name: zoemcc/LaserFileSplitter


text:
abstract type DCSource <: InverterComponent end

@def dcsource_ports begin
    state_input = Vector{Symbol}()
    inner_input = Vector{Int64}()
end

"""
Parameters of a Fixed DC Source that returns a fixed DC voltage
# Conmutable structor
```julia
FixedDCSource(voltage)
```

# Arguments
* `voltage`::Float64 : Fixed DC voltage to the converter

"""
mutable struct FixedDCSource <: DCSource
    voltage::Float64
    n_states::Int64
    states::Vector{Symbol}
    ports::Ports

        function FixedDCSource(voltage::Float64)

            @dcsource_ports

            new(voltage,
                0,
                Vector{Symbol}(),
                Ports(state_input, inner_input))
        end
end


text:
using Cairo
using ComplexVisual
@ComplexVisual.import_huge

const CVG = try using ComplexVisualGtk; ComplexVisualGtk; catch nothing; end

const fontface = cv_fontface("sans-serif")

function create_scene(trafo,
        domain, codomain; cut_test=nothing, gap=80,
        axis_label_style=cv_black → fontface → cv_fontsize(20), padding=30) # {{{
    layout = CV_StateLayout(CV_2DLayout(), CV_CyclicValue(2))
    layout = cv_do_lr_layout(cv_add(layout, trafo, domain, codomain), gap)

    setup = cv_setup_cycle_state(CV_LRSetupChain(layout))
    setup = cv_setup_lr_painters(setup; cut_test,
        parallel_hlines_style=cv_op_source → cv_antialias_best →
            cv_linewidth(3) → cv_black,
        parallel_vlines_style=cv_op_source → cv_antialias_best →
            cv_linewidth(3) → cv_white,
        parallel_hlines=cv_parallel_lines(2.0+0.0im),
        parallel_vlines= cv_parallel_lines(0.0+2.0im))
    setup = cv_setup_lr_axis(setup; label_style=axis_label_style)
    setup = cv_setup_lr_border(setup)
    padding > 0 && cv_add_padding!(setup.layout, padding)
    setup = cv_setup_domain_codomain_scene(setup)
    cv_scene_lr_start(setup; state_start=2)
    return setup.layout
end

trafo = exp
domain   = CV_Math2DCanvas(-4.0 + 4.0im, 4.0 - 4.0im, 50)
codomain = CV_Math2DCanvas(-4.0 + 4.0im, 4.0 - 4.0im, 50)

scene = create_scene(trafo, domain, codomain)

cv_save_image(cv_get_can_layout(scene), "./Exp02.png")

if CVG !== nothing
    handler = CVG.cvg_visualize(scene)
    CVG.cvg_wait_for_destroy(handler.window)

    CVG.cvg_close(handler); CVG.cv_destroy(scene)
end

# vim:syn=julia:cc=79:fdm=marker:sw=4:ts=4:


text:
struct QuadratureRule{T,N}
    order::Int
    nodes::Vector{T}
    weights::Vector{T}
end

function QuadratureRule(order, nodes::AbstractVector{T}, weights::AbstractVector{T}, TT=T) where {T}
    @assert length(nodes) == length(weights)
    QuadratureRule{TT, length(nodes)}(order, nodes, weights)
end

"""
Integrate a function `f(x)` over the interval [0,1] using the quadrature.
"""
function (quad::QuadratureRule)(f::Function)
    mapreduce(i -> quad.weights[i] * f(quad.nodes[i]), +, eachindex(quad))
end

nnodes(::QuadratureRule{T,N}) where {T,N} = N
nodes(quad::QuadratureRule) = quad.nodes
order(quad::QuadratureRule) = quad.order
weights(quad::QuadratureRule) = quad.weights

Base.eachindex(quad::QuadratureRule) = eachindex(quad.nodes, quad.weights)

Base.hash(quad::QuadratureRule, h::UInt) = hash(quad.order, hash(quad.nodes, hash(quad.weights, hash(:QuadratureRule, h))))

Base.:(==)(quad1::QuadratureRule, quad2::QuadratureRule) = (
                   quad1.order == quad2.order
                && quad1.nodes == quad2.nodes
                && quad1.weights == quad2.weights)

Base.isequal(quad1::QuadratureRule{T1,N1}, quad2::QuadratureRule{T2,N2}) where {T1,N1,T2,N2} = (
                quad1 == quad2
                && T1 == T2
                && N1 == N2)

Base.isapprox(quad1::QuadratureRule, quad2::QuadratureRule; kwargs...) = (
                            quad1.order == quad2.order
                && isapprox(quad1.nodes,   quad2.nodes;   kwargs...)
                && isapprox(quad1.weights, quad2.weights; kwargs...))

Base.eltype(::QuadratureRule{T}) where {T} = T


text:
# Autogenerated wrapper script for VOLK_jll for x86_64-w64-mingw32-cxx11
export libvolk, list_cpu_features, volk_config_info, volk_profile

using boost_jll
JLLWrappers.@generate_wrapper_header("VOLK")
JLLWrappers.@declare_library_product(libvolk, "libvolk.dll")
JLLWrappers.@declare_executable_product(list_cpu_features)
JLLWrappers.@declare_executable_product(volk_config_info)
JLLWrappers.@declare_executable_product(volk_profile)
function __init__()
    JLLWrappers.@generate_init_header(boost_jll)
    JLLWrappers.@init_library_product(
        libvolk,
        "bin\\libvolk.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        list_cpu_features,
        "bin\\list_cpu_features.exe",
    )

    JLLWrappers.@init_executable_product(
        volk_config_info,
        "bin\\volk-config-info.exe",
    )

    JLLWrappers.@init_executable_product(
        volk_profile,
        "bin\\volk_profile.exe",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
struct LowRankTAME_Return{T} <: returnType
	matchScore::Union{Int,Vector{Int}}
	motifCounts::NTuple{2,Union{Int,Vector{Int}}}
	matching::Union{Dict{Int,Int},Vector{Int}}
	embedding::Tuple{Matrix{T},Matrix{T}}
	profile::Union{Nothing,Vector{Tuple{String,Dict{String,Union{Array{Float64,1},Array{Array{Float64,1},1}}}}}}
end
#=------------------------------------------------------------------------------
              Routines for searching over alpha/beta parameters
------------------------------------------------------------------------------=#
function LowRankTAME_param_search(A::ThirdOrderSymTensor,B::ThirdOrderSymTensor;
						          iter::Int = 15,tol::Float64=1e-6,
						          U_0::Array{Float64,2} = ones(A.n,1),
						          V_0::Array{Float64,2} = ones(B.n,1),
						          alphas::Array{F,1}=[.5,1.0],
						          betas::Array{F,1} =[1000.0,100.0,10.0,1.0,0.0,0.1,0.01,0.001],
						          kwargs...) where {F <: AbstractFloat}
	A_motifs = size(A.indices,2)
	B_motifs = size(B.indices,2)
    max_triangle_match = min(A_motifs,B_motifs)

    total_triangles = size(A.indices,2) + size(B.indices,2)
    best_TAME_PP_tris = -1
	best_matching = Dict{Int,Int}()

	m = A.n
	n = B.n

	best_TAME_PP_U = ones(m,1)
	best_TAME_PP_V = ones(n,1)

    for α in alphas
        for β in betas

			U, V, triangle_count,matching =
				LowRankTAME(A,B,U_0,V_0,β,iter,tol,α;kwargs...)

            if triangle_count > best_TAME_PP_tris
                best_TAME_PP_tris = triangle_count
				best_matching = matching
				best_TAME_PP_U = copy(U)
				best_TAME_PP_V = copy(V)

            end
			println("α:$(α) -- β:$(β) -- tri_match:$(triangle_count) -- max_tris:$(max_triangle_match) -- best tri match:$best_TAME_PP_tris")
        end

    end

	return LowRankTAME_Return(best_TAME_PP_tris,(A_motifs,B_motifs),best_matching,
	                          (best_TAME_PP_U, best_TAME_PP_V), nothing)

end

function LowRankTAME_param_search_profiled(A::ThirdOrderSymTensor,B::ThirdOrderSymTensor;
                           iter::Int = 15,tol::Float64=1e-6,
						   alphas::Array{F,1}=[.5,1.0],
						   betas::Array{F,1} =[1000.0,100.0,10.0,1.0,0.0,0.1,0.01,0.001],
						   kwargs...) where {F <: AbstractFloat}
	A_motifs = size(A.indices,2)
	B_motifs = size(B.indices,2)
    max_triangle_match = min(A_motifs,B_motifs)

    total_triangles = size(A.indices,2) + size(B.indices,2)
	best_TAME_PP_tris = -1
	best_matching = Dict{Int,Int}()

	best_TAME_PP_U = ones(A.n,1)
	best_TAME_PP_V = ones(B.n,1)

	experiment_profiles = Array{Tuple{String,Dict{String,Union{Array{Float64,1},Array{Array{Float64,1},1}}}},1}(undef,0)

    for α in alphas
        for β in betas

			U, V, triangle_count,matching, experiment_profile =
				LowRankTAME_profiled(A,B,ones(A.n,1),ones(B.n,1), β,iter,tol,α;kwargs...)

			push!(experiment_profiles,("α:$(α)_β:$(β)",experiment_profile))

            if triangle_count > best_TAME_PP_tris
				best_TAME_PP_tris = triangle_count
				best_matching = matching
				best_TAME_PP_U = copy(U)
				best_TAME_PP_V = copy(V)

            end
			println("α:$(α) -- β:$(β) -- tri_match:$(best_TAME_PP_tris) -- max_tris $(max_triangle_match)")
        end

    end

	return LowRankTAME_Return(best_TAME_PP_tris,(A_motifs,B_motifs),best_matching,
							 (best_TAME_PP_U, best_TAME_PP_V), experiment_profiles)

end

function LowRankTAME_param_search(A,B;
						          iter::Int = 15,tol::Float64=1e-6,
						          U_0::Array{Float64,2} = ones(A.n,1),
						          V_0::Array{Float64,2} = ones(B.n,1),
						          alphas::Array{F,1}=[.5,1.0],
						          betas::Array{F,1} =[1000.0,100.0,10.0,1.0,0.0,0.1,0.01,0.001],
						          kwargs...) where {F <: AbstractFloat}
	A_motifs = size(A.indices,2)
	B_motifs = size(B.indices,2)
    max_motif_match = min(A_motifs,B_motifs)

    best_TAME_PP_motifs = -1
	best_matching = Dict{Int,Int}()

	m = A.n
	n = B.n

	best_TAME_PP_U = ones(m,1)
	best_TAME_PP_V = ones(n,1)

    for α in alphas
        for β in betas

			U, V, motifs_matched, matching =
				LowRankTAME(A,B,U_0,V_0,β,iter,tol,α;kwargs...)

            if  motifs_matched > best_TAME_PP_motifs
                best_TAME_PP_motifs = motifs_matched
				best_matching = matching
				best_TAME_PP_U = copy(U)
				best_TAME_PP_V = copy(V)

            end
			println("α:$α -- β:$β -- tri_match:$motifs_matched -- max_tris:$max_motif_match -- best tri match:$best_TAME_PP_motifs")
        end

    end

	return LowRankTAME_Return(best_TAME_PP_motifs,(A_motifs,B_motifs),best_matching,
	                          (best_TAME_PP_U, best_TAME_PP_V), nothing)

end

function LowRankTAME_param_search_profiled(A,B;
											iter::Int = 15,tol::Float64=1e-6,
											alphas::Array{F,1}=[.5,1.0],
											betas::Array{F,1} =[100.0,10.0,1.0,0.0],
											kwargs...) where {F <: AbstractFloat}
	A_motifs = size(A.indices,2)
	B_motifs = size(B.indices,2)

	max_motif_match = min(A_motifs,B_motifs)

	best_TAME_PP_motif = -1
	best_matching = Dict{Int,Int}()

	best_TAME_PP_U = ones(A.n,1)
	best_TAME_PP_V = ones(B.n,1)

	experiment_profiles = Array{Tuple{String,Dict{String,Union{Array{Float64,1},Array{Array{Float64,1},1}}}},1}(undef,0)

    for α in alphas
        for β in betas

			U, V, motifs_matched, matching, experiment_profile =
				LowRankTAME_profiled(A,B,ones(A.n,1),ones(B.n,1), β,iter,tol,α;kwargs...)

			push!(experiment_profiles,("α:$(α)_β:$(β)",experiment_profile))

            if motifs_matched > best_TAME_PP_motif
				best_TAME_PP_motif = motifs_matched
				best_matching = matching
				best_TAME_PP_U = copy(U)
				best_TAME_PP_V = copy(V)

            end
			println("α:$α -- β:$β -- motif_match:$best_TAME_PP_motif -- max_motifs $max_motif_match")
        end

    end

	return LowRankTAME_Return(best_TAME_PP_motif,(A_motifs,B_motifs),best_matching,
							 (best_TAME_PP_U, best_TAME_PP_V), experiment_profiles)

end

#=------------------------------------------------------------------------------
             		    Spectral Relaxation Routines
------------------------------------------------------------------------------=#

#runs TAME, but reduces down to lowest rank form first
function LowRankTAME(A::ThirdOrderSymTensor, B::ThirdOrderSymTensor,W::Array{F,2},
                       β::F, max_iter::Int,tol::F,α::F;
					   max_rank::Int = minimum((A.n,B.n)),kwargs...) where {F <:AbstractFloat}

    dimension = minimum((A.n,B.n))

	(U_k,S,VT),t = @timed svd(W)
	singular_indexes = [i for i in 1:minimum((max_rank,length(C_S))) if S[i] > S[1]*eps(Float64)*dimension]

	U = U_k[:,singular_indexes]
	V = VT[:,singular_indexes]*diagm(S[singular_indexes])

	return LowRankTAME(A,B,U,V,β,max_iter,tol,α;kwargs...)
end

"""------------------------------------------------------------------------------
  The low rank implementation of TAME, computes the terms use the mixed property
  generalized to the tensor case.

  Inputs
  ------
  * A,B - (ThirdOrderSymTensor):
	The tensors to align against one another.
  * 'U_0', 'V_0' - (Array{Float,2}):
	The low rank components of the starting iteration X = UV'. Iterates are
	normalized before the iterations begin.
  * β - (Float):
	The shift to use on the iterations.
  * α -(float):
	The mixing parameter for combining in the starting iterates.
  * 'max_iter' - (Int):
	The maximum number of iterations to run.
  * tol - (Float):
	The tolerence to solve the algorithm to. Computes the tolerance by measuring
	the absolute value of the difference between the computed eigenvalues.
  * 'max_rank' - (Int):
	Specify the maximum rank of each of the iterates. Default makes it so that
	only singular values small enough to be considered zero are truncated.
  * 'update_user' - (Int):
	Specifies how frequently output messages should be printed to the user.
	Default is -1 which means no output messages are printed, else if the value
	is k, then the kth iterate will print out an update.
  * 'no_matching' - (Bool):
	Specifies whether or not to run the matching and scoring portions of the
	algorithm. Useful if only the iterates are desired.
  * 'low_rank_matching' - (Bool):
	Specifies whether or not to run the low rank matching procedure from [1].
	This is useful when speed is needed, but may lead to regressions in the
	matching performance.

  Output
  ------
  * 'best_U', 'best_V' - (Array{Float,2})
    Returns the components to the iteration which matched the most triangles.
  * 'best_triangle_count' - (Int)
    The maximum number of triangles matched.
  * 'best_matching' - (Dict{Int,Int})
	The matching computed between the two graphs, maps the vertices of A to the
	vertices of B.

  Citation
  --------
  [1] - H. Nassar, N. Veldt, S. Mohammadi, A. Grama, and D. F. Gleich,
		“Low rank spectral network alignment,” in Proceedings of the 2018
		World Wide Web Conference, 2018, pp. 619–628
------------------------------------------------------------------------------"""
function LowRankTAME(A::ThirdOrderSymTensor, B::ThirdOrderSymTensor,
                          U_0::Array{F,2},V_0::Array{F,2}, β::F, max_iter::Int,tol::F,α::F;
						  max_rank::Int = minimum((A.n,B.n)),update_user::Int=-1,
						  no_matching=false,low_rank_matching=false,kwargs...) where {F <:AbstractFloat}

	@assert size(U_0,2) == size(V_0,2)

	dimension = minimum((A.n,B.n))

	best_triangle_count::Int = -1
	best_matching = Dict{Int,Int}()
    best_x = zeros(size(U_0,1),size(V_0,1))
    best_index = -1
	triangles = -1
	gaped_triangles = -1

	normalization_factor = sqrt(tr((V_0'*V_0)*(U_0'*U_0)))
	U_0 ./= sqrt(normalization_factor)
	V_0 ./= sqrt(normalization_factor)

	U_k = copy(U_0)
	V_k = copy(V_0)

	best_U::Array{F,2} = copy(U_k)
	best_V::Array{F,2} = copy(U_k)

    lambda = Inf

    for i in 1:max_iter

		A_comps, B_comps = get_kron_contract_comps(A,B,U_k,V_k)
		lam = tr((B_comps'*V_k)*(U_k'*A_comps))

		if α != 1.0 && β != 0.0
			U_temp = hcat(sqrt(α) * A_comps, sqrt(α * β) * U_k, sqrt(1-α) * U_0)
			V_temp = hcat(sqrt(α) * B_comps, sqrt(α * β) * V_k, sqrt(1-α) * V_0)
		elseif α != 1.0
			U_temp = hcat(sqrt(α)*A_comps, sqrt(1-α)*U_0)
			V_temp = hcat(sqrt(α)*B_comps, sqrt(1-α)*V_0)
		elseif β != 0.0
			U_temp = hcat(A_comps, sqrt(β) * U_k)
			V_temp = hcat(B_comps, sqrt(β) * V_k)
		else
			U_temp = A_comps
			V_temp = B_comps
		end

		A_Q,A_R = qr(U_temp)
		B_Q,B_R = qr(V_temp)

		core = A_R*B_R'
		C_U,C_S,C_Vt = svd(core)
		singular_indexes= [i for i in 1:1:minimum((max_rank,length(C_S))) if C_S[i] > C_S[1]*eps(Float64)*dimension]

		U_k_1 = A_Q*C_U[:,singular_indexes]
		V_k_1 = B_Q*(C_Vt[:,singular_indexes]*diagm(C_S[singular_indexes]))

		normalization_factor = sqrt(tr((V_k_1'*V_k_1)*(U_k_1'*U_k_1)))

		U_k_1 ./= sqrt(normalization_factor)
		V_k_1 ./= sqrt(normalization_factor)

		#TODO: need to redo or remove this
		Y, Z = get_kron_contract_comps(A,B,U_k_1,V_k_1)

		lam = tr((Y'*U_k_1)*(V_k_1'*Z))

		if !no_matching
			#evaluate the matchings
			if low_rank_matching
				triangles, gaped_triangles, matching = TAME_score(A,B,U_k_1,V_k_1;kwargs...)
			else
				triangles, gaped_triangles, matching = TAME_score(A,B,U_k_1*V_k_1';kwargs...)
			end

			if triangles > best_triangle_count
				best_matching = matching
				best_triangle_count  = triangles
				best_U = copy(U_k_1)
				best_V = copy(V_k_1)
			end
		end

		if update_user != -1 && i % update_user == 0
			println("λ_$i: $(lam) -- rank:$(length(singular_indexes)) -- tris:$(triangles) -- gaped_t:$(gaped_triangles)")
		end

        if abs(lam - lambda) < tol
    		return best_U, best_V, best_triangle_count, best_matching
        end

		#get the low rank factorization for the next one
		U_k = copy(U_k_1)
		V_k = copy(V_k_1)

		lambda = lam

    end

	return best_U, best_V, best_triangle_count, best_matching
end

"""------------------------------------------------------------------------------
  The low rank implementation of TAME, computes the terms use the mixed property
  generalized to the tensor case.

  Inputs
  ------
  * A,B - (ThirdOrderSymTensor):
	The tensors to align against one another.
  * 'U_0', 'V_0' - (Array{Float,2}):
	The low rank components of the starting iteration X = UV'. Iterates are
	normalized before the iterations begin.
  * β - (Float):
	The shift to use on the iterations.
  * α -(float):
	The mixing parameter for combining in the starting iterates.
  * 'max_iter' - (Int):
	The maximum number of iterations to run.
  * tol - (Float):
	The tolerence to solve the algorithm to. Computes the tolerance by measuring
	the absolute value of the difference between the computed eigenvalues.
  * 'max_rank' - (Int):
	Specify the maximum rank of each of the iterates. Default makes it so that
	only singular values small enough to be considered zero are truncated.
  * 'update_user' - (Int):
	Specifies how frequently output messages should be printed to the user.
	Default is -1 which means no output messages are printed, else if the value
	is k, then the kth iterate will print out an update.
  * 'no_matching' - (Bool):
	Specifies whether or not to run the matching and scoring portions of the
	algorithm. Useful if only the iterates are desired.
  * 'low_rank_matching' - (Bool):
	Specifies whether or not to run the low rank matching procedure from [1].
	This is useful when speed is needed, but may lead to regressions in the
	matching performance.

  Output
  ------
	* 'best_U', 'best_V' - (Array{Float,2})
	  Returns the components to the iteration which matched the most triangles.
	* 'best_triangle_count' - (Int)
	  The maximum number of triangles matched.
	* 'best_matching' - (Dict{Int,Int})
	  The matching computed between the two graphs, maps the vertices of A to the
	  vertices of B.
	* 'experiment_profile' - (Dict{String,Union{Array{F,1},Array{Array{F,1},1}}}):
	  The experiment profile computed, keys for the experiment data collected are
	  as follows.
	  +'ranks' - The ranks of each iterate X_k.
	  +'contraction_timings' - Time taken to compute each contraction.
	  +'svd_timings' - Time taken to compute the svd.
	  +'qr_timings' - Time taken to compute the QR factorations.
	  +'matched_tris' - The number of triangles matched by each iterate.
	  +'sing_vals' - The singular values of each iterate X_k.
	  +'matching_times' - Time taken to solve the matchings.
	  +'soring_timings' - Time taken to score each matching.

  Citation
  --------
	[1] - H. Nassar, N. Veldt, S. Mohammadi, A. Grama, and D. F. Gleich,
		  “Low rank spectral network alignment,” in Proceedings of the 2018
		  World Wide Web Conference, 2018, pp. 619–628
------------------------------------------------------------------------------"""
function LowRankTAME_profiled(A::ThirdOrderSymTensor, B::ThirdOrderSymTensor,
                          U_0::Array{F,2},V_0::Array{F,2}, β::F, max_iter::Int,tol::F,α::F;
						  max_rank::Int = minimum((A.n,B.n)),update_user::Int=-1,
						  no_matching::Bool=false,low_rank_matching::Bool=false,kwargs...) where {F <:AbstractFloat}

	@assert size(U_0,2) == size(V_0,2)

	dimension = minimum((A.n,B.n))

	best_triangle_count::Int = -1
	best_matching = Dict{Int,Int}()
    best_x = zeros(size(U_0,1),size(V_0,1))
    best_index = -1
	triangles = -1
	gaped_triangles = -1

	normalization_factor = sqrt(tr((V_0'*V_0)*(U_0'*U_0)))
	U_0 ./= sqrt(normalization_factor)
	V_0 ./= sqrt(normalization_factor)

	experiment_profile = Dict{String,Union{Array{F,1},Array{Array{F,1},1}}}(
		"ranks"=>Array{Float64,1}(undef,0),
		"contraction_timings"=>Array{Float64,1}(undef,0),
		"svd_timings"=>Array{Float64,1}(undef,0),
		"qr_timings"=>Array{Float64,1}(undef,0),
		"matched_tris"=>Array{Float64,1}(undef,0),
		"sing_vals"=>Array{Array{Float64,1},1}(undef,0),
		"matching_timings"=>Array{Float64,1}(undef,0),
		"scoring_timings"=>Array{Float64,1}(undef,0)
	)

	U_k = copy(U_0)
	V_k = copy(V_0)

	best_U::Array{F,2} = copy(U_k)
	best_V::Array{F,2} = copy(U_k)

    lambda = Inf

    for i in 1:max_iter

		(A_comps, B_comps),t = @timed get_kron_contract_comps(A,B,U_k,V_k)
		push!(experiment_profile["contraction_timings"],t)

		lambda_k_1 = tr((B_comps'*V_k)*(U_k'*A_comps))

		if α != 1.0 && β != 0.0
			U_temp = hcat(sqrt(α) * A_comps, sqrt(α * β) * U_k, sqrt(1-α) * U_0)
			V_temp = hcat(sqrt(α) * B_comps, sqrt(α * β) * V_k, sqrt(1-α) * V_0)
		elseif α != 1.0
			U_temp = hcat(sqrt(α)*A_comps, sqrt(1-α)*U_0)
			V_temp = hcat(sqrt(α)*B_comps, sqrt(1-α)*V_0)
		elseif β != 0.0
			U_temp = hcat(A_comps, sqrt(β) * U_k)
			V_temp = hcat(B_comps, sqrt(β) * V_k)
		else
			U_temp = A_comps
			V_temp = B_comps
		end

		(A_Q,A_R),t_A = @timed qr(U_temp)
		(B_Q,B_R),t_B = @timed qr(V_temp)
		push!(experiment_profile["qr_timings"],t_A + t_B)

		core = A_R*B_R'
		(C_U,C_S::Array{Float64,1},C_Vt),t = @timed svd(core)
		push!(experiment_profile["svd_timings"],t)

		singular_indexes= [i for i in 1:1:minimum((max_rank,length(C_S))) if C_S[i] > C_S[1]*eps(Float64)*dimension]
		push!(experiment_profile["sing_vals"],C_S)
		push!(experiment_profile["ranks"],float(length(singular_indexes)))

		U_k_1 = A_Q*C_U[:,singular_indexes]
		V_k_1 = B_Q*(C_Vt[:,singular_indexes]*diagm(C_S[singular_indexes]))

		normalization_factor = sqrt(tr((V_k_1'*V_k_1)*(U_k_1'*U_k_1)))

		U_k_1 ./= sqrt(normalization_factor)
		V_k_1 ./= sqrt(normalization_factor)

		if !no_matching
			#evaluate the matchings
			if low_rank_matching
				triangles, gaped_tris, matching, matching_time, scoring_time = TAME_score(A,B,U_k_1,V_k_1;return_timings=returnTimings(),kwargs...)
			else
				triangles, gaped_tris, matching, matching_time, scoring_time = TAME_score(A,B,U_k_1*V_k_1';return_timings=returnTimings(),kwargs...)
			end

			push!(experiment_profile["matched_tris"],float(triangles))
			push!(experiment_profile["matching_timings"],float(matching_time))
			push!(experiment_profile["scoring_timings"], float(scoring_time))

			if triangles > best_triangle_count
				best_matching = matching
				best_triangle_count  = triangles
				best_U = copy(U_k_1)
				best_V = copy(V_k_1)
			end
		end

		if update_user != -1 && i % update_user == 0
			println("λ_$i: $(lambda_k_1) -- rank:$(length(singular_indexes)) -- tris:$(triangles) -- gaped_t:$(gaped_tris)")
		end

		if abs(lambda_k_1 - lambda) < tol || i >= max_iter
			#=
			triangles,_= TAME_score(A,B,sparse(best_U*best_V');return_timings=returnTimings())
			if triangles > best_triangle_count
				best_triangle_count = triangles
			end
			=#
			return best_U, best_V, best_triangle_count,best_matching, experiment_profile
		end
		#get the low rank factorization for the next one
		U_k = copy(U_k_1)
		V_k = copy(V_k_1)

		lambda = lambda_k_1
	end

	return best_U, best_V, best_triangle_count,best_matching, experiment_profile
end

function LowRankTAME(A::SymTensorUnweighted{M}, B::SymTensorUnweighted{M},
                          U_0::Array{F,2},V_0::Array{F,2}, β::F, max_iter::Int,tol::F,α::F;
						  max_rank::Int = minimum((A.n,B.n)),update_user::Int=-1,
						  no_matching::Bool=false,low_rank_matching::Bool=false,kwargs...) where {M <:Motif, F <:AbstractFloat}

	@assert size(U_0,2) == size(V_0,2)

	dimension = minimum((A.n,B.n))

	best_matched_motif_count::Int = -1
	best_matching = Dict{Int,Int}()

	matched_motifs = -1
	missed_motifs = -1

	U_0 ./= norm(U_0)
	V_0 ./= norm(V_0)

	U_k = copy(U_0)
	V_k = copy(V_0)
	X_k = U_0*V_0'

	if α != 1.0
		X_0_with_alpha = ((1-α)*U_0)*V_0'
	end

	best_U::Array{F,2} = copy(U_k)
	best_V::Array{F,2} = copy(U_k)

    lambda = Inf

    for i in 1:max_iter

		#use pure low rank version if won't use too much memory
		if binomial(size(U_k,2) + A.order-2, A.order-1) <= 10*dimension

			(A_comps,B_comps),t = @timed get_kron_contract_comps(A,B,U_k,V_k)

			lambda_k_1 = tr((B_comps'*V_k)*(U_k'*A_comps))

			if α != 1.0 && β != 0.0
				U_temp = hcat(sqrt(α) * A_comps, sqrt(α * β) * U_k, sqrt(1-α) * U_0)
				V_temp = hcat(sqrt(α) * B_comps, sqrt(α * β) * V_k, sqrt(1-α) * V_0)
			elseif α != 1.0
				U_temp = hcat(sqrt(α)*A_comps, sqrt(1-α)*U_0)
				V_temp = hcat(sqrt(α)*B_comps, sqrt(1-α)*V_0)
			elseif β != 0.0
				U_temp = hcat(A_comps, sqrt(β) * U_k)
				V_temp = hcat(B_comps, sqrt(β) * V_k)
			else
				U_temp = A_comps
				V_temp = B_comps
			end

			A_Q,A_R = qr(U_temp)
			B_Q,B_R = qr(V_temp)

			core = A_R*B_R'

			C_U,C_S::Array{Float64,1},C_Vt = svd(core)

			singular_indexes= [i for i in 1:1:minimum((max_rank,length(C_S))) if C_S[i] > C_S[1]*eps(Float64)*dimension]

			U_k_1 = A_Q*C_U[:,singular_indexes]
			V_k_1 = B_Q*(C_Vt[:,singular_indexes]*diagm(C_S[singular_indexes]))

			normalization_factor = sqrt(tr((V_k_1'*V_k_1)*(U_k_1'*U_k_1)))

			U_k_1 ./= sqrt(normalization_factor)
			V_k_1 ./= sqrt(normalization_factor)

			X_k_1 = U_k_1*V_k_1'
		else

			#use accumulation parameter version if too much memory will be used

			X_k_1 = get_kron_contract_comps_with_accumulation_param(A,B,U_k,V_k)

			lambda_k_1 = tr((U_k'*X_k_1)*V_k)
						# equivilant to <X,UV'> but uses less memory

			if β != 0.0
				X_k_1 .+= β*X_k
			end

			if α != 1.0
				X_k_1 .*= α
				X_k_1 .+= X_0_with_alpha # has (1-α) inco
			end

			X_k_1 ./=norm(X_k_1)

			(U, S, V) = svd(X_k_1)
			singular_indexes= [i for i in 1:length(S) if S[i] > S[1]*eps(Float64)*dimension]

			U_k_1 = U[:,singular_indexes]*diagm(S[singular_indexes])
			V_k_1 = V[:,singular_indexes]

		end

		if !no_matching
			#evaluate the matchings
			if low_rank_matching
				matched_motifs, missed_motifs, matching  = TAME_score(A,B,U_k_1,V_k_1;kwargs...)
			else
				matched_motifs, missed_motifs, matching = TAME_score(A,B,X_k_1;kwargs...)
			end

			if matched_motifs > best_matched_motif_count
				best_matching = matching
				best_matched_motif_count = matched_motifs
				best_U = copy(U_k_1)
				best_V = copy(V_k_1)
			end
		end

		if update_user != -1 && i % update_user == 0
			println("λ_$i: $(lambda_k_1) -- rank:$(length(singular_indexes)) -- matched motifs:$(matched_motifs) -- missed motifs:$(missed_motifs)")
		end

		if abs(lambda_k_1 - lambda) < tol || i >= max_iter
			return best_U, best_V, best_matched_motif_count,best_matching
		end
		#get the low rank factorization for the next one
		U_k = copy(U_k_1)
		V_k = copy(V_k_1)
		X_k = copy(X_k_1)

		lambda = lambda_k_1
	end

	return best_U, best_V, best_matched_motif_count, best_matching
end

function LowRankTAME_profiled(A::SymTensorUnweighted{M}, B::SymTensorUnweighted{M},
                          U_0::Array{F,2},V_0::Array{F,2}, β::F, max_iter::Int,tol::F,α::F;
						  max_rank::Int = minimum((A.n,B.n)),update_user::Int=-1,
						  no_matching::Bool=false,low_rank_matching::Bool=false,kwargs...) where {M <:Motif, F <:AbstractFloat}

	@assert size(U_0,2) == size(V_0,2)

	dimension = minimum((A.n,B.n))

	best_matched_motif_count::Int = -1
	best_matching = Dict{Int,Int}()

	matched_motifs = -1
	missed_motifs = -1

	U_0 ./= norm(U_0)
	V_0 ./= norm(V_0)

	experiment_profile = Dict{String,Union{Array{F,1},Array{Array{F,1},1}}}(
		"ranks"=>Array{Float64,1}(undef,0),
		"contraction_timings"=>Array{Float64,1}(undef,0),
		"low_rank_factoring_timings"=>Array{Float64,1}(undef,0),
		"matched_motifs"=>Array{Float64,1}(undef,0),
		"sing_vals"=>Array{Array{Float64,1},1}(undef,0),
		"matching_timings"=>Array{Float64,1}(undef,0),
		"scoring_timings"=>Array{Float64,1}(undef,0)
	)

	U_k = copy(U_0)
	V_k = copy(V_0)
	X_k = U_0*V_0'

	if α != 1.0
		X_0_with_alpha = ((1-α)*U_0)*V_0'
	end

	best_U::Array{F,2} = copy(U_k)
	best_V::Array{F,2} = copy(U_k)

    lambda = Inf

    for i in 1:max_iter

		#use pure low rank version if won't use too much memory
		if binomial(size(U_k,2) + A.order-2, A.order-1) <= 10*dimension

			(A_comps,B_comps),t = @timed LambdaTAME.get_kron_contract_comps(A,B,U_k,V_k)
			push!(experiment_profile["contraction_timings"],t)

			lambda_k_1 = tr((B_comps'*V_k)*(U_k'*A_comps))

			if α != 1.0 && β != 0.0
				U_temp = hcat(sqrt(α) * A_comps, sqrt(α * β) * U_k, sqrt(1-α) * U_0)
				V_temp = hcat(sqrt(α) * B_comps, sqrt(α * β) * V_k, sqrt(1-α) * V_0)
			elseif α != 1.0
				U_temp = hcat(sqrt(α)*A_comps, sqrt(1-α)*U_0)
				V_temp = hcat(sqrt(α)*B_comps, sqrt(1-α)*V_0)
			elseif β != 0.0
				U_temp = hcat(A_comps, sqrt(β) * U_k)
				V_temp = hcat(B_comps, sqrt(β) * V_k)
			else
				U_temp = A_comps
				V_temp = B_comps
			end

			(A_Q,A_R),t_A = @timed qr(U_temp)
			(B_Q,B_R),t_B = @timed qr(V_temp)
			low_rank_factor_t = t_A + t_B
			#push!(experiment_profile["qr_timings"],t_A + t_B)

			core,t = @timed A_R*B_R'
			low_rank_factor_t += t
			(C_U,C_S::Array{Float64,1},C_Vt),t = @timed svd(core)
			#push!(experiment_profile["svd_timings"],t)
			low_rank_factor_t += t

			push!(experiment_profile["low_rank_factoring_timings"],low_rank_factor_t)

			singular_indexes= [i for i in 1:1:minimum((max_rank,length(C_S))) if C_S[i] > C_S[1]*eps(Float64)*dimension]
			push!(experiment_profile["sing_vals"],C_S)
			push!(experiment_profile["ranks"],float(length(singular_indexes)))

			U_k_1 = A_Q*C_U[:,singular_indexes]
			V_k_1 = B_Q*(C_Vt[:,singular_indexes]*diagm(C_S[singular_indexes]))

			normalization_factor = sqrt(tr((V_k_1'*V_k_1)*(U_k_1'*U_k_1)))

			U_k_1 ./= sqrt(normalization_factor)
			V_k_1 ./= sqrt(normalization_factor)

			X_k_1 = U_k_1*V_k_1'
		else

			#use accumulation parameter version if too much memory will be used

			X_k_1,t = @timed get_kron_contract_comps_with_accumulation_param(A,B,U_k,V_k)
			push!(experiment_profile["contraction_timings"],t)

			lambda_k_1 = tr((U_k'*X_k_1)*V_k)
						# equivilant to <X,UV'> but uses less memory

			if β != 0.0
				X_k_1 .+= β*X_k
			end

			if α != 1.0
				X_k_1 .*= α
				X_k_1 .+= X_0_with_alpha # has (1-α) incorporated already
			end

			X_k_1 ./=norm(X_k_1)

			(U, S, V),low_rank_factor_t  = @timed svd(X_k_1)
			push!(experiment_profile["low_rank_factoring_timings"],low_rank_factor_t)

			singular_indexes= [i for i in 1:length(S) if S[i] > S[1]*eps(Float64)*dimension]

			push!(experiment_profile["sing_vals"],S[singular_indexes])
			push!(experiment_profile["ranks"],float(length(singular_indexes)))

			U_k_1 = U[:,singular_indexes]*diagm(S[singular_indexes])
			V_k_1 = V[:,singular_indexes]

		end

		if !no_matching
			#evaluate the matchings
			if low_rank_matching
				matched_motifs, missed_motifs, matching, matching_time, scoring_time = TAME_score(A,B,U_k_1,V_k_1;return_timings=returnTimings(),kwargs...)
			else
				matched_motifs, missed_motifs, matching, matching_time, scoring_time = TAME_score(A,B,X_k_1;return_timings=returnTimings(),kwargs...)
			end

			push!(experiment_profile["matched_motifs"],float(matched_motifs))
			push!(experiment_profile["matching_timings"],float(matching_time))
			push!(experiment_profile["scoring_timings"], float(scoring_time))

			if matched_motifs > best_matched_motif_count
				best_matching = matching
				best_matched_motif_count = matched_motifs
				best_U = copy(U_k_1)
				best_V = copy(V_k_1)
			end
		end

		if update_user != -1 && i % update_user == 0
			println("λ_$i: $(lambda_k_1) -- rank:$(length(singular_indexes)) -- matched motifs:$(matched_motifs) -- missed motifs:$(missed_motifs)")
		end

		if abs(lambda_k_1 - lambda) < tol || i >= max_iter
			return best_U, best_V, best_matched_motif_count,best_matching, experiment_profile
		end
		#get the low rank factorization for the next one
		U_k = copy(U_k_1)
		V_k = copy(V_k_1)
		X_k = copy(X_k_1)

		lambda = lambda_k_1
	end

	return best_U, best_V, best_matched_motif_count,best_matching, experiment_profile
end

max_stars_repo_name: chuckcol/LambdaTAME

max_issues_repo_name: chuckcol/LambdaTAME

max_forks_repo_name: chuckcol/LambdaTAME


text:
using Test
using CLIMA
using CLIMA.RootSolvers

@testset "RootSolvers - compact solution correctness" begin
  f(x) = x^2 - 100^2
  f′(x) = 2x
  for FT in [Float32, Float64]
    sol = find_zero(f, FT(0.0), FT(1000.0), SecantMethod(), CompactSolution())
    @test sol.converged
    @test sol.root isa FT
    @test sol.root ≈ 100

    sol = find_zero(f, FT(0.0), FT(1000.0), RegulaFalsiMethod(), CompactSolution())
    @test sol.converged
    @test sol.root isa FT
    @test sol.root ≈ 100

    sol = find_zero(f, FT(1.0), NewtonsMethodAD(), CompactSolution())
    @test sol.converged
    @test sol.root isa FT
    @test sol.root ≈ 100

    sol = find_zero(f, f′, FT(1.0), NewtonsMethod(), CompactSolution())
    @test sol.converged
    @test sol.root isa FT
    @test sol.root ≈ 100
  end
end

@testset "RootSolvers - verbose solution correctness" begin
  f(x) = x^2 - 100^2
  f′(x) = 2x
  for FT in [Float32, Float64]
    sol = find_zero(f, FT(0.0), FT(1000.0), SecantMethod(), VerboseSolution())
    @test sol.converged
    @test sol.root isa FT
    @test sol.root ≈ 100
    @test sol.err < 1e-3
    @test sol.iter_performed < 20
    @test sol.iter_performed+1 == length(sol.root_history) == length(sol.err_history)

    sol = find_zero(f, FT(0.0), FT(1000.0), RegulaFalsiMethod(), VerboseSolution())
    @test sol.converged
    @test sol.root isa FT
    @test sol.root ≈ 100
    @test sol.err < 1e-3
    @test sol.iter_performed < 20
    @test sol.iter_performed+1 == length(sol.root_history) == length(sol.err_history)

    sol = find_zero(f, FT(1.0), NewtonsMethodAD(), VerboseSolution())
    @test sol.converged
    @test sol.root isa FT
    @test sol.root ≈ 100
    @test sol.err < 1e-3
    @test sol.iter_performed < 20
    @test sol.iter_performed+1 == length(sol.root_history) == length(sol.err_history)

    sol = find_zero(f, f′, FT(1.0), NewtonsMethod(), VerboseSolution())
    @test sol.converged
    @test sol.root isa FT
    @test sol.root ≈ 100
    @test sol.err < 1e-3
    @test sol.iter_performed < 20
    @test sol.iter_performed+1 == length(sol.root_history) == length(sol.err_history)
  end
end

if CLIMA.array_type() != Array
  CLIMA.gpu_allowscalar(false)

  @testset "RootSolvers CUDA - compact solution " begin
    for FT in [Float32, Float64]
      X0 = cu(rand(FT, 5,5))
      X1 = cu(rand(FT, 5,5)) .+ 1000
      f(x) = x^2 - 100^2

      sol = RootSolvers.find_zero.(f, X0, X1, SecantMethod(), CompactSolution())
      converged = map(x->x.converged, sol)
      X_roots = map(x->x.root, sol)
      @test all(converged)
      @test eltype(X_roots) == eltype(X0)
      @test all(X_roots .≈ 100)
    end
  end

  @testset "RootSolvers CUDA - verbose solution " begin
    for FT in [Float32, Float64]
      X0 = cu(rand(FT, 5,5))
      X1 = cu(rand(FT, 5,5)) .+ 1000
      f(x) = x^2 - 100^2

      sol = RootSolvers.find_zero.(f, X0, X1, SecantMethod(), VerboseSolution())
      converged = map(x->x.converged, sol)
      X_roots = map(x->x.root, sol)
      err = map(x->x.err, sol)
      iter_performed = map(x->x.iter_performed, sol)
      @test all(converged)
      @test eltype(X_roots) == eltype(X0)
      @test all(X_roots .≈ 100)
      @test all(err .< 1e-2)
      @test all(iter_performed .< 20)
    end
  end
end

max_stars_repo_name: emassoud/CLIMA

max_issues_repo_name: emassoud/CLIMA

max_forks_repo_name: emassoud/CLIMA


text:
#Test code
#-------------------------------------------------------------------------------

using FileIO2

#No real test code yet... just run demos:

#Use alias shorthands:
const ASCIIFmt = FileIO2.ASCIIFmt
const UTF8Fmt  = FileIO2.UTF8Fmt
const UTF16Fmt = FileIO2.UTF16Fmt
const UTF32Fmt = FileIO2.TextFormat{FileIO2.UTF32Encoding}

const ASCIIHTMLFmt = FileIO2.ASCIIHTMLFmt
const UTF16HTMLFmt = FileIO2.HTMLFormat{FileIO2.UTF16Encoding}

#For extra convenience:
const TextFile = File{TextFmt}
const ASCIIFile = File{ASCIIFmt}

const HTMLFile = File{FileIO2.HTMLFmt}

samplefile(filename) =
	joinpath(FileIO2.rootpath, "test/samplefiles", filename)

sepline = "---------------------------------------------------------------------"

println("\nText clean/simple method to print out a text file:")
println(sepline)
print(read(File{FileIO2.TextFormat}(samplefile("textfile.txt")))) #This works on abstract types too
print(read(TextFile(samplefile("textfile.txt"))))
print(read(ASCIIFile(samplefile("asciifile.txt"))))
print(read(File{UTF8Fmt}(samplefile("utf8file.txt"))))

filelist = [
	#Prefered notation (no specific encoding):
	File{TextFmt}("generictextfile.txt")
	#Even better, in some instances:
	File(TextFmt, "generictextfile.txt")

	#Verbose constructor tests:
	File{FileIO2.TextFormat{FileIO2.UTF8Encoding}}("utf8file.txt")

	#Shorthand notation (alias made by user):
	ASCIIFile("asciifile.txt")
]

println("\n=====Test text file list:")
show(filelist)

filelist = [
	HTMLFile("generichtmlfile.html")
	File{UTF16HTMLFmt}("utf16htmlfile.html")
	File{ASCIIHTMLFmt}("asciihtmlfile.html")
]

println("\n=====Test HTML file list:")
show(filelist)

println("\n...Cast to simple text file list")
println("   (to process files with >>basic text<< handling code):")
show(File(FileIO2.TextFormat, filelist))

println("\n...Same thing with shorthand notation:")
println("   (to process files with >>basic text<< handling code):")
show(File(:text, filelist))

println("\nTry auto-detecting file:")
x=File("displaynogomsg.txt")

println("\nGet extensions for an html file:")
show(Extensions(HTMLFile)); println()

println("\nTest open/read/close dispatch system:")
println(sepline)

try; read(File(:bmp, "fakefile.bmp"))
catch e; @warn(e.msg)
end

module PNGReaderTest
	using FileIO2
	struct PNGReader1 <: AbstractReader{PNGFmt}; end
	struct PNGReader2 <: AbstractReader{PNGFmt}; end
	Base.open(r::Type{PNGReader1}, path::String) = #Throw error to look for working reader
		(msg = "$PNGReader1: Failed to open $path"; @warn(msg); error(msg))
	Base.open(r::Type{PNGReader2}, path::String) = @info("\nOpening $path...")
end
#using PNGReaderTest #no longer need to do "using"

open(File(:png, "fakefile.png"), read=true)

println("\nTest text reader system:")
println(sepline)

println("\nUse introspection to auto-select \"text\" reader:")
open(File(:text, samplefile("integer_data.txt")), read=true) do reader
	@show typeof(reader)
	while !eof(reader)
		data = read(reader, Int)
		println("Read: $(typeof(data))(`$data`)")
	end
end
#NOTE: Open :csv instead of :text if ',' is a desired separator

#Ensure we use FileIO2.TextReader - to ensure consisten
println("\nForce use of FileIO2.TextReader:")
open(FileIO2.TextReader, samplefile("integer_data.txt")) do reader
	while !eof(reader)
		data = read(reader, Int)
		println("Read: $(typeof(data))(`$data`)")
	end
end

:Test_Complete


text:
immutable Normal <: ContinuousUnivariateDistribution
    mean::Float64
    std::Float64
    function Normal(mu::Real, sd::Real)
    	if sd > 0.0
    		new(float64(mu), float64(sd))
    	else
    		error("std must be positive")
    	end
    end
end

Normal(mu::Real) = Normal(mu, 1.0)
Normal() = Normal(0.0, 1.0)

const Gaussian = Normal

@_jl_dist_2p Normal norm

entropy(d::Normal) = (1.0 / 2.0) * log(2.0 * pi * e * d.std^2)

insupport(d::Normal, x::Number) = isreal(x) && isfinite(x)

kurtosis(d::Normal) = 0.0

mean(d::Normal) = d.mean

median(d::Normal) = d.mean

function mgf(d::Normal, t::Real)
	m, s = d.mean, d.std
	return exp(t * m + 0.5 * s^t * t^2)
end

function cf(d::Normal, t::Real)
	m, s = d.mean, d.std
	return exp(im * t * m - 0.5 * s^t * t^2)
end

modes(d::Normal) = [d.mean]

rand(d::Normal) = d.mean + d.std * randn()

skewness(d::Normal) = 0.0

std(d::Normal) = d.std

var(d::Normal) = d.std^2


text:
module Cryptolocker
using Dates

function gen_domain(date::Dates.Date, tld="", num=8)
    domain = ""
    year, month, day = Dates.year(date), Dates.month(date), Dates.day(date)
    for _ in 1:num
        year = ((year ⊻ 8 * year) >> 11) ⊻ ((year & 0xFFFFFFF0) << 17)
        month = ((month ⊻ 4 * month) >> 25) ⊻ 16 * (month & 0xFFFFFFF8)
        day = ((day ⊻ (day << 13)) >> 19) ⊻ ((day & 0xFFFFFFFE) << 12)
        domain *= ((year ⊻ month ⊻ day) % 25) + 'a'
    end
    domain *= tld
    return domain
end

function gen_one(date=missing, tld="", num=8)
    if date === missing
        date = Dates.Date(Dates.now())
    end
    domain = gen_domain(date, tld, num)
    return domain
end

function gen_many(cnt=100)
    res = []
    iter_num = trunc(Int,cnt/24)
    remain_num = cnt%24
    date = Dates.Date(Dates.now())
    for num in 8:31
        start_date = date + Dates.Day(trunc(Int64, rand()*1000000))
        for i in 1:iter_num
            new_date = start_date + Dates.Day(i*2+1)
            res = [res..., gen_domain(new_date, "", num)]
        end
    end
    for i in 1:remain_num
        res = [res..., gen_one(date, "", i+7)]
    end
    return res
end

end  # module Cryptolocker


text:
"""
SARSA-λ learner with surprise-modulated learning rate (core)
"""

mutable struct SarsaLModulated
    nr_states::Int
    nr_actions::Int
    α::Float64
    α_eff::Float64
    γ::Float64
    λ::Float64
    initvalue::Float64
    QSA::Array{Float64, 2}
    eligTrace::Array{Float64, 2}
    rpe::Array{Float64, 1}
    use_replacing_trace::Bool
    do_debug_print::Bool
end
function SarsaLModulated(; nr_states = 10, nr_actions = 4, α = 0.1, γ = .9,
    λ = 0.9, initvalue = 0.,
    QSA = zeros(Float64, nr_states, nr_actions) .+ initvalue,
    use_replacing_trace = true,
    do_debug_print = false)

    SarsaLModulated(nr_states, nr_actions, α, α, γ, λ, initvalue, QSA,
        zeros(Float64, nr_states, nr_actions), [0.],
        use_replacing_trace, do_debug_print)
end
export SarsaLModulated
function getlearningrate!(learner::Union{SarsaLModulated, REINFORCEModulated},
                        modulatingfactor)
    learner.α_eff = learner.α * modulatingfactor
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "Pango"
version = v"1.42.4"

# Collection of sources required to build Pango
sources = [
    ArchiveSource("http://ftp.gnome.org/pub/GNOME/sources/pango/$(version.major).$(version.minor)/pango-$(version).tar.xz",
                  "1d2b74cd63e8bd41961f2f8d952355aa0f9be6002b52c8aa7699d9f5da597c9d"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/pango-*/

# Remove misleading libtool files
rm -f ${prefix}/lib/*.la

# Be able to find libffi on ppc64le
if [[ ${target} == powerpc64le* ]]; then
    export LDFLAGS="${LDFLAGS} -Wl,-rpath-link,${prefix}/lib64"
fi

./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} \
    --disable-introspection \
    --disable-gtk-doc-html

# The generated Makefile tries to build some examples in the "tests" directory,
# but this would fail for some unknown reasons.  Let's skip it.
sed -i 's/^\(SUBDIRS = .*\) tests/\1/' Makefile
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct(["libpango", "libpango-1", "libpango-1.0"], :libpango),
    LibraryProduct(["libpangocairo", "libpangocairo-1", "libpangocairo-1.0"], :libpangocairo),
    LibraryProduct(["libpangoft2", "libpangoft2-1", "libpangoft2-1.0"], :libpangoft),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("FriBidi_jll"),
    Dependency("FreeType2_jll"),
    Dependency("Glib_jll"),
    Dependency("Fontconfig_jll"),
    Dependency("HarfBuzz_jll"),
    Dependency("Cairo_jll"),
    BuildDependency("Xorg_xorgproto_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: waralex/Yggdrasil

max_issues_repo_name: waralex/Yggdrasil

max_forks_repo_name: waralex/Yggdrasil


text:
mutable struct GreedyThickThinning <: GraphSearchStrategy
    cache::ScoreComponentCache
    max_n_parents::Int
    prior::DirichletPrior

    function GreedyThickThinning(
        cache::ScoreComponentCache;
        max_n_parents::Int=3,
        prior::DirichletPrior=UniformPrior(1.0),
        )

        new(cache, max_n_parents, prior)
    end
end

function Distributions.fit(::Type{DiscreteBayesNet}, data::DataFrame, params::GreedyThickThinning;
    ncategories::Vector{Int} = map!(i->infer_number_of_instantiations(data[!,i]), Array{Int}(ncol(data)), 1:ncol(data)),
    )

    n = ncol(data)
    parent_list = map!(i->Int[], Array{Vector{Int}}(n), 1:n)
    datamat = convert(Matrix{Int}, data)'
    score_components = bayesian_score_components(parent_list, ncategories, datamat, params.prior, params.cache)

    while true
        best_diff = 0.0
        best_parent_list = parent_list
        for i in 1:n

            # 1) add an edge (j->i)
            if length(parent_list[i]) < params.max_n_parents
                for j in deleteat!(collect(1:n), parent_list[i])
                    if adding_edge_preserves_acyclicity(parent_list, j, i)
                        new_parents = sort!(push!(copy(parent_list[i]), j))
                        new_component_score = bayesian_score_component(i, new_parents, ncategories, datamat, params.prior, params.cache)
                        if new_component_score - score_components[i] > best_diff
                            best_diff = new_component_score - score_components[i]
                            best_parent_list = deepcopy(parent_list)
                            best_parent_list[i] = new_parents
                        end
                    end
                end
            end
        end

        for i in 1:n
            for (idx, j) in enumerate(parent_list[i])

                new_parents = deleteat!(copy(parent_list[i]), idx)
                new_component_score = bayesian_score_component(i, new_parents, ncategories, datamat, params.prior, params.cache)
                if new_component_score - score_components[i] > best_diff
                    best_diff = new_component_score - score_components[i]
                    best_parent_list = deepcopy(parent_list)
                    best_parent_list[i] = new_parents
                end
            end
        end

        if best_diff > 0.0
            parent_list = best_parent_list
            score_components = bayesian_score_components(parent_list, ncategories, datamat, params.prior, params.cache)
        else
            break
        end
    end

    # construct the BayesNet
    cpds = Array{DiscreteCPD}(n)
    varnames = propertynames(data)
    for i in 1:n
        name = varnames[i]
        parents = varnames[parent_list[i]]
        cpds[i] = fit(DiscreteCPD, data, name, parents, params.prior, parental_ncategories=ncategories[parent_list[i]], target_ncategories=ncategories[i])
    end
    BayesNet(cpds)
end


text:
import HTTP, JSON2
using Test
using Dash
struct TestChild
    a::Int32
    b::Int32
end
macro sym_str(s)
    return Symbol(s)
end
@testset "Front.to_dash" begin
    @test Front.to_dash((1,3)) == (1,3)
    Front.to_dash(t::TestChild) = "$(t.a)-$(t.b)"
    t = TestChild(10, 20)
    @test Front.to_dash(t) == "$(t.a)-$(t.b)"
end

@testset "Typed property conversion" begin
    Front.to_dash(t::TestChild) = "$(t.a)-$(t.b)"
    t = TestChild(10, 20)

    a_comp = html_a(t, id = "test-a")
    @test a_comp.props[:children] == "$(t.a)-$(t.b)"

    a_comp = html_a(children = t, id = "test-a")
    @test a_comp.props[:children] == "$(t.a)-$(t.b)"

    a_comp = html_a(id = "test-a") do
        return TestChild(10, 20)
    end
    @test a_comp.props[:children] == "10-20"

end

@testset "Typed callback result conversion" begin
    Front.to_dash(t::TestChild) = "$(t.a)-$(t.b)"
    t = TestChild(10, 20)

    app = dash("Test app") do
        html_div() do
            html_div(10, id = "my-id"),
            html_div(id = "my-div")
        end
    end
    callback!(app, callid"my-id.children => my-div.children") do value
        return TestChild(value, value)
    end

    @test length(app.callbacks) == 1
    @test haskey(app.callbacks, Symbol("my-div.children"))
    @test app.callbacks[Symbol("my-div.children")].func(10) == TestChild(10, 10)

    test_json = """{"output":"my-div.children","changedPropIds":["my-id.children"],"inputs":[{"id":"my-id","property":"children","value":10}]}"""
    result = Dash.process_callback(app, test_json)
    @test result[:response][:props][:children] == "10-10"

    app = dash("Test app") do
        html_div() do
            html_div(10, id = "my-id"),
            html_div(id = "my_div"),
            html_div(id = "my_div2")
        end
    end
    callback!(app, callid"my-id.children => my_div.children, my_div2.children") do value
        return TestChild(value, value), TestChild(value + 5, value + 5)
    end

    @test length(app.callbacks) == 1
    @test haskey(app.callbacks, Symbol("..my_div.children...my_div2.children.."))
    @test app.callbacks[Symbol("..my_div.children...my_div2.children..")].func(10) == (TestChild(10, 10), TestChild(15,15))

    test_json = """{"output":"..my_div.children...my_div2.children..","changedPropIds":["my-id.children"],"inputs":[{"id":"my-id","property":"children","value":10}]}"""
    result = Dash.process_callback(app, test_json)

    @test result[:response][:my_div][:children] == "10-10"
    @test result[:response][:my_div2][:children] == "15-15"

end

@testset "Front.from_dash" begin
    @test Front.from_dash(Any, (1,3)) == (1,3)
    function Front.from_dash(::Type{TestChild}, t)
        TestChild(parse.(Int32, split(t, "-"))...)
    end

    @test Front.from_dash(TestChild, "10-20") == TestChild(10, 20)
end


text:
# TODO Move this to GPLikelihoods.jl
@doc raw"""
    NegBinomialLikelihood(r::Real)

Likelihood with a Negative Binomial distribution:
```math
    p(y|f, r) = {y + r - 1 \choose y}(1-\sigma(f))^r \sigma^y(f).
```

## Arguments
- `r::Real`, number of failures until the experiment is stopped
"""
struct NegBinomialLikelihood{Tl<:AbstractLink,Tr<:Real} <: AbstractLikelihood
    invlink::Tl
    r::Tr
end

NegBinomialLikelihood(r::Real) = NegBinomialLikelihood(LogisticLink(), r)

(lik::NegBinomialLikelihood)(f::Real) = NegativeBinomial(lik.r, 1 - lik.invlink(f))

function (lik::NegBinomialLikelihood)(f::AbstractVector{<:Real})
    return Product(lik.(f))
end

aux_field(::NegBinomialLikelihood, Ω) = getproperty(Ω, :ω)

function init_aux_variables(rng::AbstractRNG, ::NegBinomialLikelihood, n::Int)
    return TupleVector((; ω=rand(rng, PolyaGamma(1, 0), n)))
end

function init_aux_posterior(T::DataType, lik::NegBinomialLikelihood, n::Int)
    return For(TupleVector(; y=zeros(Int, n), c=zeros(T, n))) do φ
        NTDist(PolyaGamma(φ.y + lik.r, φ.c)) # Distributions uses a different parametrization
    end
end

function aux_full_conditional(lik::NegBinomialLikelihood, y::Real, f::Real)
    return NTDist(PolyaGamma(y + lik.r, abs(f)))
end

function aux_posterior!(
    qΩ, ::NegBinomialLikelihood, y::AbstractVector, qf::AbstractVector{<:Normal}
)
    φ = only(qΩ.inds)
    φ.y .= y
    map!(φ.c, qf) do fᵢ
        sqrt(second_moment(fᵢ))
    end
    return qΩ
end

function auglik_potential(lik::NegBinomialLikelihood, ::Any, y::AbstractVector)
    return ((y .- lik.r) / 2,)
end

function auglik_precision(::NegBinomialLikelihood, Ω, ::AbstractVector)
    return (Ω.ω,)
end

function expected_auglik_potential(lik::NegBinomialLikelihood, qΩ, y::AbstractVector)
    return auglik_potential(lik, qΩ, y)
end

function expected_auglik_precision(::NegBinomialLikelihood, qΩ, ::AbstractVector)
    return (tvmean(qΩ).ω,)
end

negbin_logconst(y, r::Real) = loggamma(y + r) - loggamma(y + 1) - loggamma(r)
negbin_logconst(y, r::Int) = first(logabsbinomial(y + r - 1, y))

function logtilt(lik::NegBinomialLikelihood, ω::Real, y::Integer, f::Real)
    negbin_logconst(y, lik.r) - (y + lik.r) * logtwo +
        (f * (y - lik.r) - abs2(f) * ω) / 2
end

function expected_logtilt(
    lik::NegBinomialLikelihood, qωᵢ::NTDist{<:PolyaGamma}, yᵢ::Integer, qfᵢ::Normal
)
    θ = ntmean(qωᵢ)
    return negbin_logconst(yᵢ, lik.r) - (yᵢ + lik.r) * logtwo +
           (mean(qfᵢ) * (yᵢ - lik.r) - second_moment(qfᵢ) * θ.ω) / 2
end

function aux_prior(lik::NegBinomialLikelihood, y::AbstractVector{<:Integer})
    return For(y) do yᵢ
        NTDist(aux_prior(lik, yᵢ))
    end
end

aux_prior(lik::NegBinomialLikelihood, y::Int) = PolyaGamma(lik.r + y, 0)


text:
using Gridap
using Gridap.Inference
using Gridap.Fields
using Test

u(x) = x

p1 = Point(1,2)
p2 = Point(2,1)
p3 = Point(4,3)
p4 = Point(6,1)
x = [p1,p2,p3,p4]

f = FunctionField(u)
test_field(f,x,x)

Fill(f,4)


text:
#!/usr/bin/env julia

using PencilFFTs
using PencilFFTs.PencilArrays

import FFTW
using MPI
import Pkg

using OrderedCollections: OrderedDict
using TimerOutputs

using LinearAlgebra
using Printf
using Profile
using Random

TimerOutputs.enable_debug_timings(PencilFFTs)
TimerOutputs.enable_debug_timings(PencilArrays)
TimerOutputs.enable_debug_timings(Transpositions)

FFTW.set_num_threads(1)

const PROFILE = false
const PROFILE_OUTPUT = "profile.txt"
const PROFILE_DEPTH = 8

const MEASURE_GATHER = false

const DIMS_DEFAULT = "32,64,128"

const DEV_NULL = @static Sys.iswindows() ? "nul" : "/dev/null"

const SEPARATOR = string("\n", "*"^80)

const RESULTS_HEADER =
"""# The last 4 columns show timing statistics (mean/std/min/max) in milliseconds.
    #
    # The last two letters in the name of this file determine the PencilFFTs
    # parameters used in the benchmarks.
    #
    # The first letter indicates whether dimension permutations are performed:
    #
    #     P = permutations (default in PencilFFTs)
    #     N = no permutations
    #
    # The second letter indicates the MPI transposition method:
    #
    #     I = Isend/Irecv (default in PencilFFTs)
    #     A = Alltoallv
    #
    """

mutable struct TimerData
    avg :: Float64
    std :: Float64
    min :: Float64
    max :: Float64
    TimerData() = new(0, 0, Inf, -1)
end

function Base.:*(t::TimerData, v)
    t.avg *= v
    t.std *= v
    t.min *= v
    t.max *= v
    t
end

function getenv(::Type{T}, key, default = nothing) where {T}
    s = get(ENV, key, nothing)
    if s === nothing
        default
    elseif T <: AbstractString
        s
    else
        parse(T, s)
    end
end

getenv(key, default::T) where {T} = getenv(T, key, default)

function parse_params()
    dims_str = getenv("PENCILFFTS_BENCH_DIMENSIONS", DIMS_DEFAULT)
    repetitions = getenv("PENCILFFTS_BENCH_REPETITIONS", 100)
    outfile = getenv(String, "PENCILFFTS_BENCH_OUTPUT", nothing)
    (
        dims = parse_dimensions(dims_str) :: Dims{3},
        iterations = repetitions :: Int,
        outfile = outfile :: Union{Nothing,String},
    )
end

# Slab decomposition
function create_pencils(topo::MPITopology{1}, data_dims, permutation::Val{true};
                        kwargs...)
    pen1 = Pencil(topo, data_dims, (2, ); kwargs...)
    pen2 = Pencil(pen1, decomp_dims=(3, ), permute=Permutation(2, 1, 3); kwargs...)
    pen3 = Pencil(pen2, decomp_dims=(2, ), permute=Permutation(3, 2, 1); kwargs...)
    pen1, pen2, pen3
end

function create_pencils(topo::MPITopology{1}, data_dims,
                        permutation::Val{false}; kwargs...)
    pen1 = Pencil(topo, data_dims, (2, ); kwargs...)
    pen2 = Pencil(pen1, decomp_dims=(3, ); kwargs...)
    pen3 = Pencil(pen2, decomp_dims=(2, ); kwargs...)
    pen1, pen2, pen3
end

# Pencil decomposition
function create_pencils(topo::MPITopology{2}, data_dims, permutation::Val{true};
                        kwargs...)
    pen1 = Pencil(topo, data_dims, (2, 3); kwargs...)
    pen2 = Pencil(pen1, decomp_dims=(1, 3), permute=Permutation(2, 1, 3); kwargs...)
    pen3 = Pencil(pen2, decomp_dims=(1, 2), permute=Permutation(3, 2, 1); kwargs...)
    pen1, pen2, pen3
end

function create_pencils(topo::MPITopology{2}, data_dims, permutation::Val{false};
                        kwargs...)
    pen1 = Pencil(topo, data_dims, (2, 3); kwargs...)
    pen2 = Pencil(pen1, decomp_dims=(1, 3); kwargs...)
    pen3 = Pencil(pen2, decomp_dims=(1, 2); kwargs...)
    pen1, pen2, pen3
end

function benchmark_pencils(comm, proc_dims::Tuple, data_dims::Tuple;
                           iterations=1,
                           with_permutations::Val=Val(true),
                           extra_dims::Tuple=(),
                           transpose_method=Transpositions.PointToPoint(),
                          )
    topo = MPITopology(comm, proc_dims)
    M = length(proc_dims)
    @assert M in (1, 2)

    to = TimerOutput()

    pens = create_pencils(topo, data_dims, with_permutations, timer=to)

    u = map(p -> PencilArray{Float64}(undef, p, extra_dims...), pens)

    myrank = MPI.Comm_rank(comm)
    rng = MersenneTwister(42 + myrank)
    randn!(rng, u[1])
    u[1] .+= 10 * myrank
    u_orig = copy(u[1])

    transpose_m!(a, b) = transpose!(a, b, method=transpose_method)

    # Precompile functions
    transpose_m!(u[2], u[1])
    transpose_m!(u[3], u[2])
    transpose_m!(u[2], u[3])
    transpose_m!(u[1], u[2])
    gather(u[2])

    reset_timer!(to)

    for it = 1:iterations
        transpose_m!(u[2], u[1])
        transpose_m!(u[3], u[2])
        transpose_m!(u[2], u[3])
        transpose_m!(u[1], u[2])
        MEASURE_GATHER && gather(u[2])
    end

    @assert u[1] == u_orig

    if myrank == 0
        println("\n",
            """
            Processes:               $proc_dims
            Data dimensions:         $data_dims $(isempty(extra_dims) ? "" : "× $extra_dims")
            Permutations (1, 2, 3):  $(permutation.(pens))
            Transpositions:          1 -> 2 -> 3 -> 2 -> 1
            Method:                  $(transpose_method)
            """)
        println(to, SEPARATOR)
    end

    if PROFILE
        Profile.clear()
        @profile for it = 1:iterations
            transpose_m!(u[2], u[1])
            transpose_m!(u[3], u[2])
            transpose_m!(u[2], u[3])
            transpose_m!(u[1], u[2])
        end
        if myrank == 0
            open(io -> Profile.print(io, maxdepth=PROFILE_DEPTH),
                 PROFILE_OUTPUT, "w")
        end
    end

    nothing
end

function benchmark_rfft(comm, proc_dims::Tuple, data_dims::Tuple;
                        extra_dims=(),
                        iterations=1,
                        transpose_method=Transpositions.PointToPoint(),
                        permute_dims=Val(true),
                       )
    isroot = MPI.Comm_rank(comm) == 0

    to = TimerOutput()
    plan = PencilFFTPlan(data_dims, Transforms.RFFT(), proc_dims, comm,
                         extra_dims=extra_dims,
                         permute_dims=permute_dims,
                         fftw_flags=FFTW.ESTIMATE,
                         timer=to, transpose_method=transpose_method)

    if isroot
        println("\n", plan, "\nMethod: ", plan.transpose_method)
        println("Permutations: $permute_dims\n")
    end

    u = allocate_input(plan)
    v = allocate_output(plan)
    uprime = similar(u)

    randn!(u)

    # Warm-up
    mul!(v, plan, u)
    ldiv!(uprime, plan, v)

    @assert u ≈ uprime

    reset_timer!(to)

    t = TimerData()

    for it = 1:iterations
        τ = -MPI.Wtime()
        mul!(v, plan, u)
        ldiv!(u, plan, v)
        τ += MPI.Wtime()
        t.avg += τ
        t.std += τ^2
        t.min = min(τ, t.min)
        t.max = max(τ, t.max)
    end

    t.avg /= iterations
    t.std = sqrt(t.std / iterations - t.avg^2)

    t *= 1000  # in milliseconds

    events = (to["PencilFFTs mul!"], to["PencilFFTs ldiv!"])
    @assert all(TimerOutputs.ncalls.(events) .== iterations)
    t_to = sum(TimerOutputs.time.(events)) / iterations / 1e6  # in milliseconds

    if isroot
        @printf("Average time: %.8f ms (TimerOutputs) over %d repetitions\n",
                t_to, iterations)
        @printf("              %.8f ms (MPI_Wtime) ± %.8f ms \n\n",
                t.avg, t.std / 2)
        print_timers(to, iterations, transpose_method)
        println(SEPARATOR)
    end

    t
end

struct AggregatedTimes{TM}
    transpose :: TM
    mpi  :: Float64  # MPI time in μs
    fft  :: Float64  # FFTs in μs
    data :: Float64  # data copies in μs
    others :: Float64
end

function AggregatedTimes(to::TimerOutput, transpose_method)
    repetitions = TimerOutputs.ncalls(to)

    avgtime(x) = TimerOutputs.time(x) / repetitions
    avgtime(::Nothing) = 0.0

    fft = avgtime(to["FFT"])

    tf = TimerOutputs.flatten(to)
    data = avgtime(tf["copy_permuted!"]) + avgtime(tf["copy_range!"])

    mpi = if transpose_method === Transpositions.PointToPoint()
        t = avgtime(to["MPI.Waitall!"])
        if haskey(tf, "wait receive")  # this will be false in serial mode
            t += avgtime(tf["wait receive"])
        end
        t
    elseif transpose_method === Transpositions.Alltoallv()
        avgtime(tf["MPI.Alltoallv!"]) + avgtime(to["MPI.Waitall!"])
    end

    others = 0.0
    if haskey(to, "normalise")  # normalisation of inverse transform
        others += avgtime(to["normalise"])
    end

    let scale = 1e-6  # convert to ms
        mpi *= scale / 2   # 2 transposes per iteration
        fft *= scale / 3   # 3 FFTs per iteration
        data *= scale / 2
        others *= scale
    end

    AggregatedTimes(transpose_method, mpi, fft, data, others)
end

# 2 transpositions + 3 FFTs
time_total(t::AggregatedTimes) = 2 * (t.mpi + t.data) + 3 * t.fft + t.others

function Base.show(io::IO, t::AggregatedTimes)
    maybe_newline = ""
    for p in (string(t.transpose) => t.mpi,
              "FFT" => t.fft, "(un)pack" => t.data, "others" => t.others)
        @printf io "%s  Average %-10s = %.6f ms" maybe_newline p.first p.second
        maybe_newline = "\n"
    end
    io
end

function print_timers(to::TimerOutput, iterations, transpose_method)
    println(to, "\n")

    @assert TimerOutputs.ncalls(to["PencilFFTs mul!"]) == iterations

    t_fw = AggregatedTimes(to["PencilFFTs mul!"], transpose_method)
    t_bw = AggregatedTimes(to["PencilFFTs ldiv!"], transpose_method)

    println("Forward transforms\n", t_fw)
    println("\nBackward transforms\n", t_bw)

    t_all_measured = sum(time_total.((t_fw, t_bw)))  # in milliseconds

    # Actual time taken by parallel FFTs.
    t_all = TimerOutputs.tottime(to) / 1e6 / iterations

    # Fraction the elapsed time that is not included in t_all_measured.
    t_missing = t_all - t_all_measured
    percent_missing = (1 - t_all_measured / t_all) * 100

    @printf("\nTotal from timers: %.4f ms/iteration (%.4f ms / %.2f%% missing)\n",
            t_all_measured, t_missing, percent_missing)

    nothing
end

function parse_dimensions(arg::AbstractString) :: Dims{3}
    ints = try
        sp = split(arg, ',')
        parse.(Int, sp)
    catch e
        error("Could not parse dimensions from '$arg'")
    end
    if length(ints) == 1
        N = ints[1]
        (N, N, N)
    elseif length(ints) == 3
        ntuple(n -> ints[n], Val(3))
    else
        error("Incorrect number of dimensions in '$ints'")
    end
end

make_index(::Val{true}) = 1
make_index(::Val{false}) = 2

make_index(::Transpositions.PointToPoint) = 1
make_index(::Transpositions.Alltoallv) = 2

make_index(stuff...) = CartesianIndex(make_index.(stuff))

function run_benchmarks(params)
    comm = MPI.COMM_WORLD
    Nproc = MPI.Comm_size(comm)
    myrank = MPI.Comm_rank(comm)

    dims = params.dims
    iterations = params.iterations
    outfile = params.outfile

    if myrank == 0
        @info "Global dimensions: $dims"
        @info "Repetitions:       $iterations"
    end

    # Let MPI_Dims_create choose the decomposition.
    proc_dims = let pdims = zeros(Int, 2)
        MPI.Dims_create!(Nproc, pdims)
        pdims[1], pdims[2]
    end

    transpose_methods = (Transpositions.PointToPoint(),
                         Transpositions.Alltoallv())
    permutes = (Val(true), Val(false))
    timings = Array{TimerData}(undef, 2, 2)

    map(Iterators.product(transpose_methods, permutes)) do (method, permute)
        I = make_index(permute, method)
        timings[I] = benchmark_rfft(
            comm, proc_dims, dims;
            iterations = iterations, permute_dims = permute, transpose_method = method,
        )
    end

    columns = OrderedDict{String,Union{Int,Float64}}(
        "(1) Nx" => dims[1],
        "(2) Ny" => dims[2],
        "(3) Nz" => dims[3],
        "(4) num_procs" => Nproc,
        "(5) P1" => proc_dims[1],
        "(6) P2" => proc_dims[2],
        "(7) repetitions" => iterations,
    )

    cases = (
        :PI => make_index(Val(true), Transpositions.PointToPoint()),
        :PA => make_index(Val(true), Transpositions.Alltoallv()),
        :NI => make_index(Val(false), Transpositions.PointToPoint()),
        :NA => make_index(Val(false), Transpositions.Alltoallv()),
    )

    if myrank == 0 && outfile !== nothing
        for (name, ind) in cases
            a, b = splitext(outfile)
            fname = string(a, "_", name, b)
            write_results(fname, columns, timings[ind])
        end
    end

    nothing
end

function write_results(outfile, columns, t)
    @info "Writing to $outfile"
    newfile = !isfile(outfile)
    open(outfile, "a") do io
        if newfile
            print(io, RESULTS_HEADER, "#")
            n = length(columns)
            mkname(c, name) = "($(n + c)) $name"
            names = Iterators.flatten((
                keys(columns),
                (mkname(1, "mean"), mkname(2, "std"),
                 mkname(3, "min"), mkname(4, "max"))
            ))
            for name in names
                print(io, "  ", name)
            end
            println(io)
        end
        vals = Iterators.flatten((values(columns), t.avg, t.std, t.min, t.max))
        for val in vals
            print(io, "  ", val)
        end
        println(io)
    end
    nothing
end

MPI.Init()
if MPI.Comm_rank(MPI.COMM_WORLD) == 0
    Pkg.status(mode = Pkg.PKGMODE_MANIFEST)
end
params = parse_params()
run_benchmarks(params)