text:
# *********************************************************************************
# REopt, Copyright (c) 2019-2020, Alliance for Sustainable Energy, LLC.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification,
# are permitted provided that the following conditions are met:
#
# Redistributions of source code must retain the above copyright notice, this list
# of conditions and the following disclaimer.
#
# Redistributions in binary form must reproduce the above copyright notice, this
# list of conditions and the following disclaimer in the documentation and/or other
# materials provided with the distribution.
#
# Neither the name of the copyright holder nor the names of its contributors may be
# used to endorse or promote products derived from this software without specific
# prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
# IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
# OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
# OF THE POSSIBILITY OF SUCH DAMAGE.
# *********************************************************************************
using CPLEX

#=
add a time-of-export rate that is greater than retail rate for the month of January,
check to make sure that PV does NOT export unless the site load is met first for the month of January.
=#
@testset "Do not allow_simultaneous_export_import" begin
    model = Model(optimizer_with_attributes(CPLEX.Optimizer, "CPX_PARAM_SCRIND" => 0))
    data = JSON.parsefile("./scenarios/monthly_rate.json")

    # create wholesale_rate with compensation in January > retail rate
    jan_rate = data["ElectricTariff"]["monthly_energy_rates"][1]
    data["ElectricTariff"]["wholesale_rate"] =
        append!(repeat([jan_rate + 0.1], 31 * 24), repeat([0.0], 8760 - 31*24))
    data["ElectricTariff"]["monthly_demand_rates"] = repeat([0], 12)
    data["ElectricUtility"] = Dict("allow_simultaneous_export_import" => false)

    s = Scenario(data)
    inputs = REoptInputs(s)
    results = run_reopt(model, inputs)

    @test all(x == 0.0 for (i,x) in enumerate(results["ElectricUtility"]["year_one_to_load_series_kw"][1:744])
              if results["PV"]["year_one_to_grid_series_kw"][i] > 0)
end

@testset "Solar and Storage" begin
    model = Model(optimizer_with_attributes(CPLEX.Optimizer, "CPX_PARAM_SCRIND" => 0))
    results = run_reopt(model, "./scenarios/pv_storage.json")

    @test results["PV"]["size_kw"] ≈ 217 atol=1
    @test results["Financial"]["lcc_us_dollars"] ≈ 1.23887e7 rtol=1e-5
    @test results["Storage"]["size_kw"] ≈ 56 atol=1
    @test results["Storage"]["size_kwh"] ≈ 79 atol=1
end

@testset "Minimize Unserved Load" begin
    m = Model(optimizer_with_attributes(CPLEX.Optimizer, "CPX_PARAM_SCRIND" => 0))
    results = run_reopt(m, "./scenarios/outage.json")

    @test results["expected_outage_cost"] ≈ 0
    @test sum(results["unserved_load_per_outage"]) ≈ 0
    @test value(m[:binMGTechUsed]["Generator"]) == 1
    @test value(m[:binMGTechUsed]["PV"]) == 0
    @test value(m[:binMGStorageUsed]) == 1
    @test results["Financial"]["lcc_us_dollars"] ≈ 7.3681609e7 atol=5e4

    #=
    Scenario with $0/kWh VoLL, 12x169 hour outages, 1kW load/hour, and min_resil_timesteps = 168
    - should meet 168 kWh in each outage such that the total unserved load is 12 kWh
    =#
    m = Model(optimizer_with_attributes(CPLEX.Optimizer, "CPX_PARAM_SCRIND" => 0))
    results = run_reopt(m, "./scenarios/nogridcost_minresilhours.json")
    @test sum(results["unserved_load_per_outage"]) ≈ 12

    # testing dvUnserved load, which would output 100 kWh for this scenario before output fix
    m = Model(optimizer_with_attributes(CPLEX.Optimizer, "CPX_PARAM_SCRIND" => 0))
    results = run_reopt(m, "./scenarios/nogridcost_multiscenario.json")
    @test sum(results["unserved_load_per_outage"]) ≈ 60
end

@testset "Multiple Sites" begin
    m = Model(optimizer_with_attributes(CPLEX.Optimizer, "CPX_PARAM_SCRIND" => 0))
    ps = [
        REoptInputs("./scenarios/pv_storage.json"),
        REoptInputs("./scenarios/monthly_rate.json"),
    ];
    results = run_reopt(m, ps)
    @test results[3]["Financial"]["lcc_us_dollars"] + results[10]["Financial"]["lcc_us_dollars"] ≈ 1.23887e7 + 437169.0 rtol=1e-5
end

@testset "LinDistFlow" begin
    m = Model(optimizer_with_attributes(CPLEX.Optimizer, "CPX_PARAM_SCRIND" => 0))
    ps = [
        REoptInputs("./scenarios/pv_storage.json"),
        REoptInputs("./scenarios/monthly_rate.json"),
    ];
    # make dummy REoptInputs with fixed systems for other nodes in network?
    load_nodes = Dict(
        "3" => [],
        "10" => [],
    )
    ldf_inputs = LinDistFlow.Inputs(
        "data/car10lines.dss",
        "0",
        "data/car10linecodes.dss";
        Pload=load_nodes,
        Qload=load_nodes,
        Sbase=1e6,
        Vbase=12.5e3,
        v0 = 1.00,
        v_uplim = 1.05,
        v_lolim = 0.95,
        Ntimesteps = 8760
    );
    build_reopt!(m, ps)
    LinDistFlow.build_ldf!(m, ldf_inputs, ps)
    add_objective!(m, ps)
    optimize!(m)

    results = reopt_results(m, ps)
    @test results[10]["Financial"]["lcc_us_dollars"] + results[3]["Financial"]["lcc_us_dollars"] ≈ 1.23887e7 + 437169.0 rtol=1e-5
    P0 = value.(m[:Pⱼ]["0",:]).data * ldf_inputs.Sbase / 1e3;  # converting to kW
    TotalGridPurchases = value.(m[:dvGridPurchase_3]).data + value.(m[:dvGridPurchase_10]).data;
    @test maximum(TotalGridPurchases) ≈ maximum(P0) rtol = 1e-5  # lossless model
end

## equivalent REopt Lite API Post for test 2:
#   NOTE have to hack in API levelization_factor to get LCC within 5e-5 (Mosel tol)
# {"Scenario": {
#     "Site": {
#         "longitude": -118.1164613,
#         "latitude": 34.5794343,
#         "roof_squarefeet": 5000.0,
#         "land_acres": 1.0,
#     "PV": {
#         "macrs_bonus_pct": 0.4,
#         "installed_cost_us_dollars_per_kw": 2000.0,
#         "tilt": 34.579,
#         "degradation_pct": 0.005,
#         "macrs_option_years": 5,
#         "federal_itc_pct": 0.3,
#         "module_type": 0,
#         "array_type": 1,
#         "om_cost_us_dollars_per_kw": 16.0,
#         "macrs_itc_reduction": 0.5,
#         "azimuth": 180.0,
#         "federal_rebate_us_dollars_per_kw": 350.0,
#         "dc_ac_ratio": 1.1
#     },
#     "LoadProfile": {
#         "doe_reference_name": "RetailStore",
#         "annual_kwh": 10000000.0,
#         "city": "LosAngeles"
#     },
#     "Storage": {
#         "total_rebate_us_dollars_per_kw": 100.0,
#         "macrs_option_years": 5,
#         "can_grid_charge": true,
#         "macrs_bonus_pct": 0.4,
#         "macrs_itc_reduction": 0.5,
#         "total_itc_pct": 0,
#         "installed_cost_us_dollars_per_kw": 1000.0,
#         "installed_cost_us_dollars_per_kwh": 500.0,
#         "replace_cost_us_dollars_per_kw": 460.0,
#         "replace_cost_us_dollars_per_kwh": 230.0
#     },
#     "ElectricTariff": {
#         "urdb_label": "5ed6c1a15457a3367add15ae"
#     },
#     "Financial": {
#         "escalation_pct": 0.026,
#         "offtaker_discount_pct": 0.081,
#         "owner_discount_pct": 0.081,
#         "analysis_years": 20,
#         "offtaker_tax_pct": 0.4,
#         "owner_tax_pct": 0.4,
#         "om_cost_escalation_pct": 0.025
#     }
# }}}

max_stars_repo_name: brtietz/REoptLite

max_issues_repo_name: brtietz/REoptLite

max_forks_repo_name: brtietz/REoptLite


text:
function wavenumbers_mesh(ω::T, k_effs::Vector{Complex{T}}, medium::Medium{T}, species::Vector{Specie{T}};
        dim = 2,
        tol::T = 1e-6,
        mesh_refine::T = T(0.4),
        inner_optimizer = LBFGS(),
        verbose::Bool = false,
        max_Rek::T = maximum(real(k1) for k1 in k_effs),
        min_Rek::T = minimum(real(k1) for k1 in k_effs),
        max_Imk::T = maximum(imag(k1) for k1 in k_effs),
        min_Imk::T = minimum(imag(k1) for k1 in k_effs[2:end]),
        kws...) where T<:Number

    low_tol = max(1e-4, tol)*minimum( (k1 == k2) ? Inf : abs(k1-k2) for k1 in k_effs, k2 in k_effs) # a tolerance used for a first pass with time_limit

    # the dispersion equation is given by: `dispersion(k1,k2) = 0` where k_eff = k1 + im*k2.
    dispersion = dispersion_function(ω, medium, species; tol = low_tol, dim=dim, kws...)

    kx = LinRange(min_Rek, max_Rek, Int(round(length(k_effs)/(2*mesh_refine)))) # tree shape makes this division by 2 natural
    ky = LinRange(min_Imk, max_Imk, Int(round(length(k_effs)/(2*mesh_refine))))
    k_mesh = [[x,y] for x in kx, y in ky][:]

    k_vecs = [[real(keff),imag(keff)] for keff in k_effs]

    # make slightly bigger box to be constrained within
    min_Rek = (min_Rek < zero(T)) ? min_Rek*T(1.1) : min_Rek*T(0.9)
    min_Imk = (min_Imk < zero(T)) ? min_Imk*T(1.1) : min_Imk*T(0.9)
    max_Rek = (max_Rek > zero(T)) ? max_Rek*T(1.001) : max_Rek*T(0.999)
    max_Imk = (max_Imk > zero(T)) ? max_Imk*T(1.001) : max_Imk*T(0.999)

    # Find all wavenumbers
    lower = [min_Rek, min_Imk]
    upper = [max_Rek, max_Imk]
    new_ks = [
        optimize(dispersion, lower, upper, kvec,
            Fminbox(inner_optimizer), Optim.Options(x_tol=low_tol, g_tol = low_tol^3)).minimizer
    for kvec in k_mesh]
    deleteat!(new_ks, findall(dispersion.(new_ks) .> low_tol))
    new_ks = reduce_kvecs(new_ks, low_tol/10)

    # only keep targets which are not already in k_vecs
    new_ks = [
            (findmin([norm(h - kvec) for kvec in k_vecs])[1] > low_tol) ? h : [zero(T),-one(T)]
    for h in new_ks]
    new_ks = reduce_kvecs(new_ks, low_tol)

    # Here we refine the new roots
    new_ks = map(new_ks) do k_vec
        # res = optimize(dispersion, k_vec; g_tol = tol^2.0, f_tol = tol^4.0, x_tol=tol)
        res = optimize(dispersion, lower, upper, k_vec, Fminbox(inner_optimizer),
            Optim.Options(g_tol = tol^3.0, x_tol=tol))
        if res.minimum > T(20)*tol
            [zero(T),-one(T)]
        else
            res.minimizer
        end
    end

    new_ks = [
            (findmin([norm(h - kvec) for kvec in k_vecs])[1] > 10*tol) ? h : [zero(T),-one(T)]
    for h in new_ks]
    new_ks = reduce_kvecs(new_ks, tol)

    if verbose println("New roots from mesh refiner:",new_ks) end

    # group together wavenumbers which are closer than tol
    k_vecs = reduce_kvecs([new_ks; k_vecs], tol)

    # Finally delete unphysical waves, including waves travelling backwards with almost no attenuation. This only is important in the limit of very low frequency or very weak scatterers.
    deleteat!(k_vecs, findall([k_vec[2] < -low_tol for k_vec in k_vecs]))
    deleteat!(k_vecs, findall([-low_tol < abs(k_vec[2])/k_vec[1] < zero(T) for k_vec in k_vecs]))
    # deleteat!(k_vecs, find(k_vec[2] < tol && k_vec[1] < tol for k_vec in k_vecs))

    k_effs = [k_vec[1] + k_vec[2]*im for k_vec in k_vecs]
    k_effs = sort(k_effs, by=imag)

    return k_effs
end

# Alternative grouping
# all_inds = collect(eachindex(k_vecs))
# k_vecs = map(k_vecs) do k_vec
#     nk = norm(k_vec)*sqrt(tol)
#     ind_ins = find(norm(ks - k_vec) < nk for ks in k_vecs[all_inds])
#     inds = all_inds[ind_ins]
#     deleteat!(all_inds,ind_ins)
#     isempty(inds) ? [zero(T),-one(T)] :  mean(k_vecs[inds])
# end

# alternative grouping methods

    # digs = 3 - Int(round(log(10,sqrt(tol)))) # number of digits to round to check if equal
    # k_vecs = map(groupby(k_vec -> round.(k_vec,digs), k_vecs)) do g
    #     mean(g)
    # end

# function detMM!(F,x)
#     F[1] = abs(det(MM(x[1]+im*x[2])))
# end

# Alternative solvers
# res = nlsolve(detMM!,initial_k_eff; iterations = 10000, factor=2.0)
# k_eff_nl = res.zero[1] + im*res.zero[2]
# lower = [0.,0.]; upper = [T(2)*k0,k0]
# result = optimize(dispersion, initial_k0, lower, upper; g_tol = tol^2.0, f_tol = tol^4.0)

ext: jl


text:
"""
NARP fileformat
See EUROCONTROL NEST Manual for NARP fileformat description

Example Narp file format
8214
AGGA;AUKI/GWAUNARU'U;-8.6983333333;160.6783333333;0;0;0;0
AGGE;BALALAE;-6.9833333333;155.8850000000;0;0;0;0
AGGH;HONIARA/HENDERSON;-9.4316666667;160.0533333333;0;0;0;0
AGGJ;AVU AVU;-9.8683333333;160.4100000000;0;0;0;0
AGGK;KIRA KIRA;-10.4500000000;161.8966666667;0;0;0;0
AGGL;GRACIOSA BAY;-10.7233333333;165.7800000000;0;0;0;0
AGGM;MUNDA;-8.3250000000;157.2633333333;0;0;0;0
AGGN;GIZO/NUSATAPE;-8.1016666667;156.8350000000;0;0;0;0
AGGP;PARASI;-9.6416666667;161.4252777778;0;0;0;82
"""
module Narp

export read

using Format
using CSV
using Dates
using DataFrames

const fileformat = Dict(1=>String, 2=>String, 3=>Float64, 4=>Float64,
5=>String, 6=>String, 7=>String, 8=>String)

const header_format = ["AIRPORT_ID", "AIRPORT_NAME", "LAT_DEG", "LON_DEG",
"TIS", "TRS", "TAXITIME", "ALTITUDE_FL"]

function read(file)
    df = CSV.read(file, delim=";", types=fileformat, header=header_format,
    copycols=true, datarow=2)
    reformat!(df)
    return df
end

function reformat!(df)
    df[:,:TEMP] = string2int.(df[:,:TIS])
    select!(df, Not(:TIS))
    df[:,:TIS] =  df[:,:TEMP]
    select!(df, Not(:TEMP))

    df[:,:TEMP] = string2int.(df[:,:TRS])
    select!(df, Not(:TRS))
    df[:,:TRS] =  df[:,:TEMP]
    select!(df, Not(:TEMP))

    df[:,:TEMP] = string2int.(df[:,:TAXITIME])
    select!(df, Not(:TAXITIME))
    df[:,:TAXITIME] =  df[:,:TEMP]
    select!(df, Not(:TEMP))

    df[:,:TEMP] = string2int.(df[:,:ALTITUDE_FL])
    select!(df, Not(:ALTITUDE_FL))
    df[:,:ALTITUDE_FL] =  df[:,:TEMP]
    select!(df, Not(:TEMP))
end

function string2int(str::String)
    maybeint = tryparse(Int, str)
    if maybeint == nothing
        return missing
    elseif str == ""
        return missing
    else
        return maybeint
    end
end

end # module


text:
using Optim

λ = 400:1000
θ = [5.]
pol = 0.5
beam = PlaneWave(λ, θ; p=pol)

# Refractive indices of incident (0) and substrate (2)
incident = RIdb.air(beam.λ)
emergent = RIdb.silicon(beam.λ)

# Define the RI model to use
layers = [
    LayerTMMO(incident),
    ModelFit(:looyenga; N=(ninc=incident, nhost=emergent)),
    ModelFit(:looyenga; N=(ninc=incident, nhost=emergent)),
    LayerTMMO(emergent),
]

# Set the order of the layers (build the system) to fit
order = [
    1 # incident medium
    repeat([2, 3], 9) # top DBR of the MC
    4 # substrate
]

# Absolute reflectance spectrum to fit stored in Utils
Rexp = SpectraDB.bragg_spectrum(beam.λ)

options = Optim.Options(
    g_abstol=1e-8, g_reltol=1e-8, iterations=10^5, show_trace=true, store_trace=false,
);

# Seeds for each ModelFit layer defined above (without alpha)
seed = [
    [69, 0.54], # layers[2]
    [119.0, 0.77], # layers[3]
]

solOptim = fit_tmm_optics(
    Reflectance(Rexp), seed, beam, layers;
    order=order, options=options, alg=SAMIN(), lb=0.5.*seed, ub=1.5.*seed,
)

# Seeds for each ModelFit layer defined above plus the alpha (with alpha)
seed2 = [
    [69, 0.54], # layers[2]
    [119.0, 0.77], # layers[3]
    [0.995], # alpha
]

solOptim2 = fit_tmm_optics(
    Reflectance(Rexp), seed2, beam, layers;
    order=order, options=options, alg=SAMIN(), lb=0.5.*seed2, ub=1.5.*seed2, alpha=true,
)


text:
using DrWatson
@quickactivate "EPExperiments"

using EPExperiments
using AnnealedIS

using Distributions
using Random
using StatsFuns: logsumexp
using LinearAlgebra: dot
using Logging
using LoggingExtras
using StatsPlots

const LOG_FILENAME = "out.log"

function is_sampling(logjoint, qi, qb, num_samples)
    qi_samples = rand(qi, num_samples)
    qb_samples = rand(qb, num_samples)

    i_b_samples = map(t -> (i₀ = t[1], β = t[2]), zip(qi_samples, qb_samples))
    is_weights = logjoint.(i_b_samples) - (
        logpdf.(qi, qi_samples) .+ logpdf.(qb, qb_samples)
    )
    normalised_is_weights = is_weights .- logsumexp(is_weights)

    return normalised_is_weights, qi_samples, qb_samples
end

@with_kw struct ISGroundTruthConfig
    num_samples::Int = Int(1e6)
    qi::Distributions.Distribution = truncated(Normal(10, 10), 0, 1_000)
    qb::Distributions.Distribution = truncated(Normal(1.5, 1.5), 0, Inf)
    cost_fn_name::Symbol = :cost_fn_default
    data_filename::String = "neg_binomial_data.jld2"
    seed::Int = 1234
end

DrWatson.allignore(::ISGroundTruthConfig) = (:data_filename,)

function calculate_ground_truth(is_config)
    @unpack num_samples, qi, qb, cost_fn_name, data_filename, seed = is_config
    cost_fn = COST_FUNCTIONS[cost_fn_name]

    Random.seed!(seed)

    result_folder = mkpath(joinpath(
        datadir("is_ground_truth"),
        savename(is_config)
    ))
    logger = TeeLogger(
        ConsoleLogger(),
        FileLogger(joinpath(result_folder, LOG_FILENAME))
    )

    input_data = load(datadir("sims", is_config.data_filename))
    @unpack tmax, total_population, i_0_true, β_true, γ = input_data["params"]
    Y = input_data["Y"]

    tm = bayes_sir(Y, total_population, γ, cost_fn).gamma2
    vi = Turing.VarInfo(tm)
    joint_dens = AnnealedIS.make_log_joint_density(tm, vi)

    @time begin
    log_weights, qi_samples, qb_samples = is_sampling(
        joint_dens, qi, qb, num_samples
    )
    end

    costs = cost_fn.(qb_samples, γ)

    ess = 1 / exp(logsumexp(2 * log_weights))
    target_weights = log_weights .+ log.(costs)
    ess_target = exp(2 * logsumexp(target_weights) - logsumexp(2 * target_weights))

    exp_estimate = dot(exp.(log_weights), costs)
    with_logger(logger) do
        @info "Proposals:" qi qb
        @info "Number of samples: $num_samples"
        @info "ESS p(x|y): $ess"
        @info "ESS p(x|y)f(x): $ess_target"
        @info "Expectation estimate: $exp_estimate"
    end

    @tagsave(joinpath(result_folder, "results.jld2"), Dict(
        "log_weights" => log_weights,
        "qi_samples" => qi_samples,
        "qb_samples" => qb_samples,
        "is_config" => is_config,
        "expectation_estimate" => exp_estimate
    ))

    savefig(
        plot_is_samples(qb_samples, log_weights),
        joinpath(result_folder, "beta_samples.png")
    )
end

calculate_ground_truth(ISGroundTruthConfig(
    num_samples=Int(1e8),
    cost_fn_name=:cost_fn_shift20
))


text:
#--------------------------------------------------------------
# MultilevelTT package
# example code for https://arxiv.org/abs/1802.09062
#--------------------------------------------------------------
# Interface to AMEn implementation in ttpy (Python TT-Toolbox),
# requires additional Python packages to be installed
module TTPy

# To install the required python packages via julia:
#
# using Conda, PyCall
# Conda.add("numpy")
# Conda.add("scipy")
# Conda.add("cython")
# @pyimport pip
# args = String[]
# push!(args, "install")
# push!(args, "--user")
# push!(args, "ttpy")
# pip.main(args)

using TT, PyCall

import TT.Tensor, TT.TensorMatrix

export pyvec, pymat, Tensor, TensorMatrix, amen

@pyimport tt
@pyimport tt.amen as amen_pkg

pyvec(t::Tensor) = tt.vector[:from_list](t)
pymat(t::TensorMatrix) = tt.matrix[:from_list](t)

Tensor(x::PyObject) = Tensor(tt.vector[:to_list](x))
TensorMatrix(x::PyObject) = TensorMatrix(tt.matrix[:to_list](x))

amen(A::TensorMatrix, b::Tensor, x::Tensor, ɛ::Float64,
    kickrank::Integer = 4, nswp::Integer = 20) =
  Tensor(amen_pkg.amen_solve(pymat(A), pyvec(b), pyvec(x), ɛ, kickrank, nswp))

end


text:
bandwidths(M::MulAdd) = min.(_bnds(M), prodbandwidths(M.A,M.B))
similar(M::MulAdd{<:DiagonalLayout,<:AbstractBandedLayout}, ::Type{T}, axes::NTuple{2,OneTo{Int}}) where T =
    BandedMatrix{T}(undef, axes, bandwidths(M))
similar(M::MulAdd{<:AbstractBandedLayout,<:AbstractBandedLayout}, ::Type{T}, axes::NTuple{2,OneTo{Int}}) where T =
    BandedMatrix{T}(undef, axes, bandwidths(M))
similar(M::MulAdd{<:AbstractBandedLayout,<:DiagonalLayout}, ::Type{T}, axes::NTuple{2,OneTo{Int}}) where T =
    BandedMatrix{T}(undef, axes, bandwidths(M))
similar(M::MulAdd{<:SymmetricLayout{<:AbstractBandedLayout},<:AbstractBandedLayout}, ::Type{T}, axes::NTuple{2,OneTo{Int}}) where T =
    BandedMatrix{T}(undef, axes, bandwidths(M))
similar(M::MulAdd{<:HermitianLayout{<:AbstractBandedLayout},<:AbstractBandedLayout}, ::Type{T}, axes::NTuple{2,OneTo{Int}}) where T =
    BandedMatrix{T}(undef, axes, bandwidths(M))
similar(M::MulAdd{<:TriangularLayout{uplo,trans,<:AbstractBandedLayout},<:AbstractBandedLayout}, ::Type{T}, axes::NTuple{2,OneTo{Int}}) where {uplo,trans,T} =
    BandedMatrix{T}(undef, axes, bandwidths(M))
##
# BLAS routines
##

# make copy to make sure always works

banded_gbmv!(tA, α, A, x, β, y) =
    BLAS.gbmv!(tA, size(A,1), bandwidth(A,1), bandwidth(A,2),
                α, bandeddata(A), x, β, y)

@inline function _banded_gbmv!(tA, α, A, x, β, y)
    if x ≡ y
        banded_gbmv!(tA, α, A, copy(x), β, y)
    else
        banded_gbmv!(tA, α, A, x, β, y)
    end
end

function _banded_muladd!(α::T, A, x::AbstractVector, β, y) where T
    m, n = size(A)
    (length(y) ≠ m || length(x) ≠ n) && throw(DimensionMismatch("*"))
    l, u = bandwidths(A)
    if -l > u # no bands
        _fill_lmul!(β, y)
    elseif l < 0
        _banded_muladd!(α, view(A, :, 1-l:n), view(x, 1-l:n), β, y)
    elseif u < 0
        y[1:-u] .= zero(T)
        _banded_muladd!(α, view(A, 1-u:m, :), x, β, view(y, 1-u:m))
        y
    else
        _banded_gbmv!('N', α, A, x, β, y)
    end
end

materialize!(M::BlasMatMulVecAdd{<:BandedColumnMajor,<:AbstractStridedLayout,<:AbstractStridedLayout,T}) where T<:BlasFloat =
    _banded_muladd!(M.α, M.A, M.B, M.β, M.C)

function materialize!(M::BlasMatMulVecAdd{<:BandedRowMajor,<:AbstractStridedLayout,<:AbstractStridedLayout,T}) where T<:BlasFloat
    α, A, x, β, y = M.α, M.A, M.B, M.β, M.C
    At = transpose(A)
    m, n = size(A)
    (length(y) ≠ m || length(x) ≠ n) && throw(DimensionMismatch("*"))
    l, u = bandwidths(A)
    if -l > u # no bands
      _fill_lmul!(β, y)
    elseif l < 0
      materialize!(MulAdd(α, transpose(view(At, 1-l:n, :,)), view(x, 1-l:n), β, y))
    elseif u < 0
      y[1:-u] .= zero(T)
      materialize!(MulAdd(α, transpose(view(At, :, 1-u:m)), x, β, view(y, 1-u:m)))
      y
    else
      _banded_gbmv!('T', α, At, x, β, y)
    end
end

function materialize!(M::BlasMatMulVecAdd{<:ConjLayout{<:BandedRowMajor},<:AbstractStridedLayout,<:AbstractStridedLayout,T}) where T<:BlasComplex
    α, A, x, β, y = M.α, M.A, M.B, M.β, M.C
    Ac = A'
    m, n = size(A)
    (length(y) ≠ m || length(x) ≠ n) && throw(DimensionMismatch("*"))
    l, u = bandwidths(A)
    if -l > u # no bands
        _fill_lmul!(β, y)
    elseif l < 0
        materialize!(MulAdd(α, view(Ac, 1-l:n, :,)', view(x, 1-l:n), β, y))
    elseif u < 0
        y[1:-u] .= zero(T)
        materialize!(MulAdd(α, view(Ac, :, 1-u:m)', x, β, view(y, 1-u:m)))
        y
    else
    _banded_gbmv!('C', α, Ac, x, β, y)
    end
end

##
# Non-BLAS mat vec
##

@inline function materialize!(M::MatMulVecAdd{<:AbstractBandedLayout})
    α,A,B,β,C = M.α,M.A,M.B,M.β,M.C
    _fill_lmul!(β, C)
    @inbounds for j = 1:size(A,2), k = colrange(A,j)
        C[k] += α*inbands_getindex(A,k,j)*B[j]
    end
    C
end

@inline function materialize!(M::MatMulVecAdd{<:BandedRowMajor})
    α,At,B,β,C = M.α,M.A,M.B,M.β,M.C
    A = transpose(At)
    _fill_lmul!(β, C)

    @inbounds for j = 1:size(A,2), k = colrange(A,j)
        C[j] +=  α*transpose(inbands_getindex(A,k,j))*B[k]
    end
    C
end

@inline function materialize!(M::MatMulVecAdd{<:ConjLayout{<:BandedRowMajor}})
    α,Ac,B,β,C = M.α,M.A,M.B,M.β,M.C
    A = Ac'
    _fill_lmul!(β, C)
    @inbounds for j = 1:size(A,2), k = colrange(A,j)
        C[j] += α*inbands_getindex(A,k,j)'*B[k]
    end
    C
end

####
# Matrix * Matrix
####
function banded_mul!(C::AbstractMatrix{T}, A::AbstractMatrix, B::AbstractMatrix) where T
    Am, An = size(A)
    Bm, Bn = size(B)
    Al, Au = bandwidths(A)
    Bl, Bu = bandwidths(B)
    Cl,Cu = prodbandwidths(A, B)
    size(C) == (Am,Bn) || throw(DimensionMismatch())
    C̃l,C̃u = bandwidths(C)
    # set extra bands to zero
    @inbounds for j = 1:Bn
        for k = max(j - C̃u,1):min(j-Cu-1,Am)
            inbands_setindex!(C,zero(T),k,j)
        end
        for k = max(j + Cl+1,1):min(j+C̃l,Am)
            inbands_setindex!(C,zero(T),k,j)
        end
    end

    @inbounds for j = 1:Bn, k = max(j-Cu, 1):min(j+Cl, Am)
        νmin = max(1,k-Al,j-Bu)
        νmax = min(An,k+Au,j+Bl)

        tmp = zero(T)
        for ν=νmin:νmax
            tmp = tmp + inbands_getindex(A,k,ν) * inbands_getindex(B,ν,j)
        end
        setindex!(C,tmp,k,j)
    end
    C
end

const ConjOrBandedLayout = Union{AbstractBandedLayout,ConjLayout{<:AbstractBandedLayout}}
const ConjOrBandedColumnMajor = Union{<:BandedColumnMajor,ConjLayout{<:BandedColumnMajor}}

function _banded_muladd!(α::T, A, B::AbstractMatrix, β, C) where T
    Am, An = size(A)
    Bm, Bn = size(B)
    if An != Bm || size(C, 1) != Am || size(C, 2) != Bn
        throw(DimensionMismatch("*"))
    end

    Al, Au = bandwidths(A)
    Bl, Bu = bandwidths(B)

    gbmm!('N', 'N', α, A, B, β, C)

    C
end

materialize!(M::BlasMatMulMatAdd{<:AbstractBandedLayout,<:AbstractBandedLayout,<:BandedColumnMajor}) =
    materialize!(MulAdd(M.α, convert(DefaultBandedMatrix,M.A), convert(DefaultBandedMatrix,M.B), M.β, M.C))

materialize!(M::BlasMatMulMatAdd{<:BandedColumnMajor,<:BandedColumnMajor,<:BandedColumnMajor}) =
    _banded_muladd!(M.α, M.A, M.B, M.β, M.C)

# function generally_banded_matmatmul!(C::AbstractMatrix{T}, tA::Val, tB::Val, A::AbstractMatrix{U}, B::AbstractMatrix{V}) where {T, U, V}
#     Am, An = _size(tA, A)
#     Bm, Bn = _size(tB, B)
#     if An != Bm || size(C, 1) != Am || size(C, 2) != Bn
#         throw(DimensionMismatch("*"))
#     end
#     # TODO: checkbandmatch
#
#     Al, Au = _bandwidths(tA, A)
#     Bl, Bu = _bandwidths(tB, B)
#
#     if (-Al > Au) || (-Bl > Bu)   # A or B has empty bands
#         fill!(C, zero(T))
#     elseif Al < 0
#         C[max(1,Bn+Al-1):Am, :] .= zero(T)
#         banded_matmatmul!(C, tA, tB, _view(tA, A, :, 1-Al:An), _view(tB, B, 1-Al:An, :))
#     elseif Au < 0
#         C[1:-Au,:] .= zero(T)
#         banded_matmatmul!(view(C, 1-Au:Am,:), tA, tB, _view(tA, A, 1-Au:Am,:), B)
#     elseif Bl < 0
#         C[:, 1:-Bl] .= zero(T)
#         banded_matmatmul!(view(C, :, 1-Bl:Bn), tA, tB, A, _view(tB, B, :, 1-Bl:Bn))
#     elseif Bu < 0
#         C[:, max(1,Am+Bu-1):Bn] .= zero(T)
#         banded_matmatmul!(C, tA, tB, _view(tA, A, :, 1-Bu:Bm), _view(tB, B, 1-Bu:Bm, :))
#     else
#         positively_banded_matmatmul!(C, tA, tB, A, B)
#     end
#     C
# end

###
# Special Fill Diagonal
####

function materialize!(M::MatMulMatAdd{<:DiagonalLayout{<:AbstractFillLayout},<:AbstractBandedLayout})
    M.C .= (M.α * getindex_value(M.A.diag)) .* M.B .+ M.β .* M.C
    M.C
end

function materialize!(M::MatMulMatAdd{<:AbstractBandedLayout,<:DiagonalLayout{<:AbstractFillLayout}})
    M.C .= (M.α * getindex_value(M.B.diag)) .* M.A .+ M.β .* M.C
    M.C
end


text:
"""Validate that an argument tuple is the right length and has elements of type `argtype`."""
function validate_tupled_argument(arg, argtype, argname)
    length(arg) == 3        || throw(ArgumentError("length($argname) must be 3."))
    all(isa.(arg, argtype)) || throw(ArgumentError("$argname=$arg must contain $argtype s."))
    all(arg .> 0)           || throw(ArgumentError("Elements of $argname=$arg must be > 0!"))
    return nothing
end

function validate_grid_size_and_length(sz, len, halo, x, y, z)
    validate_tupled_argument(sz, Integer, "size")
    validate_tupled_argument(halo, Integer, "halo")

    # Find domain endpoints or domain length, depending on user input:
    if !isnothing(len) # the user has specified a length!
        (!isnothing(x) || !isnothing(y) || !isnothing(z)) &&
            throw(ArgumentError("Cannot specify both length and x, y, z keyword arguments."))

        validate_tupled_argument(len, Number, "length")

        Lx, Ly, Lz = len

        # An "oceanic" default domain
        x = (0, Lx)
        y = (0, Ly)
        z = (-Lz, 0)

    else # isnothing(length) === true implies that user has not specified a length
        (isnothing(x) || isnothing(y) || isnothing(z)) &&
            throw(ArgumentError("Must supply length or x, y, z keyword arguments."))

        function coord2xyz(c)
            c == 1 && return "x"
            c == 2 && return "y"
            c == 3 && return "z"
        end

        for (i, c) in enumerate((x, y, z))
            name = coord2xyz(i)
            length(c) == 2       || throw(ArgumentError("$name length($c) must be 2."))
            all(isa.(c, Number)) || throw(ArgumentError("$name=$c should contain numbers."))
            c[2] >= c[1]         || throw(ArgumentError("$name=$c should be an increasing interval."))
        end

        Lx = x[2] - x[1]
        Ly = y[2] - y[1]
        Lz = z[2] - z[1]
    end
    return Lx, Ly, Lz, x, y, z
end


text:
module PseudoGenomes

# package code goes here

export substitute_snps_in_fasta_file
export check_variant_alleles_in_fasta_file
export read_snp_positions_from_columns
export read_vcf_columns
export make_cg_bedfile

include(Pkg.dir("PseudoGenomes","src","pseudogenome-snp.jl"))

export mask_fasta_file_with_strain_SV_deletions
export get_strain_SVs

include(Pkg.dir("PseudoGenomes","src","pseudogenome-deletions.jl"))

export make_diploid
include(Pkg.dir("PseudoGenomes","src","pseudogenome-misc.jl"))
end # module


text:
module FED

using Printf: @printf,@sprintf
using SparseArrays: SparseMatrixCSC,spzeros
using QuantumLattices: FockOperator,Operators,FOperator,MAJORANA,CREATION,rank
using QuantumLattices.Prerequisites: Float
using QuantumLattices.Mathematics.VectorSpaces: VectorSpace

import QuantumLattices: kind,id,dimension,matrix
import QuantumLattices.Mathematics.AlgebraOverFields: idtype

export FockBasis,GBasis

"""
    FockBasis{I<:Unsigned} <: VectorSpace{I}

Abatract type for all kinds of Fock basis.
"""
abstract type FockBasis{I<:Unsigned} <: VectorSpace{I} end
kind(bs::FockBasis)=bs|>typeof|>kind
idtype(bs::FockBasis)=bs|>typeof|>idtype
function Base.show(io::IO,bs::FockBasis)
    @printf io "%s:\n" repr(bs)
    for i=1:dimension(bs)
        @printf io "  %s\n" string(bs[i],base=2)
    end
end
Base.repr(bs::FockBasis)=@sprintf "%s(%s)" nameof(typeof(bs)) join(id(bs),",")

"""
    GBasis{I}(nstate::Int) where I<:Unsigned

Generic fermionic/hard-core-bosonic basis, which uses no conserved quantities.
"""
struct GBasis{I<:Unsigned} <: FockBasis{I}
    nstate::Int
end
dimension(bs::GBasis)=2^bs.nstate
Base.getindex(bs::GBasis,i::Int)=i-1
Base.searchsortedfirst(bs::GBasis,i::Int)=i+1
kind(::Type{<:GBasis})=:g
idtype(::Type{<:GBasis})=Int
id(bs::GBasis)=bs.nstate

"""
    matrix(operator::FockOperator,basis::FockBasis) -> SparseMatrixCSC{valtype(operator),Int}
    matrix(opts::Operators,basis::FockBasis,dtype::Type{<:Number}=Complex{Float}) -> SparseMatrixCSC{dtype,Int}

Get the CSC-formed sparse matrix representation of an operator / a set of operators.
"""
function matrix(operator::FockOperator,basis::FockBasis)
    @assert rank(operator)%2==0 "matrix error: two bases are needed for odd ranked Fock operators."
    ndata,data,indices,indptr=1,zeros(valtype(operator),dimension(basis)),zeros(Int,dimension(basis)),zeros(Int,dimension(basis)+1)
    eye,temp,flag=one(basis|>eltype),zeros(basis|>eltype,rank(operator)+1),true
    seqs=NTuple{rank(operator),Int}(seq-1 for seq in operator.id.seqs)
    nambus=NTuple{rank(operator),Bool}((oid.index.nambu==MAJORANA ? error("matrix error: majorana operator not supported.") : oid.index.nambu==CREATION) for oid in operator.id)
    for i=1:dimension(basis)
        flag=true
        indptr[i]=ndata
        temp[1]=basis[i]
        for j=1:rank(operator)
            (temp[j]&eye<<seqs[j]>0)==nambus[j] && (flag=false; break)
            temp[j+1]=nambus[j] ? temp[j]|eye<<seqs[j] : temp[j]&~(eye<<seqs[j])
        end
        if flag
            nsign=0
            isa(operator,FOperator) && for j=1:rank(operator) for k=0:seqs[j]-1 (temp[j]&eye<<k>0 && (nsign+=1)) end end
            indices[ndata]=searchsortedfirst(basis,temp[end])
            data[ndata]=(-1)^nsign*operator.value
            ndata+=1
        end
    end
    indptr[end]=ndata+1
    return SparseMatrixCSC(dimension(basis),dimension(basis),indptr,indices[1:ndata],data[1:ndata])
end
function matrix(opts::Operators,basis::FockBasis,dtype::Type{<:Number}=Complex{Float})
    result=spzeros(dtype,dimension(basis),dimension(basis))
    for opt in values(opts) result=result+matrix(opt,basis) end
    return result
end

end  # module


text:
struct PiezoSystem
    amp::SimpleAmplifier
    piezo::Piezo
    pid::PidController
end

amplifier(ps::PiezoSystem) = ps.amp
piezo(ps::PiezoSystem) = ps.piezo
pid_controller(ps::PiezoSystem) = ps.pid
validate_command(ps::PiezoSystem, command) = validate_command(amp(ps), command)

PIEZOSYSTEMS = Dict{String, Any}()

function retrieve_piezosystem(rig_name::String)
    if !haskey(PIEZOSYSTEMS, rig_name)
        error("Unknown rig")
    end
    return PIEZOSYSTEMS[rig_name]
end

function register_piezosystem(ps::PiezoSystem, rig_name::String)
    if haskey(PIEZOSYSTEMS, rig_name)
        error("The rig $rig_name already has a registered PiezoSystem.")
    end
    PIEZOSYSTEMS[rig_name] = ps
end


text:
using CUSPARSE
using CUDArt
using Base.Test

m = 25
n = 35
k = 10
blockdim = 5

#############
# test_gtsv #
#############

function test_gtsv!(elty)
    dl = rand(elty,m-1)
    du = rand(elty,m-1)
    d = rand(elty,m)
    B = rand(elty,m,n)
    d_dl = CudaArray(vcat([0],dl))
    d_du = CudaArray(vcat(du,[0]))
    d_d = CudaArray(d)
    d_B = CudaArray(B)
    d_B = CUSPARSE.gtsv!(d_dl,d_d,d_du,d_B)
    C = diagm(d,0) + diagm(du,1) + diagm(dl,-1)
    h_B = to_host(d_B)
    @test_approx_eq(h_B, C\B)
end

function test_gtsv(elty)
    dl = rand(elty,m-1)
    du = rand(elty,m-1)
    d = rand(elty,m)
    B = rand(elty,m,n)
    d_dl = CudaArray(vcat([0],dl))
    d_du = CudaArray(vcat(du,[0]))
    d_d = CudaArray(d)
    d_B = CudaArray(B)
    d_C = CUSPARSE.gtsv(d_dl,d_d,d_du,d_B)
    C = diagm(d,0) + diagm(du,1) + diagm(dl,-1)
    h_C = to_host(d_C)
    @test_approx_eq(h_C, C\B)
end

#####################
# test_gtsv_nopivot #
#####################

function test_gtsv_nopivot!(elty)
    dl = rand(elty,m-1)
    du = rand(elty,m-1)
    d = rand(elty,m)
    B = rand(elty,m,n)
    d_dl = CudaArray(vcat([0],dl))
    d_du = CudaArray(vcat(du,[0]))
    d_d = CudaArray(d)
    d_B = CudaArray(B)
    d_B = CUSPARSE.gtsv_nopivot!(d_dl,d_d,d_du,d_B)
    C = diagm(d,0) + diagm(du,1) + diagm(dl,-1)
    h_B = to_host(d_B)
    @test_approx_eq(h_B, C\B)
end

function test_gtsv_nopivot(elty)
    dl = rand(elty,m-1)
    du = rand(elty,m-1)
    d = rand(elty,m)
    B = rand(elty,m,n)
    d_dl = CudaArray(vcat([0],dl))
    d_du = CudaArray(vcat(du,[0]))
    d_d = CudaArray(d)
    d_B = CudaArray(B)
    d_C = CUSPARSE.gtsv_nopivot(d_dl,d_d,d_du,d_B)
    C = diagm(d,0) + diagm(du,1) + diagm(dl,-1)
    h_C = to_host(d_C)
    @test_approx_eq(h_C, C\B)
end

#########################
# test_gtsvStridedBatch #
#########################

function test_gtsvStridedBatch!(elty)
    dla = rand(elty,m-1)
    dua = rand(elty,m-1)
    da = rand(elty,m)
    dlb = rand(elty,m-1)
    dub = rand(elty,m-1)
    db = rand(elty,m)
    xa = rand(elty,m)
    xb = rand(elty,m)
    d_dl = CudaArray(vcat([0],dla,[0],dlb))
    d_du = CudaArray(vcat(dua,[0],dub,[0]))
    d_d = CudaArray(vcat(da,db))
    d_x = CudaArray(vcat(xa,xb))
    d_x = CUSPARSE.gtsvStridedBatch!(d_dl,d_d,d_du,d_x,2,m)
    Ca = diagm(da,0) + diagm(dua,1) + diagm(dla,-1)
    Cb = diagm(db,0) + diagm(dub,1) + diagm(dlb,-1)
    h_x = to_host(d_x)
    @test_approx_eq(h_x[1:m], Ca\xa)
    @test_approx_eq(h_x[m+1:2*m], Cb\xb)
end

function test_gtsvStridedBatch(elty)
    dla = rand(elty,m-1)
    dua = rand(elty,m-1)
    da = rand(elty,m)
    dlb = rand(elty,m-1)
    dub = rand(elty,m-1)
    db = rand(elty,m)
    xa = rand(elty,m)
    xb = rand(elty,m)
    d_dl = CudaArray(vcat([0],dla,[0],dlb))
    d_du = CudaArray(vcat(dua,[0],dub,[0]))
    d_d = CudaArray(vcat(da,db))
    d_x = CudaArray(vcat(xa,xb))
    d_y = CUSPARSE.gtsvStridedBatch(d_dl,d_d,d_du,d_x,2,m)
    Ca = diagm(da,0) + diagm(dua,1) + diagm(dla,-1)
    Cb = diagm(db,0) + diagm(dub,1) + diagm(dlb,-1)
    h_y = to_host(d_y)
    @test_approx_eq(h_y[1:m], Ca\xa)
    @test_approx_eq(h_y[m+1:2*m], Cb\xb)
end

types = [Float32,Float64,Complex64,Complex128]
for elty in types
    tic()
    test_gtsv!(elty)
    println("gtsv! took ", toq(), " for ", elty)
    tic()
    test_gtsv(elty)
    println("gtsv took ", toq(), " for ", elty)
    tic()
    test_gtsv_nopivot!(elty)
    println("gtsv_nopivot! took ", toq(), " for ", elty)
    tic()
    test_gtsv_nopivot(elty)
    println("gtsv_nopivot took ", toq(), " for ", elty)
    tic()
    test_gtsvStridedBatch!(elty)
    println("gtsvStridedBatch! took ", toq(), " for ", elty)
    tic()
    test_gtsvStridedBatch(elty)
    println("gtsvStridedBatch took ", toq(), " for ", elty)
end


text:
"sparse_diag_mat_mul" => ("SparseMatMul", "libSparseMatMul", "sparse_sparse_mat_mul", false),
"sparse_assembler" => ("SparseAccumulate", "libSparseAccumulator", "", false),
"sparse_mat_mul" => ("SparseMatMul", "libSparseMatMul", "sparse_sparse_mat_mul", false),
"sparse_least_square" => ("SparseLeastSquare", "libSparseLeastSquare", "sparse_least_square", true),
"sparse_indexing" => ("SparseIndexing", "libSparseIndexing", "sparse_indexing", true),
"diag_sparse_mat_mul" => ("SparseMatMul", "libSparseMatMul", "sparse_sparse_mat_mul", false),
"sparse_solver" => ("SparseSolver", "libSparseSolver", "sparse_solver", true),
"sinkhorn_knopp" => ("OT/SinkhornKnopp", "libSinkhornKnopp", "sinkhorn_knopp", true),
"sparse_concate" => ("SparseConcate", "libSparseConcate", "sparse_concate", true),
"sparse_scatter_update" => ("SparseScatterUpdate", "libSparseScatterUpdate", "sparse_scatter_update", true),
"sparse_factorization" => ("SparseFactorizationSolve", "libfactorization", "sparse_factorization", false),
"sparse_solve" => ("SparseFactorizationSolve", "libSolve", "solve", true)


text:
export SolverCorrRelu

struct SolverCorrRelu{T<:Real} <: AbstractNextBasisSolver
    maxiter::Integer
    xtol::T
    ftol::T
end

function SolverCorrRelu(;
    maxiter::Integer=100,
    ftol::T=1e-5,
    xtol::T=1e-5) where T
    SolverCorrRelu{T}(maxiter, xtol, ftol)
end

"""
Iteratively minimize a linearization. Greedy or Frank-Wolfe
"""
function (self::SolverCorrRelu{T})(prob::Problem{L, R, CS},
                   g::Vector{T},
                   uinit::Vector{T}=zeros(sizep(prob))
                   ) where {T<:Real, L, R<:Relu{T}, CS<:L2Constraint}
    X = prob.X
    u = zeros(uinit)
    obj = calc(R, X*u)'*g                                             #X
    for iter in 1:self.maxiter
        active = find(X*u.≥0)
        # printfmtln("{:2d}) {:+6.4f} {} {:4.2f} 𝒜{}/{}",
                    # iter, obj, round.(Int, 100*u), normcs(prob, u), length(active), length(g))
        ui = -X[active, :]'*g[active]
        ui = scaletonorm(prob.Ω, ui)
        unext = (iter*u + 2*ui)/(iter+2)
        objnext = calc(R, X*unext)'*g                                  #X
        if norm(unext-u) ≤ self.xtol || abs(objnext-obj) ≤ self.ftol
            break
        end
        u = unext
        obj = objnext
    end
    u, obj
end


text:
function JuliogToVerilog(ex::Expr, ad::Dict{Symbol, Arrow})
    global arrowDict = ad
    return JVmodule(ex)
end

function JVmodule(ex::Expr)
    println(ex)
    # Convert function header to a module header
    verilog = "module"
    l = length(ex.args[1].args)
    if l == 1
        warn("Conversion to Verilog: Module creation with no input/output signals")
    end
    verilog = verilog * " " * String(ex.args[1].args[1]) * "("
    for i = 2:(l-1)
        verilog = verilog * String(ex.args[1].args[i]) * ", "
    end
    verilog = verilog * String(ex.args[1].args[l]) * ") ;\n\n"

    # Convert function block into verilog syntax
    vblock = ex.args[2]
    verilog = verilog * JVblock(vblock)

    # End the module
    verilog = verilog * "\nendmodule"

    # return the completed verilog syntax String
    return verilog
end

# Convert a block of Juliog syntax to Verilog syntax
# Returns a string
# The module will have no parameters, as those will have all been removed during parameterize
function JVblock(ex::Expr)
    verilog = ""
    l = length(ex.args)

    for i = 1:l
        outstr = JVblockhelper(ex.args[i])
        verilog = verilog * outstr * "\n"
    end
    return verilog
end

function JVblockhelper(ex::Expr)
    h = ex.head
    println(ex)
    if     h == :(=)
        return JVequals(ex)
    elseif h == :(:=)
        return JVassignment(ex)
    elseif h == :macrocall
        return JVmacro(ex)
    else
        error("Hit an unexplored expression head $(h) in expr:\n$(ex)")
    end
end

# Must be creating a wire
# MOON Accomodate matrix notation here
function JVequals(ex::Expr)
    global arrowDict

    # Get the name of the new wire or reg
    if isa(ex.args[1], Expr)
        # Is a ref
        name = ex.args[1].args[1]
    else
        name = ex.args[1]
    end

    # Get the arrow info
    wt = arrowDict[name].WireType
    reg = arrowDict[name].Reg

    # Get the bit indices and string
    bl = arrowDict[name].BitLeft
    br = arrowDict[name].BitRight
    index = "[" * dec(bl) * ":" * dec(br) * "]"

    if     wt == 0 && !reg
        return "wire " * index * " " * String(name) * " ;"
    elseif wt == 0 && reg
        return "reg " * index * " " * String(name) * " ;"
    elseif wt == 1
        return "input " * index * " " * String(name) * " ;"
    elseif wt == 2 && !reg
        return "output " * index * " " * String(name) * " ;"
    elseif wt == 0 && reg
        return "output reg " * index * " " * String(name) * " ;"
    end
end

function JVassignment(ex::Expr)
    verilog = "assign "
    if isa(ex.args[1], Expr) && ex.args[1].head == :ref
        verilog = verilog * String(ex.args[1].args[1])
        # MOON account for Matrix notation here
        if isa(ex.args[1].args[2], Expr)
            # Multiple Bits refererence
            if length(ex.args[1].args[2].args) == 3
                verilog = verilog * "[" * dec(ex.args[1].args[2].args[1]) * ":" * dec(ex.args[1].args[2].args[3]) * "] = "
            else
                verilog = verilog * "[" * dec(ex.args[1].args[2].args[1]) * ":" * dec(ex.args[1].args[2].args[2]) * "] = "
            end
        else
            # Bit reference
            verilog = verilog * "[" * dec(ex.args[1].args[2]) * "] = "
        end
    else
        # Is just a simple symbol then
        verilog =  verilog * String(ex.args[1]) * " = "
    end

    rhs = ex.args[2]
    if isa(rhs, Expr)
        verilog = verilog * JVrhs(rhs)
    else
        # Is a simple symbol or int
        if isa(rhs, Int)
            verilog = verilog * dec(rhs)
        elseif isa(rhs, Char)
            verilog = verilog * string(rhs)
        else
            verilog = verilog * String(rhs)
        end
    end

    return verilog * " ;"
end

# MOON Accomodate matrix notation
# An expression needing to be converted to verilog syntax
function JVrhs(ex::Expr)
    h = ex.head
    if     h == :vcat
        return JVrhsvcat(ex)
    elseif h == :ref
        return JVrhsref(ex)
    elseif h == :call
        return JVrhscall(ex)
    elseif h == :if
        return JVrhsif(ex)
    else
        error("inside JVrhs, hit unexplored expression head $(h) in expr:\n$(ex)")
    end
end

function JVrhsvcat(ex::Expr)
    println(ex)
    verilog = "{"
    l = length(ex.args)
    bitcount = 0
    bitstr = ""
    for i = 1:l
        arg = ex.args[i]

        if     isa(arg, Int)
            # accumulate bit string until something other than int
            bitcount = bitcount + 1
            bitstr = bitstr * dec(arg)
        else
            if verilog != "{"
                verilog = verilog * ", "
            end

            if bitcount != 0
                verilog = verilog * dec(bitcount) * "b" * bitstr
                bitcount = 0
                bitstr = ""
            end

            if     isa(arg, Expr)
                # Only refs and bit repeats are allowed
                if     arg.head == :ref
                    verilog = verilog * JVrhsref(arg)
                elseif arg.head == :call
                    verilog = verilog * JVrhscall(arg)
                else
                    error("Inside JVrhsvcat, argument $(arg) is unaccepted syntax:\n$(ex)")
                end
            elseif isa(arg, Char)
                verilog = verilog * string(arg)
            else
                # Is a symbol
                verilog = verilog * String(arg)
            end
        end
        println(bitcount)
        println(bitstr)
        println(verilog)
    end
    # In case the concatenation ends with bit padding
    if bitcount != 0
        verilog = verilog * dec(bitcount) * "b" * bitstr
    end
    # Add the closing bracket and return
    return verilog * "}"
end

# MOON Accomodate matrix notation here
function JVrhsref(ex::Expr)
    verilog = String(ex.args[1])
    if isa(ex.args[2], Expr)
        # Multiple Bits refererence
        if length(ex.args[2].args) == 3
            verilog = verilog * "[" * dec(ex.args[2].args[1]) * ":" * dec(ex.args[2].args[3]) * "]"
        else
            verilog = verilog * "[" * dec(ex.args[2].args[1]) * ":" * dec(ex.args[2].args[2]) * "]"
        end
    else
        # Bit reference
        verilog = verilog * "[" * dec(ex.args[2]) * "]"
    end
    return verilog
end

function JVrhscall(ex::Expr)
    l = length(ex.args)

    # Check to see if this call is just a bitrepeat
    if l == 2 && Symbol(ex.args[1]) != :~
        # Must be a bit repeat
        verilog = "(" * dec(ex.args[2]) * "){"
        if isa(ex.args[1], Expr)
            verilog = verilog * JVrhs(ex.args[1]) * "}"
            return verilog
        else
            verilog = verilog * String(Symbol(ex.args[1])) * "}"
            return verilog
        end
    end

    argstrs = Array{String}(l)
    for i = 1:l
        arg = ex.args[i]
        if     isa(arg, Expr) && arg.head == :call
            argstrs[i] = "(" * JVrhscall(arg) * ")"
        elseif isa(arg, Expr)
            argstrs[i] = JVrhs(arg)
        elseif isa(arg, Int)
            argstrs[i] = dec(arg)
        elseif isa(arg, Char)
            argstrs[i] = string(arg)
        else
            argstrs[i] = String(Symbol(arg))
        end
    end

    if     l == 3
        return argstrs[2] * " " * argstrs[1] * " " * argstrs[3]
    elseif l == 2
        return argstrs[1] * argstrs[2]
    else
        error("In JVrhscall, Hit a :call head which has $(l) arguments. Throwing a fit:\n$(ex)")
    end
end

function JVrhsif(ex::Expr)
    # reduce condition to a single bit if acceptable
    condition = ex.args[1]
    if     isa(condition, Expr)
        verilog = "(" * JVrhsref(condition) * ") ? "
    elseif isa(condition, Symbol)
        verilog = "(" * String(condition) * ") ? "
    else
        # Must be an array
        verilog = "(" * condition[1]
        if isa(condition[2], Expr)
            verilog = verilog * JVrhsref(condition[2]) * ") ? "
        else
            verilog = verilog * String(condition[2]) * ") ? "
        end
    end

    # Must have 3 args after sanitization step
    arg = ex.args[1]
    if     isa(arg, Expr)
        argstr = "(" * JVrhs(arg) * ")"
    elseif isa(arg, Int)
        argstr = dec(arg)
    elseif isa(arg, Char)
        argstr = string(arg)
    else
        argstr = String(arg)
    end
    verilog = verilog * argstr * " : "

    arg = ex.args[2]
    if     isa(arg, Expr)
        argstr = "(" * JVrhs(arg) * ")"
    elseif isa(arg, Int)
        argstr = dec(arg)
    elseif isa(arg, Char)
        argstr = string(arg)
    else
        argstr = String(arg)
    end
    verilog = verilog * argstr

    return verilog
end

# TODO registers
function JVmacro(ex::Expr)
    mc = ex.args[1]
    if     mc == Symbol("@async")
        return JVasync(ex)
    elseif mc == Symbol("@reg")
        error("Hit an unfinished macrocall: $(mc)")
    elseif mc == Symbol("@posedge")
        error("Hit an unfinished macrocall: $(mc)")
    elseif mc == Symbol("@negedge")
        error("Hit an unfinished macrocall: $(mc)")
    elseif mc == Symbol("@delay")
        error("Hit an unfinished macrocall: $(mc)")
    elseif mc == Symbol("@block")
        return JVmacroblock(ex)
    elseif mc == Symbol("@verilog")
        return JVmacroverilog(ex)
    elseif mc == Symbol("@Julia")
        return JVmacrojulia(ex)
    else
        error("Hit an unexplored macrocall $(mc) in expr:\n$(ex)")
    end
end

function JVasync(ex::Expr)
    verilog = "always @* begin\n"

    l = length(ex.args[2].args)
    for i = 1:l
        verilog = verilog * JVasynchelper(ex.args[2].args[i], 1)
    end

    verilog = verilog * "end"
    return verilog
end

function JVasynchelper(ex::Expr, depth::Int)
    if isa(ex.args[1], Expr) && ex.args[1].head == :if
        verilog = ""
        # TODO if statements inside async
    else
        if ex.args[1].head == :ref
            verilog = String(ex.args[1].args[1])
            # MOON account for Matrix notation here
            if isa(ex.args[1].args[2], Expr)
                # Multiple Bits refererence
                if length(ex.args[1].args[2].args) == 3
                    verilog = verilog * "[" * dec(ex.args[1].args[2].args[1]) * ":" * dec(ex.args[1].args[2].args[3]) * "] = "
                else
                    verilog = verilog * "[" * dec(ex.args[1].args[2].args[1]) * ":" * dec(ex.args[1].args[2].args[2]) * "] = "
                end
            else
                # Bit reference
                verilog = verilog * "[" * dec(ex.args[1].args[2]) * "] = "
            end
        else
            # Is just a simple symbol then
            verilog =  ex.args[1] * " = "
        end

        rhs = ex.args[2]
        if isa(rhs, Expr)
            verilog = verilog * JVrhs(rhs)
        else
            # Is a simple symbol or int
            if isa(rhs, Int)
                verilog = verilog * dec(rhs)
            elseif isa(rhs, Char)
                verilog = verilog * string(rhs)
            else
                verilog = verilog * String(rhs)
            end
        end

        for i = 1:depth
            verilog = "\t" * verilog
        end

        return verilog * " ;\n"
    end
end

# TODO Verilog doesn't natively support recursion, so need to remove it
function JVmacroblock(ex::Expr)

    error("JVmacroblock is unfinished")
end

function JVmacrojulia(ex::Expr)
    verilog = "/* Pass Through Verilog */\n"
    verilog = verilog * ex.args[2]
    return verilog
end

function JVmacrojulia(ex::Expr)
    verilog = "/* Pass Through Julia\n"
    verilog = verilog * ex.args[2]
    verilog = verilog * "\n*/\n"
    return verilog
end

max_stars_repo_name: trevortheblack/Juliog

max_issues_repo_name: trevortheblack/Juliog

max_forks_repo_name: trevortheblack/Juliog


text:
export matrixtoMaple,
    xstoMaple,
    kstoMaple,
    systemtoMaple,
    WsystemtoMaple,
    toMaple

function matrixtoMaple(io, M, name)
    write(io, "\n\n$(name) := Matrix$(size(M)):\n\n")
    for i in CartesianIndices(M)
        if M[i] != 0
            write(io, "$(name)[$(i[1]),$(i[2])] := $(M[i]):\n")
        end
    end
end

macro matrixtoMaple(io, M)
    return quote
        matrixtoMaple($(esc(io)), $(esc(M)), $(esc((string(M)))))
    end
end

function xstoMaple(io, stoichiometricsources)
    xs = (1:(size(stoichiometricsources, 1)))[nonzeroslicesof(stoichiometricsources)(1)]
    write(io, "\n\nnxs := $(size(xs, 1)):\n")
    write(io, "xs := [seq(x[i], i = [")
    for x in xs[1:(end - 1)]
        write(io, "$x, ")
    end
    write(io, "$(xs[end])])];\n\n")
    write(io, "depenvars := [seq(x[i], i = [")
    for x in xs[2:(end - 2)]
        write(io, "$x, ")
    end
    write(io, "$(xs[end - 1])])];\n\n")
end

function kstoMaple(io, stoichiometricsources)
    ks = (1:(size(stoichiometricsources, 2)))[nonzeroslicesof(stoichiometricsources)(2)]
    write(io, "\n\nnks := $(size(ks, 1)):\n")
    write(io, "ks := [seq(k[i], i = [")
    for k in ks[1:(end - 1)]
        write(io, "$k, ")
    end
    write(io, "$(ks[end])])];\n\n")
end

function systemtoMaple(io)
    write(io, "v := Velocities(Y,xs):\n")
    write(io, "digK := Matrix(nks):\n")
    write(io, "for i to nks do digK[i,i] := ks[i] end do:\n")
    write(io, "S := N.digK.v:\n")
    write(io, "Seq:= equfy(S);\n")
end

function WsystemtoMaple(io, nts, W)
    write(io, "\n\nSw := copy(S):\n")
    write(io, "Wx := (W.(Vector[column](xs))) - Vector[column]([seq(T[i], i = 1 .. $(nts))]):\n")
    for (i, p) in enumerate(findpivotsof(W)(1))
        write(io, "Sw[$p] := Wx[$i]:\n")
    end
    write(io, "Sweq:= equfy(Sw):\n")
    write(io, "J := VectorCalculus[Jacobian](Sw, xs):\n")
    write(io, "DJ := (-1)^(Rank(N))*Determinant(J):\n")
end

## It needs Y and E be defined in Maple
function convexparamtoMaple(io)
    write(io, "\n\ndigL := Matrix(convert(E.(Vector[column]([seq(lambda[i], i=1..LinearAlgebra[ColumnDimension](E))])), Vector[row]), shape = diagonal):\n")
    write(io, "digH := DiagonalMatrix([seq(h[i], i = 1..LinearAlgebra[ColumnDimension](LinearAlgebra[Transpose](Y)))]):\n")
    write(io, "Jconv := N.digL.LinearAlgebra[Transpose](Y).digH:\n")
end

function toMaple(net, nxs, file::String)
    S = stoichiometriccoeffs(net, nxs) ## S may have zero rows and cols
    N = stoichiometricmatrix(S) ## Has no zero row or col
    nts, W = conservativelaws(N)
    Y = kineticorder(S)
    E = raysof(cone_positivenullspace(N))
    open(file, "w") do io
        write(io, "read(\"ModelsMatricies/Procs.mpl\"):\n")
        @matrixtoMaple io Y
        @matrixtoMaple io N
        @matrixtoMaple io W
        @matrixtoMaple io E
        xstoMaple(io, stoichiometricsources(S))
        kstoMaple(io, stoichiometricsources(S))
        systemtoMaple(io)
        WsystemtoMaple(io, nts, W)
        convexparamtoMaple(io)
    end
end


text:
using NumericalMethodsforEngineers, Test

a = [1. 2. -2.; 2. 5. -4.; -2. -4. 5.]

println("\n|a[1, 1]| = $(det(a[1, 1]))\n")
println("|a[1:2, 1:2]| = $(det(a[1:2, 1:2])))\n")
println("|a| = $(det(a))")
println("\ncholesky(a):")
cholesky(a) |> display

a = [16. 4. 8.; 4. 5. -4.; 8. -4. 22.]

println("\n|a[1, 1]| = $(det(a[1, 1]))\n")
println("|a[1:2, 1:2]| = $(det(a[1:2, 1:2])))\n")
println("|a| = $(det(a))")

println("\ncholesky(a) - upper:")
upper = cholesky(a).U
upper |> display

println("\ncholesky(a) - lower:")
lower = cholesky(a).L
lower |> display

println("\nlower * upper:")
lower * upper |> display

@test a == lower * upper


text:
module NCP_Exact

include("../Misc.jl")
include("NCP_MC.jl")
include("NCP_VB.jl")

using .Misc
import .NCP_MC, .NCP_VB

using Distributions, SpecialFunctions, Einsum, Combinatorics
import Base.Iterators: product

export allocations, full, generate
export log_marginal, alloc_dist, EP_dist
export log_posterior_T, log_PT

allocations(X::Array{ℜ,N}, H::Vararg{Ƶ}) where {ℜ<:Real, Ƶ<:Integer, N} = Channel(ctype=Array{ℜ,N+length(H)}) do c
    I, R = size(X), prod(H)
    for partition = product(map(Xᵢ -> combinations(1:(Xᵢ+R-1),R-1), X)...)
        push!(c,reshape([sᵢ for (i,pᵢ) ∈ enumerate(partition)
                        for sᵢ = diff([0, pᵢ..., X[i]+R]) .- 1],H...,I...))
    end
end

full(X::Array{ℜ,N}, T::ℜ) where {ℜ<:Real, N} = Channel(ctype=Array{ℜ,N}) do c
    miss = findall(isnan,X)
    M = length(miss)
    Tₘ = Int.(T - nansum(X))
    for partition ∈ combinations(1:(Tₘ+M-1),M-1)
        X_full = copy(X)
        X_full[miss] .= diff([0, partition..., Tₘ+M]) .- 1
        push!(c,X_full)
    end
end

function generate(T::Ƶ, R::Ƶ, I::Vararg{Ƶ}; a::ℜ=1.0) where {ℜ<:Real, Ƶ<:Integer}
    α_I = (a/R) ./ I
    α_R = a/R

    θ_R = rand(Dirichlet(R,α_R))
    θ_I = [rand(Dirichlet(Iₙ,α_Iₙ),R) for (Iₙ,α_Iₙ) in zip(I,α_I)]

    D_R = rand(Categorical(θ_R),T)
    D_I = [rand(Categorical(θ_Iₙ[:,r])) for r ∈ D_R, θ_Iₙ ∈ θ_I]

    return cooccurrences(D_I,I), (θ_R, θ_I...)
end

function log_marginal(X::Array{ℜ,N}, R::Ƶ; a=1.0, b=a/nansum(X)) where {ℜ<:Real, Ƶ<:Integer, N}
    I = size(X)
    α_I = (a/R) ./ I
    α_R = a/R

    S_R, S₊ = Array{ℜ}(undef,R), sum(X)
    S_I = [Array{ℜ}(undef,R,Iₙ) for Iₙ in I]

    log_PX = -Inf
    log_C = a*log(b) - (a+S₊)*log(b + 1) + lgamma(a + S₊) - lgamma(a)

    for S ∈ allocations(X,R)
        log_PS = log_C - sum(lgamma, S .+ 1)
        for (n, Iₙ) ∈ enumerate(I)
            S_I[n] .= sum(S,(2:n...,n+2:N+1...))
            log_PS += sum(lbeta(α_I[n] .+ S_I[n];dims=2)) - R*lbeta(α_I[n],Iₙ)
        end
        S_R .= sum(S,(2:N+1...,))
        log_PS += lbeta(α_R .+ S_R) - lbeta(α_R,R)
        log_PX = logsumexp([log_PX,log_PS])
    end
    return log_PX
end

function alloc_dist(X::Array{ℜ,N}, R::Ƶ; a=1.0, b=a/nansum(X)) where {ℜ<:Real, Ƶ<:Integer, N}
    I = size(X)
    α_I = (a/R) ./ I
    α_R = a/R

    S_R, S₊ = Array{ℜ}(undef,R), sum(X)
    S_I = [Array{ℜ}(undef,R,Iₙ) for Iₙ in I]

    nS = prod(x -> binomial(Ƶ(x+R-1),R-1), X)
    log_P = zeros(nS)

    log_C = a*log(b) - (a+S₊)*log(b + 1) + lgamma(a + S₊) - lgamma(a)

    for (s, S) ∈ enumerate(allocations(X,R))
        log_P[s] = log_C - sum(lgamma, S .+ 1)
        for (n, Iₙ) ∈ enumerate(I)
            S_I[n] .= sum(S,(2:n...,n+2:N+1...))
            log_P[s] += sum(lbeta(α_I[n] .+ S_I[n];dims=2)) - R*lbeta(α_I[n],Iₙ)
        end
        S_R .= sum(S,(2:N+1...,))
        log_P[s] += lbeta(α_R .+ S_R) - lbeta(α_R,R)
    end
    return log_P
end

function EP_dist(X::Array{ℜ,N}, R::Ƶ) where {ℜ<:Real, Ƶ<:Integer, N}
    I = size(X)

    S_R, S₊ = Array{ℜ}(undef,R), sum(X)
    S_I = [Array{ℜ}(undef,R,Iₙ) for Iₙ in I]

    nS = prod(x -> binomial(Ƶ(x+R-1),R-1), X)
    dₑₚ = zeros(nS)

    for (s, S) ∈ enumerate(allocations(X,R))
        for (n, Iₙ) ∈ enumerate(I)
            S_I[n] .= sum(S,(2:n...,n+2:N+1...))
            dₑₚ[s] += sum(!iszero, S_I[n])
        end
        S_R .= sum(S,(2:N+1...,))
        dₑₚ[s] += (N-1) * sum(!iszero, S_R) + !iszero(S₊)
    end
    return dₑₚ
end

function log_posterior_T(X::Array{ℜ,N}, T::ℜ , R::Ƶ; a=1.0, b=a/nansum(X),
    smc::Bool=false, elbo::Bool=false, particles::Int=1000, EPOCHS::Int=100, M::Int=100) where {ℜ<:Real, Ƶ<:Integer, N}
    if smc
        log_PX = [NCP_MC.smc_weight(X_full,R,particles; a=a, b=b)[1] for X_full ∈ full(X,T)]
        return logsumexp(log_PX)
    elseif elbo
        log_PX = nanmax([NCP_VB.standard_VB(X_full,R; a=a, b=b, EPOCHS=EPOCHS)[1][end] for m ∈ 1:M, X_full ∈ full(X,T)],1)
        return logsumexp(log_PX)
    else
        log_PX = [log_marginal(X_full,R; a=a, b=b) for X_full ∈ full(X,T)]
        return logsumexp(log_PX)
    end
end

end

max_stars_repo_name: mehmetburakkurutmaz/bam

max_issues_repo_name: mehmetburakkurutmaz/bam

max_forks_repo_name: atcemgil/bam


text:
@testset "half disk" begin
    ϕ = π/6  # start angle
    θ = ϕ-π/2  # outward normal to diameter
    c = 0.01 .* [cos(θ), sin(θ)]  # center of disk
    r = 2  # radius
    ∆ϕ = π  # half disk
    b = Sector(c, r, ϕ, ∆ϕ)  # note sector includes origin

    ∆vₛ = r * [cos(ϕ), sin(ϕ)]  # start of arc
    ∆vₑ = r * [cos(ϕ+∆ϕ), sin(ϕ+∆ϕ)]  # end of arc
    ∆vs = (∆vₛ, ∆vₑ)

    @test b == deepcopy(b)
    @test b ≈ Sector(c, r, ϕ+2π, ∆ϕ)
    @test b ≈ Sector(c, r, ϕ-2π, ∆ϕ)
    @test b ≈ Sector(c, r, ϕ+∆ϕ, -∆ϕ)
    @test b ≈ Sector(c, r, ϕ+∆ϕ-2π, -∆ϕ)
    @test hash(b) == hash(deepcopy(b))
    @test [0,0] ∈ b  # origin
    @test c ∈ b  # center
    @test (s = one⁻; s.*(∆vₛ+c) ∈ b && s.*(∆vₑ+c) ∈ b)  # slightly inside of end points of arc
    @test (s = one⁺; s.*(∆vₑ+c) ∉ b && s.*(∆vₑ+c) ∉ b)  # slightly outside of end points of arc

    @test all([(xn = surfpt_nearby(s.*c,b); xn ≈ (c, [cos(θ), sin(θ)])) for s = (one⁻⁻, 1, one⁺⁺)])  # handle point at center properly

    @test all([(p₀ = ∆vₛ/2 + c; p = s.*p₀; (x,n) = surfpt_nearby(p,b); atan(x[2]-c[2],x[1]-c[1]) ≈ ϕ && n ≈ [cos(θ), sin(θ)]) for s = (one⁻⁻,1,one⁺⁺)])  # around start side
    @test all([(p₀ = ∆vₑ/2 + c; p = s.*p₀; (x,n) = surfpt_nearby(p,b); atan(x[2]-c[2],x[1]-c[1]) ≈ ϕ+∆ϕ-2π && n ≈ [cos(θ), sin(θ)]) for s = (one⁻⁻,1,one⁺⁺)])  # around end side
    @test all([(α = ϕ+0.7∆ϕ; d̂ = [cos(α),sin(α)]; p = s.*r.*d̂ + c; surfpt_nearby(p,b) ≈ (r.*d̂.+c, d̂)) for s = (one⁻⁻,1,one⁺⁺)])  # around arc

    @test all([(p₀ = ∆v + c; p = one⁺⁺.*p₀; surfpt_nearby(p,b) ≈ (p₀, normalize(p-p₀))) for ∆v = ∆vs])  # outside end points of arc
    @test (p = ∆vₛ + c; surfpt_nearby(p,b) ≈ (∆vₛ+c, normalize([cos(ϕ),sin(ϕ)]+[cos(θ),sin(θ)])))   # at start point of arc
    @test (p = ∆vₑ + c; surfpt_nearby(p,b) ≈ (∆vₑ+c, normalize([cos(ϕ+∆ϕ),sin(ϕ+∆ϕ)]+[cos(θ),sin(θ)])))   # at end point of arc

    @test bounds(b) ≈ ([c[1]-r, c[2]+r*sin(ϕ+∆ϕ)], [c[1]+r*cos(ϕ), c[2]+r])
    @test checkbounds(b)

    @test (∆ = rand(2); translate(b,∆) ≈ Sector(c+∆, r, ϕ, ∆ϕ))
end  # @testset "half disk"

@testset "Sector, acute angle" begin
    ϕ = π/6  # start angle
    θ₁ = ϕ-π/2  # outward normal to start edge
    ∆ϕ = π/6  # half disk
    θ₂ = (ϕ+∆ϕ) + π/2  # outward normal to end edge

    ϕ₀ = ϕ + ∆ϕ/2
    c = -0.01 .* [cos(ϕ₀), sin(ϕ₀)]  # center
    r = 2  # radius
    b = Sector(c, r, ϕ, ∆ϕ)  # note sector includes origin

    ∆vₛ = r * [cos(ϕ), sin(ϕ)]  # start of arc
    ∆vₑ = r * [cos(ϕ+∆ϕ), sin(ϕ+∆ϕ)]  # end of arc
    ∆vs = (∆vₛ, ∆vₑ)

    @test b == deepcopy(b)
    @test b ≈ Sector(c, r, ϕ+2π, ∆ϕ)
    @test b ≈ Sector(c, r, ϕ-2π, ∆ϕ)
    @test b ≈ Sector(c, r, ϕ+∆ϕ, -∆ϕ)
    @test b ≈ Sector(c, r, ϕ+∆ϕ-2π, -∆ϕ)
    @test hash(b) == hash(deepcopy(b))
    @test [0,0] ∈ b  # origin
    @test c ∈ b  # center
    @test (s = one⁻; s.*(∆vₛ+c) ∈ b && s.*(∆vₑ+c) ∈ b)  # slightly outside of end points of arc
    @test (s = one⁺; s.*(∆vₑ+c) ∉ b && s.*(∆vₑ+c) ∉ b)  # slightly inside of end points of arc

    @test surfpt_nearby(c,b) ≈ (c, -[cos(ϕ₀), sin(ϕ₀)])  # handle point at center properly
    @test all([(p = c+0.1r.*[cos(θ),sin(θ)]; surfpt_nearby(p,b) ≈ (c, normalize(p-c))) for θ = (ϕ+π, ϕ+∆ϕ+π, ϕ₀+π, ϕ+π+0.1∆ϕ, ϕ+∆ϕ+π-0.1∆ϕ)])  # outside center

    @test all([(p₀ = ∆vₛ/2 + c; p = s.*p₀; (x,n) = surfpt_nearby(p,b); atan(x[2]-c[2],x[1]-c[1]) ≈ ϕ && n ≈ [cos(θ₁), sin(θ₁)]) for s = (one⁻⁻,1,one⁺⁺)])  # around start side
    @test all([(p₀ = ∆vₑ/2 + c; p = s.*p₀; (x,n) = surfpt_nearby(p,b); atan(x[2]-c[2],x[1]-c[1]) ≈ ϕ+∆ϕ && n ≈ [cos(θ₂), sin(θ₂)]) for s = (one⁻⁻,1,one⁺⁺)])  # around end side
    @test all([(α = ϕ+0.7∆ϕ; d̂ = [cos(α),sin(α)]; p = s.*r.*d̂ + c; surfpt_nearby(p,b) ≈ (r.*d̂.+c, d̂)) for s = (one⁻⁻,1,one⁺⁺)])  # around arc

    @test all([(p₀ = ∆v + c; p = one⁺⁺.*p₀; surfpt_nearby(p,b) ≈ (p₀, normalize(p-p₀))) for ∆v = ∆vs])  # outside end points of arc
    @test (p = ∆vₛ + c; surfpt_nearby(p,b) ≈ (∆vₛ+c, normalize([cos(ϕ),sin(ϕ)]+[cos(θ₁),sin(θ₁)])))   # at start point of arc
    @test (p = ∆vₑ + c; surfpt_nearby(p,b) ≈ (∆vₑ+c, normalize([cos(ϕ+∆ϕ),sin(ϕ+∆ϕ)]+[cos(θ₂),sin(θ₂)])))   # at end point of arc

    @test bounds(b) ≈ (c, [c[1]+r*cos(ϕ), c[2]+r*sin(ϕ+∆ϕ)])
    @test checkbounds(b)

    @test (∆ = rand(2); translate(b,∆) ≈ Sector(c+∆, r, ϕ, ∆ϕ))
end  # @testset "Sector, acute angle"

@testset "Half pipe" begin
    ϕ = π/6  # start angle
    θ = ϕ-π/2  # outward normal to diameter
    c = 0.01 .* [cos(θ), sin(θ)]  # center of disk
    r = 2  # radius
    ∆ϕ = π  # half disk
    h2 = 1.1

    b = SectoralPrism([c; 0], r, ϕ, ∆ϕ, 2h2)

    ∆vₛ = r * [cos(ϕ), sin(ϕ)]  # start of arc
    ∆vₑ = r * [cos(ϕ+∆ϕ), sin(ϕ+∆ϕ)]  # end of arc
    ∆vs = (∆vₛ, ∆vₑ)

    @test b == deepcopy(b)
    @test b ≈ SectoralPrism([c; 0], r, ϕ+2π, ∆ϕ, 2h2)
    @test b ≈ SectoralPrism([c; 0], r, ϕ-2π, ∆ϕ, 2h2)
    @test b ≈ SectoralPrism([c; 0], r, ϕ+∆ϕ, -∆ϕ, 2h2)
    @test b ≈ SectoralPrism([c; 0], r, ϕ+∆ϕ-2π, -∆ϕ, 2h2)
    @test hash(b) == hash(deepcopy(b))
    @test [0,0,0] ∈ b  # center
    @test [c; 0] ∈ b  # center
    @test all([(s = one⁻; [s.*(∆v+c);t*h2] ∈ b) for t = (-one⁻,0,one⁻), ∆v = ∆vs])  # slightly inside of end points of arc
    @test all([(s = one⁺; [s.*(∆v+c);t*h2] ∉ b) for t = (-one⁻,0,one⁻), ∆v = ∆vs])  # slightly outside of end points of arc

    @test all([(xn = surfpt_nearby([s.*c;t*h2/2], b); xn ≈ ([c;t*h2/2], [cos(θ),sin(θ),0])) for s = (one⁻⁻, 1, one⁺⁺), t = (-1,0,1)])  # handle point at center properly

    @test all([(p₀ = ∆vₛ/2 + c; p = [s.*p₀;t*h2]; (x,n) = surfpt_nearby(p,b); atan(x[2]-c[2],x[1]-c[1])≈ϕ && x[3]==t*h2 && n ≈ [cos(θ),sin(θ),0]) for s = (one⁻⁻,1,one⁺⁺), t = (-0.5,0,0.5)])  # around start side, within height
    @test all([(s = 1; p₀ = ∆vₛ/2 + c; p = [s.*p₀;t*h2]; (x,n) = surfpt_nearby(p,b); atan(x[2]-c[2],x[1]-c[1])≈ϕ && x[3]==t*h2 && n ≈ normalize([cos(θ),sin(θ),t])) for t = (-1,1)])  # at start side, at base
    @test all([(s = one⁻⁻; p₀ = ∆vₛ/2 + c; p = [s.*p₀;t*h2]; (x,n) = surfpt_nearby(p,b); x==p && n==[0,0,t]) for t = (-1,1)])  # inside start side, at base
    @test all([(s = one⁺⁺; p₀ = ∆vₛ/2 + c; p = [s.*p₀;t*h2]; (x,n) = surfpt_nearby(p,b); atan(x[2]-c[2],x[1]-c[1])≈ϕ && x[3]≈t*h2 && n ≈ normalize([cos(θ),sin(θ),0])) for t = (-1,1)])  # outside start side, at base
    @test all([(s = one⁺⁺; p₀ = ∆vₛ/2 + c; p = [s.*p₀;t*h2]; (x,n) = surfpt_nearby(p,b); atan(x[2]-c[2],x[1]-c[1])≈ϕ && x[3]≈sign(t)*h2 && n ≈ normalize(p-x)) for t = (-one⁺⁺,one⁺⁺)])  # outside start side, outside base

    @test all([(p₀ = ∆vₑ/2 + c; p = [s.*p₀;t*h2]; (x,n) = surfpt_nearby(p,b); atan(x[2]-c[2],x[1]-c[1])≈ϕ+∆ϕ-2π && x[3]==t*h2 && n ≈ [cos(θ),sin(θ),0]) for s = (one⁻⁻,1,one⁺⁺), t = (-0.5,0,0.5)])  # around end side, within height
    @test all([(s = 1; p₀ = ∆vₑ/2 + c; p = [s.*p₀;t*h2]; (x,n) = surfpt_nearby(p,b); atan(x[2]-c[2],x[1]-c[1])≈ϕ+∆ϕ-2π && x[3]==t*h2 && n ≈ normalize([cos(θ),sin(θ),t])) for t = (-1,1)])  # at end side, at base
    @test all([(s = one⁻⁻; p₀ = ∆vₑ/2 + c; p = [s.*p₀;t*h2]; (x,n) = surfpt_nearby(p,b); x==p && n==[0,0,t]) for t = (-1,1)])  # inside end side, at base
    @test all([(s = one⁺⁺; p₀ = ∆vₑ/2 + c; p = [s.*p₀;t*h2]; (x,n) = surfpt_nearby(p,b); atan(x[2]-c[2],x[1]-c[1])≈ϕ+∆ϕ-2π && x[3]≈t*h2 && n ≈ normalize([cos(θ),sin(θ),0])) for t = (-1,1)])  # outside end side, at base
    @test all([(s = one⁺⁺; p₀ = ∆vₑ/2 + c; p = [s.*p₀;t*h2]; (x,n) = surfpt_nearby(p,b); atan(x[2]-c[2],x[1]-c[1])≈ϕ+∆ϕ-2π && x[3]≈sign(t)*h2 && n ≈ normalize(p-x)) for t = (-one⁺⁺,one⁺⁺)])  # outside end side, outside base

    @test all([(α = ϕ+0.7∆ϕ; d̂ = [cos(α),sin(α)]; p = [s.*r.*d̂ + c; t*h2]; surfpt_nearby(p,b) ≈ ([r.*d̂.+c; t*h2], [d̂;0])) for s = (one⁻⁻,1,one⁺⁺), t = (-0.5,0,0.5)])  # around arc, within height
    @test all([(s = 1; α = ϕ+0.7∆ϕ; d̂ = [cos(α),sin(α)]; p = [s.*r.*d̂ + c; t*h2]; surfpt_nearby(p,b) ≈ ([r.*d̂.+c; t*h2], normalize([d̂;t]))) for t = (-1,1)])  # at arc, at base
    @test all([(s = one⁻⁻; α = ϕ+0.7∆ϕ; d̂ = [cos(α),sin(α)]; p = [s.*r.*d̂ + c; t*h2]; surfpt_nearby(p,b) ≈ (p, [0,0,t])) for t = (-1,1)])  # inside arc, at base
    @test all([(s = one⁺⁺; α = ϕ+0.7∆ϕ; d̂ = [cos(α),sin(α)]; p = [s.*r.*d̂ + c; t*h2]; surfpt_nearby(p,b) ≈ ([r.*d̂.+c; t*h2], normalize([d̂;0]))) for t = (-1,1)])  # outside arc, at base
    @test all([(s = one⁺⁺; α = ϕ+0.7∆ϕ; d̂ = [cos(α),sin(α)]; p = [s.*r.*d̂ + c; t*h2]; (x,n) = surfpt_nearby(p,b); (x,n) ≈ ([r.*d̂.+c; sign(t)*h2], normalize(p-x))) for t = (-one⁺⁺,one⁺⁺)])  # outside arc, outside base

    @test all([(p = [∆vₛ + c; t*h2]; surfpt_nearby(p,b) ≈ (p, normalize(normalize([cos(ϕ),sin(ϕ),0]+[cos(θ),sin(θ),0])+[0,0,t]))) for t = (-1,1)])  # at start point of arc
    @test all([(p = [∆vₑ + c; t*h2]; surfpt_nearby(p,b) ≈ (p, normalize(normalize([cos(ϕ+∆ϕ),sin(ϕ+∆ϕ),0]+[cos(θ),sin(θ),0])+[0,0,t]))) for t = (-1,1)])  # at start point of arc
    @test all([(p₀ = [∆v + c; t*h2]; p = one⁺⁺.*p₀; surfpt_nearby(p,b) ≈ (p₀, normalize(p-p₀))) for ∆v = ∆vs, t = (-1,1)])  # outside end points of arc

    @test bounds(b) ≈ ([c[1]-r, c[2]+r*sin(ϕ+∆ϕ), -h2], [c[1]+r*cos(ϕ), c[2]+r, h2])
    @test checkbounds(b)

    @test (∆ = rand(3); translate(b,∆) ≈ SectoralPrism([c;0]+∆, r, ϕ, ∆ϕ, 2h2))
end  # @testset "Half pipe"

@testset "Half pipe, rotated" begin
    ϕ = π/2  # start angle
    ∆ϕ = π
    r = 2  # radius
    h2 = 1.1
    θ = π/6  # tilt angle measured from z-axis towards x-axis
    ax = [sin(θ), 0, cos(θ)]

    b = SectoralPrism([0,0,0], r, ϕ, ∆ϕ, 2h2, ax)

    @test bounds(b) ≈ ([-r*cos(θ)-h2*sin(θ),-r,-h2*cos(θ)], [h2*sin(θ),r,r*sin(θ)+h2*cos(θ)])
    @test checkbounds(b)
end  # Half pipe, rotated


text:
"""

```
radial_distribution(
  sys::System{T},
  trajfile::String;
  dmax::Float64=10.,
  nbins::Int=100
) where T
```

Computes the g(r). If `T=Point3D` the system will be read as three-dimensional.

"""
function radial_distribution(
  sys::System{T},
  trajfile::String;
  dmax::Float64=10.,
  nbins::Int=100
) where T

  function shell_volume(r,Δr,T::DataType)
    if T == Point2D
      return π*((r+Δr)^2 - r^2)
    elseif T == Point3D
      return (4π/3)*((r+Δr)^3 - r^3)
    else
      error("T must be Point2D or Point3D")
    end
  end

  # Initialize histogram
  @unpack n, sides = sys
  x = zeros(T,n)
  density = n / prod(sides)

  Δr = dmax/nbins
  r = [ (i-1/2)*Δr for i in 1:nbins ]
  g = zeros(nbins)

  # Reading coordinates, computing distances, and adding to histogram
  nframes = 0
  open(trajfile) do file
    while !(eof(file))
      readline(file) # number of atoms, ignore
      readline(file) # title, ignore
      nframes += 1
      # read atomic coordinates
      for i in 1:n
        atom_data = split(readline(file))
        x[i] = T(ntuple(i -> parse(Float64,atom_data[i+1]), dim(T)))
      end
      # Compute distances and add to histogram
      for i in 1:n-1
        for j in i+1:n
          dx = image(x[j] - x[i],sides)
          idist = ceil(Int,nbins*(norm(dx)/dmax))
          if idist <= nbins
            g[idist] += 1.
          end
        end
      end
    end
  end
  println("Number of frames read = ", nframes)

  # Normalizing by the number of frames and number of atoms
  # Normalizing by the number of atoms expected in an area or volume
  for i in 1:nbins
    g[i] = g[i] / (nframes*((n-1)/2))
    g[i] = g[i] / (density*shell_volume(r[i],Δr,T))
  end

  return r, g
end


text:
function multiindexreverse(n,m)
  A = zeros(n,2m)
  @turbo for i=1:n÷2, j = 1:m
    a = 1.23
    b = 3.21
    A[i, j] = a
    A[i, m+j] = b
    r = n + 1 - i
    A[r, j] = a
    A[r, m+j] = b
  end
  return A
end
function multiindexreverse_ref(n,m)
  A = zeros(n,2m)
  @turbo for i=1:n÷2, j = 1:m
    a = 1.23
    b = 3.21
    A[i, j] = a
    A[i, m+j] = b
    r = n + 1 - i
    A[r, j] = a
    A[r, m+j] = b
  end
  return A
end

@testset "Multiple indices" begin
  for m in 1:10, n in 1:20
    @test multiindexreverse(n, m) == multiindexreverse_ref(n, m)
  end
end


text:
println("hello")
exit()


text:
include("../src/utilities.jl")

@testset begin "Interpolation"
    md = [1.,2.,3.,4.,5.]
    inc = [1.,2.,3.,4.,1.]
    tvd = [1.,2.,3.,4.,5.]
    id = [1.,1.,1.,2.,1.]

    well = Wellbore(md, inc, tvd, id)

    property1 = [1.,1.,2.,3.,4.,1.]
    property2 = [1.,1.,1.,1.,2.,1.]

    @test interpolate(well.md, property1, 4.5) == 2.5

    points = [1.5, 4.5]
    @test interpolate_all(well, [property1, property2], points) == hcat([1.5, 2.5], [1, 1.5])
end


text:
module CompositeTraits

using BinaryTraits
using BinaryTraits.Prefix: Has, Is, Not, Can
using Test

# possible traits
@trait Move
@trait CarryPassenger
@trait FourWheels
@trait Engine

# assignments
struct Acura end
@assign Acura with Can{Move}, Can{CarryPassenger}, Has{FourWheels}, Has{Engine}

struct Tricycle end
@assign Tricycle with Can{Move}, Can{CarryPassenger}

# composite
@trait Car with Can{Move}, Can{CarryPassenger}, Has{FourWheels}, Has{Engine}

function test()
    @test trait(Car, Acura) == Is{Car}()
    @test trait(Car, Tricycle) == Not{Car}()
end

end # module

using Test
@testset "Composite Traits" begin
    import .CompositeTraits
    CompositeTraits.test()
end


text:
"""
$(TYPEDEF)

Direction observation information of a `Point3` variable.
"""
mutable struct PriorPoint3{T} <: IncrementalInference.AbstractPrior where {T <: IIF.SamplableBelief}
  Z::T
  # W::Array{Float64,1} # TODO, deprecate the weight parameter
  PriorPoint3{T}() where T = new()
  PriorPoint3{T}(dist::T) where {T <: IIF.SamplableBelief} = new{T}(dist)
end
PriorPoint3(z::T) where {T <: IIF.SamplableBelief} = PriorPoint3{T}(z)

function getSample(p3::PriorPoint3)
  return (rand(p3.Z),)
end

## Serialized packing type and converters

"""
$(TYPEDEF)
"""
mutable struct PackedPriorPoint3  <: IncrementalInference.PackedInferenceType
    str::String
    PackedPriorPoint3() = new()
    PackedPriorPoint3(x::String) = new(x)
end

function convert(::Type{PriorPoint3}, d::PackedPriorPoint3)
  # Cov = reshapeVec2Mat(d.vecCov, d.dimc)
  distr = extractdistribution(d.str)
  return PriorPoint3{typeof(distr)}(distr)
end
function convert(::Type{PackedPriorPoint3}, d::PriorPoint3)
  # v2 = d.mv.Σ.mat[:];
  return PackedPriorPoint3(string(d.Z))
end


text: include("runbench.jl")


text:
# personalization.jl
## https://developer.spotify.com/documentation/web-api/reference/personalization/

## https://developer.spotify.com/documentation/web-api/reference/personalization/get-users-top-artists-and-tracks/
@doc """
# Get a User's Top Tracks
**Summary**: Get the current user's top trakcs based on calculated affinity

`offset` _Optional_: The index of the first tracks to return. Default 0.\n
`limit` _Optional_: The number of tracks to return. Default 20. \n
`time_range` _Optional_: Over what time frame the affinities are computer.
Valid Options:
* `long` : Calculated from several years of data including all new data as it becomes available
* `medium` : Approx. the last 6 months of data.
* `short` : Approx. the last 4 weeks

[Reference](https://developer.spotify.com/documentation/web-api/reference/personalization/get-users-top-artists-and-tracks/)
""" ->
function top_tracks(offset=0, limit=20, time_range="medium")
    return spotify_request("me/top/tracks?offset=$offset&limit=$limit&time_range=($time_range)_term")

end

@doc """
# Get a User's Top Artists
**Summary**: Get the current user's top artists based on calculated affinity

`offset` _Optional_: The index of the first artist to return. Default 0.\n
`limit` _Optional_: The number of artists to return. Default 20. \n
`time_range` _Optional_: Over what time frame the affinities are computer.
Valid Options:
* `long` : Calculated from several years of data including all new data as it becomes available
* `medium` : Approx. the last 6 months of data.
* `short` : Approx. the last 4 weeks

[Reference](https://developer.spotify.com/documentation/web-api/reference/personalization/get-users-top-artists-and-tracks/)
""" ->
function top_artists(offset=0, limit=20, time_range="medium")
    return spotify_request("me/top/artists?offset=$offset&limit=$limit&time_range=($time_range)_term")
end


text:
"""
Main module for `Documenter.jl` -- a documentation generation package for Julia.

Two functions are exported from this module for public use:

- [`makedocs`](@ref). Generates documentation from docstrings and templated markdown files.
- [`deploydocs`](@ref). Deploys generated documentation from *Travis-CI* to *GitHub Pages*.

$(EXPORTS)

"""
module Documenter

using DocStringExtensions
import Base64: base64decode
import Pkg

# Submodules
# ----------

include("Utilities/Utilities.jl")
include("DocSystem.jl")
include("Formats.jl")
include("Anchors.jl")
include("Documents.jl")
include("Builder.jl")
include("Expanders.jl")
include("CrossReferences.jl")
include("DocTests.jl")
include("DocChecks.jl")
include("Writers/Writers.jl")
include("Deps.jl")

import .Utilities: Selectors

# User Interface.
# ---------------

export Deps, makedocs, deploydocs, hide

"""
    makedocs(
        root    = "<current-directory>",
        source  = "src",
        build   = "build",
        clean   = true,
        doctest = true,
        modules = Module[],
        repo    = "",
    )

Combines markdown files and inline docstrings into an interlinked document.
In most cases [`makedocs`](@ref) should be run from a `make.jl` file:

```julia
using Documenter
makedocs(
    # keywords...
)
```

which is then run from the command line with:

```sh
\$ julia make.jl
```

The folder structure that [`makedocs`](@ref) expects looks like:

    docs/
        build/
        src/
        make.jl

# Keywords

**`root`** is the directory from which `makedocs` should run. When run from a `make.jl` file
this keyword does not need to be set. It is, for the most part, needed when repeatedly
running `makedocs` from the Julia REPL like so:

    julia> makedocs(root = Pkg.dir("MyPackage", "docs"))

**`source`** is the directory, relative to `root`, where the markdown source files are read
from. By convention this folder is called `src`. Note that any non-markdown files stored
in `source` are copied over to the build directory when [`makedocs`](@ref) is run.

**`build`** is the directory, relative to `root`, into which generated files and folders are
written when [`makedocs`](@ref) is run. The name of the build directory is, by convention,
called `build`, though, like with `source`, users are free to change this to anything else
to better suit their project needs.

**`clean`** tells [`makedocs`](@ref) whether to remove all the content from the `build`
folder prior to generating new content from `source`. By default this is set to `true`.

**`doctest`** instructs [`makedocs`](@ref) on whether to try to test Julia code blocks
that are encountered in the generated document. By default this keyword is set to `true`.
Doctesting should only ever be disabled when initially setting up a newly developed package
where the developer is just trying to get their package and documentation structure correct.
After that, it's encouraged to always make sure that documentation examples are runnable and
produce the expected results. See the [Doctests](@ref) manual section for details about
running doctests.

**`modules`** specifies a vector of modules that should be documented in `source`. If any
inline docstrings from those modules are seen to be missing from the generated content then
a warning will be printed during execution of [`makedocs`](@ref). By default no modules are
passed to `modules` and so no warnings will appear. This setting can be used as an indicator
of the "coverage" of the generated documentation.
For example Documenter's `make.jl` file contains:

```julia
makedocs(
    modules = [Documenter],
    # ...
)
```

and so any docstring from the module `Documenter` that is not spliced into the generated
documentation in `build` will raise a warning.

**`repo`** specifies a template for the "link to source" feature. If you are
using GitHub, this is automatically generated from the remote. If you are using
a different host, you can use this option to tell Documenter how URLs should be
generated. The following placeholders will be replaced with the respective
value of the generated link:

  - `{commit}` Git branch or tag name, or commit hash
  - `{path}` Path to the file in the repository
  - `{line}` Line (or range of lines) in the source file

For example if you are using GitLab.com, you could use

```julia
makedocs(repo = \"https://gitlab.com/user/project/blob/{commit}{path}#{line}\")
```

# Experimental keywords

In addition to standard arguments there is a set of non-finalized experimental keyword
arguments. The behaviour of these may change or they may be removed without deprecation
when a minor version changes (i.e. except in patch releases).

**`checkdocs`** instructs [`makedocs`](@ref) to check whether all names within the modules
defined in the `modules` keyword that have a docstring attached have the docstring also
listed in the manual (e.g. there's a `@docs` blocks with that docstring). Possible values
are `:all` (check all names) and `:exports` (check only exported names). The default value
is `:none`, in which case no checks are performed. If `strict` is also enabled then the
build will fail if any missing docstrings are encountered.

**`linkcheck`** -- if set to `true` [`makedocs`](@ref) uses `curl` to check the status codes
of external-pointing links, to make sure that they are up-to-date. The links and their
status codes are printed to the standard output. If `strict` is also enabled then the build
will fail if there are any broken (400+ status code) links. Default: `false`.

**`linkcheck_ignore`** allows certain URLs to be ignored in `linkcheck`. The values should
be a list of strings (which get matched exactly) or `Regex` objects. By default nothing is
ignored.

**`strict`** -- [`makedocs`](@ref) fails the build right before rendering if it encountered
any errors with the document in the previous build phases.

## Output formats

**`format`** allows the output format to be specified. Possible values are `:html` (default),
`:latex` and `:markdown`.

Documenter is designed to support multiple output formats. By default it is creates a set of
HTML files, but the output format can be controlled with the `format` keyword. The different
output formats may require additional keywords to be specified. The keywords for the default
HTML output are documented at the [`Writers.HTMLWriter`](@ref) module.

Documenter also has (experimental) support for Markdown and LaTeX / PDF outputs. See the
[Other outputs](@ref) for more information.

!!! warning

    The Markdown and LaTeX output formats will be moved to a separate package in future
    versions of Documenter. Automatic documentation deployments should not rely on it unless
    they fix Documenter to a minor version.

# See Also

A guide detailing how to document a package using Documenter's [`makedocs`](@ref) is provided
in the [Usage](@ref) section of the manual.
"""
function makedocs(; debug = false, args...)
    document = Documents.Document(; args...)
    cd(document.user.root) do
        Selectors.dispatch(Builder.DocumentPipeline, document)
    end
    debug ? document : nothing
end

"""
$(SIGNATURES)

Allows a page to be hidden in the navigation menu. It will only show up if it happens to be
the current page. The hidden page will still be present in the linear page list that can be
accessed via the previous and next page links. The title of the hidden page can be overriden
using the `=>` operator as usual.

# Usage

```julia
makedocs(
    ...,
    pages = [
        ...,
        hide("page1.md"),
        hide("Title" => "page2.md")
    ]
)
```
"""
hide(page::Pair) = (false, page.first, page.second, [])
hide(page::AbstractString) = (false, nothing, page, [])

"""
$(SIGNATURES)

Allows a subsection of pages to be hidden from the navigation menu. `root` will be linked
to in the navigation menu, with the title determined as usual. `children` should be a list
of pages (note that it **can not** be hierarchical).

# Usage

```julia
makedocs(
    ...,
    pages = [
        ...,
        hide("Hidden section" => "hidden_index.md", [
            "hidden1.md",
            "Hidden 2" => "hidden2.md"
        ]),
        hide("hidden_index.md", [...])
    ]
)
```
"""
hide(root::Pair, children) = (true, root.first, root.second, map(hide, children))
hide(root::AbstractString, children) = (true, nothing, root, map(hide, children))

"""
    deploydocs(
        root   = "<current-directory>",
        target = "build",
        repo   = "<required>",
        branch = "gh-pages",
        deps   = nothing | <Function>,
        make   = nothing | <Function>,
        devbranch = "master",
        devurl = "dev",
        versions = ["stable" => "v^", "v#.#", devurl => devurl]
    )

Converts markdown files generated by [`makedocs`](@ref) to HTML and pushes them to `repo`.
This function should be called from within a package's `docs/make.jl` file after the call to
[`makedocs`](@ref), like so

```julia
using Documenter, PACKAGE_NAME
makedocs(
    # options...
)
deploydocs(
    repo = "github.com/..."
)
```

When building the docs for a tag (i.e. a release) the documentation is deployed to
a directory with the tag name (i.e. `vX.Y.Z`) and to the `stable` directory.
Otherwise the docs are deployed to the directory determined by the `devurl` argument.

# Required keyword arguments

**`repo`** is the remote repository where generated HTML content should be pushed to. Do not
specify any protocol - "https://" or "git@" should not be present. This keyword *must*
be set and will throw an error when left undefined. For example this package uses the
following `repo` value:

```julia
repo = "github.com/JuliaDocs/Documenter.jl.git"
```

# Optional keyword arguments

**`root`** has the same purpose as the `root` keyword for [`makedocs`](@ref).

**`target`** is the directory, relative to `root`, where generated content that should be
deployed to `gh-pages` is written to. written to. It should generally be the same as
[`makedocs`](@ref)'s `build` and defaults to `"build"`.

**`branch`** is the branch where the generated documentation is pushed. If the branch does
not exist, a new orphaned branch is created automatically. It defaults to `"gh-pages"`.

**`deps`** is the function used to install any additional dependencies needed to build the
documentation. By default nothing is installed.

It can be used e.g. for a Markdown build. The following example installed the `pygments` and
`mkdocs` Python packages using the [`Deps.pip`](@ref) function:

```julia
deps = Deps.pip("pygments", "mkdocs")
```

**`make`** is the function used to specify an additonal build phase. By default, nothing gets
executed.

**`devbranch`** is the branch that "tracks" the in-development version of the  generated
documentation. By default this value is set to `"master"`.

**`devurl`** the folder that in-development version of the docs will be deployed.
Defaults to `"dev"`.

**`versions`** determines content and order of the resulting version selector in
the generated html. The following entries are valied in the `versions` vector:
 - `"v#"`: includes links to the latest documentation for each major release cycle
   (i.e. `v2.0`, `v1.1`).
 - `"v#.#"`: includes links to the latest documentation for each minor release cycle
   (i.e. `v2.0`, `v1.1`, `v1.0`, `v0.1`).
 - `"v#.#.#"`: includes links to all released versions.
 - `"v^"`: includes a link to the docs for the maximum version
   (i.e. a link `vX.Y` pointing to `vX.Y.Z` for highest `X`, `Y`, `Z`, respectively).
 - A pair, e.g. `"first" => "second"`, which will put `"first"` in the selector,
   and generate a url from which `"second"` can be accessed.
   The second argument can be `"v^"`, to point to the maximum version docs
   (as in e.g. `"stable" => "v^"`).

# See Also

The [Hosting Documentation](@ref) section of the manual provides a step-by-step guide to
using the [`deploydocs`](@ref) function to automatically generate docs and push them to
GitHub.
"""
function deploydocs(;
        root   = Utilities.currentdir(),
        target = "build",
        dirname = "",

        repo   = error("no 'repo' keyword provided."),
        branch = "gh-pages",
        latest::Union{String,Nothing} = nothing, # deprecated

        osname::Union{String,Nothing} = nothing, # deprecated
        julia::Union{String,Nothing} = nothing, # deprecated

        deps   = nothing,
        make   = nothing,

        devbranch = "master",
        devurl = "dev",
        versions = ["stable" => "v^", "v#.#", devurl => devurl]
    )
    # deprecation of latest kwarg (renamed to devbranch)
    if latest !== nothing
        Base.depwarn("The `latest` keyword argument has been renamed to `devbranch`.")
        devbranch = latest
        @info("setting `devbranch` to `$(devbranch)`.")
    end
    # deprecation/removal of `julia` and `osname` kwargs
    if julia !== nothing
        Base.depwarn("the `julia` keyword argument to `Documenter.deploydocs` is " *
            "removed. Use Travis Build Stages for determining from where to deploy instead. " *
            "See the section about Hosting in the Documenter manual for more details.", :deploydocs)
        @info("skipping docs deployment.")
        return
    end
    if osname !== nothing
        Base.depwarn("the `osname` keyword argument to `Documenter.deploydocs` is " *
            "removed. Use Travis Build Stages for determining from where to deploy instead. " *
            "See the section about Hosting in the Documenter manual for more details.", :deploydocs)
        @info("skipping docs deployment.")
        return
    end

    # Get environment variables.
    documenter_key      = get(ENV, "DOCUMENTER_KEY",       "")
    travis_branch       = get(ENV, "TRAVIS_BRANCH",        "")
    travis_pull_request = get(ENV, "TRAVIS_PULL_REQUEST",  "")
    travis_repo_slug    = get(ENV, "TRAVIS_REPO_SLUG",     "")
    travis_tag          = get(ENV, "TRAVIS_TAG",           "")

    # Other variables.
    sha = cd(root) do
        # We'll make sure we run the git commands in the source directory (root), in case
        # the working directory has been changed (e.g. if the makedocs' build argument is
        # outside root).
        try
            readchomp(`git rev-parse --short HEAD`)
        catch
            # git rev-parse will throw an error and return code 128 if it is not being
            # run in a git repository, which will make run/readchomp throw an exception.
            # We'll assume that if readchomp fails it is due to this and set the sha
            # variable accordingly.
            "(not-git-repo)"
        end
    end

    # Sanity checks
    if !isempty(travis_repo_slug) && !occursin(travis_repo_slug, repo)
        @warn("repo $repo does not match $travis_repo_slug")
    end

    # When should a deploy be attempted?
    should_deploy =
        occursin(travis_repo_slug, repo) &&
        travis_pull_request == "false"   &&
        (
            travis_branch == devbranch ||
            travis_tag    != ""
        )

    # check that the tag is valid
    if should_deploy && !isempty(travis_tag) && !occursin(Base.VERSION_REGEX, travis_tag)
        @warn("tag `$(travis_tag)` is not a valid VersionNumber")
        should_deploy = false
    end

    # check DOCUMENTER_KEY only if the branch, Julia version etc. check out
    if should_deploy && isempty(documenter_key)
        @warn("""
            DOCUMENTER_KEY environment variable missing, unable to deploy.
              Note that in Documenter v0.9.0 old deprecated authentication methods were removed.
              DOCUMENTER_KEY is now the only option. See the documentation for more information.""")
        should_deploy = false
    end

    if get(ENV, "DOCUMENTER_DEBUG", "") == "true"
        Utilities.debug("TRAVIS_REPO_SLUG       = \"$travis_repo_slug\"")
        Utilities.debug("  should occur in \"$repo\" (kwarg: repo)")
        Utilities.debug("TRAVIS_PULL_REQUEST    = \"$travis_pull_request\"")
        Utilities.debug("  deploying if equal to \"false\"")
        Utilities.debug("TRAVIS_BRANCH          = \"$travis_branch\"")
        Utilities.debug("TRAVIS_TAG             = \"$travis_tag\"")
        Utilities.debug("  deploying if branch equal to \"$devbranch\" (kwarg: devbranch) or tag is set")
        Utilities.debug("git commit SHA         = $sha")
        Utilities.debug("DOCUMENTER_KEY exists  = $(!isempty(documenter_key))")
        Utilities.debug("should_deploy          = $should_deploy")
    end

    if should_deploy
        # Add local bin path if needed.
        Deps.updatepath!()
        # Install dependencies when applicable.
        if deps !== nothing
            Utilities.log("installing dependencies.")
            deps()
        end
        # Change to the root directory and try to deploy the docs.
        cd(root) do
            Utilities.log("setting up target directory.")
            isdir(target) || mkpath(target)
            # Run extra build steps defined in `make` if required.
            if make !== nothing
                Utilities.log("running extra build steps.")
                make()
            end
            Utilities.log("pushing new documentation to remote: $repo:$branch.")
            mktempdir() do temp
                git_push(
                    root, temp, repo;
                    branch=branch, dirname=dirname, target=target,
                    tag=travis_tag, key=documenter_key, sha=sha,
                    devurl = devurl, versions = versions,
                )
            end
        end
    else
        Utilities.log("skipping docs deployment.")
        if get(ENV, "DOCUMENTER_DEBUG", "") != "true"
              Utilities.log("You can set DOCUMENTER_DEBUG to \"true\" in Travis to see more information.")
        end
    end
end

"""
    git_push(
        root, tmp, repo;
        branch="gh-pages", dirname="", target="site", tag="", key="", sha="", devurl="dev"
    )

Handles pushing changes to the remote documentation branch.
When `tag` is empty the docs are deployed to the `devurl` directory,
and when building docs for a tag they are deployed to a `vX.Y.Z` directory.
"""
function git_push(
        root, temp, repo;
        branch="gh-pages", dirname="", target="site", tag="", key="", sha="", devurl="dev",
        versions
    )
    dirname = isempty(dirname) ? temp : joinpath(temp, dirname)
    isdir(dirname) || mkpath(dirname)

    keyfile = abspath(joinpath(root, ".documenter"))
    target_dir = abspath(target)

    # The upstream URL to which we push new content and the ssh decryption commands.
    upstream = "git@$(replace(repo, "github.com/" => "github.com:"))"

    write(keyfile, String(base64decode(key)))
    chmod(keyfile, 0o600)

    try
        # Use a custom SSH config file to avoid overwriting the default user config.
        withfile(joinpath(homedir(), ".ssh", "config"),
            """
            Host github.com
                StrictHostKeyChecking no
                HostName github.com
                IdentityFile $keyfile
            """
        ) do
            cd(temp) do
                # Setup git.
                run(`git init`)
                run(`git config user.name "autodocs"`)
                run(`git config user.email "autodocs"`)

                # Fetch from remote and checkout the branch.
                run(`git remote add upstream $upstream`)
                run(`git fetch upstream`)

                try
                    run(`git checkout -b $branch upstream/$branch`)
                catch e
                    Utilities.log("Checking out $branch failed with error: $e")
                    Utilities.log("Creating a new local $branch branch.")
                    run(`git checkout --orphan $branch`)
                    run(`git commit --allow-empty -m "Initial empty commit for docs"`)
                end

                # Copy docs to `devurl`, or `stable`, `<release>`, and `<version>` directories.
                if isempty(tag)
                    devurl_dir = joinpath(dirname, devurl)
                    gitrm_copy(target_dir, devurl_dir)
                    Writers.HTMLWriter.generate_siteinfo_file(devurl_dir, devurl)
                    # symlink "latest" to devurl to preserve links (remove in some future release)
                    if devurl != "latest"
                        rm(joinpath(dirname, "latest"); recursive = true, force = true)
                        @warn(string("creating symlink from `latest` to `$(devurl)` for backwards ",
                            "compatibility with old links. In future Documenter versions this symlink ",
                            "will not be created. Please update any links that point to `latest`."))
                        cd(dirname) do; rm_and_add_symlink(devurl, "latest"); end
                    end
                else
                    tagged_dir = joinpath(dirname, tag)
                    gitrm_copy(target_dir, tagged_dir)
                    Writers.HTMLWriter.generate_siteinfo_file(tagged_dir, tag)
                end

                # Expand the users `versions` vector
                entries, symlinks = Writers.HTMLWriter.expand_versions(dirname, versions)

                # Create the versions.js file containing a list of `entries`.
                # This must always happen after the folder copying.
                Writers.HTMLWriter.generate_version_file(joinpath(dirname, "versions.js"), entries)

                # generate the symlinks, make sure we don't overwrite devurl
                cd(dirname) do
                    for kv in symlinks
                        i = findfirst(x -> x.first == devurl, symlinks)
                        if i === nothing
                            rm_and_add_symlink(kv.second, kv.first)
                        else
                            throw(ArgumentError(string("link `$(kv)` cannot overwrite ",
                                "`devurl = $(devurl)` with the same name.")))
                        end
                    end
                end

                # Add, commit, and push the docs to the remote.
                run(`git add -A .`)
                if !success(`git diff --cached --exit-code`)
                    run(`git commit -m "build based on $sha"`)
                    run(`git push -q upstream HEAD:$branch`)
                else
                    Utilities.log("New docs identical to the old -- not committing nor pushing.")
                end
            end
        end
    finally
        # Remove the unencrypted private key.
        isfile(keyfile) && rm(keyfile)
    end
end

function rm_and_add_symlink(target, link)
    if ispath(link)
        @warn "removing `$(link)` and linking `$(link)` to `$(target)`."
        rm(link; force = true, recursive = true)
    end
    symlink(target, link)
end

"""
    gitrm_copy(src, dst)

Uses `git rm -r` to remove `dst` and then copies `src` to `dst`. Assumes that the working
directory is within the git repository of `dst` is when the function is called.

This is to get around [#507](https://github.com/JuliaDocs/Documenter.jl/issues/507) on
filesystems that are case-insensitive (e.g. on OS X, Windows). Without doing a `git rm`
first, `git add -A` will not detect case changes in filenames.
"""
function gitrm_copy(src, dst)
    # --ignore-unmatch so that we wouldn't get errors if dst does not exist
    run(`git rm -rf --ignore-unmatch $(dst)`)
    cp(src, dst; force=true)
end

function withfile(func, file::AbstractString, contents::AbstractString)
    hasfile = isfile(file)
    original = hasfile ? read(file, String) : ""
    open(file, "w") do stream
        print(stream, contents)
        flush(stream) # Make sure file is written before continuing.
    end
    try
        func()
    finally
        if hasfile
            open(file, "w") do stream
                print(stream, original)
            end
        else
            rm(file)
        end
    end
end

function getenv(regex::Regex)
    for (key, value) in ENV
        occursin(regex, key) && return value
    end
    error("could not find key/iv pair.")
end

end # module