text:
module PlotChain
using Plots

function temp_at_index(stats_dict)
    temps = map(dict -> dict["temp"], stats_dict)
    reps = map(dict -> dict["reps"], stats_dict)
    temp_index = vcat(map(temprep -> fill(temprep[1], temprep[2]), zip(temps, reps))...)
    return temp_index
end

function plot_temps(stats_dict, filename="")
    temps = map(dict -> dict["temp"], stats_dict)
    x = map(d -> d[2]["reps"]*d[1], enumerate(stats_dict))
    plot(x, temps)
    if filename != ""
        savefig(filename)
    end
end

function plot_swaps(stats_dict, filename="")
    x_values = vcat(map(dict -> map(values -> values[1] + (dict[1] - 1)*dict[2]["reps"], dict[2]["changed-value"]), enumerate(stats))...)

    histogram(x_values, bins=length(stats_dict))
end

function plot_score(stats_dict, filename="")
    if ndims(stats_dict) == 1
        x_values = vcat(map(dict -> map(values -> Int(values[1] + (dict[1] - 1)*dict[2]["reps"]), dict[2]["changed-value"]), enumerate(stats_dict))...)
        func_values = vcat(map(dict -> map(values -> values[2], dict["changed-value"]), stats_dict)...)
        tempAtIndex = temp_at_index(stats_dict)[x_values]
        scatter(x_values, func_values, m = cgrad([:red, :black, :green, :white, :blue]), msw = 0, zcolor = tempAtIndex)
        gui()
        tempAtIndex
    else
        for i = 1:size(stats_dict, 1)
            x_values = vcat(map(dict -> map(values -> Int(values[1] + (dict[1] - 1)*dict[2]["reps"]), dict[2]["changed-value"]), enumerate(stats_dict[i,:]))...)
            func_values = vcat(map(dict -> map(values -> values[2], dict["changed-value"]), stats_dict[i,:])...)
            tempAtIndex = temp_at_index(stats_dict[i,:])[x_values]
            scatter!(x_values, func_values, m = cgrad([:red, :black, :green, :white, :blue]), msw = 0, zcolor = tempAtIndex)
            print(i)
        end
        gui()
    end
end

function plot_min(stats_dict, filename="")
    x_values = vcat(map(dict -> map(values -> values[1] + (dict[1] - 1)*dict[2]["reps"], dict[2]["new-min"]), enumerate(stats_dict))...)
    func_values = vcat(map(dict -> map(values -> values[2], dict["new-min"]), stats_dict)...)
    plot(x_values, func_values)
end

end

max_stars_repo_name: jacob-ruth/pt-optimization

max_issues_repo_name: jacob-ruth/pt-optimization

max_forks_repo_name: jacob-ruth/pt-optimization


text:
using Documenter
using Literate, JSON
using COBREXA

# some settings
dev_docs_folder = "dev"
pages_branch = "gh-pages"

# This must match the repo slug on github!
github_repo_slug = ENV["CI_PROJECT_NAMESPACE"] * "/" * ENV["CI_PROJECT_NAME"]

# generate notebooks
notebooks_path = joinpath(@__DIR__, "src", "notebooks")
notebooks_basenames = filter(x -> endswith(x, ".jl"), readdir(notebooks_path))
@info "base names:" notebooks_basenames
notebooks = joinpath.(notebooks_path, notebooks_basenames)
notebooks_outdir = joinpath(@__DIR__, "src", "notebooks")

for notebook in notebooks
    Literate.markdown(
        notebook,
        notebooks_outdir;
        repo_root_url = "https://github.com/$github_repo_slug/blob/master",
        nbviewer_root_url = "https://nbviewer.jupyter.org/github/$github_repo_slug/blob/gh-pages/$dev_docs_folder",
        binder_root_url = "https://mybinder.org/v2/gh/$github_repo_slug/$pages_branch?filepath=$dev_docs_folder",
    )
    Literate.notebook(notebook, notebooks_outdir)
end

# generate index.md from .template and the quickstart in README.md
readme = open(f -> read(f, String), joinpath(@__DIR__, "..", "README.md"))
quickstart =
    match(r"<!--quickstart_begin-->\n([^\0]*)<!--quickstart_end-->", readme).captures[1]
acks = match(
    r"<!--acknowledgements_begin-->\n([^\0]*)<!--acknowledgements_end-->",
    readme,
).captures[1]
ack_logos =
    match(r"<!--ack_logos_begin-->\n([^\0]*)<!--ack_logos_end-->", readme).captures[1]
index_md = open(f -> read(f, String), joinpath(@__DIR__, "src", "index.md.template"))
index_md = replace(index_md, "<!--insert_quickstart-->\n" => quickstart)
index_md = replace(index_md, "<!--insert_acknowledgements-->\n" => acks)
index_md = replace(index_md, "<!--insert_ack_logos-->\n" => ack_logos)
open(f -> write(f, index_md), joinpath(@__DIR__, "src", "index.md"), "w")

# copy the contribution guide
cp(
    joinpath(@__DIR__, "..", ".github", "CONTRIBUTING.md"),
    joinpath(@__DIR__, "src", "howToContribute.md"),
    force = true,
)

# a helper for sourcing the documentation files from directories
find_mds(path) =
    joinpath.(
        Ref(path),
        filter(x -> endswith(x, ".md"), readdir(joinpath(@__DIR__, "src", path))),
    )

# Documenter tries to guess the repo slug from git remote URL but that doesn't
# work really well here, this is the only fallback. If this breaks, "Edit on
# GitHub" links will stop working. (See Documenter.jl source in
# src/Utilities/Utilities.jl, in November 2021 it was around line 500) -mk
ENV["TRAVIS_REPO_SLUG"] = github_repo_slug

# build the docs
makedocs(
    modules = [COBREXA],
    clean = false,
    sitename = "COBREXA.jl",
    format = Documenter.HTML(
        # Use clean URLs, unless built as a "local" build
        prettyurls = !("local" in ARGS),
        assets = ["assets/favicon.ico"],
        highlights = ["yaml"],
    ),
    authors = "The developers of COBREXA.jl",
    linkcheck = !("skiplinks" in ARGS),
    pages = [
        "Home" => "index.md",
        "User guide" => [
            "Quickstart tutorials" =>
                vcat("All tutorials" => "tutorials.md", find_mds("tutorials")),
            "Advanced tutorials" =>
                vcat("All advanced tutorials" => "advanced.md", find_mds("advanced")),
            "Examples and notebooks" =>
                vcat("All notebooks" => "notebooks.md", find_mds("notebooks")),
        ],
        "Types and functions" => vcat("Contents" => "functions.md", find_mds("functions")),
        "How to contribute" => "howToContribute.md",
    ],
)

# remove the workaround (this would cause deploydocs() to get confused and try
# to deploy the travis way)
delete!(ENV, "TRAVIS_REPO_SLUG")

# replace the "edit this" links for the generated documentation
function replace_in_doc(filename, replacement)
    contents = open(f -> read(f, String), joinpath(@__DIR__, "build", filename))
    contents = replace(contents, replacement)
    open(f -> write(f, contents), joinpath(@__DIR__, "build", filename), "w")
end

replace_in_doc("index.html", "blob/master/docs/src/index.md" => "")
replace_in_doc(
    joinpath("howToContribute", "index.html"),
    "blob/master/docs/src/howToContribute.md" => "blob/master/.github/CONTRIBUTING.md",
)

# clean up notebooks -- we do not need to deploy all the stuff that was
# generated in the process
#
# extra fun: failing programs (such as plotting libraries) may generate core
# dumps that contain the dumped environment strings, which in turn contain
# github auth tokens. These certainly need to be avoided.
notebooks_names = [n[begin:end-3] for n in notebooks_basenames]
ipynb_names = notebooks_names .* ".ipynb"
notebooks_allowed_files = vcat("index.html", ipynb_names)
@info "allowed files:" notebooks_allowed_files
for (root, dirs, files) in walkdir(joinpath(@__DIR__, "build", "notebooks"))
    for f in files
        if !(f in notebooks_allowed_files)
            @info "removing notebook build artifact `$(joinpath(root, f))'"
            rm(joinpath(root, f))
        end
    end
end

# also remove the index template
rm(joinpath(@__DIR__, "build", "index.md.template"))

# Binder actually has 1.6.2 kernel (seen in October 2021), but for whatever
# reason it's called julia-1.1. Don't ask me.
for ipynb in joinpath.(@__DIR__, "build", "notebooks", ipynb_names)
    @info "changing julia version to 1.1 in `$ipynb'"
    js = JSON.parsefile(ipynb)
    js["metadata"]["kernelspec"]["name"] = "julia-1.1"
    js["metadata"]["kernelspec"]["display_name"] = "Julia 1.1.0"
    js["metadata"]["language_info"]["version"] = "1.1.0"
    open(f -> JSON.print(f, js), ipynb, "w")
end

# deploy the result
deploydocs(
    repo = "github.com/$github_repo_slug.git",
    target = "build",
    branch = pages_branch,
    devbranch = "develop",
)

ext: jl


text:
### A Pluto.jl notebook ###
# v0.12.10

using Markdown
using InteractiveUtils

# ╔═╡ e077ab16-22e5-11eb-15ed-5f0654bf7328
using Pkg, DrWatson

# ╔═╡ e75a9a06-22e5-11eb-2e36-8d4f62b830ed
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using CSV, CategoricalArrays, Plots, DataFrames
end

# ╔═╡ ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
md"## Listing4.29"

# ╔═╡ fa68607e-22e5-11eb-0558-c9a4d9f77426
begin
	df = CSV.read(datadir("companyData.csv"), DataFrame)
	mktCap = reshape(df.MarketCap, 5, 3)
	years  = levels(df.Year)
end

# ╔═╡ 4829a41e-2523-11eb-3554-d9499860583f
areaplot(years, mktCap,
	c=[:blue :red :green], labels=["A" "B" "C"],
	xlims=(minimum(years),maximum(years)), ylims=(0,6.5),
	legend=:topleft, xlabel="Years", ylabel="MarketCap")

# ╔═╡ 475ff888-22e6-11eb-2354-09f8f40a8e12
md"## End of listing4.29"

# ╔═╡ Cell order:
# ╟─ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
# ╠═e077ab16-22e5-11eb-15ed-5f0654bf7328
# ╠═e75a9a06-22e5-11eb-2e36-8d4f62b830ed
# ╠═fa68607e-22e5-11eb-0558-c9a4d9f77426
# ╠═4829a41e-2523-11eb-3554-d9499860583f
# ╟─475ff888-22e6-11eb-2354-09f8f40a8e12


text:
using SymbolicUtils, Test
using SymbolicUtils: polynormalize, Term
@testset "polyform" begin
    @syms a b c d
    @test polynormalize(a * (b + -1 * c) + -1 * (b * a + -1 * c * a)) == 0
    @eqtest polynormalize(sin(a+b)+sin(c+d)) == sin(a+b) + sin(c+d)
    @eqtest simplify(polynormalize(sin((a+b)^2)^2)) == simplify(sin(a^2+2*(b*a)+b^2)^2)
    @test simplify(polynormalize(sin((a+b)^2)^2 + cos((a+b)^2)^2)) == 1

    # cleanup rules
    @test polynormalize(Term{Number}(identity, 0)) == 0
    @test polynormalize(Term{Number}(one, 0)) == 1
    @test polynormalize(Term{Number}(zero, 0)) == 0
    @test polynormalize(identity(a * b) - b * a) == 0
    @test polynormalize(a * b - b * a) == 0
end


text:
#=
  Singular Value Decomposition
  Condition Number of Matrix
  Pseudoinverse of Matrix
  Eigen Decomposition
  LU Decomposition
  P-LU Decomposition
  Check Symmetric-Positive-Definite
  Cholesky Decomposition
  QR Decomposition
=#

function singularValueDecomposition(A)
  # accepts matrix A
  U = eigenvectors(A * transpose(A))
  V = eigenvectors(transpose(A) * A)
  #determine E
  v = cond(A)
  E = zeroesMatrix( , )
  for i = 1:length(v)
    E[i,i] = v[i]
  end
  println(U)
  print(E)
  print(transpose(V))

function condfromSVD(E)
  x = 0
  for i = length(E)[1] : -1 : 1
    if E[i,i] == 0
      continue
    else
      x = E[i,i]
      break
    end
  end
  return E[1,1] / x

function nullityfromSVD(E)
  x = 0
  for i = length(E)[1] : -1 : 1
    if E[i,i] == 0
      x + 1
    else
      return x
    end
  end

  function inversefromSVD(U,C,B,v)
    println(V)
    D = copy(E)
    for i = 1:length(v)
      D[i,i] = 1/E[i,i]
    end
    print(D)
    print(transpose(U))

function eigenDecomposition(A)
  # accepts
  B = eigenvectors(A)
  v = eigenvalues(A)
  for i = 1:length(v)
    #initialize matrix of zeroes size length(v) by length(v)
    C[i,i] = v[i]
  end
  println(B)
  print(C)
  print(inverse(B))

function luDecomposition(A)
  # admits square matrix
  B = copy(A)
  x = 0
  a = size(A)[1]
  #rows in A
  #construct identity matrix C of equal dimenstions

  for i=1:a #or use here
    for j = i+1 : a # and here
      x = B[i,i] / B[i,j]
      B[j,:] -= (x * B[i,:])
      C[i,j] = x

      C[i,j] = #sub value
      A[j,:] - #sub value

function luSolve(A,b)
  luDecomposition(A)
  y = solve(L,c)
  x = solve(U,y)
  return x

function pluSOlve(A,b)
  pluDecomposition(A)
  c = transpose(P) * b
  y = solve(L,c)
  x = Solve(U,x)
  return x

function checkSPD(A)
  # accepts symmetric matrix
  #via Sylvesters Criterion
  for i = 1:length(A[:,1])
    if determinant(A[i,i]) > 0
      conitnue
    else
      return false
    end
  end
  return true

function CholeskyDecomposition(A)
  # accepts SPD matrix A
  if checkSPD(A) = true
    continue
  else
    #try catch error here
  end

function QRDecomposition(A)
  # accepts ___ mxn
  B = copy(A)

  for i = 2:n
    a = B[:,i]
    for j = 1:i-1
      a - (dot(a,Q[:,j]) * Q[:,j])
    end
    Q[:,i] = normalize(a)
  end
  #TODO: R

max_stars_repo_name: K-P-S/Julia-Scientific-Computing

max_issues_repo_name: K-P-S/Julia-Scientific-Computing

max_forks_repo_name: K-P-S/Julia-Scientific-Computing


text:
module CtxMacro

using Mixtape
using MacroTools

foo(x) = x^5
bar(x) = x^10
apply(f, x1, x2::Val{T}) where {T} = f(x1, T)

function f(x)
    g = x < 5 ? foo : bar
    return g(2)
end

@ctx (false, true, false) struct MyMix end
allow(ctx::MyMix, m::Module) = m == CtxMacro

swap(e) = e
function swap(e::Expr)
    new = MacroTools.postwalk(e) do s
        isexpr(s, :call) || return s
        s.args[1] == Base.literal_pow || return s
        return Expr(:call, apply, Base.:(*), s.args[3:end]...)
    end
    return new
end

function transform(::MyMix, b)
    for (v, st) in b
        replace!(b, v, swap(st))
    end
    return b
end

Mixtape.@load_call_interface()
display(call(MyMix(), f, 3))

end # module


text:
# Test variable initializers
@testset "Variable Initializers" begin
    # initialize models
    m = InfiniteModel()
    @infinite_parameter(m, pars[1:2] in [0, 1], supports = [0, 1])
    @infinite_parameter(m, par in [0, 1], supports = [0, 1])
    @variable(m, x >= 0, Infinite(par), Int)
    @variable(m, y == 2, Infinite(par, pars), Bin, start = (p, ps) -> p + sum(ps))
    @variable(m, x0, Point(x, 0))
    @variable(m, 0 <= y0 <= 1, Point(y, 0, [0, 0]), Int)
    @variable(m, 0 <= z <= 1, Bin)
    @variable(m, w == 1, Int, start = 1)
    dx = @deriv(x, par)
    dy = @deriv(y, par)
    set_lower_bound(dx, 0)
    set_start_value_function(dy, (p, ps) -> p + sum(ps))
    var = build_variable(error, y, Dict{Int, Float64}(2 => 0))
    yrv = add_variable(m, var)
    data = DiscreteMeasureData(par, [0.5, 0.5], [0, 1])
    @constraint(m, c1, x + z - 2 <= 0)
    @constraint(m, c2, measure(x + y, data) - w == 0)
    @constraint(m, c3, x0 + y0 == 5)
    tm = optimizer_model(m)
    IOTO.set_parameter_supports(tm, m)
    # test transcribe_finite_variables!
    @testset "transcribe_finite_variables!" begin
        @test isa(IOTO.transcribe_finite_variables!(tm, m), Nothing)
        @test length(transcription_data(tm).finvar_mappings) == 2
        @test transcription_variable(tm, z) isa VariableRef
        @test transcription_variable(tm, w) isa VariableRef
        @test name(transcription_variable(tm, z)) == name(z)
        @test name(transcription_variable(tm, w)) == name(w)
        @test has_lower_bound(transcription_variable(tm, z))
        @test has_upper_bound(transcription_variable(tm, z))
        @test is_binary(transcription_variable(tm, z))
        @test is_fixed(transcription_variable(tm, w))
        @test is_integer(transcription_variable(tm, w))
        @test start_value(transcription_variable(tm, w)) == 1.
    end
    # test _format_infinite_info
    @testset "_format_infinite_info" begin
        var = InfiniteOpt._core_variable_object(y)
        @test !IOTO._format_infinite_info(var, [0., 0.5, 0.75]).has_lb
        @test !IOTO._format_infinite_info(var, [0., 0.5, 0.75]).has_ub
        @test IOTO._format_infinite_info(var, [0., 0.5, 0.75]).has_fix
        @test IOTO._format_infinite_info(var, [0., 0.5, 0.75]).fixed_value == 2
        @test IOTO._format_infinite_info(var, [0., 0.5, 0.75]).has_start
        @test IOTO._format_infinite_info(var, [0., 0.5, 0.75]).start == 1.25
        @test IOTO._format_infinite_info(var, [0., 0.5, 0.75]).binary
        @test !IOTO._format_infinite_info(var, [0., 0.5, 0.75]).integer
        var = InfiniteOpt._core_variable_object(x)
        @test IOTO._format_infinite_info(var, [0.]).has_lb
        @test IOTO._format_infinite_info(var, [0.]).lower_bound == 0
        @test !IOTO._format_infinite_info(var, [0.]).has_ub
        @test !IOTO._format_infinite_info(var, [0.]).has_fix
        @test !IOTO._format_infinite_info(var, [0.]).has_start
        @test isnan(IOTO._format_infinite_info(var, [0.]).start)
        @test IOTO._format_infinite_info(var, [0.]).integer
    end
    # test transcribe_infinite_variables!
    @testset "transcribe_infinite_variables!" begin
        @test isa(IOTO.transcribe_infinite_variables!(tm, m), Nothing)
        @test length(transcription_data(tm).infvar_mappings) == 2
        @test transcription_variable(tm, x) isa Vector{VariableRef}
        @test transcription_variable(tm, y) isa Vector{VariableRef}
        @test name(transcription_variable(tm, x)[1]) == "x(support: 1)"
        @test name(transcription_variable(tm, y)[3]) == "y(support: 3)"
        @test has_lower_bound(transcription_variable(tm, x)[1])
        @test is_binary(transcription_variable(tm, y)[2])
        @test is_fixed(transcription_variable(tm, y)[4])
        @test is_integer(transcription_variable(tm, x)[2])
        @test sort!(start_value.(transcription_variable(tm, y))) == [0., 1, 2, 3]
        @test supports(x) == [(0,), (1,)]
        @test length(supports(y)) == 4
    end
    # test _format_derivative_info
    @testset "_format_derivative_info" begin
        der = InfiniteOpt._core_variable_object(dy)
        @test !IOTO._format_derivative_info(der, [0., 0.5, 0.75]).has_lb
        @test !IOTO._format_derivative_info(der, [0., 0.5, 0.75]).has_ub
        @test !IOTO._format_derivative_info(der, [0., 0.5, 0.75]).has_fix
        @test IOTO._format_derivative_info(der, [0., 0.5, 0.75]).has_start
        @test IOTO._format_derivative_info(der, [0., 0.5, 0.75]).start == 1.25
        @test !IOTO._format_derivative_info(der, [0., 0.5, 0.75]).binary
        @test !IOTO._format_derivative_info(der, [0., 0.5, 0.75]).integer
        der = InfiniteOpt._core_variable_object(dx)
        @test IOTO._format_derivative_info(der, [0.]).has_lb
        @test IOTO._format_derivative_info(der, [0.]).lower_bound == 0
        @test !IOTO._format_derivative_info(der, [0.]).has_ub
        @test !IOTO._format_derivative_info(der, [0.]).has_fix
        @test !IOTO._format_derivative_info(der, [0.]).has_start
        @test isnan(IOTO._format_derivative_info(der, [0.]).start)
        @test !IOTO._format_derivative_info(der, [0.]).integer
    end
    # test transcribe_derivative_variables!
    @testset "transcribe_derivative_variables!" begin
        @test isa(IOTO.transcribe_derivative_variables!(tm, m), Nothing)
        @test length(transcription_data(tm).infvar_mappings) == 4
        @test transcription_variable(tm, dx) isa Vector{VariableRef}
        @test transcription_variable(tm, dy) isa Vector{VariableRef}
        @test name(transcription_variable(tm, dx)[1]) == (Sys.iswindows() ? "d/dpar[x(par)](support: 1)" : "∂/∂par[x(par)](support: 1)")
        @test name(transcription_variable(tm, dy)[3]) == (Sys.iswindows() ? "d/dpar[y(par, pars)](support: 3)" : "∂/∂par[y(par, pars)](support: 3)")
        @test has_lower_bound(transcription_variable(tm, dx)[1])
        @test sort!(start_value.(transcription_variable(tm, dy))) == [0., 1, 2, 3]
        @test supports(dx) == [(0,), (1,)]
        @test length(supports(dy)) == 4
    end
    # test _set_semi_infinite_variable_mapping
    @testset "_set_semi_infinite_variable_mapping" begin
        var = SemiInfiniteVariable(y, Dict{Int, Float64}(1 => 0), [1, 2], [1])
        vref = GeneralVariableRef(m, -1, SemiInfiniteVariableIndex)
        @test IOTO._set_semi_infinite_variable_mapping(tm, var, vref, SemiInfiniteVariableIndex) isa Nothing
        @test transcription_variable(vref) isa Vector{VariableRef}
        @test length(transcription_data(tm).infvar_mappings) == 5
        @test IOTO.lookup_by_support(tm, y, [0., 0, 0]) == IOTO.lookup_by_support(tm, vref, [0., 0])
        @test IOTO._set_semi_infinite_variable_mapping(tm, var, vref, ParameterFunctionIndex) isa Nothing
    end
    # test transcribe_semi_infinite_variables!
    @testset "transcribe_semi_infinite_variables!" begin
        @test IOTO.transcribe_semi_infinite_variables!(tm, m) isa Nothing
        @test transcription_variable(yrv) isa Vector{VariableRef}
        @test length(transcription_data(tm).infvar_mappings) == 6
        @test IOTO.lookup_by_support(tm, y, [1., 0., 0.]) == IOTO.lookup_by_support(tm, yrv, [1., 0])
    end
    # test _update_point_info
    @testset "_update_point_info" begin
        vref = @variable(tm)
        # test lower bound update
        fix(vref, 0, force = true)
        set_lower_bound(x0, 0)
        @test isa(IOTO._update_point_info(x0, vref), Nothing)
        @test lower_bound(vref) == 0
        # test upper bound update
        fix(vref, 0, force = true)
        set_upper_bound(x0, 0)
        @test isa(IOTO._update_point_info(x0, vref), Nothing)
        @test upper_bound(vref) == 0
        # test fix update
        fix(x0, 0, force = true)
        @test isa(IOTO._update_point_info(x0, vref), Nothing)
        @test fix_value(vref) == 0
        # test binary update
        unset_integer(x0)
        set_binary(x0)
        set_integer(vref)
        @test isa(IOTO._update_point_info(x0, vref), Nothing)
        @test is_binary(vref)
        # test integer update
        unset_binary(x0)
        set_integer(x0)
        @test isa(IOTO._update_point_info(x0, vref), Nothing)
        @test is_integer(vref)
        # test start value
        set_start_value(x0, 0.5)
        @test isa(IOTO._update_point_info(x0, vref), Nothing)
        @test start_value(vref) == 0.5
        # Undo changes
        unfix(x0)
        set_lower_bound(x0, 0)
        set_start_value(x0, NaN)
    end
    # test transcribe_point_variables!
    @testset "transcribe_point_variables!" begin
        @test isa(IOTO.transcribe_point_variables!(tm, m), Nothing)
        @test length(transcription_data(tm).finvar_mappings) == 4
        @test transcription_variable(tm, x0) == IOTO.lookup_by_support(tm, x, [0.])
        @test transcription_variable(tm, y0) == IOTO.lookup_by_support(tm, y, [0., 0., 0.])
        @test name(transcription_variable(tm, x0)) == "x(support: 1)"
        @test name(transcription_variable(tm, y0))[1:end-2] == "y(support: "
        @test lower_bound(transcription_variable(tm, x0)) == 0
        @test is_integer(transcription_variable(tm, x0))
        @test lower_bound(transcription_variable(tm, y0)) == 0
        @test upper_bound(transcription_variable(tm, y0)) == 1
        @test is_integer(transcription_variable(tm, y0))
        @test start_value(transcription_variable(tm, y0)) == 0.
        @test has_lower_bound(IOTO.lookup_by_support(tm, y, [0., 0., 0.]))
    end
end

# Test Measure and Objective Transcription
@testset "Measure/Objective Transcription" begin
    # Prepare the model
    m = InfiniteModel()
    @infinite_parameter(m, pars[1:2] in [0, 1], supports = [0, 1])
    @infinite_parameter(m, par in [0, 1], supports = [0, 1])
    @variable(m, x >= 0, Infinite(par))
    @variable(m, y == 2, Infinite(par, pars), start = (p, ps) -> p + sum(ps))
    @variable(m, 0 <= z <= 1, Bin)
    @variable(m, w == 1, Int, start = 1)
    meas1 = support_sum(x + 2w, par)
    meas2 = integral(w, par)
    meas3 = integral(y^2, pars)
    meas4 = integral(pars[1], pars[1])
    tm = transcription_model(m)
    IOTO.set_parameter_supports(tm, m)
    IOTO.transcribe_finite_variables!(tm, m)
    IOTO.transcribe_infinite_variables!(tm, m)
    tx = transcription_variable(x)
    ty = transcription_variable(y)
    tz = transcription_variable(z)
    tw = transcription_variable(w)
    # test transcribe_measures!
    @testset "transcribe_measures!" begin
        @test IOTO.transcribe_measures!(tm, m) isa Nothing
        @test transcription_variable(meas1) == tx[1] + 4tw + tx[2]
        @test transcription_variable(meas2) == tw + 0
        @test transcription_variable(meas3) isa Vector{AbstractJuMPScalar}
        @test transcription_variable(meas4) isa AffExpr
        @test supports(meas1) == ()
        @test supports(meas2) == ()
        @test sort!(supports(meas3)) == [(0.,), (1., )]
    end
    # test transcribe_objective!
    @testset "transcribe_objective!" begin
        @objective(m, Min, 2z^2 - meas1)
        @test IOTO.transcribe_objective!(tm, m) isa Nothing
        @test objective_sense(tm) == MOI.MIN_SENSE
        @test objective_function(tm) == 2tz^2 - tx[1] - 4tw - tx[2]
    end
end

# Test Constraint Methods
@testset "Constraint Transcription" begin
    # model setup
    m = InfiniteModel()
    @infinite_parameter(m, pars[1:2] in [0, 1], supports = [0, 1])
    @infinite_parameter(m, par in [0, 1], supports = [0, 0.5, 1])
    @variable(m, 0 <= x <= 2, Infinite(par, pars), Int)
    @variable(m, y == 0, Infinite(par))
    @variable(m, x0 == 0, Point(x, 0, [0, 0]), Bin)
    @variable(m, yf <= 2, Point(y, 1), Int)
    @variable(m, z, Bin)
    @constraint(m, c1, x^2 + 2y <= 42)
    @constraint(m, c2, y - z^2 == 0, DomainRestrictions(par => 0))
    @constraint(m, c3, 4x - 3 <= 2, DomainRestrictions(pars => 1, par => 1))
    @constraint(m, c4, y + z == 0, DomainRestrictions(par => [0, 0.5]))
    @constraint(m, c5, 2z^2 == 0, DomainRestrictions(par => 1))
    @constraint(m, c6, [z, x] in MOI.Zeros(2))
    tm = transcription_model(m)
    IOTO.set_parameter_supports(tm, m)
    IOTO.transcribe_finite_variables!(tm, m)
    IOTO.transcribe_infinite_variables!(tm, m)
    IOTO.transcribe_point_variables!(tm, m)
    xt = transcription_variable(x)
    yt = transcription_variable(y)
    x0t = transcription_variable(x0)
    yft = transcription_variable(yf)
    zt = transcription_variable(z)
    # test _get_info_constr_from_var
    @testset "_get_info_constr_from_var" begin
        # lower bounds
        set = moi_set(InfiniteOpt._core_constraint_object(LowerBoundRef(x)))
        expected = LowerBoundRef(IOTO.lookup_by_support(tm, x, [1., 1., 1.]))
        @test IOTO._get_info_constr_from_var(tm, x, set, [1., 1., 1.]) == expected
        @test IOTO._get_info_constr_from_var(tm, x, set, [0., 0., 0.]) isa Nothing
        # upper bounds
        set = moi_set(InfiniteOpt._core_constraint_object(UpperBoundRef(x)))
        expected = UpperBoundRef(IOTO.lookup_by_support(tm, x, [0., 1., 1.]))
        @test IOTO._get_info_constr_from_var(tm, x, set, [1., 1., 0.]) == expected
        @test IOTO._get_info_constr_from_var(tm, x, set, [0., 0., 0.]) isa Nothing
        set = moi_set(InfiniteOpt._core_constraint_object(UpperBoundRef(yf)))
        expected = UpperBoundRef(yft)
        @test IOTO._get_info_constr_from_var(tm, y, set, [1., 1., 1.]) == expected
        # fix
        set = moi_set(InfiniteOpt._core_constraint_object(FixRef(y)))
        expected = FixRef(IOTO.lookup_by_support(tm, y, [0.]))
        @test IOTO._get_info_constr_from_var(tm, y, set, [1., 1., 0.]) == expected
        @test IOTO._get_info_constr_from_var(tm, y, set, [0., 0., 1.]) isa Nothing
        set = moi_set(InfiniteOpt._core_constraint_object(FixRef(x0)))
        expected = FixRef(x0t)
        @test IOTO._get_info_constr_from_var(tm, x, set, [0., 0., 0.]) == expected
        # binary
        set = moi_set(InfiniteOpt._core_constraint_object(BinaryRef(x0)))
        expected = BinaryRef(IOTO.lookup_by_support(tm, x, [0., 0., 0.]))
        @test IOTO._get_info_constr_from_var(tm, x, set, [0., 0., 0.]) == expected
        @test IOTO._get_info_constr_from_var(tm, x, set, [1., 1., 1.]) isa Nothing
        @test IOTO._get_info_constr_from_var(tm, z, set, [1., 1., 1.]) == BinaryRef(zt)
        # integer
        set = moi_set(InfiniteOpt._core_constraint_object(IntegerRef(x)))
        expected = IntegerRef(IOTO.lookup_by_support(tm, x, [0., 1., 1.]))
        @test IOTO._get_info_constr_from_var(tm, x, set, [1., 1., 0.]) == expected
        @test IOTO._get_info_constr_from_var(tm, x, set, [0., 0., 0.]) isa Nothing
        set = moi_set(InfiniteOpt._core_constraint_object(IntegerRef(yf)))
        expected = IntegerRef(yft)
        @test IOTO._get_info_constr_from_var(tm, yf, set, [1., 1., 1.]) == expected
    end
    # test _support_in_restrictions
    @testset "_support_in_restrictions" begin
        @test IOTO._support_in_restrictions([0., 0., 0.], [3], [IntervalDomain(0, 0)])
        @test IOTO._support_in_restrictions([0., 0., 0.], Int[], IntervalDomain[])
        @test IOTO._support_in_restrictions([NaN, 0., 0.], [1], [IntervalDomain(0, 1)])
        @test !IOTO._support_in_restrictions([NaN, 0., 0.], [1, 2], [IntervalDomain(1, 1), IntervalDomain(1, 1)])
        @test !IOTO._support_in_restrictions([NaN, 0., 2.], [1, 3], [IntervalDomain(1, 1), IntervalDomain(1, 1)])
    end
    # test _process_constraint
    @testset "_process_constraint" begin
        # scalar constraint
        con = constraint_object(c1)
        func = jump_function(con)
        set = moi_set(con)
        @test IOTO._process_constraint(tm, con, func, set, zeros(3), "test1") isa ConstraintRef
        @test num_constraints(tm, typeof(func), typeof(set)) == 1
        cref = constraint_by_name(tm, "test1")
        delete(tm, cref)
        # vector constraint
        con = constraint_object(c6)
        func = jump_function(con)
        set = moi_set(con)
        @test IOTO._process_constraint(tm, con, func, set, zeros(3), "test2") isa ConstraintRef
        @test num_constraints(tm, typeof(func), typeof(set)) == 1
        cref = constraint_by_name(tm, "test2")
        delete(tm, cref)
        # fallback
        @test_throws ErrorException IOTO._process_constraint(tm, :bad, func, set,
                                                             zeros(3), "bad")
    end
    # test transcribe_constraints!
    @testset "transcribe_constraints!" begin
        @test IOTO.transcribe_constraints!(tm, m) isa Nothing
        # test the info constraint transcriptions
        @test transcription_constraint(LowerBoundRef(x)) isa Vector{ConstraintRef}
        @test transcription_constraint(UpperBoundRef(x)) isa Vector{ConstraintRef}
        @test transcription_constraint(IntegerRef(x)) isa Vector{ConstraintRef}
        @test length(transcription_constraint(LowerBoundRef(x))) == 5
        @test transcription_constraint(FixRef(x0)) == FixRef(x0t)
        @test transcription_constraint(BinaryRef(x0)) == BinaryRef(x0t)
        @test transcription_constraint(FixRef(y)) == FixRef.(yt)[1:2]
        @test transcription_constraint(UpperBoundRef(yf)) == UpperBoundRef(yft)
        @test transcription_constraint(BinaryRef(z)) == BinaryRef(zt)
        # test constraint transcriptions
        @test transcription_constraint(c1) isa Vector{ConstraintRef}
        @test length(transcription_constraint(c1)) == 6
        @test constraint_object(transcription_constraint(c2)).func == yt[1] - zt^2
        xf = IOTO.lookup_by_support(tm, x, [1., 1., 1.])
        @test constraint_object(transcription_constraint(c3)).func == 4xf
        @test constraint_object(transcription_constraint(c3)).set == MOI.LessThan(5.)
        expected = [yt[1] + zt, yt[2] + zt]
        @test jump_function.(constraint_object.(transcription_constraint(c4))) == expected
        @test constraint_object(transcription_constraint(c5)).func == 2zt^2
        @test length(transcription_constraint(c6)) == 6
        @test moi_set(constraint_object(first(transcription_constraint(c6)))) == MOI.Zeros(2)
        # test the info constraint supports
        expected = [([0., 0.], 0.5), ([0., 0.], 1.), ([1., 1.], 0.), ([1., 1.], 0.5), ([1., 1.], 1.)]
        @test sort(supports(LowerBoundRef(x))) == expected
        @test sort(supports(UpperBoundRef(x))) == expected
        @test sort(supports(IntegerRef(x))) == expected
        @test supports(FixRef(x0)) == ()
        @test supports(UpperBoundRef(yf)) == ()
        @test supports(BinaryRef(z)) == ()
        # test the constraint supports
        expected = [([0., 0.], 0.), ([0., 0.], 0.5), ([0., 0.], 1.), ([1., 1.], 0.), ([1., 1.], 0.5), ([1., 1.], 1.)]
        @test sort(supports(c1)) == expected
        @test supports(c2) == (0.,)
        @test supports(c3) == ([1., 1.], 1.)
        @test supports(c4) == [(0.0,), (0.5,)]
        @test supports(c5) == ()
        @test sort(supports(c6)) == expected
    end
end

# Test transcribe_derivative_evaluations! (TODO SEE IF THIS CAN BE MADE TO WORK WITH DEPENDENDENT PARAMETER BASED DERIVATIVES)
@testset "transcribe_derivative_evaluations!" begin
    # setup
    m = InfiniteModel()
    @infinite_parameter(m, t in [0, 1], num_supports = 4,
                        derivative_method = OrthogonalCollocation(4))
    @variable(m, y, Infinite(t))
    d1 = @deriv(y, t)
    d2 = @deriv(y, t^2)
    tm = transcription_model(m)
    IOTO.set_parameter_supports(tm, m)
    IOTO.transcribe_infinite_variables!(tm, m)
    IOTO.transcribe_derivative_variables!(tm, m)
    # main test
    @test IOTO.transcribe_derivative_evaluations!(tm, m) isa Nothing
    @test num_constraints(tm, AffExpr, MOI.EqualTo{Float64}) == 18
    @test num_supports(t) == 4
    @test num_supports(t, label = All) == 10
    @test length(supports(d1)) == 4
    @test length(supports(d2, label = All)) == 10
    @test IOTO.has_internal_supports(tm)
    @test has_generative_supports(t)
end

# Test build_transcription_model!
@testset "build_transcription_model!" begin
    # initialize model
    mockoptimizer = () -> MOIU.MockOptimizer(MOIU.UniversalFallback(MOIU.Model{Float64}()),
                                             eval_objective_value=false)
    m = InfiniteModel(mockoptimizer)
    @infinite_parameter(m, par in [0, 1], supports = [0, 1])
    @infinite_parameter(m, pars[1:2] in [0, 1], supports = [0, 1])
    @variable(m, 1 >= x >= 0, Infinite(par), Int)
    @variable(m, y == 2, Infinite(par, pars), Bin, start = 0)
    @variable(m, x0, Point(x, 0))
    @variable(m, 0 <= y0 <= 1, Point(y, 0, [0, 0]), Int)
    @variable(m, 0 <= z <= 1, Bin)
    @variable(m, w == 1, Int, start = 1)
    @finite_parameter(m, fin == 0)
    f = parameter_function(a -> 0, par)
    d1 = @deriv(y, par)
    d2 = @deriv(x, par^2)
    set_upper_bound(d1, 2)
    meas1 = support_sum(x - w - d2, par)
    meas2 = support_sum(y, pars)
    @constraint(m, c1, x + par - z + d1 + f == 0,
                DomainRestrictions(pars => 0, par => 0))
    @constraint(m, c2, z + x0 >= -3)
    @constraint(m, c3, meas1 + z == 0)
    @constraint(m, c4, meas2 - 2y0 + x + fin <= 1,
                DomainRestrictions(par => [0.5, 1]))
    @constraint(m, c5, meas2 == 0)
    @constraint(m, x + y == 83)
    @constraint(m, c6, [z, w] in MOI.Zeros(2))
    @objective(m, Min, x0 + meas1)
    # test basic usage
    tm = optimizer_model(m)
    @test IOTO.build_transcription_model!(tm, m,
                                          check_support_dims = false) isa Nothing
    # test finite variables
    zt = transcription_variable(z)
    wt = transcription_variable(w)
    @test zt isa VariableRef
    @test wt isa VariableRef
    @test name(zt) == name(z)
    @test name(wt) == name(w)
    @test has_lower_bound(zt)
    @test has_upper_bound(zt)
    @test is_binary(zt)
    @test is_fixed(wt)
    @test is_integer(wt)
    @test start_value(wt) == 1.
    # test infinite variables
    @test transcription_variable(x) isa Vector{VariableRef}
    @test transcription_variable(y) isa Vector{VariableRef}
    @test name(transcription_variable(x)[1]) == "x(support: 1)"
    @test name(transcription_variable(y)[3]) == "y(support: 3)"
    @test has_lower_bound(transcription_variable(x)[1])
    @test is_binary(transcription_variable(y)[2])
    @test is_fixed(transcription_variable(y)[4])
    @test is_integer(transcription_variable(x)[2])
    @test start_value(transcription_variable(y)[1]) == 0.
    @test supports(x) == [(0.,), (1.,)]
    @test length(supports(y)) == 4
    # test point variables
    @test transcription_variable(x0) isa VariableRef
    @test transcription_variable(y0) isa VariableRef
    @test name(transcription_variable(x0)) == "x(support: 1)"
    @test name(transcription_variable(y0))[1:end-2] == "y(support: "
    @test has_lower_bound(transcription_variable(x0))
    @test is_integer(transcription_variable(x0))
    @test has_lower_bound(transcription_variable(y0))
    @test has_upper_bound(transcription_variable(y0))
    @test is_integer(transcription_variable(y0))
    @test start_value(transcription_variable(y0)) == 0.
    # test derivatives
    d1t = transcription_variable(d1)
    d2t = transcription_variable(d2)
    @test length(d1t) == 4
    @test length(d2t) == 2
    @test upper_bound(d1t[1]) == 2
    @test supports(d2) == [(0.,), (1.,)]
    # test objective
    xt = transcription_variable(tm, x)
    @test objective_function(tm) == 2xt[1] + xt[2] - 2wt - d2t[1] - d2t[2]
    @test objective_sense(tm) == MOI.MIN_SENSE
    # test constraints
    yt = transcription_variable(y)
    dt_c1 = IOTO.lookup_by_support(tm, d1, zeros(3))
    @test constraint_object(transcription_constraint(c1)).func == -zt + xt[1] + dt_c1
    @test constraint_object(transcription_constraint(c2)).func == zt + xt[1]
    expected = transcription_variable(meas2)[2] - 2 * transcription_variable(y0) + xt[2]
    @test constraint_object(transcription_constraint(c4)).func == expected
    @test constraint_object(transcription_constraint(c3)).func == xt[1] - 2wt + xt[2] + zt - d2t[1] - d2t[2]
    @test constraint_object(transcription_constraint(c6)).func == [zt, wt]
    @test transcription_constraint(c5) isa Vector{ConstraintRef}
    @test name(transcription_constraint(c2)) == "c2(support: 1)"
    @test name(transcription_constraint(c1)) == "c1(support: 1)"
    @test supports(c1) == (0., [0., 0.])
    # test info constraints
    @test transcription_constraint(LowerBoundRef(z)) == LowerBoundRef(zt)
    @test transcription_constraint(UpperBoundRef(z)) == UpperBoundRef(zt)
    @test transcription_constraint(BinaryRef(z)) == BinaryRef(zt)
    @test transcription_constraint(FixRef(w)) == FixRef(wt)
    @test transcription_constraint(IntegerRef(w)) == IntegerRef(wt)
    @test transcription_constraint(LowerBoundRef(x0)) == LowerBoundRef(xt[1])
    @test transcription_constraint(UpperBoundRef(x0)) == UpperBoundRef(xt[1])
    @test transcription_constraint(IntegerRef(x0)) == IntegerRef(xt[1])
    @test transcription_constraint(LowerBoundRef(y0)) isa ConstraintRef
    @test transcription_constraint(UpperBoundRef(y0)) isa ConstraintRef
    @test transcription_constraint(IntegerRef(y0)) isa ConstraintRef
    @test transcription_constraint(LowerBoundRef(x)) == LowerBoundRef.(xt)
    @test transcription_constraint(UpperBoundRef(x)) == UpperBoundRef.(xt)
    @test transcription_constraint(IntegerRef(x)) == IntegerRef.(xt)
    @test transcription_constraint(FixRef(y)) isa Vector{ConstraintRef}
    @test transcription_constraint(BinaryRef(y)) isa Vector{ConstraintRef}
    @test transcription_constraint(UpperBoundRef(d1)) == UpperBoundRef.(d1t)
end


text:
module ConsoleProgressMonitor

# Use README as the docstring of the module:
@doc read(joinpath(dirname(@__DIR__), "README.md"), String) ConsoleProgressMonitor

using Logging: Logging, global_logger
using ProgressMeter: Progress, finish!, update!

include("router.jl")

const default_colors = [
    :green, :blue, :magenta, :cyan, :yellow, :red, :light_black,
    :light_green, :light_blue, :light_magenta, :light_cyan,
    :light_yellow, :light_red,
]

"""
    ProgressLogger(; colors, progress_options...)

# Keyword Arguments
- `colors :: Vector{Symbol}`: a list of colors used for progress meters.
- Other keyword arguments are used for constructing `ProgressMeter.Progress`.
"""
struct ProgressLogger <: Logging.AbstractLogger
    options::NamedTuple
    colors::Vector{Symbol}
    bars::Dict{Symbol, Progress}
    lastid::typeof(Ref(:_))
end

const _noid = gensym(:_noid)

ProgressLogger(options::NamedTuple, colors::Vector{Symbol} = default_colors) =
    ProgressLogger(options, colors, Dict(), Ref(_noid))
ProgressLogger(; colors=default_colors, options...) =
    ProgressLogger((; options...), colors)

# https://docs.julialang.org/en/latest/stdlib/Logging/#AbstractLogger-interface-1

likelytoprint(p) = time() + 1e-3 - p.tlast > p.dt

pickcycle(xs, i) = xs[mod1(i, length(xs))]

somestring(x) = x
somestring(x, xs...) =
    x isa AbstractString && !isempty(x) ? x : somestring(xs...)

function Logging.handle_message(
    logger::ProgressLogger,
    level, title, _module, group, id, file, line;
    progress = nothing,
    message = nothing,
    _...
)
    n = 1000
    if progress isa Real && (progress <= 1 || isnan(progress))
        p = get!(logger.bars, id) do
            color = pickcycle(logger.colors, length(logger.bars) + 1)
            Progress(n; color=color, logger.options...)
        end

        progress = isnan(progress) ? 0.0 : progress
        counter = floor(Int, progress * n)

        desc = somestring(title, message, p.desc)
        if !endswith(desc, " ")
            desc = string(desc, ": ")
        end
        p.desc = desc

        if logger.lastid[] ∉ (id, _noid) && likelytoprint(p)
            # Switched from unfinished progress bar:
            println(p.output)
        end
        update!(p, counter)

        if p.printed
            logger.lastid[] = id
        end
    elseif progress == "done" || (progress isa Real && progress > 1)
        p = pop!(logger.bars, id, nothing)
        p === nothing && return
        finish!(p)
        logger.lastid[] = _noid
    end
    return
end

Logging.shouldlog(::ProgressLogger, level, _module, group, id) =
    true

Logging.min_enabled_level(::ProgressLogger) = Logging.LogLevel(-1)

"""
    with_progresslogger(f; options...)

Run `f` with `ProgressLogger` enabled.
"""
with_progresslogger(f; options...) =
    Logging.with_logger(f, ProgressLogger(; options...))

"""
    install_logger(; options...)
    install_logger(logger::ProgressLogger)

Install `ProgressLogger` to global logger.

Keyword arguments `options` are passed to `ProgressLogger` constructor.
"""
install_logger(; options...) = install_logger(ProgressLogger(; options...))

function install_logger(logger::ProgressLogger)
    global previous_logger
    previous_logger = global_logger(ProgressLogRouter(global_logger(), logger))
end

"""
    uninstall_logger()

Rollback the global logger to the one before last call of `install_logger`.
"""
function uninstall_logger()
    global previous_logger
    previous_logger === nothing && return
    ans = global_logger(previous_logger)
    previous_logger = nothing
    return ans
end

end # module


text:
# This file is an example for the chemfiles library
# Any copyright is dedicated to the Public Domain.
# http://creativecommons.org/publicdomain/zero/1.0/
#!/usr/bin/env julia
using Chemfiles

topology = Topology()
add_atom!(topology, Atom("H"))
add_atom!(topology, Atom("O"))
add_atom!(topology, Atom("H"))

add_bond!(topology, 0, 1)
add_bond!(topology, 2, 1)

frame = Frame()
resize!(frame, 3)
set_topology!(frame, topology)

pos = frame.positions()
pos[1, :] = [1.0, 0.0, 0.0]
pos[2, :] = [0.0, 0.0, 0.0]
pos[3, :] = [0.0, 1.0, 0.0]

add_atom!(frame, Atom("O"), [5.0, 0.0, 0.0])
add_atom!(frame, Atom("C"), [6.0, 0.0, 0.0])
add_atom!(frame, Atom("O"), [7.0, 0.0, 0.0])
add_bond!(frame, 3, 4)
add_bond!(frame, 4, 5)

set_cell!(frame, UnitCell(10, 10, 10))

Trajectory("water-co2.pdb", 'w') do trajectory
    write(trajectory, frame)
end


text:
#=== Simplex ===#
struct Simplex{P} <: AbstractSimplex
    vs::Set{P}
    hash::UInt64
    Simplex{P}(vals) where {P} = new(vals, hash(vals))
    Simplex{P}() where {P} = new(Set{P}())
end
Simplex(splx::Vector{P}) where {P} = Simplex{P}(Set(splx))
Simplex(splx::P...) where {P} = Simplex(P[splx...])

# Private methods

Base.convert(::Type{Simplex{P}}, v::Vector{P}) where {P} = Simplex{P}(0, Set(v))
Base.show(io::IO, splx::Simplex) = show(io, "σ$(collect(splx.vs))")
Base.eltype(splx::Simplex{P}) where {P} = P

# Public methods

dim(splx::Simplex) = length(splx.vs)-1

values(splx::Simplex) = collect(splx.vs)

hash(splx::Simplex) = splx.hash

==(a::Simplex, b::Simplex) = a.hash == b.hash

function faces(splx::Simplex)
    faces = typeof(splx)[]
    for i in 1:dim(splx)+1
        face = collect(splx.vs)
        deleteat!(face,i)
        push!(faces, Simplex(face))
    end
    return faces
end

function vertecies(splx::Simplex)
    vtxs = typeof(hash(splx))[]
    for v in splx.vs
        push!(vtxs, hash(Simplex(v)))
    end
    return vtxs
end

union(u::Simplex, v::Simplex) = Simplex(collect(u.vs ∪ v.vs))

# Misc. methods

function volume(S::AbstractMatrix)
    d, vc = size(S)
    @assert d == vc-1 "Number of vertexes in simplex should be dim+1"
    v0 = S[:,1]
    return abs(det(S[:,2:end] .- v0))/prod(1:d)
end
volume(s::Simplex{Int}, X::AbstractMatrix) = volume(X[:,collect(s.vs)])

# iterator

struct Simplices{T}
    itr::T
    dim::Int
end
simplices(itr::T, dim::Int=-1) where T = Simplices{T}(itr, dim)
Base.show(io::IO, splxs::Simplices{T}) where T = print(io, "Simplex Iterator", splxs.dim < 0 ? "" : " (d=$(splxs.dim))", " for $T")

# Simplicial Complex

"""Simplicial Complex Type

Create a simplicial complex with simplices with a value type `P`.
"""
mutable struct SimplicialComplex{S<:AbstractSimplex} <: AbstractComplex
    cells::Dict{Int,Vector{S}}   # cells per dimension
    order::Dict{UInt64,Int}      # total order of cells per dimension
end
Base.show(io::IO, cplx::SimplicialComplex) = print(io, "SimplicialComplex($(size(cplx)))")

#
# Constructors
#

SimplicialComplex(::Type{S}) where {S<:AbstractSimplex} = SimplicialComplex(Dict{Int,Vector{S}}(), Dict{UInt64,Int}())
(::Type{SimplicialComplex{S}})() where {S<:AbstractSimplex} = SimplicialComplex(S)

function SimplicialComplex(splxs::Simplex...)
    S = eltype(splxs)
    cplx = SimplicialComplex(S)
    added = Set{S}()
    toprocess = Vector{S}()
    for splx in splxs
        push!(toprocess, splx) # add simplex to process its faces
        while(length(toprocess) > 0)
            tmp = pop!(toprocess) # processing simplex faces
            tmp in added && continue # skip if already processed
            addsimplex!(cplx, tmp) # add simples to complex
            push!(added, tmp) # mark as processed
            for f in faces(tmp) # add simplex faces for processing
                push!(toprocess, f)
            end
        end
    end

    return cplx
end

SimplicialComplex(splxs::Vector{<:AbstractSimplex}) = SimplicialComplex(splxs...)

# ---------------
# Private Methods
# ---------------

"""Add a simplex to the complex (as is) and return a complex size, a dimension of simplex and its index in it

This function **doesn't** add missing faces of the simplex to the complex. Use `addsimplices!` function for instead.
"""
function addsimplex!(cplx::SimplicialComplex, splx::AbstractSimplex)
    d = dim(splx)
    d < 0 && return nothing
    !haskey(cplx.cells, d) && setindex!(cplx.cells, celltype(cplx)[], d)
    push!(cplx.cells[d], splx)
    cplx.order[hash(splx)] = length(cplx.cells[d]) # length(cplx.order)+1
    return splx
end

"""Add a simplex to the complex and all of its faces recursivly and return a complex size, a dimension of simplex and its index in it"""
function addsimplices!(cplx::SimplicialComplex, splx::AbstractSimplex)
    CT = celltype(cplx)
    # cache already
    added = Set{CT}()

    toprocess = CT[]
    ret = CT[]
    addedidxs = NTuple{3,Integer}[]

    # add simplex to complex
    s = addsimplex!(cplx, splx)
    push!(addedidxs, (sum(size(cplx)), dim(s), hash(s)))
    push!(ret, s)
    for f in faces(splx) # add simplex faces for processing
        push!(toprocess, f)
    end

    # process simplex faces and add missing simplexes
    while(length(toprocess) > 0)
        tmp = pop!(toprocess) # processing simplex faces
        d = dim(tmp)

        cplx[tmp] !== nothing && continue # skip if already in complex
        tmp in added && continue # skip if already processed

        s = addsimplex!(cplx, tmp) # add simples to the complex
        push!(addedidxs, (sum(size(cplx)), dim(s), hash(s)))
        push!(added, tmp) # mark as processed
        push!(ret, s)

        dim(tmp) == 0 && continue # do not create faces for 0-simplexes
        for f in faces(tmp) # add simplex faces for processing
            push!(toprocess, f)
        end
    end

    return ret
end

#
# AbstractComplex Interface
#
celltype(cplx::SimplicialComplex) = eltype(valtype(cplx.cells))

function cells(cplx::SimplicialComplex)
    CCT = valtype(cplx.cells)
    length(cplx.cells) == 0 && return CCT[] # no cell in complex
    dims = maximum(keys(cplx.cells))
    return CCT[haskey(cplx.cells, d) ? cplx.cells[d] : CCT() for d in 0:dims]
end

cells(cplx::SimplicialComplex, d::Int) = get(cplx.cells, d,  nothing)

function boundary(cplx::SimplicialComplex, idx::IX, d::Int, ::Type{PID}) where {PID, IX<:Integer}
    ch = Chain(d-1, PID, IX)
    d == 0 && return ch

    pos = one(PID)
    neg = -one(PID)

    splx = cplx[idx, d]
    splx === nothing && return ch

    sgn = true
    for face in faces(splx)
        fidx = cplx[face]
        push!(ch, (sgn ? pos : neg), fidx)
        sgn = !sgn
    end

    return ch
end

function coboundary(cplx::SimplicialComplex, idx::IX, d::Int, ::Type{PID}) where {PID, IX<:Integer}
    ch = Chain(d+1, PID, IX)
    d == 0 && return ch

    # (δϕ)(c) = ϕ(∂c)
    for c in cells(cplx, d+1)
        i = hash(c)
        for (coef, elem) in boundary(cplx, i, d+1, PID)
            elem == idx && push!(ch, coef, i)
        end
    end

    return ch
end

Base.push!(cplx::SimplicialComplex, splx::AbstractSimplex; recursive=false) =
    recursive ? addsimplices!(cplx, splx) : [addsimplex!(cplx, splx)]

Base.position(cplx::SimplicialComplex, idx::Integer, d::Int) where {C<:AbstractCell} = get(cplx.order, idx, nothing)

#
# Miscellaneous
#

Base.similar(cplx::SimplicialComplex) = SimplicialComplex(eltype(celltype(cplx)))

function Base.read(io::IO, ::Type{SimplicialComplex{S}}) where {S<:AbstractSimplex}
    splxs = S[]
    P = eltype(S())
    while !eof(io)
        l = chomp(readline(io))
        si = map(e->parse(P, e), split(l, ' '))
        push!(splxs, Simplex(si...))
    end
    return SimplicialComplex(splxs)
end

function Base.write(io::IO, cplx::SimplicialComplex)
    zsplxs = Dict{Int,UInt}()
    for s in cplx.cells[0]
        idx = cplx.order[hash(s)]
        zsplxs[first(values(s))] = idx
        write(io, "$idx")
        write(io, 0x0A)
    end
    for d in 1:dim(cplx)
        for s in cplx.cells[d]
            for (j,i) in enumerate(values(s))
                write(io, "$(zsplxs[i])")
                if j>d
                    write(io, 0x0A)
                else
                    write(io, ' ')
                end
            end
        end
    end
    return
end

#
# Complex simplex iterator
#
Base.length(splxs::Simplices{C}) where C<:AbstractComplex =
    splxs.dim < 0 ? sum(size(splxs.itr)) : size(splxs.itr, splxs.dim)

Base.eltype(iter::Simplices{C}) where C<:AbstractComplex = celltype(iter.itr)

# State (total id, dim, dim id)
function Base.iterate(iter::Simplices{C}, (tid, d, did)=(0, iter.dim, 1)) where C<:AbstractComplex
    tid >= length(iter) && return nothing
    if d < 0
        d = 0
    end
    if did > size(iter.itr, d)
        d += 1
        did = 1
    end
    return cells(iter.itr, d)[did], (tid+1, d, did+1)
end


text:
struct StandingController{M<:MomentumBasedController}
    lowlevel::M
    robotmass::Float64

    foottasks::Dict{RigidBody{Float64}, SpatialAccelerationTask}
    linmomtask::LinearMomentumRateTask
    pelvistask::AngularAccelerationTask
    jointtasks::Dict{JointID, JointAccelerationTask{Revolute{Float64}}}

    comgains::PDGains{Float64,Float64}
    pelvisgains::PDGains{Float64,Float64}
    jointgains::Dict{JointID, PDGains{Float64, Float64}}

    comref::Point3D{SVector{3, Float64}}
    jointrefs::Dict{JointID, Float64}
end

function StandingController(
        lowlevel::MomentumBasedController,
        feet::Vector{<:RigidBody},
        pelvis::RigidBody,
        nominalstate::MechanismState;
        joint_regularization::Float64 = 0.05,
        linear_momentum_weight::Float64 = 1.0,
        comgains::PDGains = PDGains(10., 2 * sqrt(10.0)),
        pelvisgains::PDGains = PDGains(20., 2 * sqrt(20.0)),
        jointgains = Dict(JointID(j) => PDGains(100.0, 20.) for j in tree_joints(lowlevel.state.mechanism)),
        comref::Point3D = center_of_mass(nominalstate) - FreeVector3D(root_frame(lowlevel.state.mechanism), 0., 0., 0.05),
        jointrefs = Dict(JointID(j) => configuration(nominalstate, j)[1] for j in tree_joints(lowlevel.state.mechanism) if joint_type(j) isa Revolute))
    mechanism = lowlevel.state.mechanism
    world = root_body(mechanism)
    worldframe = root_frame(mechanism)
    m = mass(mechanism)

    regularize!.(Ref(lowlevel), tree_joints(mechanism), joint_regularization)

    foottasks = Dict(foot => SpatialAccelerationTask(mechanism, path(mechanism, world, foot)) for foot in feet)
    addtask!.(Ref(lowlevel), collect(values(foottasks)))

    linmomtask = LinearMomentumRateTask(mechanism, centroidal_frame(lowlevel))
    addtask!(lowlevel, linmomtask, linear_momentum_weight)

    pelvistask = AngularAccelerationTask(mechanism, path(mechanism, world, pelvis))
    addtask!(lowlevel, pelvistask)

    revolutejoints = filter(j -> joint_type(j) isa Revolute, tree_joints(mechanism))
    positioncontroljoints = setdiff(revolutejoints, vcat((collect(task.path) for task in values(foottasks))...))
    jointtasks = Dict(JointID(j) => JointAccelerationTask(j) for j in positioncontroljoints)
    addtask!.(Ref(lowlevel), collect(values(jointtasks)))

    StandingController(
        lowlevel, m,
        foottasks, linmomtask, pelvistask, jointtasks,
        comgains, pelvisgains, jointgains,
        comref, jointrefs)
end

function (controller::StandingController)(τ::AbstractVector, t::Number, state::MechanismState)
    # Linear momentum control
    m = controller.robotmass
    c = center_of_mass(state)
    centroidal = centroidal_frame(controller.lowlevel)
    world_to_centroidal = Transform3D(c.frame, centroidal, -c.v)
    e = FreeVector3D(centroidal, -transform(controller.comref, world_to_centroidal).v)
    ė = FreeVector3D(centroidal, linear(transform(momentum(state), world_to_centroidal)) / m)
    c̈ = pd(controller.comgains, e, ė)
    setdesired!(controller.linmomtask, m * c̈)

    # Pelvis orientation control
    pelvistask = controller.pelvistask
    pelvisgains = controller.pelvisgains
    pelvis = target(pelvistask.path)
    Hpelvis = transform_to_root(state, pelvis)
    Tpelvis = transform(twist_wrt_world(state, pelvis), inv(Hpelvis))
    ωdesired = FreeVector3D(Tpelvis.frame, pd(pelvisgains, rotation(Hpelvis), Tpelvis.angular))
    setdesired!(pelvistask, ωdesired)

    # Joint position control
    for jointid in keys(controller.jointtasks)
        task = controller.jointtasks[jointid]
        gains = controller.jointgains[jointid]
        ref = controller.jointrefs[jointid]
        v̇desired = pd(gains, configuration(state, jointid)[1], ref, velocity(state, jointid)[1], 0.0)
        setdesired!(task, v̇desired)
    end

    controller.lowlevel(τ, t, state)
    τ
end


text:
module Mash

using ..Graphs: sequence

const map = UInt64[
    4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
    4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
    4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
    4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
    4, 0, 4, 1, 4, 4, 4, 2, 4, 4, 4, 4, 4, 4, 4, 4,
    4, 4, 4, 4, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
    4, 0, 4, 1, 4, 4, 4, 2, 4, 4, 4, 4, 4, 4, 4, 4,
    4, 4, 4, 4, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
    4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
    4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
    4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
    4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
    4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
    4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
    4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
    4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
]

const maxU64 = typemax(UInt64)

"""
    struct Minimizer
        value    :: UInt64
        position :: UInt64
    end

A minimizer is a kmer that, given a hash function that maps kmers to integers, is the minimum kmer within a given set of kmers.
The value is the result of applying the hash function to the kmer.
The position is a bitpacked integer that includes reference ID, locus, and strand
"""
struct Minimizer
    value    :: UInt64
    position :: UInt64
end

# transliteration of the invertible hash function found in minimap
"""
    hash(x::UInt64, mask::UInt64)

A transliteration of Jenkin's invertible hash function for 64 bit integers.
Bijectively maps any kmer to an integer.
"""
function hash(x::UInt64, mask::UInt64)::UInt64
    x = (~x + (x << 21)) & mask
    x = x ⊻ x >> 24
    x = (x + (x << 3) + (x << 8)) & mask
    x = x ⊻ x >> 14
    x = (x + (x << 2) + (x << 4)) & mask
    x = x ⊻ x >> 28
    x = (x + (x << 31)) & mask
    return x
end

"""
    sketch(seq::Array{UInt8}, k::Int, w::Int, id::Int)

Sketch a linear sequence into a vector of minimizers.
`k` sets the kmer size.
`w` sets the number of contiguous kmers that will be used in the window minimizer comparison.
`id` is a unique integer that corresponds to the sequence. It will be bitpacked into the minimizer position.
"""
function sketch(seq::Array{UInt8}, k::Int, w::Int, id::Int)
    (k < 0 || k > 32)  && error("k='$(k)' must be ∈ [0,32]")
    (w < 0 || w > 255) && error("w='$(w)' must be ∈ [0,255]")

    fwd :: UInt64 = 0
    rev :: UInt64 = 0

    mask  :: UInt64 = (1 << (2*k)) - 1;
    shift :: UInt64 = 2*(k-1);

    min = Minimizer(maxU64, maxU64)
    minimizer = Minimizer[]
    window = fill(Minimizer(maxU64,maxU64), w)

    l = 0
    bi, mi = 1, 1
    for (locus,nuc) in enumerate(seq)
        c = map[nuc+1] # need to offset as Julia is 1-indexed
        new = if c ≥ 4
            l = 0
            Minimizer(maxU64,maxU64)
        else
            fwd = ((fwd<<2) | c) & mask
            rev = ((rev>>2) | ((3⊻c) << shift))
            l += 1

            if l ≥ k
                pos = (UInt64(id) << 32) | (UInt64(locus) << 1)
                if fwd ≤ rev
                    Minimizer(hash(fwd,mask),pos)
                else
                    Minimizer(hash(rev,mask),pos|1)
                end
            else
                Minimizer(maxU64,maxU64)
            end
        end

        window[bi] = new
        if (l == w+k-1) && (min.value != maxU64)
            for i in (bi+1):w
                if min.value == window[i].value && min.position != window[i].position
                    push!(minimizer, window[i])
                end
            end
            for i in 1:bi
                if min.value == window[i].value && min.position != window[i].position
                    push!(minimizer, window[i])
                end
            end
        end

        if new.value < min.value
            if l ≥ w+k && min.value != maxU64
                push!(minimizer, min)
            end
            min = new
            mi = bi
        elseif bi == mi
            if (l ≥ w+k-1) && min.value != maxU64
                push!(minimizer, min)
            end

            min = Minimizer(maxU64, min.position)
            for i in (bi+1):w
                if window[i].value < min.value
                    mi  = i
                    min = window[i]
                end
            end

            for i in 1:bi
                if window[i].value < min.value
                    mi  = i
                    min = window[i]
                end
            end

            if (l ≥ w+k-1) && min.value != maxU64
                for i in (bi+1):w
                    if min.value == window[i].value && min.position != window[i].position
                        push!(minimizer, window[i])
                    end
                end
                for i in 1:bi
                    if min.value == window[i].value && min.position != window[i].position
                        push!(minimizer, window[i])
                    end
                end
            end
        end

        bi += 1
        if bi > w
            bi = 1
        end
    end

    if min.value != maxU64
        push!(minimizer, min)
    end

    return minimizer
end

tuples(iter) = ((x,y) for (i,x) in enumerate(iter) for y in iter[i:end])

"""
    distance(graphs...; k=15, w=100)

Compute the pairwise distance between all input graphs.
Distance is the set distance between minimizers.
Linear-time algorithm using hash collisions.
"""
function distance(graphs...; k=15, w=100)
    sequences = Dict(seq for graph in graphs for seq in sequence(graph))

    names = collect(keys(sequences))
    seqs = [sequences[name] for name in names]

    minimizers = Minimizer[]
    for (i,seq) in enumerate(seqs)
        append!(minimizers, sketch(Array{UInt8}(seq), k, w, i))
    end

    sort!(minimizers, by=(m)->m.value)
    distance = zeros(Float64, length(seqs), length(seqs))

    l, r = 1, 1
    while l ≤ length(minimizers)
        while r ≤ length(minimizers) && minimizers[r].value == minimizers[l].value
            r += 1
        end

        hits = sort(collect(Set(Int(m.position >> 32) for m in minimizers[l:(r-1)])))

        for (i,j) in tuples(hits)
            distance[i,j] += 1
        end

        l = r
    end

    for i ∈ 1:length(seqs)
        for j ∈ (i+1):length(seqs)
            distance[i,j] = 1 - distance[i,j]/distance[i,i]
            distance[j,i] = distance[i,j]
        end
        distance[i,i] = 0
    end

    return distance, names
end

end

max_stars_repo_name: neherlab/pangraph

max_issues_repo_name: neherlab/pangraph

max_forks_repo_name: neherlab/pangraph


text:
type Portfolio <: AbstractTimeSeries
    timestamp::Vector{Date}
    values::Matrix{Float64}
    colnames::Vector{ASCIIString}
#    books::Vector{OrderBook}
    blotters::Vector{Portfolio}

    function Portfolio(timestamp::Vector{Date},
                       values::Matrix{Float64},
                       colnames::Vector{ASCIIString},
#                       books::Vector{OrderBook},
                       blotters::Vector{Portfolio})

                       nrow, ncol = size(values, 1), size(values, 2)
                       nrow != size(timestamp, 1) ? error("values must match length of timestamp"):
                       ncol != size(colnames,1) ? error("column names must match width of array"):
                       timestamp != unique(timestamp) ? error("there are duplicate dates"):
                       ~(flipud(timestamp) == sort(timestamp) || timestamp == sort(timestamp)) ? error("dates are mangled"):
                       flipud(timestamp) == sort(timestamp) ?
                       new(flipud(timestamp), flipud(values), colnames):
                       new(timestamp, values, colnames, books, blotters)
    end
end

#function Portfolio(books::Vector{OrderBook}, blotters::Vector{Blotter})
function Portfolio(blotters::Vector{Blotter})
    #tstamps = sort(unique(vcat([b.timestamp for b in books], [bb.timestamp for bb in blotters])))
    tstamps = sort(unique(vcat([b.timestamp for b in blotters])))
#   vals  are each blotter's log return (zero when NA) plus
#         the equity and cah values
#   colnames should be each blotter's ticker, + equity and cash
end


text:
struct AngularPositionAroundPoint{M,N,S,m} <: TaskMap{M,N,S}
    position_center::SVector{m,S}
end
AngularPositionAroundPoint{M,N}(c::SVector{m,S}) where {M,N,S,m} =
    AngularPositionAroundPoint{M,N,S,dim(M)}(c)

struct AngularPositionAroundPointT{M,N,S,m} <: TaskMapT{M,N,S}
    base_map::AngularPositionAroundPoint{M,N,S,m}
end
AngularPositionAroundPointT{M,N}(c::SVector{m,S}) where {M,N,S,m} =
    AngularPositionAroundPointT{M,N,S,m}(AngularPositionAroundPoint{M,N}(c))

# S1
codomain_coord_rep(::AngularPositionAroundPoint{R2,S1}) = ChartRep()

function task_map_emb_chart(::EmbRep, ::ChartRep, xme, task_map::AngularPositionAroundPoint{R2,S1},
        ::Chart{Angleπ,S1})
    a = xme - task_map.position_center
    d = norm(a)
    if d == 0
        return SA[0.]
    else
        â = a ./ norm(a)
        return SA[atan(â[2], â[1])]
    end
end

function task_map_emb_chart(::EmbRep, ::ChartRep, xme, task_map::AngularPositionAroundPoint{R2,S1},
        C2π::Chart{Angle2π,S1})
    Cπ = Chart{Angleπ,S1}()
    xnπ = task_map_emb_chart(xme, task_map, Cπ)
    chart_transition(xnπ, Cπ, C2π)
end

# S2
codomain_coord_rep(::AngularPositionAroundPoint{R3,S2}) = EmbRep()

function task_map_emb(::EmbRep, ::EmbRep, xme, task_map::AngularPositionAroundPoint{R3,S2})
    Δx = xme - task_map.position_center
    d = norm(Δx)
    if d == 0
        return @SVector zeros(3)
    else
        return Δx/norm(Δx)
    end
end


text:
using Documenter
using MyAwesomePackage

if haskey(ENV, "DOCSARGS")
    for arg in split(ENV["DOCSARGS"])
        (arg in ARGS) || push!(ARGS, arg)
    end
end

#push!(LOAD_PATH, "../src/")
makedocs(
    modules = [MyAwesomePackage],
    sitename = "MyAwesomePackage.jl Documentation",
    pages = [
        "Index" => "index.md",
        "Another page" => "anotherPage.md",
    ],

    format = Documenter.HTML(
        #prettyurls = true
        canonical = "https://sgibbs409.github.io/MyAwesomePackage.jl/stable/",
        prettyurls = !("local" in ARGS)
    ),

    strict = !("strict=false" in ARGS),
    doctest = ("doctest=only" in ARGS) ? :only : true,
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo = "github.com/sgibbs409/MyAwesomePackage.jl.git",
    devbranch="main",
    target="build",
    push_preview = true,
    #versions = ["stable" => "v^", "v#.#", "dev" => "v^"],
    #versions = ["v#.#", "dev" => "v^"],
    #versions = ["stable" => "v^", "v#.#"]
)


text:
const DOC_PLOT_PARAMS = """
- **`title`** : text to display on the top of the plot.

- **`xlabel`** : text to display on the `x` axis of the plot.

- **`ylabel`** : text to display on the `y` axis of the plot.

- **`xscale`** : `x`-axis scale `(:identity, :ln, :log2, :log10)`,
  or scale function e.g. `x -> log10(x)`.

- **`yscale`** : `y`-axis scale.

- **`labels`** : can be used to hide the labels by setting `labels = false`.

- **`border`** : the style of the bounding box of the plot.
  Supports `:corners`, `:solid`, `:bold`, `:dashed`, `:dotted`, `:ascii`, and `:none`.

- **`compact`** : compact plot (labels), defaults to `false`.

- **`blend`** : blend colors on the underlying canvas, defaults to `true`.

- **`margin`** : number of empty characters to the left of the whole plot.

- **`padding`** : Space of the left and right of the plot between the labels and the canvas.

- **`color`** : color of the drawing.
  Can be any of `:green`, `:blue`, `:red`, `:yellow`, `:cyan`, `:magenta`, `:white`,
  `:normal`, an integer in the range `0`-`255` for Color256 palette,
  or a tuple of three integers as RGB components.

- **`colorbar_lim`** : colorbar limit, defaults to `(0, 1)`.

- **`colorbar_border`** : colorbar border, defaults to `:solid`.

- **`width`** : number of characters per row that should be used for plotting.
"""

#! format: off
const BORDER_SOLID = (
    tl = '┌',
    tr = '┐',
    bl = '└',
    br = '┘',
    t = '─',
    l = '│',
    b = '─',
    r = '│',
)
const BORDER_CORNERS = (
    tl = '┌',
    tr = '┐',
    bl = '└',
    br = '┘',
    t = ' ',
    l = ' ',
    b = ' ',
    r = ' ',
)
const BORDER_BARPLOT = (
    tl = '┌',
    tr = '┐',
    bl = '└',
    br = '┘',
    t = ' ',
    l = '┤',
    b = ' ',
    r = ' ',
)
const BORDER_BOLD = (
    tl = '┏',
    tr = '┓',
    bl = '┗',
    br = '┛',
    t = '━',
    l = '┃',
    b = '━',
    r = '┃',
)
const BORDER_NONE = (
    tl = ' ',
    tr = ' ',
    bl = ' ',
    br = ' ',
    t = ' ',
    l = ' ',
    b = ' ',
    r = ' ',
)
const BORDER_BNONE = (
    tl = Char(0x2800),
    tr = Char(0x2800),
    bl = Char(0x2800),
    br = Char(0x2800),
    t = Char(0x2800),
    l = Char(0x2800),
    b = Char(0x2800),
    r = Char(0x2800),
)
const BORDER_DASHED = (
    tl = '┌',
    tr = '┐',
    bl = '└',
    br = '┘',
    t = '╌',
    l = '┊',
    b = '╌',
    r = '┊',
)
const BORDER_DOTTED = (
    tl = '⡤',
    tr = '⢤',
    bl = '⠓',
    br = '⠚',
    t = '⠤',
    l = '⡇',
    b = '⠒',
    r = '⢸',
)
const BORDER_ASCII = (
    tl = '+',
    tr = '+',
    bl = '+',
    br = '+',
    t = '-',
    l = '|',
    b = '-',
    r = '|',
)
const BORDERMAP = (
    solid   = BORDER_SOLID,
    corners = BORDER_CORNERS,
    barplot = BORDER_BARPLOT,
    bold    = BORDER_BOLD,
    none    = BORDER_NONE,
    bnone   = BORDER_BNONE,
    dashed  = BORDER_DASHED,
    dotted  = BORDER_DOTTED,
    ascii   = BORDER_ASCII,
)
const MARKERS = (
    circle    = '⚬',
    rect      = '▫',
    diamond   = '◇',
    hexagon   = '⬡',
    cross     = '✚',
    xcross    = '✖',
    utriangle = '△',
    dtriangle = '▽',
    rtriangle = '▷',
    ltriangle = '◁',
    pentagon  = '⬠',
    star4     = '✦',
    star5     = '★',
    star6     = '✶',
    star8     = '✴',
    vline     = '|',
    hline     = '―',
    (+)       = '+',
    (x)       = '⨯',
)
const SUPERSCRIPT = Dict(
    # '.' => '‧',  # U+2027: Hyphenation Point
    # '.' => '˙',  # U+02D9: Dot Above
    # '.' => '⸳',  # U+2E33: Raised Dot
    # '.' => '⋅',  # U+22C5: Dot Operator
    # '.' => '·',  # U+00B7: Middle Dot
    '.' => '⸱',  # U+2E31: Word Separator Middle Dot
    '-' => '⁻',
    '+' => '⁺',
    '0' => '⁰',
    '1' => '¹',
    '2' => '²',
    '3' => '³',
    '4' => '⁴',
    '5' => '⁵',
    '6' => '⁶',
    '7' => '⁷',
    '8' => '⁸',
    '9' => '⁹',
)
#! format: on

const COLOR_CYCLE = [:green, :blue, :red, :magenta, :yellow, :cyan]

const FSCALES = (identity = identity, ln = log, log2 = log2, log10 = log10)  # forward
const ISCALES = (identity = identity, ln = exp, log2 = exp2, log10 = exp10)  # inverse
const BASES = (identity = nothing, ln = "ℯ", log2 = "2", log10 = "10")

const MarkerType = Union{Symbol,Char,AbstractString}
const UserColorType = Union{Integer,Symbol,NTuple{3,Integer},Nothing}  # allowed color type
const JuliaColorType = Union{Symbol,Int}  # color type for printstyled (defined in base/util.jl)
const ColorType = Union{Nothing,UInt8}  # internal UnicodePlots color type

# standard terminals seem to respect a 4:3 aspect ratio
# unix.stackexchange.com/questions/148569/standard-terminal-font-aspect-ratio
# retrocomputing.stackexchange.com/questions/5629/why-did-80x25-become-the-text-monitor-standard
const ASPECT_RATIO = 4 / 3

# default display size for the default BrailleCanvas (which has aspect ratio = 2) ==> (40, 15)
const DEFAULT_HEIGHT = Ref(15)
const DEFAULT_WIDTH = Ref(2round(Int, ASPECT_RATIO * DEFAULT_HEIGHT[]))

function char_marker(marker::MarkerType)::Char
    if marker isa Symbol
        get(MARKERS, marker, MARKERS[:circle])
    else
        length(marker) == 1 || throw(error("`marker` keyword has a non unit length"))
        marker[1]
    end
end

iterable(obj::AbstractVector) = obj
iterable(obj) = Iterators.repeated(obj)

fscale(x, s::Symbol) = FSCALES[s](x)
iscale(x, s::Symbol) = ISCALES[s](x)

# support arbitrary scale functions
fscale(x, f::Function) = f(x)
iscale(x, f::Function) = f(x)

function transform_name(tr, basename = "")
    name = string(tr isa Union{Symbol,Function} ? tr : typeof(tr))  # typeof(...) for functors
    name == "identity" && return basename
    name = occursin("#", name) ? "custom" : name
    string(basename, " [", name, "]")
end

roundable(num::Number) = isinteger(num) & (typemin(Int) <= num < typemax(Int))
compact_repr(num::Number) = repr(num, context = :compact => true)

ceil_neg_log10(x) =
    roundable(-log10(x)) ? ceil(Integer, -log10(x)) : floor(Integer, -log10(x))
round_up_tick(x, m) = (
    x == 0 ? 0 :
    (x > 0 ? ceil(x, digits = ceil_neg_log10(m)) : -floor(-x, digits = ceil_neg_log10(m)))
)
round_down_tick(x, m) = (
    x == 0 ? 0 :
    (x > 0 ? floor(x, digits = ceil_neg_log10(m)) : -ceil(-x, digits = ceil_neg_log10(m)))
)
round_up_subtick(x, m) = (
    x == 0 ? 0 :
    (
        x > 0 ? ceil(x, digits = ceil_neg_log10(m) + 1) :
        -floor(-x, digits = ceil_neg_log10(m) + 1)
    )
)
round_down_subtick(x, m) = (
    x == 0 ? 0 :
    (
        x > 0 ? floor(x, digits = ceil_neg_log10(m) + 1) :
        -ceil(-x, digits = ceil_neg_log10(m) + 1)
    )
)
float_round_log10(x::F, m) where {F<:AbstractFloat} = (
    x == 0 ? F(0) :
    (
        x > 0 ? round(x, digits = ceil_neg_log10(m) + 1)::F :
        -round(-x, digits = ceil_neg_log10(m) + 1)::F
    )
)
float_round_log10(x::Integer, m) = float_round_log10(float(x), m)
float_round_log10(x) = x > 0 ? float_round_log10(x, x) : float_round_log10(x, -x)

function superscript(s::AbstractString)
    v = collect(s)
    for (i, k) in enumerate(v)
        v[i] = get(SUPERSCRIPT, k, k)
    end
    String(v)
end

function plotting_range(xmin, xmax)
    diff = xmax - xmin
    xmax = round_up_tick(xmax, diff)
    xmin = round_down_tick(xmin, diff)
    Float64(xmin), Float64(xmax)
end

function plotting_range_narrow(xmin, xmax)
    diff = xmax - xmin
    xmax = round_up_subtick(xmax, diff)
    xmin = round_down_subtick(xmin, diff)
    Float64(xmin), Float64(xmax)
end

extend_limits(vec, limits) = extend_limits(vec, limits, :identity)

function extend_limits(vec, limits, scale::Union{Symbol,Function})
    mi, ma = map(Float64, extrema(limits))
    if mi == 0 && ma == 0
        mi, ma = map(Float64, extrema(vec))
    end
    diff = ma - mi
    if diff == 0
        ma = mi + 1
        mi = mi - 1
    end
    if string(scale) != "identity"
        fscale(mi, scale), fscale(ma, scale)
    else
        all(iszero.(limits)) ? plotting_range_narrow(mi, ma) : (mi, ma)
    end
end

sort_by_keys(dict::Dict) = sort!(collect(dict), by = x -> x[1])

function sorted_keys_values(dict::Dict; k2s = true)
    if k2s  # check and force key type to be of AbstractString type if necessary
        kt, vt = eltype(dict).types
        if !(kt <: AbstractString)
            dict = Dict(string(k) => v for (k, v) in pairs(dict))
        end
    end
    keys_vals = sort_by_keys(dict)
    first.(keys_vals), last.(keys_vals)
end

print_color(color::UserColorType, io::IO, args...) =
    printstyled(io, string(args...); color = julia_color(color))

function crayon_256_color(color::UserColorType)::ColorType
    color in (:normal, :default, :nothing, nothing) && return nothing
    ansicolor = Crayons._parse_color(color)
    if ansicolor.style == Crayons.COLORS_16
        return Crayons.val(ansicolor) % 60
    elseif ansicolor.style == Crayons.COLORS_24BIT
        return Crayons.val(Crayons.to_256_colors(ansicolor))
    end
    Crayons.val(ansicolor)
end

julia_color(color::Integer)::JuliaColorType = Int(color)
julia_color(color::Nothing)::JuliaColorType = :normal
julia_color(color::Symbol)::JuliaColorType = color
julia_color(color)::JuliaColorType = julia_color(crayon_256_color(color))

@inline function set_color!(
    colors::Array{ColorType,2},
    x::Int,
    y::Int,
    color::ColorType,
    blend::Bool,
)
    if color === nothing || colors[x, y] === nothing || !blend
        colors[x, y] = color
    else
        colors[x, y] |= color
    end
    nothing
end

out_stream_size(out_stream::Union{Nothing,IO}) =
    out_stream === nothing ? (DEFAULT_WIDTH[], DEFAULT_HEIGHT[]) : displaysize(out_stream)
out_stream_width(out_stream::Union{Nothing,IO})::Int = out_stream_size(out_stream)[1]
out_stream_height(out_stream::Union{Nothing,IO})::Int = out_stream_size(out_stream)[2]

function _handle_deprecated_symb(symb, symbols)
    if symb === nothing
        symbols
    else
        Base.depwarn(
            "The keyword `symb` is deprecated in favor of `symbols`",
            :BarplotGraphics,
        )
        [symb]
    end
end


text:
function Base.reinterpret(::typeof(reshape), ::Type{T}, a::A) where {T,S,A<:CuArray{S}}
    isbitstype(T) || throwbits(S, T, T)
    isbitstype(S) || throwbits(S, T, S)
    if sizeof(S) == sizeof(T)
        N = ndims(a)
    elseif sizeof(S) > sizeof(T)
        rem(sizeof(S), sizeof(T)) == 0 || throwintmult(S, T)
        N = ndims(a) + 1
    else
        rem(sizeof(T), sizeof(S)) == 0 || throwintmult(S, T)
        N = ndims(a) - 1
        N > -1 || throwsize0(S, T, "larger")
        axes(a, 1) == Base.OneTo(sizeof(T) ÷ sizeof(S)) || throwsize1(a, T)
    end
    paxs = axes(a)
    new_axes = if sizeof(S) > sizeof(T)
        (Base.OneTo(div(sizeof(S), sizeof(T))), paxs...)
    elseif sizeof(S) < sizeof(T)
        Base.tail(paxs)
    else
        paxs
    end
    reshape(reinterpret(T, vec(a)), new_axes)
end

@noinline function throwintmult(S::Type, T::Type)
    throw(ArgumentError("`reinterpret(reshape, T, a)` requires that one of `sizeof(T)` (got $(sizeof(T))) and `sizeof(eltype(a))` (got $(sizeof(S))) be an integer multiple of the other"))
end
@noinline function throwsize1(a::AbstractArray, T::Type)
    throw(ArgumentError("`reinterpret(reshape, $T, a)` where `eltype(a)` is $(eltype(a)) requires that `axes(a, 1)` (got $(axes(a, 1))) be equal to 1:$(sizeof(T) ÷ sizeof(eltype(a))) (from the ratio of element sizes)"))
end
@noinline function throwbits(S::Type, T::Type, U::Type)
    throw(ArgumentError("cannot reinterpret `$(S)` as `$(T)`, type `$(U)` is not a bits type"))
end
@noinline function throwsize0(S::Type, T::Type, msg)
    throw(ArgumentError("cannot reinterpret a zero-dimensional `$(S)` array to `$(T)` which is of a $msg size"))
end


text:
export mipdecoHeuristic, heuristicKnapsack

function heuristicKnapsack(mc,gc,delta)
    mt = copy(mc)
    pred = 0
    cIndex = reverse(sortperm(abs.(gc)))
    wCount = 0
	for j=1:length(cIndex)
		ind = cIndex[Int(j)]
			  if  gc[ind] > 0 && mc[ind] == 1
				mt[ind] = 0
				pred -= gc[ind]
                wCount += 1
			  elseif  gc[ind] < 0 && mc[ind] == 0
				mt[ind] = 1
				pred += gc[ind]
                wCount += 1
			  end
            if (wCount >= delta)
                break
            end
	  end
    return mt,pred
end

function  mipdecoHeuristic(mc::BitArray{1}, pInv::InverseParam, pMis;
	                out::Int=2, getNeighborhood::Function=x->trues(size(x)),
					delta::Int=32, sigma::Float64=0.25,outFile::String="")

	mc          = copy(mc)
    iter        = 0
	TermFlag	= 0
	alpha       = pInv.alpha
	maxIter     = pInv.maxIter      #  Max. no. iterations.

	## Time statistics
	timeMIP = 0.0

	## evaluate function and derivatives
    Jc,gc,Dc,Fc,Rc,timeF,timeR = jInvObj(mc,pMis,pInv)
    rId = getNeighborhood(mc)
    gc .*= rId
    J0 = Jc;
	His = zeros(maxIter,3)

	### Initialization for trust region algorithm
	outStr = @sprintf("\n %4s\t%08s\t%08s\t%08s\t%08s\t%08s\t%08s\t%06s\n",
					  	"iter", "Fc", "Rc","alpha","Jc/J0","pred","ared","delta")
    if out>=2; print(outStr); end
	if !isempty(outFile)
		f = open(outFile, "w")
		write(f, outStr)
	end

      outStr = @sprintf("%3d\t%3.2e\t%3.2e\t%3.2e\t%3.2e\t%3.2e\t%3.2e\t%3d\n",
		         iter, Fc, Rc,alpha[1],Jc/J0,0,0,delta)
      if out>=2; print(outStr); end
      if !isempty(outFile)
		  write(f, outStr)
	  end

	while (delta > 0)
	  iter  += 1
	  toMIP = @elapsed begin
          mt,pred = heuristicKnapsack(mc,gc,delta)
      end
      pred *=-1
      timeMIP += toMIP

      Jt,gt,Dt,Ft,Rt,toF,toR = jInvObj(mt,pMis,pInv)
      timeF += toF
      timeR += toR
	  ared  = Jc - Jt  	# actual reduction
	  deltaO = delta

      outStr = @sprintf("%3d\t%3.2e\t%3.2e\t%3.2e\t%3.2e\t%3.2e\t%3.2e\t%3d\n",
		         iter, Fc, Rc,alpha[1],Jc/J0,pred,ared,delta)
      if out>=2; print(outStr); end
      if !isempty(outFile)
		  write(f, outStr)
	  end

        if ared > 0
	    mc .= mt
	    (Jc,Dc,Fc,Rc) = (Jt,Dt,Ft,Rt)

	    if ared >= sigma*pred;delta *= 2;end
	      gc = gt
          rId = getNeighborhood(mc)
          gc .*= rId
        else
         delta = floor(.5*delta)
	  end
	  His[iter,:] = [Jc Jt deltaO]
	  if iter >= maxIter
          TermFlag = 1
          break
      end
	end 	#while of trust region
	println("Trust region on entire space: MisFit, Regularizer, Total $(Fc) , $(Rc)	, $(Jc)")
	println("Time to compute Data term, Regularizer, and MIP: $(timeF), $(timeR), $(timeMIP) \n")
	if out>=1
		if TermFlag==0
			println("Trust region iterated $iter times with final objective $(Jc)." )
		else
			println("Trust region reached iteration limit with final objective $(Jc).")
		end
	end
    if !isempty(outFile)
	  close(f)
  end
	return mc,Dc,TermFlag,His
end

max_stars_repo_name: JuliaInv/ConvDiffMIPDECO

max_issues_repo_name: JuliaInv/ConvDiffMIPDECO

max_forks_repo_name: JuliaInv/ConvDiffMIPDECO


text:
#TODO: Document interface to support static size
# Define maybestaticsize, maybestaticlength, and maybestaticfirstindex

struct Static{N} <: Number end
Base.@pure Static(N) = Static{N}()
@inline Static(::Val{N}) where {N} = Static{N}()
@inline Base.Val(::Static{N}) where {N} = Val{N}()

Base.CartesianIndex(I::Tuple{<:Static,Vararg}) = CartesianVIndex(I)
Base.CartesianIndex(I::Tuple{<:Integer,<:Static,Vararg}) = CartesianVIndex(I)
Base.CartesianIndex(I::Tuple{<:Integer,<:Integer,<:Static,Vararg}) = CartesianVIndex(I)
Base.CartesianIndex(I::Tuple{<:Integer,<:Integer,<:Integer,<:Static,Vararg}) = CartesianVIndex(I)
Base.CartesianIndex(I::Tuple{<:Integer,<:Integer,<:Integer,<:Integer,<:Static,Vararg}) = CartesianVIndex(I)

abstract type AbstractStaticUnitRange <: AbstractUnitRange{Int} end
# rank 2, but 3 unknowns; 5 types can express all different posibilities
# 1: all unknown, UnitRange; 2: all three known (only two must be made explicit):
struct StaticUnitRange{L,U} <: AbstractStaticUnitRange end
Base.@pure StaticUnitRange(L,U) = StaticUnitRange{L,U}() # Do I use this definition anywhere?
# Then each specifying one of the three parameters
struct StaticLowerUnitRange{L} <: AbstractStaticUnitRange
    U::Int
end
struct StaticUpperUnitRange{U} <: AbstractStaticUnitRange
    L::Int
end
struct StaticLengthUnitRange{N} <: AbstractStaticUnitRange
    L::Int
end

@inline Base.first(::StaticUnitRange{L}) where {L} = L
@inline Base.first(::StaticLowerUnitRange{L}) where {L} = L
@inline Base.first(r::StaticUpperUnitRange) = r.L
@inline Base.first(r::StaticLengthUnitRange) = r.L
@inline Base.last(::StaticUnitRange{L,U}) where {L,U} = U
@inline Base.last(r::StaticLowerUnitRange) = r.U
@inline Base.last(::StaticUpperUnitRange{U}) where {U} = U
@inline Base.last(r::StaticLengthUnitRange{N}) where {N} = vadd(r.L, N - 1)

Base.show(io::IO, r::AbstractStaticUnitRange) = print(io, "$(first(r)):$(last(r))")

@inline Base.iterate(x::AbstractStaticUnitRange) = (i = unwrap(first(x)); (i,i))

@inline Base.:(:)(::Static{L}, ::Static{U}) where {L,U} = StaticUnitRange{L,U}()
@inline Base.:(:)(::Static{L}, U::Int) where {L} = StaticLowerUnitRange{L}(U)
@inline Base.:(:)(L::Int, ::Static{U}) where {U} = StaticUpperUnitRange{U}(L)
@inline unwrap(x) = x
@inline unwrap(::Val{N}) where {N} = N
@inline unwrap(::Type{Val{N}}) where {N} = N
@inline unwrap(::Static{N}) where {N} = N
@inline unwrap(::Type{Static{N}}) where {N} = N

@inline StaticLowerUnitRange{L}(::Static{U}) where {L,U} = StaticUnitRange{L,U}()
@inline StaticUpperUnitRange{U}(::Static{L}) where {L,U} = StaticUnitRange{L,U}()
@generated StaticLengthUnitRange{N}(::Static{L}) where {L,N} = StaticUnitRange{L,L+N-1}()

@inline maybestaticsize(A, ::Val{I}) where {I} = size(A, I)
@inline maybestaticsize(A::AbstractArray{<:Any,0}, ::Val{1:2}) = (1, 1)
@inline maybestaticsize(A::AbstractVector, ::Val{1:2}) = (length(A), 1)
@inline maybestaticsize(A::AbstractMatrix, ::Val{1:2}) = size(A)
@inline maybestaticsize(A::AbstractArray, ::Val{1:2}) = (size(A,1),size(A,2))
@inline maybestaticlength(A) = length(A)

@inline maybestaticfirstindex(A::AbstractArray, ::Val{I}) where {I} = firstindex(A, I)
@inline maybestaticfirstindex(A::Array, ::Val{I}) where {I} = Static{1}()

# @inline maybestaticeachindex(A::AbstractArray) = maybestaticrange(eachindex(A))

@inline maybestaticrange(r::Base.OneTo) = StaticLowerUnitRange{1}(last(r))
@inline maybestaticrange(r::UnitRange) = r
@inline maybestaticrange(r::AbstractStaticUnitRange) = r
@inline maybestaticrange(r) = maybestaticfirst(r):maybestaticlast(r)
# @inline maybestaticaxis(A::AbstractArray, ::Val{I}) where {I} = maybestaticfirstindex(A, Val{I}()):maybestaticsize(A, Val{I}())

@inline maybestaticsize(::NTuple{N}, ::Val{1}) where {N} = Static{N}() # should we assert that i == 1?
@inline maybestaticlength(::NTuple{N}) where {N} = Static{N}()
@inline maybestaticsize(::Adjoint{T,V}, ::Val{1}) where {T,V<:AbstractVector{T}} = Static{1}()
@inline maybestaticsize(::Transpose{T,V}, ::Val{1}) where {T,V<:AbstractVector{T}} = Static{1}()
@inline maybestaticlength(B::Adjoint) = maybestaticlength(parent(B))
@inline maybestaticlength(B::Transpose) = maybestaticlength(parent(B))
@inline maybestaticsize(B::Adjoint{T,A}, ::Val{1}) where {T,A<:AbstractMatrix{T}} = maybestaticsize(parent(B), Val{2}())
@inline maybestaticsize(B::Adjoint{T,A}, ::Val{2}) where {T,A<:AbstractMatrix{T}} = maybestaticsize(parent(B), Val{1}())
@inline maybestaticsize(B::Transpose{T,A}, ::Val{1}) where {T,A<:AbstractMatrix{T}} = maybestaticsize(parent(B), Val{2}())
@inline maybestaticsize(B::Transpose{T,A}, ::Val{2}) where {T,A<:AbstractMatrix{T}} = maybestaticsize(parent(B), Val{1}())

@inline Base.iszero(::Static{0}) = true
@inline Base.iszero(::Static) = false
@generated function Base.divrem(N::Integer, ::Static{L}) where {L}
    if ispow2(L)
        quote
            $(Expr(:meta,:inline))
            d = N >>> $(intlog2(L))
            r = N & $(L-1)
            d, r
        end
    else
        quote
            $(Expr(:meta,:inline))
            divrem(N, $L)
        end
    end
end
@generated Base.divrem(::Static{N}, ::Static{D}) where {N,D} = divrem(N, D)

@inline vadd(::Static{N}, i::Number) where {N} = vadd(N, i)
@inline vadd(i::Number, ::Static{N}) where {N} = vadd(i, N)
@inline vsub(::Static{N}, i::Number) where {N} = vsub(N, i)
@inline vsub(i::Number, ::Static{N}) where {N} = vsub(i, N)

@inline vadd(::Static{M}, ::Static{N}) where {M,N} = vadd(M, N)
@inline vmul(::Static{N}, i::Number) where {N} = vmul(N, i)
@inline vmul(i::Number, ::Static{N}) where {N} = vmul(i, N)
@inline vmul(::Static{N}, i) where {N} = vmul(N, i)
@inline vmul(i, ::Static{N}) where {N} = vmul(i, N)
@inline vmul(::Static{M}, ::Static{N}) where {M,N} = vmul(M, N)
@inline vsub(::Static{M}, ::Static{N}) where {M,N} = vsub(M, N)
@inline Base.:+(::Static{N}, i) where {N} = vadd(N, i)
@inline Base.:+(i, ::Static{N}) where {N} = vadd(N, i)
@inline Base.:+(::Static{M}, ::Static{N}) where {M,N} = vadd(M, N)
@inline Base.:*(::Static{N}, i) where {N} = vmul(N, i)
@inline Base.:*(i, ::Static{N}) where {N} = vmul(N, i)
@inline Base.:*(::Static{M}, ::Static{N}) where {M,N} = vmul(M, N)
@inline Base.:-(::Static{N}, i) where {N} = vsub(N, i)
@inline Base.:-(i, ::Static{N}) where {N} = vsub(i, N)
@inline Base.:-(::Static{M}, ::Static{N}) where {M,N} = vsub(M, N)
@inline Base.checked_add(::Static{N}, i) where {N} = Base.checked_add(N, i)
@inline Base.checked_add(i, ::Static{N}) where {N} = Base.checked_add(i, N)
@generated Base.checked_add(::Static{M}, ::Static{N}) where {M,N} = Static{Base.checked_add(M, N)}()
@inline Base.checked_sub(::Static{N}, i) where {N} = Base.checked_sub(N, i)
@inline Base.checked_sub(i, ::Static{N}) where {N} = Base.checked_sub(i, N)
@generated Base.checked_sub(::Static{M}, ::Static{N}) where {M,N} = Static{Base.checked_sub(M, N)}()
@inline Base.:>>(::Static{N}, i) where {N} = N >> i
@inline Base.:>>(i, ::Static{N}) where {N} = i >> N
@inline Base.:>>(::Static{M}, ::Static{N}) where {M,N} = M >> N
@inline Base.:<<(::Static{N}, i) where {N} = N << i
@inline Base.:<<(i, ::Static{N}) where {N} = i << N
@inline Base.:<<(::Static{M}, ::Static{N}) where {M,N} = M << N
@inline Base.:>>>(::Static{N}, i) where {N} = N >>> i
@inline Base.:>>>(i, ::Static{N}) where {N} = i >>> N
@inline Base.:>>>(::Static{M}, ::Static{N}) where {M,N} = M >>> N
@inline Base.:&(::Static{N}, i) where {N} = N & i
@inline Base.:&(i, ::Static{N}) where {N} = N & i
@inline Base.:&(::Static{M}, ::Static{N}) where {M,N} = M & N
@inline Base.:>(::Static{N}, i) where {N} = N > i
@inline Base.:>(i, ::Static{N}) where {N} = i > N
@inline Base.:>(::Static{M}, ::Static{N}) where {M,N} = M > N
@inline Base.:<(::Static{N}, i) where {N} = N < i
@inline Base.:<(i, ::Static{N}) where {N} = i < N
@inline Base.:<(::Static{M}, ::Static{N}) where {M,N} = M < N
@inline Base.:(==)(::Static{M}, i) where {M} = M == i
@inline Base.:(==)(i, ::Static{M}) where {M} = M == i
@inline Base.:(==)(::Static{M}, ::Static{N}) where {M,N} = false
@inline Base.:(==)(::Static{M}, ::Static{M}) where {M} = true

@inline vadd(::Tuple{}, ::Tuple{}) = tuple()
@inline vadd(a::Tuple{I1,Vararg}, b::Tuple{}) where {I1} = a
@inline vadd(a::Tuple{}, b::Tuple{I2,Vararg}) where {I2} = b
@inline vadd(a::Tuple{I1}, b::Tuple{I2}) where {I1,I2} = (vadd(a[1],b[1]),)
@inline vadd(a::Tuple{I1,I3}, b::Tuple{I2,I4}) where {I1,I2,I3,I4} = (vadd(a[1],b[1]),vadd(a[2],b[2]),)
@inline vadd(a::Tuple{I1,Vararg}, b::Tuple{I2,Vararg}) where {I1,I2} = (vadd(a[1],b[1]),vadd(Base.tail(a),Base.tail(b))...)

function static_promote(i, j)
    i == j || throw("$i ≠ $j")
    i
end
function static_promote(::Static{M}, i) where {M}
    M == i || throw("$M ≠ $i")
    Static{M}()
end
function static_promote(i, ::Static{M}) where {M}
    M == i || throw("$M ≠ $i")
    Static{M}()
end
static_promote(::Static{M}, ::Static{N}) where {M, N} = throw("$M ≠ $N")
static_promote(::Static{M}, ::Static{M}) where {M} = Static{M}()

@generated staticp1(::Static{N}) where {N} = Static{N+1}()
@inline staticp1(N) = vadd(N, one(N))
@inline staticp1(i::Tuple{}) = tuple()
@inline staticp1(i::Tuple{I}) where {I} = @inbounds (staticp1(i[1]),)
@inline staticp1(i::Tuple{I1,I2}) where {I1,I2} = @inbounds (staticp1(i[1]), staticp1(i[2]))
@inline staticp1(i::Tuple{I1,I2,I3,Vararg}) where {I1,I2,I3} = @inbounds (staticp1(i[1]), staticp1(Base.tail(i))...)
@generated staticm1(::Static{N}) where {N} = Static{N-1}()
@inline staticm1(N) = vsub(N, one(N))
@inline staticm1(i::Tuple{}) = tuple()
@inline staticm1(i::Tuple{I}) where {I} = @inbounds (staticm1(i[1]),)
@inline staticm1(i::Tuple{I1,I2}) where {I1,I2} = @inbounds (staticm1(i[1]), staticm1(i[2]))
@inline staticm1(i::Tuple{I1,I2,I3,Vararg}) where {I1,I2,I3} = @inbounds (staticm1(i[1]), staticm1(Base.tail(i))...)
@generated staticmul(::Type{T}, ::Static{N}) where {T,N} = Static{sizeof(T) * N}()
@generated staticmul(::Type{T}, ::Val{N}) where {T,N} = Val{sizeof(T) * N}()
@inline staticmul(::Type{T}, N) where {T} = vmul(N, sizeof(T))
@inline staticmul(::Type{T}, i::Tuple{}) where {T} = tuple()
@inline staticmul(::Type{T}, i::Tuple{I}) where {T,I} = @inbounds (vmul(i[1], sizeof(T)),)
@inline staticmul(::Type{T}, i::Tuple{I1,I2}) where {T,I1,I2} = @inbounds (vmul(sizeof(T), i[1]), vmul(sizeof(T), i[2]))
@inline staticmul(::Type{T}, i::Tuple{I1,I2,I3,Vararg}) where {T,I1,I2,I3} = @inbounds (vmul(sizeof(T), i[1]), staticmul(T, Base.tail(i))...)
@inline Base.ntuple(f::F, ::Static{N}) where {F,N} = ntuple(f, Val{N}())

struct LazyP1{T}
    data::T
end
@inline staticm1(d::LazyP1) = d.data

@inline maybestaticfirst(A) = first(A)
@inline maybestaticfirst(::StaticUnitRange{L}) where {L} = Static{L}()
@inline maybestaticfirst(::StaticLowerUnitRange{L}) where {L} = Static{L}()
@inline maybestaticfirst(::Base.OneTo) = Static{1}()

@inline _maybestaticfirst(A::Tuple{}) = tuple()
@inline _maybestaticfirst(A::Tuple{I}) where {I} = (maybestaticfirst(@inbounds(A[1])),)
@inline _maybestaticfirst(A::Tuple{I1,I2}) where {I1,I2} = @inbounds (maybestaticfirst(A[1]), maybestaticfirst(A[2]))
@inline _maybestaticfirst(A::Tuple{I1,I2,I3,Vararg}) where {I1,I2,I3} = (maybestaticfirst(@inbounds A[1]), _maybestaticfirst(Base.tail(A))...)
@inline maybestaticfirst(A::CartesianIndices) = CartesianVIndex(_maybestaticfirst(A.indices))

@inline maybestaticlast(A) = last(A)
@inline maybestaticlast(::StaticUnitRange{L,U}) where {L,U} = Static{U}()
@inline maybestaticlast(::StaticUpperUnitRange{U}) where {U} = Static{U}()
@inline _maybestaticlast(A::Tuple{}) = tuple()
@inline _maybestaticlast(A::Tuple{I}) where {I} = (maybestaticlast(@inbounds(A[1])),)
@inline _maybestaticlast(A::Tuple{I1,I2}) where {I1,I2} = @inbounds (maybestaticlast(A[1]), maybestaticlast(A[2]))
@inline _maybestaticlast(A::Tuple{I1,I2,I3,Vararg}) where {I1,I2,I3} = (maybestaticlast(@inbounds A[1]), _maybestaticlast(Base.tail(A))...)
@inline maybestaticlast(A::CartesianIndices) = CartesianVIndex(_maybestaticlast(A.indices))

# Static 0
const Zero = Static{0}
@inline vsub(::Zero, i) = vsub(i)
@inline vsub(i, ::Zero) = i
@inline vsub(::Zero, i::Number) = vsub(i)
@inline vsub(i::Number, ::Zero) = i
@inline vsub(::Zero, ::Zero) = Zero()
@inline vadd(::Zero, ::Zero) = Zero()
@inline vadd(::Zero, a) = a
@inline vadd(a, ::Zero) = a
@inline vadd(::Zero, a::Number) = a
@inline vadd(a::Number, ::Zero) = a
@inline vmul(::Zero, ::Any) = Zero()
@inline vmul(::Any, ::Zero) = Zero()
@inline vmul(::Zero, ::Number) = Zero()
@inline vmul(::Number, ::Zero) = Zero()
for T ∈ [:Int,:SVec]
    @eval @inline vadd(::Zero, a::$T) = a
    @eval @inline vadd(a::$T, ::Zero) = a
    @eval @inline vsub(::Zero, a::$T) = vsub(a)
    @eval @inline vsub(a::$T, ::Zero) = a
    @eval @inline vmul(::Zero, ::$T) = Zero()
    @eval @inline vmul(::$T, ::Zero) = Zero()
end
@inline vadd(::Zero, a::_MM) = a
@inline vadd(a::_MM, ::Zero) = a
@inline vload(ptr::Ptr, ::Zero) = vload(ptr)
@inline vload(ptr::Ptr, ::_MM{W,Zero}) where {W} = vload(Val{W}(), ptr)
@inline vload(ptr::Ptr, ::_MM{W,Zero}, m::Mask) where {W} = vload(Val{W}(), ptr, m.u)
@inline vstore!(ptr::Ptr{T}, v::T, ::Zero) where {T} = vstore!(ptr, v)
@inline vnoaliasstore!(ptr::Ptr{T}, v::T, ::Zero) where {T} = vnoaliasstore!(ptr, v)
@inline vstore!(ptr::Ptr{T}, v, ::Zero) where {T} = vstore!(ptr, convert(T,v))
@inline vnoaliasstore!(ptr::Ptr{T}, v, ::Zero) where {T} = vnoaliasstore!(ptr, convert(T,v))
@inline vstore!(ptr::Ptr{T}, v::Integer, ::Zero) where {T <: Integer} = vstore!(ptr, v % T)
@inline vnoaliasstore!(ptr::Ptr{T}, v::Integer, ::Zero) where {T <: Integer} = vnoaliasstore!(ptr, v % T)
# @inline vstore!(ptr::Ptr{T}, v::T, ::Zero, m::Mask) where {T} = vstore!(ptr, v, m.u)
# @inline vnoaliasstore!(ptr::Ptr{T}, v::T, ::Zero, m::Mask) where {T} = vnoaliasstore!(ptr, v, m.u)
for V ∈ [:(NTuple{W,Core.VecElement{T}}), :(SVec{W,T})]
    @eval @inline vstore!(ptr::Ptr{T}, v::$V, ::Zero) where {W,T} = vstore!(ptr, v)
    @eval @inline vstore!(ptr::Ptr{T}, v::$V, ::_MM{W,Zero}) where {W,T} = vstore!(ptr, v)
    @eval @inline vnoaliasstore!(ptr::Ptr{T}, v::$V, ::Zero) where {W,T} = vnoaliasstore!(ptr, v)
    @eval @inline vnoaliasstore!(ptr::Ptr{T}, v::$V, ::_MM{W,Zero}) where {W,T} = vnoaliasstore!(ptr, v)
    for M ∈ [:(Mask{W}), :Unsigned]
        @eval @inline vstore!(ptr::Ptr{T}, v::$V, ::Zero, m::$M) where {W,T} = vstore!(ptr, v, m)
        @eval @inline vstore!(ptr::Ptr{T}, v::$V, ::_MM{W,Zero}, m::$M) where {W,T} = vstore!(ptr, v, m)
        @eval @inline vnoaliasstore!(ptr::Ptr{T}, v::$V, ::Zero, m::$M) where {W,T} = vnoaliasstore!(ptr, v, m)
        @eval @inline vnoaliasstore!(ptr::Ptr{T}, v::$V, ::_MM{W,Zero}, m::$M) where {W,T} = vnoaliasstore!(ptr, v, m)
    end
end

const One = Static{1}
@inline vmul(::One, a) = a
@inline vmul(a, ::One) = a
@inline vmul(::One, a::Number) = a
@inline vmul(a::Number, ::One) = a
@inline vmul(::One, ::One) = One()
@inline vmul(::One, ::Zero) = Zero()
@inline vmul(::Zero, ::One) = Zero()

for T ∈ [:Int,:SVec]
    @eval @inline vmul(::One, a::$T) = a
    @eval @inline vmul(a::$T, ::One) = a
end

struct LazyStaticMul{N,T}
    data::T
end

# @inline vmul(::Static{N}, i::I) where {N,I<:Union{Integer,_MM}} = LazyStaticMul{N,I}(i)
# @inline vmul(i::I, ::Static{N}) where {N,I<:Union{Integer,_MM}} = LazyStaticMul{N,I}(i)
@inline vmul(::Static{N}, i::I) where {N,I<:Integer} = LazyStaticMul{N,I}(i)
@inline vmul(i::I, ::Static{N}) where {N,I<:Integer} = LazyStaticMul{N,I}(i)
@inline vadd(a::LazyStaticMul{N}, b) where {N} = vadd(vmul(N, a.data), b)
@inline vadd(b, a::LazyStaticMul{N}) where {N} = vadd(b, vmul(N, a.data))
@inline vadd(a::LazyStaticMul, ::Static{0}) = a
@inline vadd(::Static{0}, a::LazyStaticMul) = a
@inline vadd(a::LazyStaticMul{N}, ::Static{M}) where {M,N} = vadd(vmul(N, a.data), M)
@inline vadd(::Static{M}, a::LazyStaticMul{N}) where {M,N} = vadd(vmul(N, a.data), M)
@inline vadd(a::LazyStaticMul{N}, b::LazyStaticMul{M}) where {M,N} = vadd(vmul(N, a.data), vmul(M, b.data))
@inline vadd(a::LazyStaticMul{N,T}, b::LazyStaticMul{N,T}) where {N,T} = LazyStaticMul{N,T}(vadd(a.data, b.data))

@inline vadd(a::LazyStaticMul{N,I}, ::_MM{W,Static{0}}) where {W,N,I<:Integer} = LazyStaticMul{N,_MM{W,I}}(_MM{W}(a.data))
@inline vadd(::_MM{W,Static{0}}, a::LazyStaticMul{N,I}) where {W,N,I<:Integer} = LazyStaticMul{N,_MM{W,I}}(_MM{W}(a.data))

# @inline extract_data(a::LazyStaticMul{N}) where {N} = vmul(a.data, N)
# @inline vload(ptr, lsm::LazyStaticMul{N}) where {N} = vload(gep(ptr, lsm.data, Val{N}()))
# @inline vload(ptr, lsm::LazyStaticMul{N,_MM{W,I}}) where {N,W,I} = vload(Val{W}(), gep(ptr, lsm.data, Val{N}()))
# @inline vload(ptr, lsm::LazyStaticMul{N,SVec{W,I}}) where {N,W,I} = vload(Val{W}(), gep(ptr, lsm.data, Val{N}()))
# @inline vload(ptr::Ptr{T}, lsm::LazyStaticMul{N,_Vec{W,I}}) where {T,N,W,I} = vload(_Vec{W,T}, gep(ptr, lsm.data, Val{N}()))

@inline vload(ptr, lsm::LazyStaticMul{N}   ) where {N} = vload(ptr, vmul(N, lsm.data)   )
@inline vload(ptr, lsm::LazyStaticMul{N}, m) where {N} = vload(ptr, vmul(N, lsm.data), m)
@inline vstore!(ptr, v, lsm::LazyStaticMul{N}   ) where {N} = vstore!(ptr, v, vmul(N, lsm.data)   )
@inline vstore!(ptr, v, lsm::LazyStaticMul{N}, m) where {N} = vstore!(ptr, v, vmul(N, lsm.data), m)
@inline vnoaliasstore!(ptr, v, lsm::LazyStaticMul{N}   ) where {N} = vnoaliasstore!(ptr, v, vmul(N, lsm.data)   )
@inline vnoaliasstore!(ptr, v, lsm::LazyStaticMul{N}, m) where {N} =! vnoaliasstore(ptr, v, vmul(N, lsm.data), m)

@inline vload(ptr, lsm::LazyStaticMul{N,<:_MM{W}}   ) where {N,W} = vload(Val{W}(), ptr, vmul(N, lsm.data.i)   )
@inline vload(ptr, lsm::LazyStaticMul{N,<:_MM{W}}, m) where {N,W} = vload(Val{W}(), ptr, vmul(N, lsm.data.i), m)
@inline vstore!(ptr, v, lsm::LazyStaticMul{N,<:_MM}   ) where {N} = vstore!(ptr, v, vmul(N, lsm.data.i)   )
@inline vstore!(ptr, v, lsm::LazyStaticMul{N,<:_MM}, m) where {N} = vstore!(ptr, v, vmul(N, lsm.data.i), m)
@inline vnoaliasstore!(ptr, v, lsm::LazyStaticMul{N,<:_MM}   ) where {N} = vnoaliasstore!(ptr, v, vmul(N, lsm.data.i)   )
@inline vnoaliasstore!(ptr, v, lsm::LazyStaticMul{N,<:_MM}, m) where {N} = vnoaliasstore!(ptr, v, vmul(N, lsm.data.i), m)
for N ∈ [1,2,4,8,6,10]
    @eval @inline vload(ptr, lsm::LazyStaticMul{$N}) = vload(gep(ptr, lsm.data, Val{$N}()))
    @eval @inline vload(ptr, lsm::LazyStaticMul{$N}, m) = vload(gep(ptr, lsm.data, Val{$N}()), m)
    @eval @inline vstore!(ptr, v, lsm::LazyStaticMul{$N}) = vstore!(gep(ptr, lsm.data, Val{$N}()), v)
    @eval @inline vstore!(ptr, v, lsm::LazyStaticMul{$N}, m) = vstore!(gep(ptr, lsm.data, Val{$N}()), v, m)

    @eval @inline vload(ptr, lsm::LazyStaticMul{$N,_MM{W,I}}) where {W,I}  = vload(Val{W}(), gep(ptr, lsm.data.i, Val{$N}()))
    @eval @inline vload(ptr, lsm::LazyStaticMul{$N,_MM{W,I}}, m) where {W,I} = vload(Val{W}(), gep(ptr, lsm.data.i, Val{$N}()), m)
    @eval @inline vstore!(ptr, v, lsm::LazyStaticMul{$N,_MM{W,I}}) where {W,I} = vstore!(gep(ptr, lsm.data.i, Val{$N}()), v)
    @eval @inline vstore!(ptr, v, lsm::LazyStaticMul{$N,_MM{W,I}}, m) where {W,I} = vstore!(gep(ptr, lsm.data.i, Val{$N}()), v, m)
end


text:
module WaveletInformation

using Reexport
@reexport using MLBase
using ScikitLearn
using Wavelets
using Combinatorics
using MultivariateStats
using GaussianMixtures

include("wavelet_decode.jl")
include("MI.jl")
include("wpca.jl")
include("scikit_compat.jl")

end