text:
module LinearRegression

using NamedArrays:length
using LinearAlgebra:length
using Distributions:length
export regress, predict_and_stats, linRegRes

using Base: Tuple, Int64
using StatsBase:eltype, isapprox, length, coefnames, push!
using Distributions
using Printf, NamedArrays
using StatsBase
using StatsModels, DataFrames

include("sweep_operator.jl")
include("utilities.jl")
include("newey_west.jl")

"""
    struct linRegRes

    Store results from the regression
"""
struct linRegRes
    extended_inverse::Matrix                # Store the extended inverse matrix
    coefs::Union{Nothing,Vector}            # Store the coefficients of the fitted model
    white_types::Union{Nothing,Vector}      # Store the type of White's covariance estimator(s) used
    hac_types::Union{Nothing,Vector}        # Store the type of White's covariance estimator(s) used
    stderrors::Union{Nothing,Vector}        # Store the standard errors for the fitted model
    white_stderrors::Union{Nothing,Vector}  # Store the standard errors modified for the White's covariance estimator
    hac_stderrors::Union{Nothing,Vector}    # Store the standard errors modified for the Newey-West covariance estimator
    t_values::Union{Nothing,Vector}         # Store the t values for the fitted model
    white_t_values::Union{Nothing,Vector}   # Store the t values modified for the White's covariance estimator
    hac_t_values::Union{Nothing,Vector}     # Store the t values modified for the Newey-West covariance estimator
    p::Int64                                # Store the number of parameters (including the intercept as a parameter)
    MSE::Union{Nothing,Float64}             # Store the Mean squared error for the fitted model
    intercept::Bool                         # Indicate if the model has an intercept
    R2::Union{Nothing,Float64}              # Store the R-squared value for the fitted model
    ADJR2::Union{Nothing,Float64}           # Store the adjusted R-squared value for the fitted model
    RMSE::Union{Nothing,Float64}            # Store the Root mean square error for the fitted model
    AIC::Union{Nothing,Float64}             # Store the Akaike information criterion for the fitted model
    σ̂²::Union{Nothing,Float64}              # Store the σ̂² for the fitted model
    p_values::Union{Nothing,Vector}         # Store the p values for the fitted model
    white_p_values::Union{Nothing,Vector}   # Store the p values modified for the White's covariance estimator
    hac_p_values::Union{Nothing,Vector}     # Store the p values modified for the Newey-West covariance estimator
    ci_up::Union{Nothing,Vector}            # Store the upper values confidence interval of the coefficients
    ci_low::Union{Nothing,Vector}           # Store the lower values confidence interval of the coefficients
    white_ci_up::Union{Nothing,Vector}      # Store the upper values confidence interval of the coefficients for White covariance estimators
    white_ci_low::Union{Nothing,Vector}     # Store the upper values confidence interval of the coefficients for White covariance estimators
    hac_ci_up::Union{Nothing,Vector}        # Store the upper values confidence interval of the coefficients for Newey-West covariance estimators
    hac_ci_low::Union{Nothing,Vector}       # Store the upper values confidence interval of the coefficients for Newey-West covariance estimators
    observations                            # Store the number of observations used in the model
    t_statistic::Union{Nothing,Float64}     # Store the t statistic
    VIF::Union{Nothing,Vector}              # Store the Varince inflation factor
    modelformula                            # Store the model formula
    dataschema                              # Store the dataschema
    updformula                              # Store the updated model formula (after the dataschema has been applied)
    alpha                                   # Store the alpha used to compute the confidence interval of the coefficients
end

"""
    function Base.show(io::IO, lr::linRegRes)

    Display information about the fitted model
"""
function Base.show(io::IO, lr::linRegRes)
    println(io, "Model definition:\t", lr.modelformula)
    println(io, "Used observations:\t", lr.observations)
    println(io, "Model statistics:")
    # Display stats when available
    if !isnothing(lr.R2) && !isnothing(lr.ADJR2)
        @printf(io, "  R²: %g\t\t\tAdjusted R²: %g\n", lr.R2, lr.ADJR2)
    elseif !isnothing(lr.R2)
        @printf(io, "  R²: %g\n", lr.R2)
    end

    if !isnothing(lr.MSE) && !isnothing(lr.RMSE)
        @printf(io, "  MSE: %g\t\t\tRMSE: %g\n", lr.MSE, lr.RMSE)
    elseif !isnothing(lr.MSE)
        @printf(io, "  MSE: %g\n", lr.MSE)
    end

    if length(lr.white_types) + length(lr.hac_types) == 0
        if !isnothing(lr.σ̂²) && !isnothing(lr.AIC)
            @printf(io, "  σ̂²: %g\t\t\tAIC: %g\n", lr.σ̂², lr.AIC)
        elseif !isnothing(lr.σ̂²)
            @printf(io, "  σ̂²: %g\n", lr.σ̂²)
        elseif !isnothing(lr.AIC)
            @printf(io, "  AIC: %g\n", lr.AIC)
        end
    end

    if !isnothing(lr.ci_low) || !isnothing(lr.ci_up)
        @printf(io, "Confidence interval: %g%%\n", (1 - lr.alpha) * 100 )
    end

    if length(lr.white_types) + length(lr.hac_types) == 0
        helper_print_table(io, "Coefficients statistics:",
            [lr.coefs, lr.stderrors, lr.t_values, lr.p_values, lr.ci_low, lr.ci_up, lr.VIF],
            ["Coefs", "Std err", "t", "Pr(>|t|)", "low ci", "high ci", "VIF"],
            lr.updformula)
    end

    if length(lr.white_types) > 0
        for (cur_i, cur_type) in enumerate(lr.white_types)
            helper_print_table(io, "White's covariance estimator ($(Base.Unicode.uppercase(string(cur_type)))):",
                [lr.coefs, lr.white_stderrors[cur_i], lr.white_t_values[cur_i], lr.white_p_values[cur_i], lr.white_ci_low[cur_i], lr.white_ci_up[cur_i] ],
                ["Coefs", "Std err", "t", "Pr(>|t|)", "low ci", "high ci"],
                lr.updformula)
        end
    end

    if length(lr.hac_types) > 0
        for (cur_i, cur_type) in enumerate(lr.hac_types)
            helper_print_table(io, "Newey-West's covariance estimator:",
                [lr.coefs, lr.hac_stderrors[cur_i], lr.hac_t_values[cur_i], lr.hac_p_values[cur_i], lr.hac_ci_low[cur_i], lr.hac_ci_up[cur_i] ],
                ["Coefs", "Std err", "t", "Pr(>|t|)", "low ci", "high ci"],
                lr.updformula)
        end
    end
end

"""
    function getVIF(x, intercept, updf, df, p)

    (internal) Calculates the VIF, Variance Inflation Factor, for a given regression.
    When the has an intercept use the simplified formula. When there is no intercept use the classical formula.
"""
function getVIF(x, intercept, updf, df, p)
    if intercept
        if p == 1
            return [0., 1.]
        end
        return vcat(0, diag(inv(cor(@view(x[:, 2:end])))))
    else
        if p == 1
            return [0.]
        end
        return diag(inv(cor(x)))
    end
end

"""
    function getSST(y, intercept)

    (internal) Calculates "total sum of squares" see link for description.
    https://en.wikipedia.org/wiki/Total_sum_of_squares
    When the mode has no intercept the SST becomes the sum of squares of y
"""
function getSST(y, intercept)
    SST = zero(eltype(y))
    if intercept
        ȳ = mean(y)
        SST = sum(abs2.(y .- ȳ))
    else
        SST = sum(abs2.(y))
    end

    return SST
end

"""
    lr_predict(xs, coefs, intercept::Bool)

    (internal) Predict the ŷ given the x(s) and the coefficients of the linear regression.
"""
function lr_predict(xs, coefs, intercept::Bool)
    if intercept
        return muladd(@view(xs[:, 2:end]), @view(coefs[2:end]), coefs[1])
    else
        return muladd(xs, coefs, zero(eltype(coefs)))
    end
end

"""
    function hasintercept(f::StatsModels.FormulaTerm)
    (internal) return a tuple with the first item being true when the formula has an intercept term, the second item being the potentially updated formula.
    If there is no intercept indicated add one.
    If the intercept is specified as absent (y ~ 0 + x) then do not change.
"""
function hasintercept(f::StatsModels.FormulaTerm)
    intercept = true
    if f.rhs isa ConstantTerm{Int64}
        intercept = convert(Bool, f.rhs.n)
        return intercept, f
    elseif f.rhs isa Tuple
        for t in f.rhs
            if t isa ConstantTerm{Int64}
                intercept = convert(Bool, t.n)
                return intercept, f
            end
        end
    end
    f = FormulaTerm(f.lhs, InterceptTerm{true}() + f.rhs)
    return intercept, f
end

"""
    function regress(f::StatsModels.FormulaTerm, df::DataFrames.DataFrame; α::Float64=0.05, req_stats=["all"], remove_missing=false, cov=[:none])

    Estimate the coefficients of the regression, given a dataset and a formula.

    The formula details are provided in the StatsModels package and the behaviour aims to be similar as what the Julia GLM package provides.
    The data shall be provided as a DataFrame without missing data.
    If remove_missing is set to true a copy of the dataframe will be made and the row with missing data will be removed.
"""
function regress(f::StatsModels.FormulaTerm, df::DataFrames.DataFrame; α::Float64=0.05, req_stats=["all"], remove_missing=false, cov=[:none])
    intercept, f= hasintercept(f)

    (α > 0. && α < 1.) || throw(ArgumentError("α must be between 0 and 1"))

    copieddf = df
    if remove_missing
        copieddf = copy(df[: , Symbol.(keys(schema(f, df).schema))])
        dropmissing!(copieddf)
    end

    needed_stats = get_needed_model_stats(req_stats)

    dataschema = schema(f, copieddf)
    updatedformula = apply_schema(f, dataschema)

    y, x = modelcols(updatedformula, copieddf)
    n, p = size(x)
    xy = [x y]

    xytxy = xy' * xy

# mandatory stats
    sse = sweep_op_full!(xytxy)[end]
    coefs = xytxy[1:p, end]
    mse = xytxy[p + 1, p + 1] / (n - p)

# optional stats
    total_scalar_stats = Set([:sse, :mse, :sst, :r2, :adjr2, :rmse, :aic, :sigma, :t_statistic, :vif])
    total_vector_stats = Set([:coefs, :stderror, :t_values, :p_values, :ci])

    scalar_stats = Dict{Symbol,Union{Nothing,Float64}}(intersect(total_scalar_stats, needed_stats) .=> nothing)
    vector_stats = Dict{Symbol,Union{Nothing,Vector}}(intersect(total_vector_stats, needed_stats) .=> nothing)

# optional stats
    if :sst in needed_stats
        scalar_stats[:sst] = getSST(y, intercept)
    end
    if :r2 in needed_stats
        scalar_stats[:r2] = 1. - (sse / scalar_stats[:sst])
    end
    if :adjr2 in needed_stats
        scalar_stats[:adjr2] = 1. - ((n - convert(Int64, intercept)) * (1. - scalar_stats[:r2])) / (n - p)
    end
    if :rmse in needed_stats
        scalar_stats[:rmse] = √mse
    end
    if :aic in needed_stats
        scalar_stats[:aic] =  n * log(sse / n) + 2p
    end
    if :sigma in needed_stats
        scalar_stats[:sigma] = mse
    end
    if :t_statistic in needed_stats
        scalar_stats[:t_statistic] = quantile(TDist(n - p), 1 - α / 2)
    end
    if :stderror in needed_stats
        vector_stats[:stderror] = sqrt.(diag(mse * @view(xytxy[1:end - 1, 1:end - 1])))
    end
    if :t_values in needed_stats
        vector_stats[:t_values] = coefs ./ vector_stats[:stderror]
    end
    if :p_values in needed_stats
        vector_stats[:p_values] = ccdf.(Ref(FDist(1., (n - p))), abs2.(vector_stats[:t_values]))
    end
    if :ci in needed_stats
        vector_stats[:ci] = vector_stats[:stderror] * scalar_stats[:t_statistic]
    end
    if :vif in needed_stats
        vector_stats[:vif] = getVIF(x, intercept, updatedformula, copieddf, p)
    end

# robust estimators stats
    needed_white, needed_hac = get_needed_robust_cov_stats(cov)
    white_types = Vector{Symbol}()
    white_stds = Vector{Vector}()
    white_t_vals = Vector{Vector}()
    white_p_vals = Vector{Vector}()
    white_ci_up = Vector{Vector}()
    white_ci_low = Vector{Vector}()

    if length(needed_white) > 0
        for t in needed_white
            if t in white_types
                continue
            end
            cur_type, cur_std = heteroscedasticity(t, x, y, coefs, intercept, n, p, xytxy)
            push!(white_types, cur_type)
            push!(white_stds, cur_std)

            if !isnothing(get(vector_stats, :t_values, nothing))
                cur_t_vals = coefs ./ cur_std
                push!(white_t_vals, cur_t_vals)
            else
                white_t_vals = nothing
            end

            if !isnothing(get(vector_stats, :p_values, nothing))
                cur_p_vals = ccdf.(Ref(FDist(1., (n - p))), abs2.(cur_t_vals))
                push!(white_p_vals, cur_p_vals)
            else
                white_p_vals = nothing
            end
            if !isnothing(get(vector_stats, :ci, nothing))
                cur_ci = cur_std * scalar_stats[:t_statistic]
                cur_ci_up = coefs .+ cur_ci
                cur_ci_low = coefs .- cur_ci
                push!(white_ci_up, cur_ci_up)
                push!(white_ci_low, cur_ci_low)
            else
                white_ci_up = nothing
                white_ci_low = nothing
            end
        end
    end

    hac_types = Vector{Symbol}()
    hac_stds = Vector{Vector}()
    hac_t_vals = Vector{Vector}()
    hac_p_vals = Vector{Vector}()
    hac_ci_up = Vector{Vector}()
    hac_ci_low = Vector{Vector}()

    if length(needed_hac) > 0
        for t in needed_hac
            if t in hac_types
                continue
            end
            cur_type, cur_std = HAC(t, x, y, coefs, intercept, n, p)
            push!(hac_types, cur_type)
            push!(hac_stds, cur_std)

            if !isnothing(get(vector_stats, :t_values, nothing))
                cur_t_vals = coefs ./ cur_std
                push!(hac_t_vals, cur_t_vals)
            else
                hac_t_vals = nothing
            end
            if !isnothing(get(vector_stats, :p_values, nothing))
                cur_p_vals = ccdf.(Ref(FDist(1., (n - p))), abs2.(cur_t_vals))
                push!(hac_p_vals, cur_p_vals)
            else
                hac_p_vals = nothing
            end
            if !isnothing(get(vector_stats, :ci, nothing))
                cur_ci = cur_std * scalar_stats[:t_statistic]
                cur_ci_up = coefs .+ cur_ci
                cur_ci_low = coefs .- cur_ci
                push!(hac_ci_up, cur_ci_up)
                push!(hac_ci_low, cur_ci_low)
            else
                hac_ci_up = nothing
                hac_ci_low = nothing
            end
        end
    end

    sres = linRegRes(xytxy, coefs,
        white_types, hac_types,
        get(vector_stats, :stderror, nothing), white_stds, hac_stds,
        get(vector_stats, :t_values, nothing), white_t_vals, hac_t_vals,
        p, mse, intercept, get(scalar_stats, :r2, nothing),
        get(scalar_stats, :adjr2, nothing), get(scalar_stats, :rmse, nothing), get(scalar_stats, :aic, nothing), get(scalar_stats, :sigma, nothing),
        get(vector_stats, :p_values, nothing), white_p_vals, hac_p_vals,
        haskey(vector_stats, :ci) ? coefs .+ vector_stats[:ci] : nothing,
        haskey(vector_stats, :ci) ? coefs .- vector_stats[:ci] : nothing,
        white_ci_up, white_ci_low,
        hac_ci_up, hac_ci_low,
        n, get(scalar_stats, :t_statistic, nothing), get(vector_stats, :vif, nothing), f, dataschema, updatedformula, α)

    return sres
end

"""
    function HAC(t::Symbol, x, y, coefs, intercept, n, p)

    (Internal) Return the relevant HAC (heteroskedasticity and autocorrelation consistent) estimator.
    In the current version only Newey-West is implemented.
"""
function HAC(t::Symbol, x, y, coefs, intercept, n, p)
    inv_xtx = inv(x' * x)
        e = y - lr_predict(x, coefs, intercept)
    xe = x .* e
    return (t, sqrt.(diag(n * inv_xtx * newey_west(xe) * inv_xtx)))
end

"""
    function heteroscedasticity(t::Symbol, x, y, coefs, intercept, n, p)

    (Internal) Compute the standard errors modified for the White's covariance estimator.
    Currently support HC0, HC1, HC2 and HC3. When :white is passed, select HC3 when the number of observation is below 250 otherwise select HC0.
"""
function heteroscedasticity(t::Symbol, x, y, coefs, intercept, n, p, xytxy)
    inv_xtx = inv(x' * x)
    XX = @view(xytxy[1:end - 1, 1:end - 1])
    e = y - lr_predict(x, coefs, intercept)
    xe = x .* e

    if t == :white && n < 250
        t = :hc3
    elseif t == :white && n >= 250
        t = :hc0
    end

    if t == :hc0
        xetxe = xe' * xe
        return (:hc0, sqrt.(diag(XX * xetxe * XX)))
    elseif t == :hc1
        scale = (n / (n - p))
        xetxe = xe' * xe
        return (:hc1, sqrt.(diag(XX * xetxe * XX .* scale)))
    elseif t == :hc2
        leverage = diag(x * inv(x'x) * x')
        scale = @.( 1. / (1. - leverage))
        xe = @.(xe .* real(sqrt(Complex(scale))))
        xetxe = xe' * xe
        return (t, sqrt.(diag(XX * xetxe * XX)))
    elseif t == :hc3
        leverage = diag(x * inv(x'x) * x')
        scale = @.( 1. / (1. - leverage)^2)
        xe = @.(xe .* real(sqrt(Complex(scale))))
        xetxe = xe' * xe
        return (t, sqrt.(diag(XX * xetxe * XX)))
    else
        throw(error("Unknown symbol ($(t)) used as the White's covariance estimator"))
    end

end

"""
    function predict_and_stats(lr::linRegRes, df::DataFrames.DataFrame, α=0.05)

    Using the estimated coefficients from the regression make predictions, and calculate related statistics.
"""
function predict_and_stats(lr::linRegRes, df::DataFrames.DataFrame; α=0.05, req_stats=["none"])

    copieddf = df[: , Symbol.(keys(schema(lr.modelformula, df).schema))]
    dropmissing!(copieddf)
    dataschema = schema(lr.modelformula, copieddf)
    updatedformula = apply_schema(lr.modelformula, dataschema)
    y, x = modelcols(updatedformula, copieddf)
    n, p = size(x)

    needed, present = get_prediction_stats(req_stats)
    needed_stats = Dict{Symbol,Vector}()
    for sym in needed
        needed_stats[sym] = zeros(length(n))
    end
    if :leverage in needed
        needed_stats[:leverage] = diag(x * inv(x'x) * x')
    end
    if :predicted in needed
        needed_stats[:predicted] = lr_predict(x, lr.coefs, lr.intercept)
    end
    if :residuals in needed
        needed_stats[:residuals] = y .- needed_stats[:predicted]
    end
    if :stdp in needed
        if isnothing(lr.σ̂²)
            throw(ArgumentError(":stdp requires that the σ̂² (:sigma) was previously calculated through the regression"))
        end
        if length(lr.white_types) + length(lr.hac_types) > 0
            println(io, "The STDP statistic that relies on Sigma^2 has been requested. At least one robust covariance have been requested indicating that the assumptions needed for Sigma^2 may not be present.")
        end
        needed_stats[:stdp] = sqrt.(needed_stats[:leverage] .* lr.σ̂²)
    end
    if :stdi in needed
        if isnothing(lr.σ̂²)
            throw(ArgumentError(":stdi requires that the σ̂² (:sigma) was previously calculated through the regression"))
        end
        if length(lr.white_types) + length(lr.hac_types) > 0
            println(io, "The STDI statistic that relies on Sigma^2 has been requested. At least one robust covariance have been requested indicating that the assumptions needed for Sigma^2 may not be present.")
        end
        needed_stats[:stdi] = sqrt.((1. .+ needed_stats[:leverage]) .* lr.σ̂²)
    end
    if :stdr in needed
        if isnothing(lr.σ̂²)
            throw(ArgumentError(":stdr requires that the σ̂² (:sigma) was previously calculated through the regression"))
        end
        if length(lr.white_types) + length(lr.hac_types) > 0
            println(io, "The STDR statistic that relies on Sigma^2 has been requested. At least one robust covariance have been requested indicating that the assumptions needed for Sigma^2 may not be present.")
        end
        needed_stats[:stdr] = sqrt.((1. .- needed_stats[:leverage]) .* lr.σ̂²)
    end
    if :student in needed
        if length(lr.white_types) + length(lr.hac_types) > 0
            println(io, "The student statistic that relies on Sigma^2 has been requested. At least one robust covariance have been requested indicating that the assumptions needed for Sigma^2 may not be present.")
        end
        needed_stats[:student] = needed_stats[:residuals] ./ needed_stats[:stdr]
    end
    if :rstudent in needed
        if length(lr.white_types) + length(lr.hac_types) > 0
            println(io, "The rstudent statistic that relies on Sigma^2 has been requested. At least one robust covariance have been requested indicating that the assumptions needed for Sigma^2 may not be present.")
        end
        needed_stats[:rstudent] = needed_stats[:student] .* real.(sqrt.(complex.((n .- p .- 1 ) ./ (n .- p .- needed_stats[:student].^2 ), 0)))
    end
    if :lcli in needed
        if length(lr.white_types) + length(lr.hac_types) > 0
            println(io, "The LCLI statistic that relies on Sigma^2 has been requested. At least one robust covariance have been requested indicating that the assumptions needed for Sigma^2 may not be present.")
        end
        needed_stats[:lcli] = needed_stats[:predicted] .- (lr.t_statistic .* needed_stats[:stdi])
    end
    if :ucli in needed
        if length(lr.white_types) + length(lr.hac_types) > 0
            println(io, "The UCLI statistic that relies on Sigma^2 has been requested. At least one robust covariance have been requested indicating that the assumptions needed for Sigma^2 may not be present.")
        end
        needed_stats[:ucli] = needed_stats[:predicted] .+ (lr.t_statistic .* needed_stats[:stdi])
    end
    if :lclp in needed
        if length(lr.white_types) + length(lr.hac_types) > 0
            println(io, "The LCLP statistic that relies on Sigma^2 has been requested. At least one robust covariance have been requested indicating that the assumptions needed for Sigma^2 may not be present.")
        end
        needed_stats[:lclp] = needed_stats[:predicted] .- (lr.t_statistic .* needed_stats[:stdp])
    end
    if :uclp in needed
        if length(lr.white_types) + length(lr.hac_types) > 0
            println(io, "The UCLP statistic that relies on Sigma^2 has been requested. At least one robust covariance have been requested indicating that the assumptions needed for Sigma^2 may not be present.")
        end
        needed_stats[:uclp] = needed_stats[:predicted] .+ (lr.t_statistic .* needed_stats[:stdp])
    end
    if :press in needed
        needed_stats[:press] = needed_stats[:residuals] ./ (1. .- needed_stats[:leverage])
    end
    if :cooksd in needed
        if length(lr.white_types) + length(lr.hac_types) > 0
            println(io, "The CooksD statistic that relies on Sigma^2 has been requested. At least one robust covariance have been requested indicating that the assumptions needed for Sigma^2 may not be present.")
        end
        needed_stats[:cooksd] = needed_stats[:stdp].^2 ./  needed_stats[:stdr].^2 .* needed_stats[:student].^2 .* (1 / lr.p)
    end

    for sym in present
        copieddf[!, sym] = needed_stats[sym]
end

    return copieddf
end

end # end of module definition


text:
include("euler/euler.jl")
using .BigIntegers: BigInteger

function compute(d::Int, n::Int)::Int
    numerator, denominator = BigInteger(1), BigInteger(1)
    a, b = d, d - 1
    result = 0
    for _ ∈ 1:n
        numerator, denominator = numerator + a * denominator, numerator + b * denominator
        if length(numerator) ≠ length(denominator) result += 1 end
    end
    return result
end

ext: jl

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
function area_balance(
    psi_container::PSIContainer,
    expression::Symbol,
    area_mapping::Dict{String, Array{PSY.Bus, 1}},
    branches,
)
    time_steps = model_time_steps(psi_container)
    remove_undef!(psi_container.expressions[expression])
    nodal_net_balance = psi_container.expressions[expression]
    constraint_bal = JuMPConstraintArray(undef, keys(area_mapping), time_steps)
    participation_assignment_up = JuMPConstraintArray(undef, keys(area_mapping), time_steps)
    participation_assignment_dn = JuMPConstraintArray(undef, keys(area_mapping), time_steps)
    assign_constraint!(psi_container, "area_dispatch_balance", constraint_bal)
    area_balance = get_variable(psi_container, ActivePowerVariable, PSY.Area)
    for (k, buses_in_area) in area_mapping
        for t in time_steps
            area_net = model_has_parameters(psi_container) ? zero(PGAE) : JuMP.AffExpr(0.0)
            for b in buses_in_area
                JuMP.add_to_expression!(area_net, nodal_net_balance[PSY.get_number(b), t])
            end
            constraint_bal[k, t] =
                JuMP.@constraint(psi_container.JuMPmodel, area_balance[k, t] == area_net)
        end
    end

    expr_up = get_expression(psi_container, :emergency_up)
    expr_dn = get_expression(psi_container, :emergency_dn)

    assign_constraint!(
        psi_container,
        "participation_assignment_up",
        participation_assignment_up,
    )
    assign_constraint!(
        psi_container,
        "participation_assignment_dn",
        participation_assignment_dn,
    )

    for area in keys(area_mapping), t in time_steps
        participation_assignment_up[area, t] =
            JuMP.@constraint(psi_container.JuMPmodel, expr_up[area, t] == 0)
        participation_assignment_dn[area, t] =
            JuMP.@constraint(psi_container.JuMPmodel, expr_dn[area, t] == 0)
    end

    return
end


text:
export transform, transform!

transform(x::Element, tm::Array{T,2}) where T<:Number = _transform(x,tm)
transform(x::Group, tm::Array{T,2}) where T<:Number = _transform(x,tm)
function _transform(x::Element, tm::Array{T,2}) where T<:Number
  newx = deepcopy(x)
  transform!(newx,tm)
  return newx
end

transform!(::Element, ::Array{T,2}) where T<:Number = nothing
function transform!(x::Array{T,1}, tm::Array{Float64,2}) where T<:Element
  for i in eachindex(x)
    transform!(x[i],tm)
  end
end
function transform!(element::Element, tm::Array{Float64,2},
                    context::Context = Context(element))
  for e in traverseTree(element)
    _transform!(e,tm,context)
  end
end
function tmscaletofirstunits(tm::Array{Float64,2},
                             element::Element,
                             context::Context)
  scale = scaletofirstunits(element,context)
  invscale = 1/scale
  scalexyz(invscale,invscale,invscale)*tm*scalexyz(scale,scale,scale)
end

_transform!(::Element, ::Array{Float64,2}, ::Context) = nothing
function _transform!(n::Node, tm::Array{Float64,2})
  n.xyz1[1:4] = tm*n.xyz1
  return nothing
end
function _transform!(n::Node, tm::Array{Float64,2}, context::Context)
  _transform!(n,tmscaletofirstunits(tm,n,context))
  return nothing
end
function _transform!(s::Segment, tm::Array{Float64,2}, context::Context)
  #assumes nodes are transformed as part of a group
  s.wxwywz.xyz[1:3] = tm[1:3,1:3]*s.wxwywz.xyz
  s.wxwywz.isdefault = false
  return nothing
end
function _transform!(x::UniformPlane, tm::Array{Float64,2}, context::Context)
  tmscale = tmscaletofirstunits(tm,x,context)
  x.corner1[1:4] = tmscale*x.corner1
  x.corner2[1:4] = tmscale*x.corner2
  x.corner3[1:4] = tmscale*x.corner3
  for i in eachindex(x.nodes)
    _transform!(x.nodes[i],tmscale)
  end
  for i in eachindex(x.holes)
    _transform!(x.holes[i],tmscale)
  end
  return nothing
end
function _transform!(x::Point, tm::Array{Float64,2})
  x.xyz1[1:4] = tm*x.xyz1
  return nothing
end
function _transform!(x::Rect, tm::Array{Float64,2})
  x.corner1[1:4] = tm*x.corner1
  x.corner2[1:4] = tm*x.corner2
  return nothing
end
function _transform!(x::Circle, tm::Array{Float64,2})
  x.xyz1[1:4] = tm*x.xyz1
  return nothing
end

# for creating Groups of elements
function transform(e::Element, transformationlist)
  elementarraygroup(e, transformationlist)
end
function transform(e::Group, transformationlist)
  newgroup = elementarraygroup(e, transformationlist)
  newterms = terms(newgroup)
  for key in keys(terms(elements(newgroup)[1]))
    newterms[key] = Array{Node}(undef,0)
  end
  for subgroup in elements(newgroup)
    for (key,value) in terms(subgroup)
      addterm!(newterms, key, value)
    end
  end
  return newgroup
end
function elementarraygroup(e::Element, transformationlist)
  newgroup = Group()
  newelements = elements(newgroup)
  for tm in transformationlist
    push!(newelements,transform(e,tm))
  end
  return newgroup
end
function addterm!(newterms::Dict{Symbol,Union{Node,Array{Node,1}}},
                           key::Symbol, value::Array{Node,1})
  append!(newterms[key],value)
  return nothing
end
function addterm!(newterms::Dict{Symbol,Union{Node,Array{Node,1}}},
                           key::Symbol, value::Node)
  push!(newterms[key],value)
  return nothing
end

"""
    transform(element, transformation_matrix)
    transform!(element, transformation_matrix)

Transform (rotate, translate, scale, etc...) a element by 4x4
[transform matrix](# http://www.cs.mtu.edu/~shene/COURSES/cs3621/NOTES/geometry/geo-tran.html)
 `tm`.

 Transform will modify the coordinates of `Element`s and the `wx`, `wy`, and `wz`
 parameters of `Segment`.  `transform` of a `Segment` will not modify its `Node`s.
 `transform` a `Group` containing the `Node`s and `Segment` instead.  Typically
 `transform` would only be applied to `Group` objects.

`transform(element, transformation_matrix_list)`

Return a `Group` of copies of `element` transformed by each matrix in
`transformation_matrix_list`.

If the `element` is a `Group`, the terms of the returned group has the same
keys as `element`, but the values are arrays of `Nodes` of all the copies of
elements.  In other words, the elements are electricaly in parallel.
"""
transform, transform!


text:
@testset "nls_solver.jl" begin

    @testset "memory allocs" begin

        n = 10
        X = Float64[1:n;]
        Y = rand(n)

        model = Gaussian_Peak()

        θ = Float64[1,0,1]

        nls = NLS_ForwardDiff_From_Model2Fit(model,X,Y)

        @test NLS_Solver.parameter_size(nls) == 3
        @test (@benchmark(NLS_Solver.eval_r($nls,$θ))).allocs == 1
        @test (@benchmark(NLS_Solver.eval_r_J($nls,$θ))).allocs == 9

    end

end


text:
using Documenter
using ClickHouse

makedocs(
    modules = [ClickHouse],
    sitename = "ClickHouse.jl",
    pages = [
        "index.md",
        "usage.md",
        "api.md",
    ]
)


text:
#=
0: Optimal solution found
-1: Time limit exceeded
-2: NAN occured
1: Iteration limit reached
2: Convergence to an infeasible point
4: Current solution estimate cannot be improved
5: Current solution estimate cannot be improved. It appears to be optimal,
but desired accuracy could not be achieved.

0: Optimal solution found
-1: Time limit exceeded
1: Maximal number of iterations exceeded
11: NaN or Inf occured
12: Number of iterations exceeded in restoration phase
16: Point is (almost) feasible, but restoration phase is called
17: Convergence to stationary point for infeasibility
18: Restoration phase cannot further improve feasibility
=#
function real_problems(problem)
    if problem["origin"] == "real" || problem["origin"] == "modelling"
        return true
    else
        return false
    end
end

function lrg_problem_func(min_num::Int64, max_num::Int64)
    function func(problem)
        regular = problem["derivative_order"] >= 2 && problem["regular"] == true && problem["constraints"]["number"] >= 1
        min_size_ok = problem["variables"]["number"] >= min_num && problem["constraints"]["number"] >= min_num
        #max_size = problem["variables"]["number"] <= max_num && problem["constraints"]["number"] <= max_num
        max_size_ok = problem["constraints"]["number"] + problem["variables"]["number"] <= max_num

        if min_size_ok && max_size_ok && regular
            return true
        else
          return false
        end
    end

    return func
end

function jump_status_conversion(info::problem_summary2; MAX_IT::Int64=3000)
  status = info.status
  if status == :Optimal
    return :Optimal
  elseif status == :Infeasible
    return :primal_infeasible
  elseif status == :Unbounded
    return :dual_infeasible
  elseif status == :UserLimit
    if info.it_count >= MAX_IT
      return :MAX_IT
    else
      return :MAX_TIME
    end
  elseif status == :Error
    if info.it_count > 0
      return :ERROR
    else
      return :INIT_ERROR
    end
  else
    return :UNKNOWN
  end
end

function convert_JuMP(results::Dict{String,problem_summary2})
    for (prb_name, info) in results
      info.status = jump_status_conversion(info)
    end

    return results
end

function status_conversion(status::Int64)
    if status == 0
      return :Optimal
    elseif status == 17 || status == 2
      return :primal_infeasible
    #elseif status ==
    #  return :dual_infeasible
    elseif status == 1
      return :MAX_IT
    elseif status == -1
      return :MAX_TIME
    else
      return :ERROR
    end
end

function parse_number(something)
  if isa(something,Number)
      return something
  else
      return -1
  end
end

function load_results(file_name::String)
    results = readdlm(file_name,' ')
    results_dict = Dict{String, problem_summary2}()
    for i = 1:size(results,1)
      problem_name = results[i,1]
      results_dict[problem_name] = problem_summary2()

      results_dict[problem_name].it_count = parse_number(results[i,4])
      results_dict[problem_name].total_time = parse_number(results[i,7])
      results_dict[problem_name].status = status_conversion(parse_number(results[i,10]))
    end
    return results_dict
end

Infty = 99999999

function alg_success(status::Symbol)
    return (status == :Optimal  || status == :primal_infeasible || status == :dual_infeasible)
end

function iteration_list(results::Dict{String, problem_summary2}; MAX_IT::Int64=3000)
  problem_list = collect(keys(results))
  its = zeros(length(problem_list))

  for i = 1:length(problem_list)
      problem_name = problem_list[i]
      info = results[problem_name]

      if alg_success(info.status) && info.it_count < MAX_IT
        if info.it_count < 1
          println("warning: non-positive iteration count for $problem_name")
          its[i] = 1
        else
          its[i] = info.it_count
        end
      else
        its[i] = Infty
      end
  end

  return its
end

function best_its(it_list_by_method::Dict)
    n = length(first(values(it_list_by_method)))
    best_its = ones(n) * Infty
    for (method_name, it_list) in it_list_by_method
        @assert( length(it_list) == n )
        for i = 1:n
          best_its[i] = min(best_its[i], it_list[i])
        end
    end

    return best_its
end

function overlap(dics::Dict{String, Dict{String,problem_summary2}})
  overlap = collect(keys(first(dics)[2]))
  for (method_name, data) in dics
    problem_list = collect(keys(data))
    overlap = intersect(problem_list, overlap)
  end

  return restrict_problems(dics, overlap)

  #@show overlap
  #=new_dic = Dict{String, Dict{String,problem_summary2}}()
  for (method_name, data) in dics
    new_dic[method_name] = Dict{String,problem_summary2}()
    for problem_name in overlap
      new_dic[method_name][problem_name] = data[problem_name]
    end
  end=#
end

function all_status(res::Dict{String, Dict{String,problem_summary2}},status::Symbol)
  method_list = keys(res)

  opt_res = Dict{String, Dict{String,problem_summary2}}()
  for method in method_list
      opt_res[method] = Dict{String,problem_summary2}()
  end

  problem_list = keys(first(res)[2])
  for problem_name in problem_list
      include_prob = true
      for method in method_list
        if res[method][problem_name].status != status
          include_prob = false
        end
      end

      if include_prob
        for method in method_list
           opt_res[method][problem_name] = res[method][problem_name]
        end
      end
  end

  return opt_res
end

function restrict_to_set(dics::Dict{String, Dict{String,problem_summary2}}, set)
    p_list = []
    for (method_name,data) in dics
      p_list = union(get_list_of(data, set), p_list)
    end
    return restrict_problems(dics, p_list)
end

function get_list_of(results::Dict{String, problem_summary2}, set)
  new_problem_list = []
  for (problem_name, info) in results
      if info.status in set
        push!(new_problem_list,problem_name)
      end
  end

  return new_problem_list
end

function restrict_problems(results, problem_list)
  new_results = Dict{String,Dict{String,problem_summary2}}()
  for (method_name, sum_data) in results
      new_results[method_name] = Dict()
      for (problem_name, info) in sum_data
        if problem_name in problem_list
          new_results[method_name][problem_name] = info
        end
      end
  end
  return new_results
end

function tot(results::Dict{String, problem_summary2}, set)
  problem_list = collect(keys(results))
  tot = 0

  for i = 1:length(problem_list)
      problem_name = problem_list[i]
      info = results[problem_name]
      if info.status in set
        tot += 1
      end
  end

  return tot
end

function tot_failures(results::Dict{String, problem_summary2};MAX_IT::Int64=3000)
  problem_list = collect(keys(results))
  tot = 0

  for i = 1:length(problem_list)
      problem_name = problem_list[i]
      info = results[problem_name]
      if alg_success(info.status) && info.it_count <= MAX_IT
        # optimal
      else
        tot += 1
      end
  end

  return tot
end

function failure_list(results::Dict{String, problem_summary2};MAX_IT=3000)
  problem_list = collect(keys(results))
  failure_list = zeros(length(problem_list))

  for i = 1:length(problem_list)
      problem_name = problem_list[i]
      info = results[problem_name]
      if alg_success(info.status) && info.it_count <= MAX_IT
        failure_list[i] = 0
      else
        failure_list[i] = 1
      end
  end

  return failure_list
end

function shared_failures(results::Dict{String, Dict{String,problem_summary2}})
  problem_list = collect(keys(first(results)[2]))
  all_fail = ones(length(problem_list))
  #@show length(all_fail)
  for (method_name, data) in results
      ls = failure_list(data)
      #@show length(ls)
      all_fail = all_fail .* ls
  end
  return sum(all_fail)
end

function failure_causes(results::Dict{String, problem_summary2})
  problem_list = collect(keys(results))
  failure_causes = Dict()

  for i = 1:length(problem_list)
      problem_name = problem_list[i]
      info = results[problem_name]
      if !alg_success(info.status) #&& info.it_count <= 3000
        if info.status in keys(failure_causes)
          failure_causes[info.status] += 1
        else
          failure_causes[info.status] = 1
        end
      end
  end

  return failure_causes
end

function list_combined_failures(results::Dict{String, Dict{String,problem_summary2}})
  problem_list = collect(keys(first(results)[2]))
  all_fail = ones(length(problem_list))
  #@show length(all_fail)
  for (method_name, data) in results
      for (problem_name,info) in data
        ls = failure_list(data)
        #@show length(ls)
        all_fail = all_fail .* ls
      end
  end
  println("shared failures = ", sum(all_fail))
end

function list_combined_success(results::Dict{String, Dict{String,problem_summary2}})
  problem_list = collect(keys(first(results)[2]))
  all_succeed = ones(length(problem_list))
  #@show length(all_fail)
  for (method_name, data) in results
      for (problem_name,info) in data
        ls = 1.0 - failure_list(data)
        #@show length(ls)
        all_succeed = all_succeed .* ls
      end
  end
  println("shared successes = ", sum(all_succeed))
end

#=
function outcomes_table(results::Dict{String, Dict{String,problem_summary2}})
    outcomes = Dict()
    for problem_name in problem_list
        key = ()
        for (method_name, data) in results
            key = info.status
            #outcomes[(info.status,] += 1
        end
    end
end
=#

function print_failure_problems(results::Dict{String, Dict{String,problem_summary2}};MAX_IT::Int64=3000)
  problem_list = collect(keys(first(results)[2]))

  println("FAILURES")
  for i = 1:length(problem_list)
      problem_name = problem_list[i]
      first_failure = true

      for (method_name, data) in results
        info = results[method_name][problem_name]
        if alg_success(info.status) && info.it_count <= MAX_IT
          # optimal
        else
          if first_failure
            print(problem_name)
            first_failure = false
          end
          print(" ", method_name)
        end
      end
      if first_failure == false
        print("\n")
      end
  end
end

function remove_errors(dics::Dict{String, Dict{String,problem_summary2}}, error_set)
    new_dic = Dict{String, Dict{String,problem_summary2}}()
    for (method_name, data) in dics
      new_dic[method_name] = Dict{String,problem_summary2}()
      for (problem_name, info) in data
        if !(info.status in error_set)
          new_dic[method_name][problem_name] = info
        end
      end
    end
    return new_dic
end

function select(dics::Dict{String, Dict{String,problem_summary2}}, select_set)
    new_dic = Dict{String, Dict{String,problem_summary2}}()
    for (method_name, data) in dics
      new_dic[method_name] = Dict{String,problem_summary2}()
      for (problem_name, info) in data
        if (info.status in select_set)
          new_dic[method_name][problem_name] = info
        end
      end
    end
    return new_dic
end

function union_results(dics::Dict{String, Dict{String,problem_summary2}}, master::Dict{String, Dict{String,problem_summary2}})
    names_list = []
    for (method_name, data) in dics
      names_list = union(collect(keys(data)),names_list)
    end

    new_dic = Dict{String, Dict{String,problem_summary2}}()
    for (method_name, data) in results
      new_dic[method_name] = Dict{String,problem_summary2}()
      for (problem_name, info) in data
        if problem_name in names_list
          new_dic[method_name][problem_name] = info
        end
      end
    end
    return new_dic
end

function show_results(results::Dict{String, Dict{String,problem_summary2}})
    problem_list = collect(keys(first(results)[2]))
    for problem_name in problem_list
      for (method_name, data) in results
         info = data[problem_name]
         println(method_name, " ", problem_name, " ", info.it_count, " ", info.status)
      end
    end
end

function compute_its_etc(overlapping_results; MAX_IT::Int64=3000)
    its = Dict{String,Array{Int64,1}}()
    println("Quartiles of iterations:")
    for (method_name, sum_data) in overlapping_results
        its[method_name] = iteration_list(sum_data, MAX_IT=MAX_IT);
        d = its[method_name]
        println(OnePhase.pd(method_name,20), " = ", OnePhase.rd(quantile(d,0.25)), OnePhase.rd(quantile(d,0.5)), OnePhase.rd(quantile(d,0.75)))
    end

    println("Times:")
    times = Dict{String,Array{Float64,1}}()
    for (method_name, sum_data) in overlapping_results
        times[method_name] = []
        for (problem_name,info) in sum_data
          if alg_success(info.status)
              t = info.total_time
          else
              t = Infty
          end

          push!(times[method_name], t);
        end
        tm = times[method_name]
        println(OnePhase.pd(method_name,20), " = ", OnePhase.rd(quantile(tm,0.25)), OnePhase.rd(quantile(tm,0.5)), OnePhase.rd(quantile(tm,0.75)))
        #println(method_name, " mean=", rd(mean(times[method_name])), " median=", rd(median(times[method_name])))
    end

    println("Quartiles of iteration ratios to best solver:")
    best = best_its(its)

    ratios = Dict()
    for (method_name, val) in its
      ratios[method_name] = its[method_name] ./ best;
      lrg = its[method_name] .>= Infty
      ratios[method_name][lrg] = Inf

      d = ratios[method_name]
      println(OnePhase.pd(method_name,20), " = ", OnePhase.rd(quantile(d,0.25)), OnePhase.rd(quantile(d,0.5)), OnePhase.rd(quantile(d,0.75)))
    end
    return its, best, ratios, times
end
####################################################################################################
####################################################################################################
#############################       DUAL VARIABLE STUFF                #############################
####################################################################################################
####################################################################################################

function compute_dual_hist(results)
    dual_hist = Array{Float64,1}()
    t = 0
    for it_hist in results
      if it_hist.t > t
        push!( dual_hist, it_hist.y_norm )
        t += 1
      end
    end
    #dual_maxes[folder_name][problem_name] = maximum(dual_hist)

    return dual_hist
end

function get_lists(summary)
    problem_list = collect(keys(first(summary)[2]))
    method_list = keys(summary)

    return problem_list, method_list
end

function get_all_optimal(summary, problem_list) # get all problems where all solvers get the optimal solution
    remove_these = Array{String,1}()

    for problem_name in problem_list
      for (method_name, sum_data) in summary_netlib
          if sum_data[problem_name].status != :Optimal
            push!(remove_these, problem_name)
          end
      end
    end

    all_optimal = Array{String,1}()

    for problem_name in problem_list
        if !(problem_name in remove_these)
          push!(all_optimal, problem_name)
        end
    end

    return all_optimal
end

function summarize_by_iteration(hist)
    sum_hist = Array{abstract_alg_history,1}()
    for i = 2:length(hist)
      if hist[i].t > hist[i-1].t
        push!(sum_hist, hist[i-1])
      end
    end
    push!(sum_hist, hist[end])

    return sum_hist
end

function aggregate_summary_by_it_for_plots(sh, field_list::Dict{String,Symbol})
    data_us = Dict()

    for (label, field) in field_list
        data_us[field] = zeros(length(sh));
    end

    for i = 1:length(sh)
        for (label, field) in field_list
          data_us[field][i] = getfield(sh[i],field);
        end
    end

    return data_us
end


text:
include("../utils.jl")
include("group_utils.jl")

@testset "Rotation action" begin
    M = Rotations(2)
    G = SpecialOrthogonal(2)
    A_left = RotationAction(Euclidean(2), G)
    A_right = RotationAction(Euclidean(2), G, RightAction())

    @test repr(A_left) == "RotationAction($(repr(Euclidean(2))), $(repr(G)), LeftAction())"
    @test repr(A_right) ==
          "RotationAction($(repr(Euclidean(2))), $(repr(G)), RightAction())"

    types_a = [Matrix{Float64}]

    types_m = [Vector{Float64}]

    @test g_manifold(A_left) == Euclidean(2)
    @test base_group(A_left) == G
    @test isa(A_left, AbstractGroupAction{LeftAction})
    @test base_manifold(G) == M

    for (i, T_A, T_M) in zip(1:length(types_a), types_a, types_m)
        angles = (0.0, π / 2, 2π / 3, π / 4)
        a_pts = convert.(T_A, [[cos(ϕ) -sin(ϕ); sin(ϕ) cos(ϕ)] for ϕ in angles])
        m_pts = convert.(T_M, [[0.0, 1.0], [-1.0, 0.0], [1.0, 1.0]])
        X_pts = convert.(T_M, [[1.0, 2.0], [3.0, 4.0], [5.0, 6.0]])

        atol = if eltype(T_M) == Float32
            2e-7
        else
            1e-15
        end
        test_action(
            A_left,
            a_pts,
            m_pts,
            X_pts;
            test_optimal_alignment=true,
            test_diff=true,
            atol=atol,
        )

        test_action(
            A_right,
            a_pts,
            m_pts,
            X_pts;
            test_optimal_alignment=true,
            test_diff=true,
            atol=atol,
        )
    end
end

@testset "Rotation around axis action" begin
    M = Circle()
    G = RealCircleGroup()
    axis = [sqrt(2) / 2, sqrt(2) / 2, 0.0]
    A = Manifolds.RotationAroundAxisAction(axis)

    types_a = [Ref(Float64)]

    types_m = [Vector{Float64}]

    @test g_manifold(A) == Euclidean(3)
    @test base_group(A) == G
    @test isa(A, AbstractGroupAction{LeftAction})
    @test base_manifold(G) == M

    for (i, T_A, T_M) in zip(1:length(types_a), types_a, types_m)
        angles = (0.0, π / 2, 2π / 3, π / 4)
        a_pts = convert.(T_A, [angles...])
        m_pts = convert.(T_M, [[0.0, 1.0, 0.0], [-1.0, 0.0, 1.0], [1.0, 1.0, -2.0]])
        X_pts = convert.(T_M, [[1.0, 2.0], [3.0, 4.0], [5.0, 6.0]])

        atol = if eltype(T_M) == Float32
            2e-7
        else
            1e-15
        end
        test_action(
            A,
            a_pts,
            m_pts,
            X_pts;
            test_optimal_alignment=false,
            test_diff=false,
            test_mutating_group=false,
            test_switch_direction=false,
            atol=atol,
        )
    end
end


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

module BaseDocs

immutable Keyword
    name :: Symbol
end
macro kw_str(text) Keyword(Symbol(text)) end

"Hello, Human."
kw"hello", kw"hi"

"""
**Welcome to Julia $(string(VERSION)).** The full manual is available at

    http://docs.julialang.org/

as well many great tutorials and learning resources:

    http://julialang.org/learning/

For help on a specific function or macro, type `?` followed
by its name, e.g. `?fft`, or `?@time`, and press enter.
"""
kw"help", kw"?", kw"julia"

"""
`using` will load the given module or package and make some of its names available for
use (see also `export`). For example:

    using Gadfly

loads the plotting package, Gadfly, so that the `plot` function can be used.

Names can be used via dot syntax, whether they are exported or not:

    Gadfly.plot(...)

If you don't want to use the packages exports directly, see also `import`. If you're not
sure, `using` is almost definitely what you want.
"""
kw"using"

"""
    import Gadfly

`import`, like `using`, will load modules and packages for use. Unlike `using`, however,
it will *not* make any `export`ed names available for use. To use Gadfly's `plot`
function after importing it, for example, you have to write:

    Gadfly.plot(...)

Import can also be used with specific names, for example

    import Gadfly: plot, render

This syntax is used when you want to extend the modules functions with new methods.
"""
kw"import"

"""
`export` is used within modules and packages to tell Julia which functions should be
made available to the user. For example:

    module Test
    export foo # foo is exported, but bar isn't
    foo(x) = x
    bar(y) = y
    end

    using Test
    foo(1) # 1
    bar(1) # Error: bar not defined
    Test.bar(1) # 1
"""
kw"export"

"""
`abstract` declares a type that cannot be instantiated, and serves only as a node in the
type graph, thereby describing sets of related concrete types: those concrete types
which are their descendants. Abstract types form the conceptual hierarchy which makes
Julia’s type system more than just a collection of object implementations. For example:

    abstract Number
    abstract Real <: Number

`abstract Number` has no supertype, whereas `abstract Real` is an abstract subtype of `Number`.
"""
kw"abstract"

"""
`module` declares a Module, which is a separate global variable workspace.  Within a
module, you can control which names from other modules are visible (via importing), and
specify which of your names are intended to be public (via exporting). For example:

    module
    import Base.show
    export MyType, foo

    type MyType
        x
    end

    bar(x) = 2x
    foo(a::MyType) = bar(a.x) + 1
    show(io, a::MyType) = print(io, "MyType \$(a.x)")
    end

Modules allow you to create top-level definitions without worrying about name conflicts
when your code is used together with somebody else’s.
"""
kw"module"

"""
`baremodule` declares a module that does not contain `using Base`
or a definition of `eval`.  It does still import `Core`.
"""
kw"baremodule"

"""
`bitstype` declares a concrete type whose data consists of plain old bits. Classic
examples of bits types are integers and floating-point values. Some example built-in
bits type declarations:

    bitstype 32 Char
    bitstype 8  Bool <: Integer

The first parameter indicates how many bits of storage the type requires. Currently,
only sizes that are multiples of 8 bits are supported. The second parameter gives the
name of the type.  The `Bool` declaration shows how a bits type can be optionally
declared to be a subtype of some supertype.
"""
kw"bitstype"

"""
`macro` defines a method to include generated code in the final body of a program. A
macro maps a tuple of arguments to a returned expression, and the resulting expression
is compiled directly rather than requiring a runtime `eval()` call. Macro arguments may
include expressions, literal values, and symbols. For example:

    macro sayhello(name)
        return :( println("Hello, ", \$name) )
    end

This macro takes one argument: `name`. When `@sayhello` is encountered, the quoted
expression is expanded to interpolate the value of the argument into the final
expression.
"""
kw"macro"

"""
`importall` imports all names exported by the specified module, as if `import` were used
individually on all of them.  For example:

    importall Distributions

As with `import`, functions imported by `importall` can be extended.
"""
kw"importall"

"""
`local` introduces a new local variable. For example:

    function foo(n)
        x = 0
        for i = 1:n
            local x
            x = i
        end
        x
    end

    julia> foo(10)
    0

Here `local x` introduces a separate `x` inside the loop, so the function returns `0`.
"""
kw"local"

"""
`global x` makes `x` in the current scope and its inner scopes refer to the global
variable of that name.   In the example below, `global` is needed so the function can
modify the global variable `z`:

    z=3
    function foo()
        global z=6
    end

    julia> foo()
    6
    julia> z
    6

Without the `global` declaration in `foo()`, a new local variable would have been
created inside foo(), and the `z` in the global scope would have remained equal to `3`.
"""
kw"global"

"""
`let` statements allocate new variable bindings each time they run. Whereas an
assignment modifies an existing value location, `let` creates new locations. This
difference is only detectable in the case of variables that outlive their scope via
closures. The `let` syntax accepts a comma-separated series of assignments and variable
names:

    let var1 = value1, var2, var3 = value3
        code
    end

The assignments are evaluated in order, with each right-hand side evaluated in the scope
before the new variable on the left-hand side has been introduced. Therefore it makes
sense to write something like `let x = x`, since the two `x` variables are distinct and
have separate storage.
"""
kw"let"

"""
`quote` creates multiple expression objects in a block without using the explicit `Expr`
constructor. For example:

    ex = quote
        x = 1
        y = 2
        x + y
    end

Unlike the other means of quoting, `:( ... )`, this form introduces `QuoteNode` elements
to the expression tree, which must be considered when directly manipulating the tree.
For other purposes, `:( ... )` and `quote .. end` blocks are treated identically.
"""
kw"quote"

"""
`'` is the conjugate transposition operator:

    julia> A = reshape(1:4, 2,2)
    2×2 Array{Int64,2}:
     1  3
     2  4

    julia> A'
    2×2 Array{Int64,2}:
     1  2
     3  4

    julia> B = A + im
    2×2 Array{Complex{Int64},2}:
     1+1im  3+1im
     2+1im  4+1im

    julia> B'
    2×2 Array{Complex{Int64},2}:
     1-1im  2-1im
     3-1im  4-1im

"""
kw"'"

"""
`.'` is the transposition operator:

    julia> A = reshape(1:4, 2,2)
    2×2 Array{Int64,2}:
     1  3
     2  4

    julia> A.'
    2×2 Array{Int64,2}:
     1  2
     3  4

    julia> B = A + im
    2×2 Array{Complex{Int64},2}:
     1+1im  3+1im
     2+1im  4+1im

    julia> B.'
    2×2 Array{Complex{Int64},2}:
     1+1im  2+1im
     3+1im  4+1im

    julia> v = [1,2,3]
    3-element Array{Int64,1}:
     1
     2
     3

    julia> v.'
    1×3 RowVector{Int64,Array{Int64,1}}:
     1  2  3

"""
kw".'"

"""
`const` is used to declare global variables which are also constant. In almost all code
(and particularly performance sensitive code) global variables should be declared
constant in this way.

    const x = 5

Note that "constant-ness" is not enforced inside containers, so if `x` is an array or
dictionary (for example) you can still add and remove elements.

Technically, you can even redefine `const` variables, although this will generate a
warning from the compiler. The only strict requirement is that the *type* of the
variable does not change, which is why `const` variables are much faster than regular
globals.
"""
kw"const"

"""
Functions are defined with the `function` keyword:

    function add(a, b)
        return a + b
    end

Or the short form notation:

    add(a, b) = a + b

The use of the `return` keyword is exactly the same as in other languages, but is often
optional. When it's not used, the last expression in the function body will be returned
by default:

    function compare(a, b)
        a == b && return "equal to"
        a < b ? "less than" : "greater than"
    end
"""
kw"function"

"""
`return` can be used function bodies to exit early and return a given value, e.g.

    function compare(a, b)
        a == b && return "equal to"
        a < b ? "less than" : "greater than"
    end

In general you can place a `return` statement anywhere within a function body, including
within deeply nested loops or conditionals, but be careful with `do` blocks. For
example:

    function test1(xs)
        for x in xs
            iseven(x) && return 2x
        end
    end

    function test2(xs)
        map(xs) do x
            iseven(x) && return 2x
            x
        end
    end

In the first example, the return breaks out of its enclosing function as soon as it hits
an even number, so `test1([5,6,7])` returns `12`.

You might expect the second example to behave the same way, but in fact the `return`
there only breaks out of the *inner* function (inside the `do` block) and gives a value
back to `map`. `test2([5,6,7])` then returns `[5,12,7]`.
"""
kw"return"

"""
`if`-`elseif`-`else` performs conditional evaluation, which allows portions of code to
be evaluated or not evaluated depending on the value of a boolean expression. Here is
the anatomy of the `if`-`elseif`-`else` conditional syntax:

    if x < y
        println("x is less than y")
    elseif x > y
        println("x is greater than y")
    else
        println("x is equal to y")
    end

If the condition expression `x < y` is true, then the corresponding block is evaluated;
otherwise the condition expression `x > y` is evaluated, and if it is true, the
corresponding block is evaluated; if neither expression is true, the `else` block is
evaluated. The `elseif` and `else` blocks are optional, and as many `elseif` blocks as
desired can be used.
"""
kw"if", kw"elseif", kw"else"

"""
`for` loops repeatedly evaluate the body of the loop by iterating over a sequence of
values. For example:

    for i in [1,4,0]
        println(i)
    end
"""
kw"for"

"""
`while` loops repeatedly evaluate a conditional expression, and continues evaluating the
body of the while loop so long as the expression remains `true`. If the condition
expression is false when the while loop is first reached, the body is never evaluated.
For example:

    while i <= 5
        println(i)
        i += 1
    end
"""
kw"while"

"""
`end` marks the conclusion of a block of expressions. In the example below, `end` marks
the conclusion of a `function`.

    function foo()
        println("hello, world")
    end

`end` marks the conclusion of all kinds of expression blocks: `module`, `type`, `begin`,
`let`, `for`, etc.

In addition, `end` may be used when indexing into an array to represent the last index
of each dimension:

    x[1:end, 2:end-1]
"""
kw"end"

"""
A `try/catch` statement allows for `Exception`s to be tested for. For example, a
customized square root function can be written to automatically call either the real or
complex square root method on demand using `Exception`s:

    f(x) = try
        sqrt(x)
    catch
        sqrt(complex(x, 0))
    end

`try/catch` statements also allow the `Exception` to be saved in a variable, e.g. `catch y`.

The `catch` clause is not strictly necessary; when omitted, the default return value is
`nothing`. The power of the `try/catch` construct lies in the ability to unwind a deeply
nested computation immediately to a much higher level in the stack of calling functions.
"""
kw"try", kw"catch"

"""
`finally` provides a way to run some code when a given block of code exits, regardless
of how it exits. For example, here is how we can guarantee that an opened file is
closed:

    f = open("file")
    try
        operate_on_file(f)
    finally
        close(f)
    end

When control leaves the `try` block (for example due to a `return`, or just finishing
normally), `close(f)` will be executed. If the `try` block exits due to an exception,
the exception will continue propagating. A `catch` block may be combined with `try` and
`finally` as well. In this case the `finally` block will run after `catch` has handled
the error.
"""
kw"finally"

"""
`break` breaks out of a loop immediately. For example

    i = 0
    while true
        i += 1
        i > 10 && break
        println(i)
    end

prints the numbers 1 to 10.
"""
kw"break"

"""
`continue` skips the rest of the current loop, then carries on looping. For example

    for i = 1:10
        iseven(i) && continue
        println(i)
    end

prints the numbers 1, 3, 5..., skipping the even numbers.
"""
kw"continue"

"""
The `do` keyword creates an anonymous function. For example

    map(1:10) do x
        2x
    end

is equivalent to `map(x->2x, 1:10)`.

Use multiple arguments like so:

    map(1:10, 11:20) do x, y
        x + y
    end
"""
kw"do"

"""
The "splat" operator, `...`, represents a sequence of arguments. For example

    add(xs...) = reduce(+, xs)

can take any number of arguments:

    add(1, 2, 3, 4, 5)

`...` can also be used to apply a function to a sequence of arguments like so:

    add([1, 2, 3]...) # 6
    add(7, 1:100..., 1000:1100...) # 111107
"""
kw"..."

"""
`;` has a similar role in Julia as in many C-like languages, and is used to delimit the
end of the previous statement. `;` is not necessary after new lines, but can be used to
separate statements on a single line or to join statements into a single expression:

    function foo()
        println("Hello, "); println("World!")
        return true
    end

    foo() = (println("Hello, World!"); true)

`;` is also used to suppress output in the REPL and similar interfaces.
"""
kw";"

"""
    x && y

Short-circuiting boolean AND.
"""
kw"&&"

"""
    x || y

Short-circuiting boolean OR.
"""
kw"||"

"""
    ccall((symbol, library) or function_pointer, ReturnType, (ArgumentType1, ...), ArgumentValue1, ...)

Call function in C-exported shared library, specified by `(function name, library)`
tuple, where each component is a string or symbol.

Note that the argument type tuple must be a literal tuple, and not a tuple-valued
variable or expression. Alternatively, `ccall` may also be used to call a function
pointer, such as one returned by `dlsym`.

Each `ArgumentValue` to the `ccall` will be converted to the corresponding
`ArgumentType`, by automatic insertion of calls to `unsafe_convert(ArgumentType,
cconvert(ArgumentType, ArgumentValue))`. (See also the documentation for each of these
functions for further details.) In most cases, this simply results in a call to
`convert(ArgumentType, ArgumentValue)`.
"""
kw"ccall"

"""
    llvmcall(IR::String, ReturnType, (ArgumentType1, ...), ArgumentValue1, ...)
    llvmcall((declarations::String, IR::String), ReturnType, (ArgumentType1, ...), ArgumentValue1, ...)

Call LLVM IR string in the first argument. Similar to an LLVM function `define` block,
arguments are available as consecutive unnamed SSA variables (%0, %1, etc.).

The optional declarations string contains external functions declarations that are
necessary for llvm to compile the IR string. Multiple declarations can be passed in by
separating them with line breaks.

Note that the argument type tuple must be a literal tuple, and not a tuple-valued
variable or expression.

Each `ArgumentValue` to `llvmcall` will be converted to the corresponding
`ArgumentType`, by automatic insertion of calls to `unsafe_convert(ArgumentType,
cconvert(ArgumentType, ArgumentValue))`. (see also the documentation for each of these
functions for further details). In most cases, this simply results in a call to
`convert(ArgumentType, ArgumentValue)`.

See `test/llvmcall.jl` for usage examples.
"""
Core.Intrinsics.llvmcall

"""
`begin...end` denotes a block of code.

    begin
        println("Hello, ")
        println("World!")
    end

Usually `begin` will not be necessary, since keywords such as `function` and `let`
implicitly begin blocks of code. See also `;`.
"""
kw"begin"

"""
At their most basic, Julia types are specified as a name and a set of fields.

    type Point
        x
        y
    end

Fields can have type restrictions, which may be parametrised:

    type Point{X}
        x::X
        y::Float64
    end

Type can also declare an abstract super type via `<:` syntax:

    type Point <: AbstractPoint
        ...

See the manual for more details, such as information on inner constructors.
"""
kw"type"

"""
Introduce a new name for an already expressible type. For example, in `base/boot.jl`,
`UInt` is type aliased to either `UInt64` or `UInt32` as appropriate for the size of
pointers on the system:

    if Int === Int64
        typealias UInt UInt64
    else
        typealias UInt UInt32
    end

For parametric types, `typealias` can be convenient for providing names in cases where
some parameter choices are fixed.  In `base` for example:

    typealias Vector{T} Array{T,1}
"""
kw"typealias"

"""
`immutable` acts in the same way as `type`, but declares that the fields of the type may
not be set after construction. See `type` and the manual for more information.
"""
kw"immutable"

"""
    @__LINE__ -> Int

`@__LINE__` expands to the line number of the call-site.
"""
kw"@__LINE__"

"""
    ans

A variable referring to the last computed value, automatically set at the interactive prompt.
"""
kw"ans"

"""
    nothing

The singleton instance of type `Void`, used by convention when there is no value to return
(as in a C `void` function). Can be converted to an empty `Nullable` value.
"""
nothing

"""
    ANY

Equivalent to `Any` for dispatch purposes, but signals the compiler to skip code
generation specialization for that field.
"""
ANY

"""
    DevNull

Used in a stream redirect to discard all data written to it. Essentially equivalent to
/dev/null on Unix or NUL on Windows. Usage:

```julia
run(pipeline(`cat test.txt`, DevNull))
```
"""
DevNull

end


text:
struct priorInfo{T1, T2, T3 <:Function}
    f::T1
    df::T2
    d2f::T3
end

function initialize_prior(range, type="Uniform")
    if type=="Uniform"
        return priorInfo(x -> 1/(range[2]-range[1]), x-> 0, x->0)
    else
        throw(ArgumentError("Prior type not supported."))
    end
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
__precompile__()

using CPLEX
include("../src/FeasibilityPump.jl")
using .FeasibilityPump
__compilefp__()
using Random

function mainBasic()
    files = "../data/" .* readdir("../data/")

    nSeed = 3

    for file in files
        fileName = split(file, ['/'])[end][1:end-4]

        println("--- "*file*" ---")
        println("- Solving with feasibility pump:")
        for i in 1:nSeed
            GC.gc()
            Random.seed!(i)
            myModel = CPLEX.Model(CPLEX.Env())
            CPLEX.read_model(myModel, file)
            CPLEX.set_param!(myModel.env, "CPX_PARAM_THREADS", 1)

            param = initializeparametersfp()
            #setparamfp!(param, "logLevel", 1)

            timeTmp = @timed begin
                st, x, z = feasibilitypump(myModel, param)
            end
            print("-Time: ")
            println(timeTmp[2])
            println("-Status: ", st)
            println("-ObjectiveValue: ", z)
            println()
            println()
        end
    end
end

mainBasic()

max_stars_repo_name: remi-garcia/FeasibilityPump

max_issues_repo_name: remi-garcia/FeasibilityPump

max_forks_repo_name: remi-garcia/FeasibilityPump


text:
### Diagonal ops

# we use Base.Diagonal(f) for diagonal operators so very little specific code is
# actually needed here.

simulate(D::DiagOp{F}) where {F<:Field} = sqrt(D) * white_noise(F)

# automatic basis conversion (and NaN-zeroing)
(*)(D::DiagOp{<:Field{B}}, f::Field) where {B} = D.diag .* B(f)
(\)(D::DiagOp{<:Field{B}}, f::Field) where {B} = nan2zero.(D.diag .\ B(f))

# broadcasting
BroadcastStyle(::StructuredMatrixStyle{<:DiagOp{F}}, ::StructuredMatrixStyle{<:DiagOp{<:ImplicitField}}) where {F<:Field} = StructuredMatrixStyle{DiagOp{F}}()
BroadcastStyle(::StructuredMatrixStyle{<:DiagOp{<:ImplicitField}}, ::StructuredMatrixStyle{<:DiagOp{F}}) where {F<:Field} = Base.Broadcast.Unknown()
function similar(bc::Broadcasted{<:StructuredMatrixStyle{<:DiagOp{F}}}, ::Type{T}) where {F<:Field,T}
    Diagonal(similar(typeof(BroadcastStyle(F)),T))
end
diag_data(D::Diagonal) = D.diag
diag_data(x) = x
function copyto!(dest::DiagOp, bc::Broadcasted{<:StructuredMatrixStyle})
    copyto!(dest.diag, map_bc_args(diag_data, bc))
    dest
end

# the generic versions of these kind of suck so we need these specializized
# versions:
(*)(x::Adjoint{<:Any,<:Field}, D::Diagonal) = (D*parent(x))'
(*)(x::Adjoint{<:Any,<:Field}, D::Diagonal, y::Field) = x*(D*y)
diag(L::DiagOp) = L.diag

# use getproperty here to ensure no basis conversion is done
getindex(D::DiagOp, s::Symbol) = Diagonal(getproperty(D.diag,s))

# for printing
size(::DiagOp{<:ImplicitField}) = ()
axes(::DiagOp{<:ImplicitField}, I) = OneTo(0)

# the generic version of this is prohibitively slow so we need this
hash(D::DiagOp, h::UInt64) = hash(D.diag, h)

### Derivative ops

# These ops just store the coordinate with respect to which a derivative is
# being taken, and each Field type F implements how to actually take the
# derivative.
#
# ∇diag{coord, covariance, prefactor} is a Field which represents the diagonal
# entries of the Fourier representation of the gradient operation, with `coord`
# specifiying the coordinate for the gradient direction (1 or 2), `covariance`
# being :covariant or :contravariant, and multiplied by a `prefactor` (generally
# 1 or -1 to keep track of if an adjoint has been taken). ∇[i] returns one of
# these wrapped in a Diagonal so that multipliying by these does the necessary
# basis conversion. ∇ (along with ∇ⁱ, and ∇ᵢ) are StaticVectors which can also
# be used with FieldVector algebra, e.g. ∇*f.
#
# Note: we are cheating a bit, because while ∇ is diagonal when a derivative is
# taken in Fourier space via FFT's, its tri-diagonal if taking a map-space
# derivative, so conceptually these should not be wrapped in a Diagonal. Its a
# bit ugly and I'd like to fix it, but since the auto-basis conversion mechanism
# is based on Diagonal, keep things this way for now. Individual field types
# just need to intercept the broadcast machinery at copyto! to implement the
# map-space derivative, e.g. see: flat_generic.jl.
#
# Also note: We define the components of vectors, including ∇, to be with respect to
# the _unnormalized_ covariant or contravariant basis vectors, hence ∇ⁱ = d/dxᵢ
# and ∇ᵢ = d/dxⁱ. This is different than with respect to the _normalized)
# covariant basis vectors, which, e.g., in spherical coordinates, gives the more
# familiar ∇ = (d/dθ, 1/sinθ d/dϕ), (but whose components are neither covariant
# nor contravariant).

struct ∇diag{coord, covariance, prefactor} <: ImplicitField{DerivBasis,Spin,Pix} end
struct ∇²diag <: ImplicitField{DerivBasis,Spin,Pix} end

# adjoint(D::Diagonal{<:∇diag}) calls conj(D.diag) and here lazily we keep track
# of that a conjugate was taken
conj(::∇diag{coord,covariance,prefactor}) where {coord,covariance,prefactor} = ∇diag{coord,covariance,-prefactor}()
-(::DiagOp{∇diag{coord,covariance,prefactor}}) where {coord,covariance,prefactor} = Diagonal(∇diag{coord,covariance,-prefactor}())

# Gradient vector which can be used with FieldVector algebra.
struct ∇Op{covariance,prefactor} <: StaticVector{2,Diagonal{Float32,∇diag{<:Any,covariance,prefactor}}} end
getindex(::∇Op{covariance,prefactor}, i::Int) where {covariance,prefactor} = Diagonal(∇diag{i,covariance,prefactor}())
const ∇ⁱ = ∇Op{:contravariant,1}()
const ∇ᵢ = ∇Op{:covariant,1}()
const ∇ = ∇ⁱ # ∇ is contravariant by default if not specified
const ∇² = Diagonal(∇²diag())

@doc doc"""
    gradhess(f)

Compute the gradient $g^i = \nabla^i f$, and the hessian, $H_j^{\,i} = \nabla_j \nabla^i f$.
"""
function gradhess(f)
    g = ∇ⁱ*f
    (g=g, H=SMatrix{2,2}([permutedims(∇ᵢ * g[1]); permutedims(∇ᵢ * g[2])]))
end

# this is not strictly true (∇[1] is generically a gradient w.r.t. the first
# coordinate, e.g. ∂θ), but this is useful shorthand to have for the flat-sky:
const ∂x = ∇[1]
const ∂y = ∇[2]

### FuncOp

# An Op which applies some arbitrary function to its argument.
# Transpose and/or inverse operations which are not specified will return an error.
@kwdef struct FuncOp <: ImplicitOp{Basis,Spin,Pix}
    op   = nothing
    opᴴ  = nothing
    op⁻¹ = nothing
    op⁻ᴴ = nothing
end
SymmetricFuncOp(;op=nothing, op⁻¹=nothing) = FuncOp(op,op,op⁻¹,op⁻¹)
*(L::FuncOp, f::Field) = L.op   != nothing ? L.op(f)   : error("op*f not implemented")
*(f::Field, L::FuncOp) = L.opᴴ  != nothing ? L.opᴴ(f)  : error("f*op not implemented")
\(L::FuncOp, f::Field) = L.op⁻¹ != nothing ? L.op⁻¹(f) : error("op\\f not implemented")
adjoint(L::FuncOp) = FuncOp(L.opᴴ,L.op,L.op⁻ᴴ,L.op⁻¹)
const IdentityOp = FuncOp(identity,identity,identity,identity)
inv(L::FuncOp) = FuncOp(L.op⁻¹,L.op⁻ᴴ,L.op,L.opᴴ)
adapt_structure(to, L::FuncOp) = FuncOp(adapt(to, fieldvalues(L))...)

### BandPassOp

# An op which applies some bandpass, like a high or low-pass filter. This object
# stores the bandpass weights, Wℓ, and each Field type F should implement
# broadcast_data(::Type{F}, ::BandPassOp) to describe how this is actually
# applied.

abstract type HarmonicBasis <: Basislike end

struct BandPass{W<:InterpolatedCℓs} <: ImplicitField{HarmonicBasis,Spin,Pix}
    Wℓ::W
end
BandPassOp(ℓ,Wℓ) = Diagonal(BandPass(InterpolatedCℓs(promote(collect(ℓ),collect(Wℓ))...)))
BandPassOp(Wℓ::InterpolatedCℓs) = Diagonal(BandPass(Wℓ))
HighPass(ℓ;Δℓ=50) = BandPassOp(0:10000, [zeros(ℓ-Δℓ); @.((cos($range(π,0,length=2Δℓ))+1)/2); ones(10001-ℓ-Δℓ)])
LowPass(ℓ;Δℓ=50) = BandPassOp(0:(ℓ+Δℓ-1), [ones(ℓ-Δℓ); @.(cos($range(0,π,length=2Δℓ))+1)/2])
MidPass(ℓmin,ℓmax;Δℓ=50)  = BandPassOp(0:(ℓmax+Δℓ-1), [zeros(ℓmin-Δℓ);  @.(cos($range(π,0,length=2Δℓ))+1)/2; ones(ℓmax-ℓmin-2Δℓ); @.((cos($range(0,π,length=2Δℓ))+1)/2)])

### ParamDependentOp

# A LinOp which depends on some parameters, θ.
# L(;θ...) recomputes the operator at a given set of parameters, but the
# operator can also be used as-is in which case it is evaluated at a fiducial θ
# (which is stored inside the operator when it is first constructed).

@doc doc"""
    ParamDependentOp(recompute_function::Function)
    ParamDependentOp(recompute_function!::Function, mem)

Creates an operator which depends on some parameters $\theta$ and can be
evaluated at various values of these parameters.

There are two forms to construct this operator. In the first form,
`recompute_function` should be a function which accepts keyword arguments for
$\theta$ and returns the operator. Each keyword must have a default value; the
operator will act as if evaluated at these defaults unless it is explicitly
evaluated at other parameters. In the second form, we can preallocate some
memory for the results `mem`, in which case `recompute_function!` should
additionally accept a single positional argument holding this memory, which
should then be assigned in-place.

Example:

```julia
Cϕ₀ = Diagonal(...) # some fixed Diagonal operator
Cϕ = ParamDependentOp((;Aϕ=1)->Aϕ*Cϕ₀) # create ParamDependentOp

Cϕ(Aϕ=1.1) * ϕ   # Cϕ(Aϕ=1.1) is equal to 1.1*Cϕ₀
Cϕ * ϕ           # Cϕ alone will act like Cϕ(Aϕ=1) because that was the default above

# a version which preallocates the memory:
Cϕ = ParamDependentOp((mem;Aϕ=1)->(@. mem = Aϕ*Cϕ₀), similar(Cϕ₀))
```

Note: if you are doing parallel work, global variables referred to in the
`recompute_function` need to be distributed to all workers. A more robust
solution is to avoid globals entirely and instead ensure all variables are
"closed" over (and hence will automatically get distributed). This will happen
by default if defining the `ParamDependentOp` inside any function, or can be
forced at the global scope by wrapping everything in a `let`-block, e.g.:

```julia
Cϕ = let Cϕ₀=Cϕ₀
    ParamDependentOp((;Aϕ=1)->Aϕ*Cϕ₀)
end
```

After executing the code above, `Cϕ` is now ready to be shipped to any workers
and will work regardless of what global variables are defined on these workers.
"""
struct ParamDependentOp{B, S, P, L<:LinOp{B,S,P}, F<:Function} <: ImplicitOp{B,S,P}
    op::L
    recompute_function::F
    parameters::Vector{Symbol}
    inplace::Bool
end
function ParamDependentOp(recompute_function::Function)
    ParamDependentOp(recompute_function(), recompute_function, get_kwarg_names(recompute_function), false)
end
function ParamDependentOp(recompute_function!::Function, mem)
    op = recompute_function!(similar(mem))
    ParamDependentOp(op, (mem=mem;θ...)->(recompute_function!(mem;θ...);mem), get_kwarg_names(recompute_function!), true)
end
function (L::ParamDependentOp)(mem=nothing;θ...)
    if (mem!=nothing)
        L.inplace || throw(ArgumentError("Can't pass preallocated memory to non-inplace ParamDependentOp."))
        mem isa typeof(L.op) || throw(ArgumentError("Preallocated memory passed to ParamDependentOp should be $(typeof(L.op)), not $(typeof(mem))."))
    end
    if depends_on(L,θ)
        dependent_θ = filter(((k,_),)->k in L.parameters, pairs(θ))
        # type annotation here for if any Core.Box'ed variables slipped into our recompute_function:
        L.recompute_function((mem==nothing ? () : (mem,))...; dependent_θ...) :: typeof(L.op)
    else
        L.op
    end
end
(L::ParamDependentOp)(θ::NamedTuple) = L(;θ...)
*(L::ParamDependentOp, f::Field) = L.op * f
\(L::ParamDependentOp, f::Field) = L.op \ f
for F in (:inv, :pinv, :sqrt, :adjoint, :Diagonal, :diag, :simulate, :zero, :one, :logdet)
    @eval $F(L::ParamDependentOp) = $F(L.op)
end
depends_on(L::ParamDependentOp, θ) = depends_on(L, keys(θ))
depends_on(L::ParamDependentOp, θ::Tuple) = any(L.parameters .∈ Ref(θ))
depends_on(L,                   θ) = false

adapt_structure(to, L::ParamDependentOp) =
    ParamDependentOp(adapt(to, L.op), adapt(to, L.recompute_function), L.parameters, L.inplace)

### LazyBinaryOp

# we use LazyBinaryOps to create new operators composed from other operators
# which don't actually evaluate anything until they've been multiplied by a
# field
struct LazyBinaryOp{F,A<:Union{LinOrAdjOp,Scalar},B<:Union{LinOrAdjOp,Scalar}} <: ImplicitOp{Basis,Spin,Pix}
    a::A
    b::B
    LazyBinaryOp(op,a::A,b::B) where {A,B} = new{op,A,B}(a,b)
end
# creating LazyBinaryOps
for op in (:+, :-, :*)
    @eval ($op)(a::ImplicitOrAdjOp,          b::ImplicitOrAdjOp)          = LazyBinaryOp($op,a,b)
    @eval ($op)(a::Union{LinOrAdjOp,Scalar}, b::ImplicitOrAdjOp)          = LazyBinaryOp($op,a,b)
    @eval ($op)(a::ImplicitOrAdjOp,          b::Union{LinOrAdjOp,Scalar}) = LazyBinaryOp($op,a,b)
    # explicit vs. lazy binary operations on Diagonals:
    @eval ($op)(D1::DiagOp{<:Field{B}},  D2::DiagOp{<:Field{B}})  where {B}     = Diagonal(broadcast($op,D1.diag,D2.diag))
    @eval ($op)(D1::DiagOp{<:Field{B1}}, D2::DiagOp{<:Field{B2}}) where {B1,B2} = LazyBinaryOp($op,D1,D2)
end
/(op::ImplicitOrAdjOp, n::Real) = LazyBinaryOp(/,op,n)
literal_pow(::typeof(^), op::ImplicitOrAdjOp, ::Val{-1}) = inv(op)
literal_pow(::typeof(^), op::ImplicitOrAdjOp, ::Val{n}) where {n} = LazyBinaryOp(^,op,n)
^(op::ImplicitOrAdjOp, n::Int) = LazyBinaryOp(^,op,n)
inv(op::ImplicitOrAdjOp) = LazyBinaryOp(^,op,-1)
-(op::ImplicitOrAdjOp) = -1 * op
# evaluating LazyBinaryOps
for op in (:+, :-)
    @eval *(lz::LazyBinaryOp{$op}, f::Field) = ($op)(lz.a * f, lz.b * f)
end
*(lz::LazyBinaryOp{/}, f::Field) = (lz.a * f) / lz.b
*(lz::LazyBinaryOp{*}, f::Field) = lz.a * (lz.b * f)
\(lz::LazyBinaryOp{*}, f::Field) = lz.b \ (lz.a \ f)
*(lz::LazyBinaryOp{^}, f::Field) = foldr((lz.b>0 ? (*) : (\)), fill(lz.a,abs(lz.b)), init=f)
adjoint(lz::LazyBinaryOp{F}) where {F} = LazyBinaryOp(F,adjoint(lz.b),adjoint(lz.a))
ud_grade(lz::LazyBinaryOp{op}, args...; kwargs...) where {op} = LazyBinaryOp(op,ud_grade(lz.a,args...;kwargs...),ud_grade(lz.b,args...;kwargs...))
adapt_structure(to, lz::LazyBinaryOp{op}) where {op} = LazyBinaryOp(op, adapt(to,lz.a), adapt(to,lz.b))

### OuterProdOp

# an operator L which represents L = M*M'
# this could also be represented by a LazyBinaryOp, but this allows us to
# do simulate(L) = M * whitenoise
struct OuterProdOp{TM<:LinOp} <: ImplicitOp{Basis,Spin,Pix}
    M::TM
end
simulate(L::OuterProdOp{<:DiagOp{F}}) where {F} = L.M * white_noise(F)
simulate(L::OuterProdOp{<:LazyBinaryOp{*}}) = L.M.a * sqrt(L.M.b) * simulate(L.M.b)
pinv(L::OuterProdOp{<:LazyBinaryOp{*}}) = OuterProdOp(pinv(L.M.a)' * pinv(L.M.b)')
*(L::OuterProdOp, f::Field) = L.M * (L.M' * f)
\(L::OuterProdOp, f::Field) = L.M' \ (L.M \ f)
adjoint(L::OuterProdOp) = L
adapt_structure(to, L::OuterProdOp) = OuterProdOp(adapt(to, L.M))


text:
const SuffixVector = Vector{Vector{DNA}}

mutable struct RNASuffixArray{A<:BioSequences.Alphabet,I<:Integer,K<:Integer}
   sai::Vector{I}
   meta::Vector{K}
   depth::SuffixVector
   species::Vector{SuffixVector}
   length::Int
   isphylo::Bool

   function RNASuffixArray( seq::BioSequences.Sequence, len::Int, metadata::K=one(K))
      sai   = SuffixArrays.suffixsort( String(seq) ).index + one(I)
      meta  = K[metadata for i in 1:length(sai)]
      depth = SuffixVector(len)
      for i in 1:len
         depth[i] = map(x->x+i-1 > length(seq) ? gap(A) : seq[x+i-1], sai)
      end
      empty = Vector{SuffixVector}()
      new( sai, meta, depth, empty, len, false )
   end

   function RNASuffixArray( maf::MAFRecord, tree::PhyloTree, len::Int, metadata::K=one(K))
      deletegaps!( maf )
      for i in 2:length(maf.species)
         @assert length(maf.species[i].sequence) == length(maf.species[1].sequence)
      end
      sai   = SuffixArrays.suffixsort( String(maf[1].sequence) ).index + one(I)
      meta  = K[metadata for i in 1:length(sai)]
      depth = SuffixVector(len)
      for i in 1:len
         depth[i] = map(x->x+i-1 > length(maf[1].sequence) ? gap(A) : maf.species[1].sequence[x+i-1], sai)
      end
      species = Vector{SuffixVector}(length(tree.order)-1)
      m = 2
      for i in 1:length(species)
         curind = tree.index[maf[m].name] - 1
         species[i] = SuffixVector(len)
         if i != curind
            for j in 1:len
               species[i][j] = [gap(A) for i in 1:length(sai)]
            end
         else
           for j in 1:len
              species[i][j] = map(x->x+j-1 > length(maf[m].sequence) ? gap(A) : maf[m].sequence[x+j-1], sai)
           end
           m < length(maf) && (m += 1)
         end
      end
      new( sai, meta, depth, species, len, true )
   end
end

function Base.:(<=)( a::SuffixVector, ai, b::SuffixVector, bi )
   for i in 1:length(a)
      if a[i][ai] < b[i][bi]
         return true
      elseif a[i][ai] > b[i][bi]
         return false
      end
   end
   true
end

function Base.push!{A,I,K}( suf::RNASuffixArray{A,I,K}, seq::BioSequences.Sequence;
                            metadata::K=convert(K,maximum(suf.meta)+1) )

   suf.isphylo && error("Cannot push! single BioSequence when isphylo is enabled!")

   # suffix sort
   const suftmp  = RNASuffixArray{A,I,K}( seq, suf.length, metadata )

   # initialize
   const newsai    = zeros(I, length(suf.sai) + length(suftmp.sai) )
   const newmeta   = zeros(K, length(suf.sai) + length(suftmp.sai) )
   const newdepth  = SuffixVector(suf.length)
   for i in 1:suf.length
      newdepth[i] = Vector{DNA}(length(suf.sai) + length(suftmp.sai))
   end

   # merge
   i,j,k = 1,1,1
   while i <= length(suf.sai) || j <= length(suftmp.sai)
      if j > length(suftmp.sai) || (i <= length(suf.sai) && # does suf have next least suffix?
                                    <=( suf.depth, i, suftmp.depth, j ))
         newsai[k]  = suf.sai[i]
         newmeta[k] = suf.meta[i]
         for d in 1:suf.length
            newdepth[d][k] = suf.depth[d][i]
         end
         i += 1
      else # suftmp has least suffix
         newsai[k]  = suftmp.sai[j]
         newmeta[k] = suftmp.meta[j]
         for d in 1:suf.length
            newdepth[d][k] = suftmp.depth[d][j]
         end
         j += 1
      end
      k += 1
   end

   # set new values of suf
   suf.sai     = newsai
   suf.meta    = newmeta
   suf.depth   = newdepth
   suf
end

function Base.push!{A,I,K}( suf::RNASuffixArray{A,I,K}, maf::MAFRecord, tree::PhyloTree;
                            metadata::K=convert(K,maximum(suf.meta)+1) )

   suf.isphylo || error("Cannot push! MAF block when isphylo isn't enabled!")
   length(suf.species) == length(tree.order)-1 || error("# species in tree not equal to suffix array!")

   # suffix sort
   const suftmp  = RNASuffixArray{A,I,K}( maf, tree, suf.length, metadata )

   # initialize
   const newsai     = zeros(I, length(suf.sai) + length(suftmp.sai) )
   const newmeta    = zeros(K, length(suf.sai) + length(suftmp.sai) )
   const newdepth   = SuffixVector(suf.length)
   const newspecies = Vector{SuffixVector}(length(tree.order)-1)

   for i in 1:suf.length
      newdepth[i] = Vector{DNA}(length(suf.sai) + length(suftmp.sai))
   end
   for i in 1:length(newspecies)
      newspecies[i] = SuffixVector(suf.length)
      for j in 1:suf.length
         newspecies[i][j] = Vector{DNA}(length(suf.sai) + length(suftmp.sai))
      end
   end

   # merge
   i,j,k = 1,1,1
   while i <= length(suf.sai) || j <= length(suftmp.sai)
      if j > length(suftmp.sai) || (i <= length(suf.sai) && # does suf have next least suffix?
                                    <=( suf.depth, i, suftmp.depth, j ))
         newsai[k]  = suf.sai[i]
         newmeta[k] = suf.meta[i]
         for d in 1:suf.length
            newdepth[d][k] = suf.depth[d][i]
         end
         for s in 1:length(suf.species)
            for d in 1:suf.length
               newspecies[s][d][k] = suf.species[s][d][i]
            end
         end
         i += 1
      else # suftmp has least suffix
         newsai[k]  = suftmp.sai[j]
         newmeta[k] = suftmp.meta[j]
         for d in 1:suf.length
            newdepth[d][k] = suftmp.depth[d][j]
         end
         for s in 1:length(suf.species)
            for d in 1:suf.length
               newspecies[s][d][k] = suftmp.species[s][d][j]
            end
         end
         j += 1
      end
      k += 1
   end

   # set new values of suf
   suf.sai     = newsai
   suf.meta    = newmeta
   suf.depth   = newdepth
   suf.species = newspecies
   suf
end

struct RNAGenomeCoord
   name::String
   offset::Int
   strand::Bool
   length::Int
end

mutable struct RNADuplexArray{A,I,K}
   fwd::RNASuffixArray{A,I,K}
   rev::RNASuffixArray{A,I,K}
   seqs::Vector{BioSequences.Sequence}
   coords::Vector{RNAGenomeCoord}
   depth::Int
   isphylo::Bool

   function RNADuplexArray( seq::BioSequences.Sequence, depth::Int;
                            seqname::String="chr", genomepos::Int=1, strand::Bool=true )
      fwd    = RNASuffixArray{A,I,K}( seq, depth )
      rev    = RNASuffixArray{A,I,K}( reverse(seq), depth )
      seqs   = BioSequences.Sequence[seq]
      coords = RNAGenomeCoord[RNAGenomeCoord(seqname, genomepos, strand, length(seq))]
      return new( fwd, rev, seqs, coords, depth, false )
   end

   function RNADuplexArray( left::BioSequences.Sequence, right::BioSequences.Sequence, depth::Int;
                            left_seqname::String="chr", right_seqname::String="chr",
                            left_genomepos::Int=1, right_genomepos::Int=1,
                            left_strand::Bool=true, right_strand::Bool=true )
      fwd    = RNASuffixArray{A,I,K}( left,  depth, convert(K, 1) )
      rev    = RNASuffixArray{A,I,K}( reverse(right), depth, convert(K, 2) )
      seqs   = BioSequences.Sequence[left, right]

      coords = RNAGenomeCoord[RNAGenomeCoord(left_seqname, left_genomepos, left_strand, length(left)),
                              RNAGenomeCoord(right_seqname, right_genomepos, right_strand, length(right))]
      return new( fwd, rev, seqs, coords, depth, false )
   end

   function RNADuplexArray( seqs::Vector{BioSequences.Sequence}, depth::Int;
                            seqnames::Vector{String}=String["chr" for i in 1:length(seqs)],
                            genomepos::Vector{Int}=ones(Int, length(seqs)),
                            strands::Vector{Bool}=trues(length(seqs)))
      length(seqs) >= 1 || error("Cannot build RNADuplexArray from seq vector of length 0!")
      length(seqs) == length(seqnames) == length(genomepos) == length(strands) || error("all RNADuplexArray input vectors must be equal size!")
      fwd    = RNASuffixArray{A,I,K}( seqs[1], depth )
      rev    = RNASuffixArray{A,I,K}( reverse(seqs[1]), depth )
      coords = RNAGenomeCoord[RNAGenomeCoord(seqnames[i], genomepos[i], strands[i], length(seqs[i])) for i in 1:length(seqs)]
      for i in 2:length(seqs)
         push!( fwd, seqs[i] )
         push!( rev, reverse(seqs[i]) )
      end
      return new( fwd, rev, seqs, coords, depth, false )
   end

   function RNADuplexArray( maf::MAFRecord, tree::PhyloTree, depth::Int )
      refmaf = maf[1]
      fwd    = RNASuffixArray{A,I,K}( maf, tree, depth )
      reverse!( maf )
      rev    = RNASuffixArray{A,I,K}( maf, tree, depth )
      seqs   = BioSequences.Sequence[refmaf.sequence]
      coords = RNAGenomeCoord[RNAGenomeCoord(refmaf.chr, refmaf.position, refmaf.strand, length(refmaf.sequence))]
      return new( fwd, rev, seqs, coords, depth, true )
   end
end

function Base.push!{A,I,K}( rda::RNADuplexArray{A,I,K}, seq::BioSequences.Sequence;
                            seqname::String="chr", genomepos::Int=1, strand::Bool=true )
   rda.isphylo && error("Cannot push single sequence to RNADuplexArray with isphylo=true!")
   push!( rda.fwd, seq )
   push!( rda.rev, reverse(seq) )
   push!( rda.seqs, seq )
   push!( rda.coords, RNAGenomeCoord(seqname, genomepos, strand, length(seq)) )
end

function Base.push!{A,I,K}( rda::RNADuplexArray{A,I,K}, maf::MAFRecord, tree::PhyloTree )
   refmaf = maf[1]
   push!( rda.fwd, maf, tree )
   reverse!( maf )
   push!( rda.rev, maf, tree )
   push!( rda.seqs, refmaf.sequence )
   push!( rda.coords, RNAGenomeCoord(refmaf.chr, refmaf.position, refmaf.strand, length(refmaf.sequence)) )
end

tochar( strand::Bool ) = strand ? '+' : '-'

function positions{N <: NucleicAcid}( nucs::Vector{N}, lo::Int, hi::Int, alpha::Tuple{Vararg{N}} )
   map( x->searchsorted(nucs, x, lo, hi, Base.Order.ForwardOrdering()), alpha )
end

function traverse{A,I,K}( dsa::RNADuplexArray{A,I,K}, foldrange::UnitRange=3:typemax(Int);
                          bulge_max::Int=zero(Int), mismatch_max::Int=zero(Int),
                          raw_output::String="", minfold::Float64=-5.0,
                          single::PhyloTree=EMPTY_TREE, paired::PhyloTree=EMPTY_TREE)

   const pairs_idx      = index(A, pairs)
   const bulges_idx     = index(A, gaps)
   const mismatches_idx = index(A, mismatches)

   const alpha          = alphabet(A)

   const duplex         = RNADuplex()
   const deprange       = 1:dsa.depth
   const energy_max     = 0.01

   const intervals      = DuplexCollection{K}()

   if raw_output != ""
      const io = open( raw_output, "w" )
      const stream = ZlibDeflateOutputStream( io )
      const output = true
   else
      const output = false
   end

   @inline function _traverse( fdepth::Int64, rdepth::Int64,
                               frange::UnitRange{Int}, rrange::UnitRange{Int},
                               bulge_n::Int64, mismatch_n::Int64,
                               from_bulge::Bool, bulge_left::Bool,
                               from_mismatch::Bool )

      # obtain position range for current depth
      if fdepth <= dsa.fwd.length && rdepth <= dsa.rev.length
         const franges = positions( dsa.fwd.depth[fdepth], first(frange), last(frange), alpha )
         const rranges = positions( dsa.rev.depth[rdepth], first(rrange), last(rrange), alpha )
      else
         pop!(duplex)
         return
      end

      # recurse through valid watson-crick pairs
      for (l,r) in pairs_idx
         if length(franges[l]) >= 1 && length(rranges[r]) >= 1
            push!( duplex, convert(RNAPair, onehot(l), onehot(r)) )
            if output
               tab_write( stream, (depth, energy(duplex)))
            end
            if fdepth in deprange && rdepth in deprange && energy(duplex) < minfold
               @inbounds for ix in franges[l], jx in rranges[r]
                     const i = dsa.fwd.sai[ix]
                     const j = dsa.rev.sai[jx]
                     const fwd_entry = dsa.coords[ dsa.fwd.meta[ix] ]
                     const rev_entry = dsa.coords[ dsa.rev.meta[jx] ]
                     const ik = revoffset(i, fwd_entry.length)
                     const jk = revoffset(j, rev_entry.length)
                     const ifirst,ilast = fwd_entry.strand ? (i,i+fdepth-1) : (ik-fdepth+1,ik)
                     const jfirst,jlast = rev_entry.strand ? (jk-rdepth+1,jk) : (j,j+rdepth-1)
                     const fg = fwd_entry.offset-1
                     const rg = rev_entry.offset-1
                     fwd_entry.name == rev_entry.name && !(((rg+jfirst) - (fg+ilast)) - 1 in foldrange) && continue
                     if dsa.isphylo
                        evo = paired == EMPTY_TREE ? EvoDuplex(duplex, single, dsa, ix, jx) :
                                                     EvoDuplex(duplex, single, paired, dsa, ix, jx)
                        score(evo) < 0.0 && continue
                        push!( intervals, DuplexInterval( Interval(fwd_entry.name, fg+ifirst, fg+ilast, tochar(fwd_entry.strand), dsa.fwd.meta[ix]),
                                                          Interval(rev_entry.name, rg+jfirst, rg+jlast, tochar(rev_entry.strand), dsa.rev.meta[jx]),
                                                          evo ) )
                     else
                        push!( intervals, DuplexInterval( Interval(fwd_entry.name, fg+ifirst, fg+ilast, tochar(fwd_entry.strand), dsa.fwd.meta[ix]),
                                                          Interval(rev_entry.name, rg+jfirst, rg+jlast, tochar(rev_entry.strand), dsa.rev.meta[jx]),
                                                          deepcopy(duplex) ) )
                     end
               end
            end
            _traverse(  fdepth + 1, rdepth + 1,
                        franges[l], rranges[r],
                        bulge_n, mismatch_n,
                        false, false, false )
         end
      end

      # recurse through bulges
      if bulge_n < bulge_max &&
         (bulge_n >= 1 ? from_bulge : !from_bulge) &&
         fdepth > 1 && rdepth > 1 && duplex.energy[end] < energy_max
         for (l,r) in bulges_idx
            if l == 0
               (from_bulge && !bulge_left) && continue # only bulge one way
               if length(rranges[r]) >= 1 && rdepth + 1 <= dsa.rev.length
                  push!( duplex, convert(RNABulge, zero(UInt8), onehot(r)) )
                  _traverse( fdepth, rdepth + 1,
                             frange, rranges[r],
                             bulge_n + 1, mismatch_n,
                             true, true, from_mismatch )
               end
            elseif r == 0
               (from_bulge && bulge_left) && continue
               if length(franges[l]) >= 1 && fdepth + 1 <= dsa.fwd.length
                   push!( duplex, convert(RNABulge, onehot(l), zero(UInt8)) )
                   _traverse( fdepth + 1, rdepth,
                              franges[l], rrange,
                              bulge_n + 1, mismatch_n,
                              true, false, from_mismatch )
               end
            end
         end
      end

      # recurse through mismatches
      if mismatch_n < mismatch_max &&
         (mismatch_n >= 1 ? from_mismatch : !from_mismatch) &&
         fdepth > 1 && rdepth > 1 &&
         !from_bulge && duplex.energy[end] < energy_max &&
         fdepth < dsa.fwd.length && rdepth < dsa.rev.length
         for (l,r) in mismatches_idx
            if length(franges[l]) >= 1 && length(franges[r]) >= 1
                push!( duplex, convert(RNAMismatch, onehot(l), onehot(r)) )
                _traverse( fdepth + 1, rdepth + 1,
                           franges[l], rranges[r],
                           bulge_n, mismatch_n + 1,
                           from_bulge, bulge_left, true )
            end
         end
      end

      pop!( duplex ) # clean up this level
      return
   end

   _traverse( one(Int), one(Int),
              1:length(dsa.fwd.sai), 1:length(dsa.rev.sai),
              zero(Int), zero(Int),
              false, false, false )

   if output
      close(stream)
      close(io)
   end

   intervals
end


text:
@testset "Types" begin
    @testset "Ground Types" begin
        @test ispassive(Clock()) == true
        @test ispassive(UInteger()) == true
        @test ispassive(SInteger()) == true
        @test ispassive(Analog()) == true
        constructors = [Clock, () -> UInteger(rand(UInt)), () -> SInteger(rand(UInt)), () -> Analog(rand(UInt))]
        @testset for (x, y) in Base.product(constructors, constructors)
            if x == y
                @test (x() == y()) == true
                @test (x() ≈ y()) == true
            else
                @test (x() == y()) == false
                @test (x() ≈ y()) == false
            end
        end
    end

    @testset "Vector Types" begin
        @test ispassive([Clock(), Clock()]) == true
        @test ([Clock(), Clock()] == [Clock(), Clock()]) == true
        @test ([Clock(), Clock()] == [Clock(), UInteger()]) == false
        @test ([Clock(), Clock()] ≈ [Clock(), Clock()]) == true
        @test ([Clock(), Clock()] ≈ [Clock(), UInteger()]) == false
    end

    @testset "Bundle Types" begin
        x = Bundle([:a => UInteger(), :b => UInteger()])
        y = Bundle([:a => UInteger(), :b => UInteger()])
        w = Bundle([:a => (UInteger(), true), :b => (UInteger(), false)])
        z = Bundle([:a => SInteger(), :b => UInteger()])
        q = Bundle([:b => UInteger(), :a => UInteger()])
        @test ispassive(x) == true
        @test (x == y) == true
        @test (x == w) == false
        @test (x == z) == false
        @test (x ≈ q) == true
        @test (x ≈ z) == false
    end
end