text:
using GenParticleFilters, Gen, Test, Logging

@gen (static) function line_step(t::Int, x::Float64, slope::Float64)
    x = x + 1
    outlier ~ bernoulli(0.1)
    y ~ normal(x * slope, outlier ? 10.0 : 1.0)
    return x
end

line_unfold = Unfold(line_step)

@gen (static) function line_model(n::Int)
    slope ~ uniform_discrete(-2, 2)
    line ~ line_unfold(n, 0, slope)
    return line
end

@load_generated_functions()

generate_line(n::Int, slope::Float64=0.) =
    choicemap([(:line => i => :y, i*slope) for i in 1:n]...)

@testset "Particle initialization" begin

@testset "Initialize with default proposal" begin
state = pf_initialize(line_model, (0,), choicemap(), 100)
@test all(-2 <= tr[:slope] <= 2 for tr in get_traces(state))
state = pf_initialize(line_model, (1,), generate_line(1), 100)
@test all(tr[:line => 1 => :y] == 0 for tr in get_traces(state))
state = pf_initialize(line_model, (10,), generate_line(10), 100)
@test all(tr[:line => 10 => :y] == 0 for tr in get_traces(state))
end

@gen line_propose(s) =
    slope ~ uniform_discrete(0, 0)
@gen outlier_propose(idxs) =
    [{:line => i => :outlier} ~ bernoulli(0.0) for i in idxs]

@testset "Initialize with custom proposal" begin
state = pf_initialize(line_model, (0,), choicemap(), line_propose, (0,), 100)
@test all(tr[:slope] == 0 for tr in get_traces(state))
state = pf_initialize(line_model, (1,), choicemap(), outlier_propose, ([1],), 100)
@test all(tr[:line => 1 => :outlier] == false for tr in get_traces(state))
state = pf_initialize(line_model, (10,), choicemap(), outlier_propose, ([10],), 100)
@test all(tr[:line => 10 => :outlier] == false for tr in get_traces(state))
end

@testset "Initialize for dynamic model structure" begin
state = pf_initialize(line_model, (0,), choicemap(), 10, true)
@test typeof(state) == ParticleFilterState{Trace}
state = pf_initialize(line_model, (0,), choicemap(), line_propose, (0,), 10, true)
@test typeof(state) == ParticleFilterState{Trace}
end

end

@testset "Particle update" begin

@testset "Update with default proposal" begin
state = pf_initialize(line_model, (0,), choicemap(), 100)
state = pf_update!(state, (10,), (UnknownChange(),), generate_line(10))
@test all(tr[:line => 10 => :y] == 0 for tr in get_traces(state))
@test all(w != 0 for w in get_log_weights(state))
end

@gen outlier_propose(tr, idxs) =
    [{:line => i => :outlier} ~ bernoulli(0.0) for i in idxs]
@gen outlier_reverse(tr, idxs) =
    [{:line => i => :outlier} ~ bernoulli(0.1) for i in idxs]

@testset "Update with custom proposal" begin
state = pf_initialize(line_model, (0,), choicemap(), 100)
state = pf_update!(state, (10,), (UnknownChange(),), generate_line(10),
                   outlier_propose, (1:10,))
@test all(tr[:line => 10 => :y] == 0 for tr in get_traces(state))
@test all(tr[:line => 10 => :outlier] == false for tr in get_traces(state))
@test all(w != 0 for w in get_log_weights(state))
end

@testset "Update with custom forward and backward proposals" begin
state = pf_initialize(line_model, (10,), generate_line(10), 100)
state = pf_update!(state, (10,), (UnknownChange(),), choicemap(),
                   outlier_propose, (1:10,), outlier_reverse, (1:10,))
@test all(tr[:line => 10 => :y] == 0 for tr in get_traces(state))
@test all(tr[:line => 10 => :outlier] == false for tr in get_traces(state))
@test all(w != 0 for w in get_log_weights(state))
end

@testset "Update with custom proposal and trace transform" begin
@gen proposal(tr, idxs) = [{(:outlier, i)} ~ bernoulli(0.0) for i in idxs]
@transform remap (p_old, _) to (p_new, _) begin
    idxs = get_args(p_old)[2]
    for i in idxs
        @copy(p_old[(:outlier, i)], p_new[:line => i => :outlier])
    end
end
state = pf_initialize(line_model, (0,), choicemap(), 100)
state = pf_update!(state, (10,), (UnknownChange(),), generate_line(10),
                   proposal, (1:10,), remap)
@test all(tr[:line => 10 => :y] == 0 for tr in get_traces(state))
@test all(tr[:line => 10 => :outlier] == false for tr in get_traces(state))
@test all(w != 0 for w in get_log_weights(state))
end

@testset "Update with custom bidirectional proposals and trace transform" begin
@gen function fwd_kernel(tr, idxs)
    if (flip ~ bernoulli(0.5))
        [{(:outlier, i)} ~ bernoulli(0.0) for i in idxs]
    else
        slope ~ uniform_discrete(0, 0)
    end
end
@gen function bwd_kernel(tr, idxs)
    if (flip ~ bernoulli(0.5))
        [{(:outlier, i)} ~ bernoulli(0.1) for i in idxs]
    else
        slope ~ uniform_discrete(-2, 2)
    end
end
@transform line_transform (p_old, q_fwd) to (p_new, q_bwd) begin
    flip = @read(q_fwd[:flip], :discrete)
    @write(q_bwd[:flip], flip, :discrete)
    if flip
        idxs = get_args(q_fwd)[2]
        for i in idxs
            @copy(q_fwd[(:outlier, i)], p_new[:line => i => :outlier])
            @copy(p_old[:line => i => :outlier], q_bwd[(:outlier, i)])
        end
    else
        @copy(q_fwd[:slope], p_new[:slope])
        @copy(p_old[:slope], q_bwd[:slope])
    end
end
is_involution!(line_transform)
state = pf_initialize(line_model, (5,), generate_line(5), 100)
state = pf_update!(state, (10,), (UnknownChange(),), choicemap(),
                   fwd_kernel, (1:10,), bwd_kernel, (1:5,),
                   line_transform, true)
@test all(tr[:line => 5 => :y] == 0 for tr in get_traces(state))
@test all(tr[:line => 5 => :outlier] == false || tr[:slope] == 0
          for tr in get_traces(state))
@test all(w != 0 for w in get_log_weights(state))
end

@testset "Update with different proposals per view" begin
state = pf_initialize(line_model, (0,), choicemap(), 100)
substate = pf_update!(state[1:50], (10,), (UnknownChange(),), generate_line(10))
@test all(tr[:line => 10 => :y] == 0 for tr in get_traces(substate))
@test all(w != 0 for w in get_log_weights(substate))
substate = pf_update!(state[51:end], (10,), (UnknownChange(),),
                      generate_line(10), outlier_propose, (10,))
@test all(tr[:line => 10 => :y] == 0 for tr in get_traces(substate))
@test all(tr[:line => 10 => :outlier] == false for tr in get_traces(substate))
@test all(w != 0 for w in get_log_weights(state))
end

end

@testset "Particle resampling" begin

@testset "Multinomial resampling" begin
# Test that new traces have the correct parents
state = pf_initialize(line_model, (10,), generate_line(10), 100)
old_traces = get_traces(state)
state = pf_multinomial_resample!(state)
new_traces = get_traces(state)
@test new_traces == old_traces[state.parents]
end

@testset "Residual resampling" begin
# Test that no resampling occurs if all weights are equal
state = pf_initialize(line_model, (0,), choicemap(), 100)
old_traces = get_traces(state)
state = pf_residual_resample!(state)
new_traces = get_traces(state)
@test new_traces == old_traces

# Test that at least the minimum number of copies are resampled
state = pf_initialize(line_model, (10,), generate_line(10), 100)
old_traces = get_traces(state)
weights = get_norm_weights(state)
min_copies = floor.(Int, weights * 100)
state = pf_residual_resample!(state)
new_traces = get_traces(state)
copies = [sum([t1 == t2 for t1 in new_traces]) for t2 in old_traces]
@test new_traces == old_traces[state.parents]
@test all(copies .>= min_copies)

# Same test but with a custom priority function
p_fn = w -> w / 2
state = pf_initialize(line_model, (10,), generate_line(10), 100)
old_traces = get_traces(state)
log_priorities = p_fn.(get_log_weights(state))
weights = exp.(log_priorities .- logsumexp(log_priorities))
min_copies = floor.(Int, weights * 100)
state = pf_residual_resample!(state; priority_fn=p_fn)
new_traces = get_traces(state)
copies = [sum([t1 == t2 for t1 in new_traces]) for t2 in old_traces]
@test new_traces == old_traces[state.parents]
@test all(copies .>= min_copies)
end

@testset "Stratified resampling" begin
# Test that no resampling occurs if all weights are equal
state = pf_initialize(line_model, (0,), choicemap(), 100)
old_traces = get_traces(state)
state = pf_stratified_resample!(state)
new_traces = get_traces(state)
@test new_traces == old_traces

# Test that the highest weight particle has the right number of copies
state = pf_initialize(line_model, (10,), generate_line(10), 100)
old_traces = get_traces(state)
weights = get_norm_weights(state)
max_weight, max_idx = findmax(weights)
min_copies = floor.(Int, max_weight * 100)
state = pf_stratified_resample!(state; sort_particles=true)
new_traces = get_traces(state)
copies = sum([tr == old_traces[max_idx] for tr in new_traces])
@test new_traces == old_traces[state.parents]
@test copies >= min_copies

# Same test but with a custom priority function
p_fn = w -> w / 2
state = pf_initialize(line_model, (10,), generate_line(10), 100)
old_traces = get_traces(state)
log_priorities = p_fn.(get_log_weights(state))
weights = exp.(log_priorities .- logsumexp(log_priorities))
max_weight, max_idx = findmax(weights)
min_copies = floor.(Int, max_weight * 100)
state = pf_stratified_resample!(state; sort_particles=true, priority_fn=p_fn)
new_traces = get_traces(state)
copies = sum([tr == old_traces[max_idx] for tr in new_traces])
@test new_traces == old_traces[state.parents]
@test copies >= min_copies
end

end

@testset "Particle rejuvenation" begin

@testset "Move-reweight kernel" begin
out_addr = :line => 1 => :outlier
observations = choicemap((:line => 1 => :y, 0))
trace, _ = generate(line_model, (1,), observations)
slope, outlier = trace[:slope], trace[out_addr]

# Test selection variant
expected_w(out_old, out_new, slope) =
    logpdf(normal, 0, slope, out_new ? 10. : 1.) -
    logpdf(normal, 0, slope, out_old ? 10. : 1.)
trs_ws = [move_reweight(trace, select(out_addr)) for i in 1:100]
@test all(w ≈ expected_w(outlier, tr[out_addr], slope) for (tr, w) in trs_ws)

# Test proposal variant
expected_w(out_old, out_new, slope) =
    logpdf(bernoulli, out_new, 0.1) - logpdf(bernoulli, out_old, 0.1) +
    logpdf(normal, 0, slope, out_new ? 10. : 1.) -
    logpdf(normal, 0, slope, out_old ? 10. : 1.) +
    (out_old == out_new ? 0.0 :
        logpdf(bernoulli, out_old, 0.9) - logpdf(bernoulli, out_old, 0.1))
@gen outlier_propose(tr, idx) = {:line => idx => :outlier} ~ bernoulli(0.9)
trs_ws = [move_reweight(trace, outlier_propose, (1,)) for i in 1:100]
@test all(w ≈ expected_w(outlier, tr[out_addr], slope) for (tr, w) in trs_ws)
end

@testset "Move-accept rejuvenation" begin
# Log which particles were rejuvenated
buffer = IOBuffer()
logger = SimpleLogger(buffer, Logging.Debug)
state = pf_initialize(line_model, (10,), generate_line(10, 1.), 100)
old_traces = get_traces(state)
with_logger(logger) do
    pf_move_accept!(state, metropolis_hastings, (select(:slope),), 1)
end

# Extract acceptances from debug log
lines = split(String(take!(buffer)), "\n")
lines = filter(s -> occursin("Accepted: ", s), lines)
accepts = [match(r".*Accepted: (\w+).*", l).captures[1] for l in lines]
accepts = parse.(Bool, accepts)

# Check that only traces that were accepted are rejuvenated
new_traces = get_traces(state)
@test all(a ? t1 !== t2 : t1 === t2
          for (a, t1, t2) in zip(accepts, old_traces, new_traces))
end

@testset "Move-reweight rejuvenation" begin
# Log which particles were rejuvenated
buffer = IOBuffer()
logger = SimpleLogger(buffer, Logging.Debug)
state = pf_initialize(line_model, (10,), generate_line(10, 1.), 100)
old_weights = copy(get_log_weights(state))
with_logger(logger) do
    pf_move_reweight!(state, move_reweight, (select(:slope),), 1)
end
new_weights = copy(get_log_weights(state))

# Extract relative weights from debug log
lines = split(String(take!(buffer)), "\n")
lines = filter(s -> occursin("Rel. Weight: ", s), lines)
rel_weights = [match(r".*Rel\. Weight: (.+)\s*", l).captures[1] for l in lines]
rel_weights = parse.(Float64, rel_weights)

# Check that weights are adjusted accordingly
@test all(isapprox.(new_weights, old_weights .+ rel_weights; atol=1e-3))
end

@testset "Rejuvenation on separate views" begin
# Log which particles were rejuvenated
buffer = IOBuffer()
logger = SimpleLogger(buffer, Logging.Debug)
state = pf_initialize(line_model, (10,), generate_line(10, 1.), 100)
old_traces = get_traces(state)[1:50]
old_weights = get_log_weights(state)[51:end]

with_logger(logger) do
    pf_move_accept!(state[1:50], metropolis_hastings, (select(:slope),), 1)
    pf_move_reweight!(state[51:end], move_reweight, (select(:slope),), 1)
end

# Extract acceptances and relative weights from debug log
lines = split(String(take!(buffer)), "\n")
a_lines = filter(s -> occursin("Accepted: ", s), lines)
accepts = [match(r".*Accepted: (\w+).*", l).captures[1] for l in a_lines]
accepts = parse.(Bool, accepts)
r_lines = filter(s -> occursin("Rel. Weight: ", s), lines)
rel_weights = [match(r".*Rel\. Weight: (.+)\s*", l).captures[1] for l in r_lines]
rel_weights = parse.(Float64, rel_weights)

# Check that only traces that were accepted are rejuvenated
new_traces = get_traces(state)[1:50]
@test all(a ? t1 !== t2 : t1 === t2
          for (a, t1, t2) in zip(accepts, old_traces, new_traces))
# Check that weights are adjusted accordingly
new_weights = get_log_weights(state)[51:end]
@test all(isapprox.(new_weights, old_weights .+ rel_weights; atol=1e-3))

end

end

@testset "Utility functions" begin
@gen model() = x ~ normal(0, 0)
state = pf_initialize(model, (), choicemap(), 100)

@test mean(state, :x) == 0
@test mean(state) == 0
@test var(state, :x) == 0
@test var(state) == 0

weights = get_norm_weights(state)
@test sum(weights) ≈ 1.0
ess = get_ess(state)
@test ess ≈ sum(weights)^2 / sum(weights .^ 2)
end


text:
using PackageCompiler

using Cbc
using Clp
using Geodesy
using JSON
using JSONSchema
using JuMP
using MathOptInterface
using ProgressBars

pkg = [:Cbc, :Clp, :Geodesy, :JSON, :JSONSchema, :JuMP, :MathOptInterface, :ProgressBars]

@info "Building system image..."
create_sysimage(pkg, sysimage_path = "build/sysimage.so")

ext: jl

max_stars_repo_name: ANL-CEEESA/RELOG

max_issues_repo_name: ANL-CEEESA/RELOG

max_forks_repo_name: ANL-CEEESA/RELOG


text:
module SpringCollab2020SafeRespawnCrumble

using ..Ahorn, Maple

@mapdef Entity "SpringCollab2020/safeRespawnCrumble" SafeRespawnCrumble(x::Integer, y::Integer, width::Integer=Maple.defaultBlockWidth)

const placements = Ahorn.PlacementDict(
    "Safe Respawn Crumble (Spring Collab 2020)" => Ahorn.EntityPlacement(
        SafeRespawnCrumble,
        "rectangle"
    )
)

Ahorn.minimumSize(entity::SafeRespawnCrumble) = 8, 0
Ahorn.resizable(entity::SafeRespawnCrumble) = true, false

function Ahorn.selection(entity::SafeRespawnCrumble)
    x, y = Ahorn.position(entity)
    width = Int(get(entity.data, "width", 8))

    return Ahorn.Rectangle(x, y, width, 8)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::SafeRespawnCrumble, room::Maple.Room)
    texture = "objects/SpringCollab2020/safeRespawnCrumble/tile"

    # Values need to be system specific integer
    x = Int(get(entity.data, "x", 0))
    y = Int(get(entity.data, "y", 0))

    width = Int(get(entity.data, "width", 8))
    tilesWidth = div(width, 8)

    Ahorn.Cairo.save(ctx)

    Ahorn.rectangle(ctx, 0, 0, width, 8)
    Ahorn.clip(ctx)

    for i in 0:ceil(Int, tilesWidth)
        Ahorn.drawImage(ctx, texture, 8 * i, 0, 0, 0, 8, 8)
    end

    Ahorn.restore(ctx)
end

end


text:
import Base.LinAlg: axpy!

export svdl

"""
Matrix of the form
 d1          a_1
    d2       a_2
      ...    ...
         d_l a_l
             d_l+1 e_l+1
                   ...  e_k-1
                        d_k
"""
type BrokenArrowBidiagonal{T} <: AbstractMatrix{T}
    dv::Vector{T}
    av::Vector{T}
    ev::Vector{T}
end

Base.size(B::BrokenArrowBidiagonal) = (n=length(B.dv); (n, n))
function Base.size(B::BrokenArrowBidiagonal, n::Int)
    if n==1 || n==2
        return length(B.dv)
    else
        throw(ArgumentError("invalid dimension $n"))
    end
end

function Base.getindex{T}(B::BrokenArrowBidiagonal{T}, i::Int, j::Int)
    n = size(B, 1)
    k = length(B.av)
    if !(1 ≤ i ≤ n && 1 ≤ j ≤ n)
        throw(BoundsError())
    end

    if i == j
        return B.dv[i]
    elseif i ≤ k && j == k+1
        return B.av[i]
    elseif i > k && j == i+1
        return B.ev[i-k]
    else
        return zero(T)
    end
end

function Base.full{T}(B::BrokenArrowBidiagonal{T})
    n = size(B, 1)
    k = length(B.av)
    M = zeros(T, n, n)
    for i=1:n
        M[i, i] = B.dv[i]
    end
    for i=1:k
        M[i,k+1] = B.av[i]
    end
    for i=k+1:n-1
        M[i, i+1] = B.ev[i-k]
    end
    M
end

Base.svdfact(B::BrokenArrowBidiagonal) = svdfact(full(B)) #XXX This can be much faster

"""
Partial factorization object which is an approximation of a matrix

    A ≈ P * [B 0; 0 β] * Q
"""
type PartialFactorization{T, Tr} <: Factorization{T}
    P :: Matrix{T}
    Q :: Matrix{T}
    B :: AbstractMatrix{Tr}
    β :: Tr
end

"""
Compute some singular values (and optionally vectors) using Golub-Kahan-Lanczos
bidiagonalization \cite{Golub1965} with thick restarting \cite{Wu2000}.

# Inputs

- `A` : The matrix or matrixlike object whose singular values are desired
- `l` : The number of singular values requested.
        Default: 6

# Keyword inputs

- `v0` : The starting guess vector in the domain of `A`.
         The length of `q` should be the number of columns in `A`.
         Default: A random unit vector.
- `k` : The maximum number of Lanczos vectors to compute before restarting.
        Default: `2*l`
- `j` : The number of vectors to keep at the end of the restart.
        Default: `l`. We don't recommend j < l.
- `maxiter`: Maximum number of iterations to run
             Default: `minimum(size(A))`
- `verbose`: Whether to print information at each iteration
             Default: false
- `tol`    : Maximum absolute error in each desired singular value.
             Default: `√eps()`
- `reltol` : Maximum error in each desired singular value relative to the
             estimated norm of the input matrix.  Default: `√eps()`
- `restart`: Which restarting algorithm to use. Valid choices are:
             - `:ritz`: Thick restart with Ritz values [Wu2000]. Default
             - `:harmonic`: Restart with harmonic Ritz values [Baglama2005]
- `doplot` : Plot a history of the Ritz value convergence. Requires the
             [UnicodePlots.jl](https://github.com/Evizero/UnicodePlots.jl.git)
             package to be installed.
             Default: `false`
- `vecs`   : Return singular vectors also.
             - `:both`: Both left and right singular vectors are returned
             - `:left`: Only the left singular vectors are returned
             - `:right`: Only the right singular vectors are returned
             - `:none`: No singular vectors are returned (Default)
- `dolock` : If `true`, locks converged Ritz values, removing them from the
             Krylov subspace being searched in the next macroiteration.
             Default: `false`

# Output

- `Σ`: A list of the desired singular values if `vecs == :none` (the default),
    otherwise returns an `SVD` object with the desired singular vectors filled in
- `L`: The computed partial factorizations of A
- `Bs`: The history of the computed projected matrices
- `ritzvalhist`: Ritz values computed at each iteration
- `convhist`: Convergence data

# Implementation notes

The implementation of thick restarting follows closely that of SLEPc as
described in [Hernandez2008]. Thick restarting can be turned off by setting `k
= maxiter`, but most of the time this is not desirable.

The singular vectors are computed directly by forming the Ritz vectors from the
product of the Lanczos vectors `L.P`/`L.Q` and the singular vectors of `L.B`.
Additional accuracy in the singular triples can be obtained using inverse
iteration.

# References

```bibtex
@article{Golub1965,
    author = {Golub, G. and Kahan, W.},
    doi = {10.1137/0702016},
    journal = {Journal of the Society for Industrial and Applied Mathematics
        Series B Numerical Analysis},
    volume = 2,
    number = 2,
    pages = {205--224},
    title = {Calculating the Singular Values and Pseudo-Inverse of a Matrix},
    year = 1965
}

@article{Wu2000,
    author = {Wu, Kesheng and Simon, Horst},
    journal = {SIAM Journal on Matrix Analysis and Applications},
    number = 2,
    pages = {602--616},
    title = {Thick-Restart {L}anczos Method for Large Symmetric Eigenvalue Problems},
    volume = 22,
    year = 2000
}

@article{Baglama2005,
    author = {Baglama, James and Reichel, Lothar},
    doi = {10.1137/04060593X},
    journal = {SIAM Journal on Scientific Computing},
    number = 1,
    pages = {19--42},
    title = {Augmented Implicitly Restarted {L}anczos Bidiagonalization Methods},
    volume = 27,
    year = 2005
}

@article{Hernandez2008,
    author = {Hern\'{a}ndez, Vicente and Rom\'{a}n, Jos\'{e} E and Tom\'{a}s,
    Andr\'{e}s},
    journal = {Electronic Transactions on Numerical Analysis},
    pages = {68--85},
    title = {A Robust and Efficient Parallel {SVD} Solver based on Restarted
        {L}anczos Bidiagonalization},
    url = {http://etna.mcs.kent.edu/volumes/2001-2010/vol31/abstract.php?vol=31\&pages=68-85},
    volume = 31,
    year = 2008
}
```

"""
function svdl(A, l::Int=min(6, size(A,1)); tol::Real=√eps(), k::Int=2l,
    maxiter::Int=minimum(size(A)), method::Symbol=:ritz, kwargs...
    )
    history = ConvergenceHistory()
    history[:tol] = tol
    reserve!(BitArray, history,:conv, maxiter)
    reserve!(history,[:ritz,:resnorm], maxiter, l)
    Bs_type = (method == :ritz) ? BrokenArrowBidiagonal : UpperTriangular
    reserve!(Bs_type, history,:Bs, maxiter)
    reserve!(history,:betas, maxiter)
    X, L = svdl_method!(history, A, l; k=k, tol=tol, maxiter=maxiter, method=method, kwargs...)
    X, L, history
end

#Method implementation
function svdl_method!(log::ConvergenceHistory, A, l::Int=min(6, size(A,1)); k::Int=2l,
    j::Int=l, v0::AbstractVector = Vector{eltype(A)}(randn(size(A, 2))) |> x->scale!(x, inv(norm(x))),
    maxiter::Int=minimum(size(A)), tol::Real=√eps(), reltol::Real=√eps(),
    verbose::Bool=false, method::Symbol=:ritz, doplot::Bool=false, vecs=:none, dolock::Bool=false)
#function svdl(A, l::Int=6, k::Int=2l,
#    j::Int=l, v0::AbstractVector = Vector{eltype(A)}(randn(size(A, 2))) |> x->scale!(x, inv(norm(x))),
#    maxiter::Int=minimum(size(A)), tol::Real=√eps(), reltol::Real=√eps(),
#    verbose::Bool=false, method::Symbol=:ritz, doplot::Bool=false, vecs=:none, dolock::Bool=false)

    T0 = time_ns()
    @assert k>1
    L = build(log, A, v0, k)

    iter = 0
    local F
    for iter in 1:maxiter
        nextiter!(log)

        #@assert size(L.B) == (k, k)
        F = svdfact(L.B)::LinAlg.SVD{eltype(v0), typeof(real(one(eltype(v0)))), Matrix{eltype(v0)}}
        iter==1 && @assert eltype(F)==eltype(v0)
        if method == :ritz
            thickrestart!(A, L, F, j)
        elseif method == :harmonic
            harmonicrestart!(A, L, F, j)
        else
            throw(ArgumentError("Unknown restart method $method"))
        end
        extend!(log, A, L, k)
        if verbose
            elapsedtime = round((time_ns()-T0)*1e-9, 3)
            info("Iteration $iter: $elapsedtime seconds")
        end

        conv = isconverged(L, F, l, tol, reltol, log, verbose)

        push!(log, :conv, conv)
        push!(log, :ritz, F[:S][1:k])
        push!(log, :Bs, deepcopy(L.B))
        push!(log, :betas, L.β)

        #Lock
        if method == :ritz && dolock
            for i in eachindex(conv)
                if conv[i]
                    L.B.av[i] = 0
                end
            end
        end
        all(conv) && (setconv(log, true); break)
    end
    shrink!(log)

    #Compute singular vectors as necessary and return them in the output
    values = F[:S][1:l]
    m, n = size(A)

    leftvecs = if vecs == :left || vecs == :both
        L.P*view(F[:U], :, 1:l)
    else
        zeros(eltype(v0), m, 0)
    end

    rightvecs = if vecs == :right || vecs == :both
        (view(L.Q, :, 1:size(L.Q,2)-1)*view(F[:V], :, 1:l))'
    else
        zeros(eltype(v0), 0, n)
    end

    if vecs == :none
        values, L
    else
        LinAlg.SVD(leftvecs, values, rightvecs), L
    end
end

"""
Determine if any singular values in a partial factorization have converged.

# Inputs

- `L` : A `PartialFactorization` computed by an iterative method such as `svdl`
- `F` : A `SVD` factorization computed for `L.B`
- `k` : Number of singular values to check
- `tol`: Absolute tolerance for a Ritz value to be considered converged
- `reltol`: Relative tolerance for a Ritz value to be considered converged
- `verbose`: If `true`, prints out all the results of convergence tests.
             Default: `false`.

# Implementation notes

This convergence test routine uses a variety of different tests.

1. The crudest estimate of the error bound is a simple error bound which dates
back at least to Wilkinson's classic book [Wilkinson1965:Ch.3 §53 p.170]

```bibtex
@book{Wilkinson1965,
    address = {Oxford, UK},
    author = {J H Wilkinson},
    publisher = {Oxford},
    title = {The Algebraic Eigenvalue Problem},
    year = 1965
}
```

2. When the Ritz values become sufficiently well-separated, more refined estimates
can be derived from the Rayleigh-Ritz properties of the Krylov process, as
described in [Wilkinson1965:Ch.3 §54-55 p.173, Yamamoto1980, Ortega1990]

```bibtex
@book{Ortega1990,
    address = {Philadelphia, PA},
    author = {Ortega, James M},
    doi = {10.1137/1.9781611971323},
    edition = {2},
    publisher = {SIAM},
    series = {Classics in Applied Mathematics},
    title = {Numerical Analysis: A Second Course},
    url = {http://epubs.siam.org/doi/book/10.1137/1.9781611971323},
    year = 1990
}

@article{Yamamoto1980,
    author = {Yamamoto, Tetsuro},
    doi = {10.1007/BF01396059},
    journal = {Numerische Mathematik},
    number = 2,
    pages = {189--199},
    title = {Error bounds for computed eigenvalues and eigenvectors},
    volume = {34},
    year = 1980
}

@article{Geurts1982,
author = {Geurts, A J},
doi = {10.1007/BF01399313},
journal = {Numerische Mathematik},
month = feb,
number = {1},
pages = {85--96},
title = {A contribution to the theory of condition},
volume = {39},
year = {1982}
}

@article{Deif1989,
author = {Deif, A.},
doi = {10.1007/BF01396348},
journal = {Numerische Mathematik},
month = jun,
number = {6},
pages = {625--626},
title = {A relative backward perturbation theorem for the eigenvalue problem},
volume = {56},
year = {1989}
}

```
"""
function isconverged(L::PartialFactorization, F::Base.LinAlg.SVD, k::Int, tol::Real,
        reltol::Real, log::ConvergenceHistory, verbose::Bool=false
        )

    @assert tol ≥ 0

    σ = F[:S][1:k]
    Δσ= L.β*abs(F[:U][end, 1:k])

    #Best available eigenvalue bounds
    δσ = copy(Δσ)

    #Update better error bounds from Rayleigh-Ritz
    #In theory these only apply if we know all the values
    #But so long as the gap between the converged Ritz values and all the
    #others is larger than the smallest empirical spectral gap d then we're fine.
    #Can also get some eigenvector statistics!!
    if k > 1
        d = Inf
        for i in eachindex(σ), j=1:i-1
            d = min(d, abs(σ[i] - σ[j]))
        end
        verbose && println("Smallest empirical spectral gap: ", d)
        #Chatelein 1993 - normwise backward error associated with approximate invariant subspace
        #Use the largest singular (Ritz) value to estimate the 2-norm of the matrix
        verbose && println("Normwise backward error associated with subspace: ", L.β/σ[1])
        for i in eachindex(Δσ)
            α = Δσ[i]

            #Simple error bound
            if 2α ≤ d
                #Rayleigh-Ritz bounds
                x = abs(α/(d-α)*√(1+(α/(d-α))^2))
                verbose && println("Rayleigh-Ritz error bound on eigenvector: $x")
                #δσ[i] = min(δσ[i], x)

                y = α^2/d #[Wilkinson:Ch.3 Appendix (4), p.188]
                verbose && println("Rayleigh-Ritz error bound on eigenvalue: $y")
                δσ[i] = min(δσ[i], y)
            end

	    verbose && println("Ritz value ", i, ": ", σ[i], " ± ", signif(δσ[i], 3))

            #Estimate of the normwise backward error [Deif 1989]
            ##Use the largest singular (Ritz) value to estimate the 2-norm of the matrix
            verbose && println("Normwise backward error estimate: ", α/σ[1])

            #Geurts, 1982 - componentwise backward error also known
            #A. J. Geurts, (1982), A contribution to the theory of condition,
            #Numer.  Math., #39, 85-96.
        end
    end

    #Estimate condition number and see if two-sided reorthog is needed
    if verbose && (F[:S][1]/F[:S][end]) > 1/√eps()
        warn("Two-sided reorthogonalization should be used but is not implemented")
    end

    push!(log, :resnorm, δσ[1:k])
    conv = (δσ[1:k] .< max(tol, reltol*σ[1]))::BitVector
end

#Hernandez2008
function build{T}(log::ConvergenceHistory, A, q::AbstractVector{T}, k::Int)
    m, n = size(A)
    Tr = typeof(real(one(T)))
    β = norm(q)
    scale!(q, inv(β))
    p = A*q
    α = convert(Tr, norm(p))
    scale!(p, inv(α))
    extend!(log, A, PartialFactorization(
        reshape(p, m, 1), reshape(q, n, 1), Bidiagonal([α], Tr[], true), β
        ), k)
end

"""
Thick restart (with ordinary Ritz values)

# Reference

[Hernandez2008]
"""
function thickrestart!{T,Tr}(A, L::PartialFactorization{T,Tr},
        F::Base.LinAlg.SVD{Tr,Tr}, l::Int)

    k = size(F[:V], 1)
    m, n = size(A)
    #@assert size(L.P) == (m, k)
    #@assert size(L.Q) == (n, k+1)

    Q = view(L.Q, :,1:k)*view(F[:V], :,1:l)
    L.Q = [Q view(L.Q, :, k+1)]
    #Be pedantic about ensuring normalization
    #L.Q = qr(L.Q)[1]
    #@assert all([norm(L.Q[:,i]) ≈ 1 for i=1:size(L.Q,2)])

    f = A*view(L.Q, :, l+1)
    ρ = L.β * reshape(F[:U][end, 1:l], l)
    L.P = view(L.P, :, 1:k)*view(F[:U], :, 1:l)

    #@assert ρ[i] ≈ f⋅L.P[:, i]
    f -= L.P*ρ
    α = convert(Tr, norm(f))
    scale!(f, inv(α))
    L.P = [L.P f]

    g = A'f - α*L.Q[:, end]
    L.β = β = norm(g)
    L.B = BrokenArrowBidiagonal([F[:S][1:l]; α], ρ, typeof(β)[])
    #@assert size(L.P, 2) == size(L.Q, 2) == size(L.B, 2)
    L
end

"""
Thick restart with harmonic Ritz values

# Reference

[Baglama2005] - note that they have P and Q swapped relative to our notation,
which follows that of [Hernandez2008]
"""
function harmonicrestart!{T,Tr}(A, L::PartialFactorization{T,Tr},
        F::Base.LinAlg.SVD{Tr,Tr}, k::Int)

    m = size(L.B, 1)::Int
    #@assert size(L.P,2)==m==size(L.Q,2)-1

    F0 = F# svdfact(L.B)
    ρ = L.β*F0[:U][end,:] #Residuals of singular values

    #Construct broken arrow matrix
    if VERSION < v"0.5.0-"
        BA = [diagm(F0[:S]) ρ']
    else #slicing behavior changed in 0.5
        BA = [diagm(F0[:S]) ρ]
    end
    F2 = svdfact!(BA, thin=false)

    #Take k largest triplets
    Σ = (F2[:S]::Vector{Tr})[1:k]
    U = F0[:U]*view(F2[:U],:,1:k)
    M = eye(T, m+1)
    M[1:m, 1:m] = F0[:V]::Matrix{T}
    M = M * F2[:V]
    Mend = M[end, 1:k]
    #Compute scaled residual from the harmonic Ritz problem
    r0 = zeros(Tr, m)
    r0[end] = 1
    isa(L.B, Bidiagonal) && @assert L.B.isupper
    r = try
        #(L.B\r0)
        A_ldiv_B!(L.B, r0)
    catch exc
        if isa(exc, LinAlg.LAPACKException) ||
                isa(exc, LinAlg.SingularException) #B\r is singular

            pinv(full(L.B))*r0
        else rethrow(exc) end
    end::Vector{Tr}
    scale!(r, L.β)
    M::Matrix{T} = view(M,1:m, :) + r*view(M,m+1:m+1,:)

    M2 = zeros(T, m+1, k+1)
    M2[1:m, 1:k] = M[:,1:k]
    M2[1:m, k+1] = -r
    M2[m+1, k+1] = 1
    Q, R = qr(M2)

    Q = L.Q*Q
    P = L.P*view(U,:,1:k)

    if VERSION < v"0.5.0-"
        R = view(R,1:k+1,1:k) + view(R,:,k+1)*Mend
    else
        R = view(R,1:k+1,1:k) + view(R,:,k+1)*Mend'
    end

    f = A*view(Q,:,k+1)
    f -= P*(P'f)
    α = convert(Tr, norm(f))
    scale!(f, inv(α))
    P = [P f]
    B = UpperTriangular{Tr,Matrix{Tr}}([(Diagonal(Σ)*triu(R')); zeros(Tr,1,k) α])
    #@assert size(P, 2) == size(B, 1) == size(Q, 2)
    g = A'f
    q = view(Q,:,k+1)
    #g-= (g⋅q)*q
    axpy!(-(g⋅q), q, g)
    β = convert(Tr, norm(g))
    scale!(g, inv(β))
    #@assert size(P, 2) == size(Q, 2) == size(B, 2)
    L.P = P
    L.Q = Q
    L.B = B
    L.β = β
    L
end

"""
Extend a PartialFactorization L using GKL bidiagonalization with k extra pairs
of Lanczos vectors

# Input

- `A`: matrix or linear map generating the Lanczos vectors

- `L`: `PartialFactorization` object

- `orthleft::Bool`: whether or not to orthogonalize left Lanczos vectors

- `orthright::Bool`: whether or not to orthogonalize right Lanczos vectors

- `α::Real`: criterion for doing a second reorthogonalization. Default: 1/√2

# Implementation notes

The implementation mostly follows the description in [Simon2000,Hernandez2008].

The reorthogonalization method used is using double classical Gram-Schmidt
full reorthogonalization. As explained in the numerical analysis literature by
Kahan, Golub, Rutishauser, and others in the 1970s, double classical
Gram-Schmidt reorthogonalization always suffices to keep vectors orthogonal to
within machine precision. As described in [Bjorck2015], α is a threshold
for determinining when the second orthogonalization is necessary. -log10(α) is
the number of (decimal) digits lost due to cancellation. Common choices are
α=0.1 [Rutishauser] and α=1/√2 [Daniel1976] (our default).

In most situations it suffices to orthogonalize either the left vectors or the
right vectors, except when the matrix norm exceeds `1/√eps(eltype(A))`, in
which case it will be necessary to orthogonalize both sets of vectors. See
[Simon2000].

@book{Bjorck2015,
author = {Bj{\"{o}}rck, {\AA}ke},
doi = {10.1007/978-3-319-05089-8},
publisher = {Springer},
series = {Texts in Applied Mathematics},
title = {Numerical Methods in Matrix Computations},
year = {2015}
}

@article{Simon2000,
        author = {Simon, Horst D. and Zha, Hongyuan},
        doi = {10.1137/S1064827597327309},
        journal = {SIAM Journal on Scientific Computing},
        number = {6},
        pages = {2257--2274},
        title = {Low-Rank Matrix Approximation Using the {Lanczos} Bidiagonalization Process with Applications},
        volume = {21},
        year = {2000}
}

@article{Daniel1976,
author = {Daniel, J. W. and Gragg, W. B. and Kaufman, L. and Stewart, G. W.},
doi = {10.1090/S0025-5718-1976-0431641-8},
journal = {Mathematics of Computation},
number = {136},
pages = {772--795},
title = {Reorthogonalization and stable algorithms for updating the {Gram-Schmidt QR} factorization},
volume = {30},
year = {1976}
}

"""
function extend!{T,Tr}(
    log::ConvergenceHistory, A, L::PartialFactorization{T, Tr}, k::Int,
    orthleft::Bool=false, orthright::Bool=true, α::Real = 1/√2
    )

    l = size(L.B, 2)::Int-1
    p = L.P[:,l+1]

    m, n = size(A)
    q = zeros(T, n)
    #@assert l+1 == size(L.B, 1) == size(L.B, 2) == size(L.P, 2)

    if isa(L.B, UpperTriangular) #Cannot be appended
        B = zeros(Tr, k, k)
        B[1:size(L.B,1), 1:size(L.B,2)] = L.B
        L.B = UpperTriangular(B)
        #@assert size(L.B) == (k, k)
    end

    β = L.β
    for j=l+1:k
        log.mtvps+=1
        Ac_mul_B!(q, A, p) #q = A'p

        if orthright #Orthogonalize right Lanczos vector
            #Do double classical Gram-Schmidt reorthogonalization
            oldqnorm = norm(q)
            q -= L.Q*(L.Q'q)
            if norm(q) ≤ α * oldqnorm
                q -= L.Q*(L.Q'q)
            end
        end

        β = norm(q)
        scale!(q, inv(β))

        L.Q = [L.Q q]
        j==k && break

        log.mvps+=1
        #p = A*q - β*p
        A_mul_B!(p, A, q)
        Base.LinAlg.axpy!(-β, view(L.P, :, j), p)

        if orthleft #Orthogonalize left Lanczos vector
            #Do double classical Gram-Schmidt reorthogonalization
            oldpnorm = norm(p)
            p -= L.P*(L.P'p)
            if norm(p) ≤ α * oldpnorm
                p -= L.P*(L.P'p)
            end
        end

        α = norm(p)
        scale!(p, inv(α))
        if isa(L.B, Bidiagonal) || isa(L.B, BrokenArrowBidiagonal)
            push!(L.B.dv, α)
            push!(L.B.ev, β)
        else
            L.B[j+1, j+1] = α
            L.B[j  , j+1] = β
        end
        L.P = [L.P p]
    end
    L.β = β
    L
end

let
    B = BrokenArrowBidiagonal([1, 2, 3], [1, 2], Int[])
    @assert full(B) == [1 0 1; 0 2 2; 0 0 3]
end


text:
_settled = ContinuousCallback(
    (u, t, integrator) -> u[2],
    (integrator) -> terminate!(integrator),
    initialize=(c,u,t,integrator) -> u[2] != 0 || terminate!(integrator)
)

"""
    transform(prob::CauchyProblem) -> DifferentialEquations.ODEProblem

Transform `prob` into an ODE problem in terms of ϕ. The ODE problem is set up to terminate automatically
(`.retcode === :Terminated`) when the steady state is reached.

See also: [`DifferentialEquations`](https://diffeq.sciml.ai/stable/)
"""
function transform(prob::CauchyProblem)
    u0 = @SVector [prob.b, prob.d_dϕb]
    ϕb = float(prob.ϕb)
    ODEProblem(transform(prob.eq), u0, (ϕb, typemax(ϕb)), callback=_settled)
end

monotonicity(odeprob::ODEProblem)::Int = sign(odeprob.u0[2])

function _integrate(odeprob::ODEProblem; limit=nothing)
    ode_maxiters = 1000

    if !isnothing(limit)
        past_limit = DiscreteCallback(
            let direction=monotonicity(odeprob)
                (u, t, integrator) -> direction*u[1] > direction*limit
            end,
            (integrator) -> terminate!(integrator)
        )
        return OrdinaryDiffEq.solve(odeprob, RadauIIA5(),
                                    callback=past_limit, verbose=false, maxiters=ode_maxiters)
    end

    return OrdinaryDiffEq.solve(odeprob, RadauIIA5(), verbose=false, maxiters=ode_maxiters)
end

function _integrate(prob::CauchyProblem; limit=nothing)
    _integrate(transform(prob), limit=limit)
end

function solve(prob::CauchyProblem)

    @argcheck isindomain(prob.eq, prob.b) DomainError(prob.b, "prob.b not valid for the given equation")

    odesol = _integrate(prob)

    if odesol.retcode !== :Terminated
        throw(SolvingError("could not find a solution to the problem"))
    end

    return Solution(odesol, prob.eq, 0)
end


text:
# Autogenerated wrapper script for utf8proc_jll for armv6l-linux-musleabihf
export libutf8proc

JLLWrappers.@generate_wrapper_header("utf8proc")
JLLWrappers.@declare_library_product(libutf8proc, "libutf8proc.so.2")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libutf8proc,
        "lib/libutf8proc.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module SemialgebraicSets

using Random

using MultivariatePolynomials
const MP = MultivariatePolynomials

const APL = AbstractPolynomialLike

export AbstractSemialgebraicSet, AbstractBasicSemialgebraicSet, AbstractAlgebraicSet
export FullSpace, AlgebraicSet, BasicSemialgebraicSet, addequality!, addinequality!

# Semialgebraic set described by polynomials with coefficients in T
abstract type AbstractSemialgebraicSet end

abstract type AbstractBasicSemialgebraicSet <: AbstractSemialgebraicSet end
abstract type AbstractAlgebraicSet <: AbstractBasicSemialgebraicSet end

addinequality!(S::AbstractAlgebraicSet, p) = throw(ArgumentError("Cannot add inequality to an algebraic set"))

struct FullSpace <: AbstractAlgebraicSet
end
function Base.show(io::IO, ::FullSpace)
    print(io, "R^n")
end
nequalities(::FullSpace) = 0
equalities(::FullSpace) = []
MP.changecoefficienttype(S::Type{FullSpace}, T::Type) = S

function MP.changecoefficienttype(set::AbstractSemialgebraicSet, T::Type)
    return convert(MP.changecoefficienttype(typeof(set), T), set)
end

Base.intersect(S::AbstractAlgebraicSet, T::FullSpace) = S
Base.intersect(S::FullSpace, T::AbstractAlgebraicSet) = T
Base.intersect(S::FullSpace, T::FullSpace) = S

# If `intersect(S, T)` is not implemented, this method will `StackOverflow`.
Base.intersect(S::AbstractSemialgebraicSet, T::AbstractSemialgebraicSet, args...; kws...) = intersect(intersect(S, T), args...; kws...)
# The keywords are only used when transforming `Element`
# into `BasicSemialgebraicSet`.
Base.intersect(set::AbstractSemialgebraicSet; kws...) = set

include("groebner.jl")
include("ideal.jl")
include("solve.jl")
include("variety.jl")
include("basic.jl")

include("fix.jl")
include("macro.jl")

end # module


text:
@test @inferred(Octavian.cld_fast(Octavian.StaticInt(6), Octavian.StaticInt(2))) == @inferred(Octavian.cld_fast(11, Octavian.StaticInt(4))) == @inferred(Octavian.cld_fast(Octavian.StaticInt(7), 3)) == @inferred(Octavian.cld_fast(9, 4)) == 3
@test @inferred(Octavian.cld_fast(Octavian.StaticInt(7), Octavian.StaticInt(2))) == @inferred(Octavian.cld_fast(13, Octavian.StaticInt(4))) == @inferred(Octavian.cld_fast(Octavian.StaticInt(8), 2)) == @inferred(Octavian.cld_fast(15, 4)) == 4


text:
## Make MNIST Graph
using MLDatasets
train_x, train_y = MNIST.traindata()
test_x,  test_y  = MNIST.testdata()
##
m = size(train_x,1)*size(train_x,2)
n = size(train_x,3)
X = Float64.(reshape(train_x, m, n))
## Build the 5-NN graph
using NearestNeighbors
@time T = KDTree(X)
##
k = 4
idxs = knn(T, X, k+1)[1]
##
ei = Int[]
ej = Int[]
for i=1:n
  for j=idxs[i]
    if i != j
      push!(ei,i)
      push!(ej,j)
    end
  end
end
##
using SparseArrays
A = sparse(ei,ej,1,n,n)
A = max.(A,A')
##
using DelimitedFiles
writedlm("mnist-train-$k.edges", hcat(findnz(A)[1:2]...))


text:
#=

        Creates a matrix of raw noise.
        Thanks to Cormullion for this example!

=#

using Luxor, Colors, Images, ImageFiltering, ImageView

m = @imagematrix begin
    background("black")
    t = Tiler(400, 400, 400, 400)
    for (pos, n) in t
        sethue(HSB(360 * noise(t.currentrow * 0.009, t.currentcol * 0.009), 0.8, 0.8))
        box(pos, 1, 1, :fillstroke)
    end
end 400 400

function convertmatrixtocolors(m)
    return convert.(Colors.RGBA, m)
end

img = convertmatrixtocolors(m)
imshow(img)

max_stars_repo_name: TheCedarPrince/NeuriViz

max_issues_repo_name: TheCedarPrince/NeuriViz

max_forks_repo_name: TheCedarPrince/NeuriViz


text:
export get_diet

"""
    get_diet(your_age, your_sex, your_weight; kwargs)

Return a personalized list of foods and their recommended intake amount in grams per day to fulfill all your nutritional requirements with minimum calories. `your_age` is in years, `your_sex` can be `male` or `female`, `your_weight` is in kilograms.

## kwargs:

* `only_groups`: Array with list of food group IDs. Only include foods in in these categories. Defaults to `[100, 500, 900, 1000, 1100, 1200, 1300, 1500, 1600, 1700, 2000]`. Here is a list of all IDs and their corresponding food category:
  1. 0100: Dairy and Egg Products
  1. 0200: Spices and Herbs
  1. 0300: Baby Foods
  1. 0400: Fats and Oils
  1. 0500: Poultry Products
  1. 0600: Soups, Sauces, and Gravies
  1. 0700: Sausages and Luncheon Meats
  1. 0800: Breakfast Cereals
  1. 0900: Fruits and Fruit Juices
  1. 1000: Pork Products
  1. 1100: Vegetables and Vegetable Products
  1. 1200: Nut and Seed Products
  1. 1300: Beef Products
  1. 1400: Beverages
  1. 1500: Finfish and Shellfish Products
  1. 1600: Legumes and Legume Products
  1. 1700: Lamb, Veal, and Game Products
  1. 1800: Baked Products
  1. 1900: Sweets
  1. 2000: Cereal Grains and Pasta
  1. 2100: Fast Foods
  1. 2200: Meals, Entrees, and Side Dishes
  1. 2500: Snacks
  1. 3500: American Indian/Alaska Native Foods
  1. 3600: Restaurant Foods
* `exclude_food`: Array with list of keywords. Exclude foods that have the such words in their description. Defaults to `["pickle", "pickled", "w/ added", "with added", "sweetened", "canned", "with artificial", "w/ artificial"]`.
* 'excluded_nutrients'=[]: Provide a list of IDs of nutrients not to be considered as a constraint to be fulfilled. By default the list is empty, and the following nutrients and their minimum and maximum recommened daily intakes are taken as constrains.
    * nutrient                : nutrient id
    1. Calcium                : 301
    1. Carbohydrate           : 205
    1. Choline                : 421
    1. Copper                 : 312
    1. Fluoride               : 313
    1. Folate                 : 417
    1. Histidine              : 512
    1. Iron                   : 303
    1. Isoleucine             : 503
    1. Leucine                : 504
    1. Lysine                 : 505
    1. Magnesium              : 304
    1. Manganese              : 315
    1. Methionine+cysteine    : 506
    1. Niacin                 : 406
    1. Pantothenic Acid       : 410
    1. Phenylalanine+tyrosine : 508
    1. Phosphorus             : 305
    1. Potassium              : 306
    1. Protein                : 203
    1. Riboflavin             : 405
    1. Selenium               : 317
    1. Sodium                 : 307
    1. Thiamin                : 404
    1. Threonine              : 502
    1. Total Fiber            : 291
    1. Tryptophan             : 501
    1. Valine                 : 510
    1. Vitamin A              : 320
    1. Vitamin B12            : 418
    1. Vitamin B6             : 415
    1. Vitamin C              : 401
    1. Vitamin D              : 328
    1. Vitamin E              : 323
    1. Vitamin K              : 430
    1. Zinc                   : 309
"""
function get_diet(your_age, your_sex, your_weight; only_groups=[100, 500, 900, 1000, 1100, 1200, 1300, 1500, 1600, 1700, 2000], exclude_food = ["pickle", "pickled", "w/ added", "with added", "sweetened", "canned", "with artificial", "w/ artificial"], excluded_nutrients=[])

	println("\nQuerying the database...\n")
	nutrs, dri, groups, foodids, foodnames, nfoods, nutamounts, calories, dri_ids, mins, maxs = loaddata(your_age, your_sex, your_weight; only_groups=only_groups);

  dri_ids, mins, maxs, nutamounts = exclude_nutrients(dri_ids, mins, maxs, nutamounts, [204, 255])
  if length(excluded_nutrients) > 0
    dri_ids, mins, maxs, nutamounts = exclude_nutrients(dri_ids, mins, maxs, nutamounts, excluded_nutrients)
  end
  exclude_indices, exclude_ids = return_excluded_food_indices(nutrs, exclude_food, foodids)

	println("\nOptimizing...\n")
  m, status, objval, yval = optimize(mins, maxs, nutamounts, calories, dri_ids, nfoods, exclude_indices);

  println("\nOptimization status: $status \n")
  if isnan(objval)
    return
  end
  println("Calorie intake: $(round(Int, objval)) \n")

  println("\n\nYour recommended diet per day:\n")

	toeatindices = findall(a-> a>0, yval);
	toeatnames = foodnames[toeatindices];
	toeatids = foodids[toeatindices];
	toeactamounts = yval[toeatindices];
	foodlongname = find_longdes(nutrs, foodids);
  toeatlongnames = foodlongname[toeatindices];
  diet_ids = Int[]
  diet_names = String[]
  diet_amounts = Float64[]
  for (index, nameamount) in enumerate(zip(toeatnames, toeactamounts, toeatlongnames))
    # amounts are in 100 grams. So the * 100 here converts the recommended amounts to grams
    amountingrams = round(nameamount[2] * 100, digits=3)
    if amountingrams > 1
      push!(diet_ids, toeatids[index])
      push!(diet_names, nameamount[3])
      push!(diet_amounts, amountingrams)
    end
  end

  print_diet(diet_names, diet_amounts)

  return diet_ids, diet_names, diet_amounts, m, status, objval, yval
end

function print_diet(diet_names, diet_amounts)
  println("")
  for i in 1:length(diet_names)
    println(i, ". ", diet_names[i], ": ", diet_amounts[i], " grams")
  end
end


text:
module ReversibleTuringMachine

using NiLang

export RTM, Quadruple, quadruples
export run!, is_deterministic, is_reversible
export BLANK, SLASH, LEFT, RIGHT, STAY

const BLANK = -10000000
const SLASH = -10000001
const LEFT = -10000002
const RIGHT = -10000003
const STAY = -10000004

struct RTM{QT,RT}
    qs::QT
    qf::QT
    rules::RT
end

struct Quadruple{QT, ST}
    q1::QT
    s1::ST
    s2::ST
    q2::QT
end

function quadruples(rules)
    [Quadruple(rules[i,:]...) for i=1:size(rules, 1)]
end

Base.show(io::IO, ::MIME"text/plain", state::Quadruple) = show(io, state)
function Base.show(io::IO, state::Quadruple)
    q1, s1, s2, q2 = state.q1, state.s1, state.s2, state.q2
    print(io, q1, ' ')
    if s1 == SLASH  # move rule
        print(io, "/ ")
        if s2 == RIGHT
            print(io, "→ ")
        elseif s2 == LEFT
            print(io, "← ")
        elseif s2 == STAY
            print(io, "↔ ")
        else
            error("wrong argument for SLASH rule $s2.")
        end
    elseif s1 >= 0 || s1 == BLANK # symbol rule
        print(io, s1 == BLANK ? '_' : s1, ' ')
        print(io, s2 == BLANK ? '_' : s2, ' ')
    end
    print(io, q2)
end

function is_deterministic(rtm::RTM)
    n = length(rtm.rules)
    for i = 1:n
        for j = i+1:n
            ri = rtm.rules[i]
            rj = rtm.rules[j]
            if ri.q1 == rj.q1 && ri.s1 == rj.s1
                return false
            elseif ri.q1 == rj.q1 && (ri.s1 == SLASH || rj.s1 == SLASH)
                return false
            end
        end
    end
    true
end

function is_reversible(rtm::RTM)
    n = length(rtm.rules)
    for i = 1:n
        for j = i+1:n
            ri = rtm.rules[i]
            rj = rtm.rules[j]
            if ri.q2 == rj.q2 && ri.s2 == rj.s2
                return false
            elseif ri.q2 == rj.q2 && (ri.s1 == SLASH || rj.s1 == SLASH)
                return false
            end
        end
    end
    true
end

"""
    run!(rtm::RTM, tape, loc, pc; visualize=false)

Run a reversible turing machine `rtm` on `tape`.
`loc` is the location of reading head,
`pc` is the program counter used to search rules, normally, one can just set it to `1`.
"""
@i function run!(rtm::RTM, tape, loc, pc; visualize=false)
    q ← rtm.qs
    while (q != rtm.qf, q != rtm.qs)
        # state q, at loc
        if (visualize, ~)
            @safe loc0 = loc
            @safe s0 = tape[loc]
            @safe q0 = q
            instr(q, tape[loc], rtm.rules[pc], loc)
            @safe if loc != loc0 || s0 != tape[loc] || q != q0
                viz_tape(tape, q, loc, rtm.rules[pc])
            end
        else
            instr(q, tape[loc], rtm.rules[pc], loc)
        end
        pc += identity(1)
        if (pc == length(rtm.rules)+1, pc == 1)
            pc -= identity(length(rtm.rules))
        end
    end
    q → rtm.qf
end

"""
    instr(q, s, cmd::Quadruple, loc)

For turing machine in state `q`, at location `loc`, reading symbol `s`,
"Try" to run an instruction `cmd`.
"""
@i function instr(q, s, cmd::Quadruple, loc)
    if (q == cmd.q1 && s == cmd.s1, q == cmd.q2 && s == cmd.s2)
        # Symbol rule
        q += cmd.q2 - cmd.q1
        s += cmd.s2 - cmd.s1
    elseif (q == cmd.q1 && cmd.s1 == SLASH, q == cmd.q2 && cmd.s1 == SLASH)
        # Shift rule
        q += cmd.q2 - cmd.q1
        if (cmd.s2 == RIGHT, ~)
            loc += identity(1)
        elseif (cmd.s2 == LEFT, ~)
            loc -= identity(1)
        end
    end
end

function viz_tape(tape, q, loc, comment="")
    print(' '^(4*(loc-1)+1))
    print(" $q")
    println()
    print(' '^(4*(loc-1)+1))
    print(" ↓")
    if comment != ""
        print(" ($comment)")
    end
    println()
    print('|')
    for elem in tape
        s = string(elem)[1]
        print(" $s |")
    end
    println()
    println("—"^(4*length(tape)+1))
end

end # module


text:
module ContinuousTransformations

import Base.rand

using ArgCheck
using AutoHashEquals
import Distributions:
    ContinuousUnivariateDistribution, ContinuousMultivariateDistribution, logpdf
using DocStringExtensions
using LinearAlgebra
using MacroTools
using Parameters
using StatsFuns
using Lazy
# using Unrolled

export
    # intervals
    AbstractInterval, RealLine, ℝ, PositiveRay, ℝ⁺, NegativeRay, ℝ⁻,
    Segment, width,
    # general API
    domain, image, isincreasing, transform, inverse, logjac,
    transform_and_logjac, @increment_logjac, ContinuousTransformation,
    # univariate transformations
    UnivariateTransformation,
    Affine, IDENTITY,
    Negation, NEGATION, Logistic, LOGISTIC, RealCircle, REALCIRCLE, Exp, EXP,
    Logit, LOGIT, InvRealCircle, INVREALCIRCLE, Log, LOG,
    affine_bridge, default_transformation, bridge,
    # multivariate transformations NOTE API may change
    UnitVector, CorrelationCholeskyFactor, lkj_correlation_cholesky_logpdf,
    # grouped transformations
    GroupedTransformation, get_transformation, ArrayTransformation,
    TransformationTuple,
    # wrapped transformations
    TransformationWrapper, TransformLogLikelihood, get_loglikelihood,
    TransformDistribution, get_distribution, logpdf_in_domain, logpdf_in_image,
    # utilities
    ungrouping_map, all_finite

import Base: in, length, size, ∘, show, getindex, intersect, extrema, minimum,
    maximum, isfinite, isinf, isapprox

import Statistics: middle

# utilities

"""
    $SIGNATURES

Define an `isapprox` method, comparing the given fields in type `T`.
"""
macro define_isapprox(T, fields...)
    body = foldl((a, b) -> :($(a) && $(b)),
                 :(isapprox(x.$(field), y.$(field); rtol=rtol, atol=atol))
                 for field in fields)
    quote
        function Base.isapprox(x::$T, y::$T; rtol::Real=√eps(), atol::Real=0)
            $(body)
        end
    end
end

"""
    $SIGNATURES

Define a singleton type with the given name and supertype (specified
as `name <: supertype`), and a constant which defaults to the name in
uppercase.
"""
macro define_singleton(name_and_supertype, constant = nothing)
    @capture name_and_supertype name_ <: supertype_
    if constant == nothing
        constant = Symbol(uppercase(string(name)))
    end
    quote
        Core.@__doc__ struct $(esc(name)) <: $(esc(supertype)) end
        Core.@__doc__ const $(esc(constant)) = $(esc(name))()
    end
end

"""
    _fma(x, y, z)

Placeholder for `Base.fma` until
https://github.com/JuliaDiff/ReverseDiff.jl/issues/86 is fixed.
"""
_fma(x, y, z) = x*y+z

"""
    $SIGNATURES

Define [`transform`](@ref) and [`logjac`](@ref) using
[`transform_and_logjac`](@ref).
"""
macro define_from_transform_and_logjac(T)
    S = esc(T)
    quote
        (t::$S)(x) = transform_and_logjac(t, x)[1]
        ContinuousTransformations.logjac(t::$S, x) = transform_and_logjac(t, x)[2]
    end
end

# abstract interface for transformations

"""
    logjac(t, x)

The log of the determinant of the Jacobian of `t` at `x`.
```
"""
function logjac end

"""
    inverse(t, x)

Return ``t⁻¹(x)``.

    inverse(t)

Return the transformation ``t⁻¹``.
"""
inverse(t, x) = inverse(t)(x)

"""
    domain(transformation)

Return the domain of the transformation.
"""
function domain end

"""
    image(transformation)

Return the image of the transformation.
"""
function image end

"""
$TYPEDEF

Continuous bijection ``D ⊂ ℝ^n→ I ⊂ ℝ^n`` or ``D ⊂ ℝ → I ⊂ ℝ``.
"""
abstract type ContinuousTransformation <: Function end

"""
    transform(t, x)

Return ``t(x)``. NOTE: this equivalent to the callable syntax `t(x)`, which is
also available.
"""
transform(t::ContinuousTransformation, x) = t(x)

"""
    y, logjac = $SIGNATURES

Return the transformed value and the determinant of the log Jacobian.

Equivalent to `t(x), logjac(t, x)`, but may be faster because of reused
components.
"""
function transform_and_logjac(t::ContinuousTransformation, x)
    transform(t, x), logjac(t, x)
end

"""
    $SIGNATURES

Return `transformation(x)` while incrementing `logjac` with the log Jacobian
determinant.
"""
macro increment_logjac(logjac, transformation, x)
    quote
        y, lj = transform_and_logjac($(esc(transformation)), $(esc(x)))
        $(esc(logjac)) += lj
        y
    end
end

"""
    $SIGNATURES

Return `transformation(x)` while incrementing `logjac` with the log Jacobian
determinant.
"""
macro increment_logjac(logjac, transformation_x)
    @capture transformation_x transformation_(x_)
    quote
        @increment_logjac $(esc(logjac)) $(esc(transformation)) $(esc(x))
    end
end

"""
    rhs_string(transformation, term)

Return the formula representing the hand side of the `transformation`, with
`term` as the argument.
"""
function rhs_string end

transformation_string(t, x = "x") = x * " ↦ " * rhs_string(t, x)

# this is necessary because `display` calls this method, which is defined for
# Function and would match
show(io::IO, ::MIME"text/plain", t::ContinuousTransformation) =
    println(io, transformation_string(t, "x"))

show(io::IO, t::ContinuousTransformation) =
    print(io, transformation_string(t, "x"))

"""
    bridge(dom, img, [transformation])

Return a transformation that maps `dom` to `img`.

The `transformation` argument may be used to specify a particular transformation
family, otherwise `default_transformation` is used.
"""
function bridge end

# intervals

"""
$TYPEDEF

Abstract supertype for all univariate intervals. It is not specified whether
they are open or closed.
"""
abstract type AbstractInterval end

isinf(x::AbstractInterval) = !isfinite(x)

isapprox(::AbstractInterval, ::AbstractInterval; rtol=√eps(), atol=0) = false

extrema(x::AbstractInterval) = minimum(x), maximum(x)

"""
    RealLine()

The real line. Use the constant [`ℝ`](@ref).
"""
struct RealLine <: AbstractInterval end

isapprox(::RealLine, ::RealLine; rtol=√eps(), atol=0) = true

show(io::IO, ::RealLine) = print(io, "ℝ")

"A constant for the real line."
const ℝ = RealLine()

in(x::Real, ::RealLine) = true

minimum(::RealLine) = -Inf
maximum(::RealLine) = Inf

isfinite(::RealLine) = false

"""
    PositiveRay(left)

The real numbers above `left`. See [`ℝ⁺`](@ref).
"""
@auto_hash_equals struct PositiveRay{T <: Real} <: AbstractInterval
    left::T
    function PositiveRay{T}(left::T) where T
        @argcheck isfinite(left) "Need finite endpoint."
        new(left)
    end
end

PositiveRay(left::T) where {T} = PositiveRay{T}(left)

in(x::Real, ray::PositiveRay) = ray.left ≤ x
minimum(ray::PositiveRay) = ray.left
maximum(::PositiveRay{T}) where T = T(Inf)
isfinite(::PositiveRay) = false

"The positive real numbers."
const ℝ⁺ = PositiveRay(0.0)

@define_isapprox PositiveRay left

"""
    NegativeRay(right)

The real numbers below `right`. See [`ℝ⁻`](@ref).
"""
@auto_hash_equals struct NegativeRay{T <: Real} <: AbstractInterval
    right::T
    function NegativeRay{T}(right::T) where T
        @argcheck isfinite(right) "Need finite endpoint."
        new(right)
    end
end

NegativeRay(right::T) where {T} = NegativeRay{T}(right)

in(x::Real, ray::NegativeRay) = x ≤ ray.right
minimum(::NegativeRay{T}) where T = -T(Inf)
maximum(ray::NegativeRay) = ray.right
isfinite(::NegativeRay) = false

"The negative real numbers."
const ℝ⁻ = NegativeRay(0.0)

@define_isapprox NegativeRay right

"""
    Segment(left, right)

The real numbers between `left` and `right`, with
``-∞ < \\text{left} < \\text{right} < ∞`` enforced.
"""
@auto_hash_equals struct Segment{T <: Real} <: AbstractInterval
    left::T
    right::T
    function Segment{T}(left::T, right::T) where T
        @argcheck isfinite(left) && isfinite(right) "Need finite endpoints."
        @argcheck left < right "Need strictly increasing endpoints."
        new(left, right)
    end
end

Segment(left::T, right::T) where {T <: Real} = Segment{T}(left, right)

Segment(left::Real, right::Real) = Segment(promote(left, right)...)

in(x::Real, s::Segment) = s.left ≤ x ≤ s.right
minimum(s::Segment) = s.left
maximum(s::Segment) = s.right
isfinite(::Segment) = true

@define_isapprox Segment left right

"""
    $SIGNATURES

Width of a finite interval.
"""
width(s::Segment) = s.right - s.left

middle(s::Segment) = middle(s.left, s.right)

# intersections

## general fallback method. define specific methods with the following
## argument ordering Segment < PositiveRay < NegativeRay < RealLine < all
intersect(a::AbstractInterval, b::AbstractInterval) = intersect(b, a)

intersect(a::RealLine, b::AbstractInterval) = b

"""
    $SIGNATURES

Helper function for forming a segment when possible. Internal, not exported.
"""
@inline function _maybe_segment(a, b)
    # NOTE Decided not to represent the empty interval, as it has no use in
    # the context of this package. Best to throw an error as soon as possible.
    a < b ? Segment(a, b) :
        throw(ArgumentError("intersection of intervals is empty"))
end

intersect(a::Segment, b::Segment) =
    _maybe_segment(max(a.left, b.left), min(a.right, b.right))

intersect(a::Segment, b::PositiveRay) =
    _maybe_segment(max(b.left, a.left), a.right)

intersect(a::Segment, b::NegativeRay) =
    _maybe_segment(a.left, min(a.right, b.right))

intersect(a::PositiveRay, b::PositiveRay) = PositiveRay(max(a.left, b.left))

intersect(a::PositiveRay, b::NegativeRay) = _maybe_segment(a.left, b.right)

intersect(a::NegativeRay, b::NegativeRay) = NegativeRay(min(a.right, b.right))

# univariate transformations

"""
$TYPEDEF

Univariate monotone transformation, either *increasing* or *decreasing* on the
whole domain (thus, a bijection).
"""
abstract type UnivariateTransformation <: ContinuousTransformation end

length(::UnivariateTransformation) = 1

size(::UnivariateTransformation) = ()

"""
    isincreasing(transformation)

Return `true` (`false`), when the transformation is monotonically increasing
(decreasing).
"""
function isincreasing end

"""
$TYPEDEF

Trait that is useful for domain and image calculations. See [`RRStable`](@ref).
"""
abstract type RRStability end

"""
$TYPEDEF

Trait that indicates that a univariate transformation

1. maps ``ℝ`` to ``ℝ``,

2. supports mapping intervals, and

3. maps subtypes of [`AbstractInterval`](@ref) to the same type.
"""
struct RRStable <: RRStability end

"""
$TYPEDEF

Trait that indicates that a univariate transformation is *not* [`RRStable`](@ref).
"""
struct NotRRStable <: RRStability end

"""
    $SIGNATURES

Return either the trait [`RRStable`](@ref) and `NotRRStable`.
"""
RR_stability(::UnivariateTransformation) = NotRRStable()

∘(::RRStability, ::RRStability) = NotRRStable()

∘(::RRStable, ::RRStable) = RRStable()

# Affine and Negation
#
# These play a special role since map subtypes of `AbstractInterval` to
# themselves, and thus are useful for composition.

"""
    Affine(α, β)

Mapping ``ℝ → ℝ`` using ``x ↦ α⋅x + β``.

``α > 0`` is enforced, see [`Negation`](@ref).
"""
@auto_hash_equals struct Affine{T <: Real} <: UnivariateTransformation
    α::T
    β::T
    function Affine{T}(α::T, β::T) where T
        @argcheck α > 0 DomainError
        new(α, β)
    end
end

Affine(α::T, β::T) where T = Affine{T}(α, β)
Affine(α, β) = Affine(promote(α, β)...)

"Identity (as an affine transformation)."
const IDENTITY = Affine(1, 0)

domain(::Affine) = ℝ
image(::Affine) = ℝ
(t::Affine)(x) = _fma(x, t.α, t.β)
logjac(t::Affine, x) = log(abs(t.α))
inverse(t::Affine) = Affine(1/t.α, -t.β/t.α)
isincreasing(t::Affine) = true

RR_stability(::Affine) = RRStable()

(t::Affine)(x::Segment) = Segment(t(x.left), t(x.right))
(t::Affine)(x::PositiveRay) = PositiveRay(t(x.left))
(t::Affine)(x::NegativeRay) = NegativeRay(t(x.right))
(t::Affine)(::RealLine) = ℝ

function rhs_string(t::Affine, term)
    @unpack α, β = t
    α == 1 || (term = "$(round(α; sigdigits = 4))⋅" * term)
    β == 0 || (term = term * " + $(round(β, sigdigits = 4))")
    term
end

"""
    Negation()

Mapping ``ℝ → ℝ`` using ``x ↦ -x``.
"""
@define_singleton Negation <: UnivariateTransformation

domain(::Negation) = ℝ
image(::Negation) = ℝ
(::Negation)(x) = -x
logjac(::Negation, x) = zero(x)
inverse(::Negation) = NEGATION
isincreasing(::Negation) = false

RR_stability(::Negation) = RRStable()

(::Negation)(x::Segment) = Segment(-x.right, -x.left)
(::Negation)(x::PositiveRay) = NegativeRay(-x.left)
(::Negation)(x::NegativeRay) = PositiveRay(-x.right)
(::Negation)(::RealLine) = ℝ

rhs_string(::Negation, term) = "-" * term

# transformations from ℝ to subsets

"""
    Logistic()

Mapping ``ℝ → (0,1)`` using ``x ↦ 1/(1+\\exp(-x))``.
"""
@define_singleton Logistic <: UnivariateTransformation

domain(::Logistic) = ℝ
image(::Logistic) = Segment(0, 1)
(t::Logistic)(x) = logistic(x)
logjac(::Logistic, x) = -(log1pexp(x)+log1pexp(-x))
inverse(::Logistic) = LOGIT
isincreasing(::Logistic) = true

rhs_string(::Logistic, term) = "logistic($term)"

"""
    RealCircle()

Mapping ``ℝ → (-1,1)`` using ``x ↦ x/√(1+x^2)``.
"""
@define_singleton RealCircle <: UnivariateTransformation

domain(::RealCircle) = ℝ
image(::RealCircle) = Segment(-1, 1)
(t::RealCircle)(x) = isinf(x) ? sign(x) : x/√(1+x^2)
logjac(::RealCircle, x) = -1.5*log1psq(x)
inverse(::RealCircle) = INVREALCIRCLE # x/√(1-x^2)
isincreasing(::RealCircle) = true

rhs_string(::RealCircle, term) = "realcircle($term)"

"""
    Exp()

Mapping ``ℝ → ℝ⁺`` using ``x ↦ \\exp(x)``.
"""
@define_singleton Exp <: UnivariateTransformation

domain(::Exp) = ℝ
image(::Exp) = ℝ⁺
(t::Exp)(x) = exp(x)
logjac(::Exp, x) = x
inverse(::Exp) = LOG
isincreasing(::Exp) = true

rhs_string(::Exp, term) = "exp($term)"

# transformations to ℝ from subsets

"""
    Logit()

Mapping ``(0,1) → ℝ`` using ``x ↦ \\log(x/(1-x))``.
"""
@define_singleton Logit <: UnivariateTransformation

domain(::Logit) = Segment(0, 1)
image(::Logit) = ℝ
(t::Logit)(x) = logit(x)
logjac(::Logit, x) = -log(x)-log1p(-x)
inverse(::Logit) = LOGISTIC
isincreasing(::Logit) = true

rhs_string(::Logit, term) = "logit($term)"

"""
    InvRealCircle()

Mapping ``(-1,1) → ℝ`` using ``x ↦ x/√(1-x^2)``.
"""
@define_singleton InvRealCircle <: UnivariateTransformation

domain(::InvRealCircle) = Segment(-1, 1)
image(::InvRealCircle) = ℝ
(t::InvRealCircle)(x) = x/√(1-x^2)
logjac(::InvRealCircle, x) = -1.5*log1p(-abs2(x))
inverse(::InvRealCircle) = REALCIRCLE
isincreasing(::InvRealCircle) = true

rhs_string(::InvRealCircle, term) = "realcircle⁻¹($term)"

"""
    Log()

Mapping ``ℝ → ℝ⁺`` using ``x ↦ \\exp(x)``.
"""
@define_singleton Log <: UnivariateTransformation

domain(::Log) = ℝ⁺
image(::Log) = ℝ
(t::Log)(x) = log(x)
logjac(::Log, x) = -log(x)
inverse(::Log) = EXP
isincreasing(::Log) = true

rhs_string(::Log, term) = "log($term)"

# composed transformations

"""
    ComposedTransformation(f, g)

Compose two univariate transformations, resulting in the mapping
``f∘g``, or ``x ↦ f(g(x))`.

Use the `∘` operator for construction.
"""
struct ComposedTransformation{Tf <: UnivariateTransformation,
                              Tg <: UnivariateTransformation} <:
                                  UnivariateTransformation
    f::Tf
    g::Tg
end

RR_stability(t::ComposedTransformation) = RR_stability(t.f) ∘ RR_stability(t.g)

(c::ComposedTransformation)(x) = c.f(c.g(x))

rhs_string(t::ComposedTransformation, term) =
    rhs_string(t.f, rhs_string(t.g, term))

# composed image and domain uses traits, effectively allowing calculation for
# composing with Affine and Negation

"""
    composed_image(f_RR_stability, g_RR_stability, f, g)
"""
composed_image(::RRStable, ::RRStability, f, g) = f(image(g))

composed_image(::RRStability, ::RRStable, f, g) = image(f)

composed_image(::RRStable, ::RRStable, f, g) = ℝ

"""
    composed_domain(f_RR_stability, g_RR_stability, f, g)
"""
composed_domain(::RRStable, ::RRStability, f, g) = domain(g)

composed_domain(::RRStability, ::RRStable, f, g) = inverse(g)(domain(f))

composed_domain(::RRStable, ::RRStable, f, g) = ℝ

image(t::ComposedTransformation) =
    composed_image(RR_stability(t.f), RR_stability(t.g), t.f, t.g)

domain(t::ComposedTransformation) =
    composed_domain(RR_stability(t.f), RR_stability(t.g), t.f, t.g)

function logjac(t::ComposedTransformation, x)
    @unpack f, g = t
    y = g(x)
    log_g′x = logjac(g, x)
    log_f′y = logjac(f, y)
    log_f′y + log_g′x
end

inverse(t::ComposedTransformation, x) = inverse(t.g, inverse(t.f, x))

inverse(t::ComposedTransformation) =
    ComposedTransformation(inverse(t.g), inverse(t.f))

isincreasing(c::ComposedTransformation) = isincreasing(c.f) == isincreasing(c.g)

∘(f::UnivariateTransformation, g::UnivariateTransformation) =
    ComposedTransformation(f, g)

∘(f::Affine, g::Affine) = Affine(f.α*g.α, _fma(f.α, g.β, f.β))

# calculated transformations

"""
    affine_bridge(interval1, interval1)

Return an affine transformation between two intervals of the same type.
"""
function affine_bridge(x::Segment, y::Segment)
    α = width(y) / width(x)
    β = middle(y) - middle(x) * α
    Affine(α, β)
end

affine_bridge(::RealLine, ::RealLine) = IDENTITY

affine_bridge(x::PositiveRay, y::PositiveRay) = Affine(1, y.left - x.left)
affine_bridge(x::NegativeRay, y::NegativeRay) = Affine(1, y.right - x.right)
affine_bridge(x::PositiveRay, y::NegativeRay) =
    Affine(1, y.right + x.left) ∘ NEGATION
affine_bridge(x::NegativeRay, y::PositiveRay) =
    Affine(1, y.left + x.right) ∘ NEGATION

"""
    default_transformation(dom, img)

Return a transformation from `dom` that can be mapped to `img` using
`affine_bridge`.
"""
default_transformation(::RealLine, ::Segment) = LOGISTIC
default_transformation(::RealLine, ::PositiveRay) = EXP
default_transformation(::RealLine, ::NegativeRay) = EXP

default_transformation(::Segment, ::RealLine) = LOGIT
default_transformation(::PositiveRay, ::RealLine) = LOG
default_transformation(::NegativeRay, ::RealLine) = LOG

bridge(dom::RealLine, img, t = default_transformation(dom, img)) =
    affine_bridge(image(t), img) ∘ t

bridge(dom, img::RealLine, t = default_transformation(dom, img)) =
    t ∘ affine_bridge(dom, domain(t))

bridge(::RealLine, ::RealLine) = IDENTITY

# Special constants

"""
Hyperbolic tangent transformation.

An affine stretch of LOGISTIC to (-1, 1).
"""
const TANH = bridge(ℝ, Segment(-1.0, 1.0), LOGISTIC)

# UnitVector

"""
    UnitVector(n)

Transform `n-1` real numbers to a unit vector of length `n`, under the Euclidean
norm.
"""
struct UnitVector <: ContinuousTransformation
    n::Int
    function UnitVector(n::Int)
        @argcheck n ≥ 1 "Dimension should be positive."
        new(n)
    end
end

length(t::UnitVector) = t.n - 1

rhs_string(t::UnitVector, term) = "UnitVector($(t.n))($term)"

function transform_and_logjac!(ys::AbstractVector{T},
                               t::UnitVector, xs::AbstractVector) where T
    @unpack n = t
    @argcheck length(xs) == n - 1
    @argcheck length(ys) == n
    remainder = one(T)
    logjac = zero(T)
    for (i, x) in enumerate(xs)
        z = @increment_logjac logjac TANH(x)
        ys[i] = z * √remainder
        logjac += log(remainder) / 2
        remainder -= abs2(z) * remainder
    end
    ys[end] = √remainder
    ys, logjac
end

transform_and_logjac(t::UnitVector, zs::AbstractVector{T}) where T =
    transform_and_logjac!(Vector{T}(undef, t.n), t, zs)

function inverse!(xs::AbstractVector{T},
                  t::UnitVector, ys::AbstractVector) where T
    @unpack n = t
    @argcheck length(xs) == n - 1
    @argcheck length(ys) == n
    remainder = one(T)
    for i in 1:(n - 1)
        y = ys[i]
        z = y / √remainder
        xs[i] = inverse(TANH, z)
        remainder -= abs2(y)
    end
    xs
end

inverse(t::UnitVector, ys::AbstractVector{T}) where T =
    inverse!(Vector{T}(undef, t.n - 1), t, ys)

@define_from_transform_and_logjac UnitVector

# Cholesky factor of a correlation matrix

"""
    CorrelationCholeskyFactor(n)

Cholesky factor of a correlation matrix of size `n`.
"""
struct CorrelationCholeskyFactor <: ContinuousTransformation
    n::Int
    function CorrelationCholeskyFactor(n)
        @argcheck n ≥ 1 "Dimension should be positive."
        new(n)
    end
end

"""
    $SIGNATURES

Number of elements in a triangle of an ``n×n`` square matrix, including the
diagonal.
"""
triangle_length(n) = n * (n-1) ÷ 2

length(t::CorrelationCholeskyFactor) = triangle_length(t.n)

rhs_string(t::CorrelationCholeskyFactor, term) =
    "CorrelationCholeskyFactor($(t.n))($term)"

function transform_and_logjac(t::CorrelationCholeskyFactor,
                              x::AbstractVector{T}) where T
    @unpack n = t
    @argcheck length(x) == length(t)
    logjac = zero(T)
    cumulative_index = 0
    L = zeros(T, n, n)
    for i in 1:n
        k = i - 1
        _, lj = transform_and_logjac!(@view(L[i, 1:i]), UnitVector(i),
                                      @view(x[cumulative_index .+ (1:k)]))
        logjac += lj
        cumulative_index += k
    end
    LowerTriangular(L), logjac
end

@define_from_transform_and_logjac CorrelationCholeskyFactor

function inverse(t::CorrelationCholeskyFactor,
                 L::LowerTriangular{T}) where T
    @unpack n = t
    @argcheck size(L) == (n, n)
    x = Vector{T}(undef, triangle_length(n))
    cumulative_index = 0
    for i in 1:n
        k = i - 1
        inverse!(@view(x[cumulative_index .+ (1:k)]), UnitVector(i),
                 @view(L[i, 1:i]))
        cumulative_index += k
    end
    x
end

# helper function for Cholesky decompositions

"""
    $SIGNATURES

Log PDF of the LKJ distribution (Lewandowski et al 2009) for correlation
matrices.

A correlation matrix ``Ω=LL'`` has the density ``|Ω|^(η-1)``. However, it is
usually not necessary to construct ``Ω``, so this function is parametrized in
terms of a Cholesky decomposition `L`.

Note that this function **does not check if `L` yields a valid correlation
matrix**. Use [`CorrelationCholeskyFactor`](@ref) for generating valid
arguments.

Valid values are ``η > 0``. When ``η > 1``, the distribution is unimodal at the
identity, while ``0 < η < 1`` has a trough. ``η = 2`` is recommended as a vague
prior.

When ``η = 1``, the density is uniform. Note however that the function should
still be invoked, because of the Jacobian correction of the transformation.
"""
function lkj_correlation_cholesky_logpdf(L::Union{LowerTriangular,
                                                  UpperTriangular}, η)
    @argcheck η > 0
    z = diag(L)
    n = size(L, 1)
    sum(log.(z) .* ((n:-1:1) .+ 2*(η-1))) + log(2) * n
end

# grouped transformations

"""
$TYPEDEF

Abstract type for grouped transformations.

A grouped transformation takes a vector, and transforms contiguous blocks of
elements to some output type, determined by the specific transformation type.

All subtypes support

- `length`: return the length of the vector that can be used as an argument

- callable object for the transformation

- `logjac`, and `inverse`,

- `domain` and `image`, which may have specific interpretation for their result
  types depending on the concrete subtype.
"""
abstract type GroupedTransformation <: ContinuousTransformation end

"""
    $SIGNATURES

Return the transformation from a wrapper object.
"""
get_transformation(d) = d.transformation

# array transformations

"""
    ArrayTransformation(transformation, dims)
    ArrayTransformation(transformation, dims...)

Apply transformation to a vector, returning an array of the given dimensions.

[`domain`](@ref), [`image`](@ref), and [`isincreasing`](@ref) return the
corresponding values for the underlying transformation.
"""
struct ArrayTransformation{T <: UnivariateTransformation,
                           D} <: GroupedTransformation
    transformation::T
    function ArrayTransformation(transformation::T,
                                 dims::Tuple{Vararg{Int64, N}}) where {T,N}
        @argcheck all(dims .> 0) "Invalid dimensions."
        new{T,dims}(transformation)
    end
end

ArrayTransformation(transformation, dims::Int...) =
    ArrayTransformation(transformation, dims)

size(t::ArrayTransformation{T,D}) where {T,D} = D

length(t::ArrayTransformation) = prod(size(t))

function transformation_string(t::ArrayTransformation, term)
    s = size(t)
    s_print = length(s) == 1 ? s[1] : s
    "$(transformation_string(t.transformation, term)) for $(s_print) elements"
end

(t::ArrayTransformation)(x) = reshape((t.transformation).(x), size(t))

function logjac(t::ArrayTransformation, x)
    @argcheck length(x) == length(t) DimensionMismatch()
    lj(x) = logjac(t.transformation, x)
    sum(lj, x)
end

inverse(t::ArrayTransformation, x) =
    reshape(inverse.(t.transformation, x), size(t))

@forward ArrayTransformation.transformation isincreasing, image, domain

# transformation tuple

"""
    TransformationTuple(transformations::Tuple)
    TransformationTuple(transformations...)

A tuple of [`ContinuousTransformation`](@ref)s. Given a vector of matching
`length`, each takes as many reals as needed, and returns the result as a tuple.
"""
@auto_hash_equals struct
    TransformationTuple{T <: Tuple{Vararg{ContinuousTransformation}}} <:
        GroupedTransformation
    transformations::T
end

TransformationTuple(ts::ContinuousTransformation...) = TransformationTuple(ts)

function transformation_string(t::TransformationTuple, x)
    s = IOBuffer()
    print(s, "TransformationTuple")
    for (tt, ix) in zip(t.transformations,
                        transformation_indexes(t.transformations))
        print(s, "\n    ")
        print(s, transformation_string(tt, "$(x)[$(ix)]"))
    end
    String(take!(s))
end

length(t::TransformationTuple) = sum(length(t) for t in t.transformations)

domain(t::TransformationTuple) = domain.(t.transformations)

image(t::TransformationTuple) = image.(t.transformations)

getindex(t::TransformationTuple, ix) = t.transformations[ix]

function next_indexes(acc, t)
    l = length(t)
    acc + l, acc .+ (1:l)
end

next_indexes(acc, t::UnivariateTransformation) = acc+1, acc+1

#=@unroll=# function transformation_indexes(ts)
    acc = 0
    result = ()
    #=@unroll=# for t in ts
        acc, ix = next_indexes(acc, t)
        result = (result..., ix)
    end
    result
end

function (t::TransformationTuple)(x)
    ts = t.transformations
    map((t,ix) -> t(x[ix]), ts, transformation_indexes(ts))
end

function logjac(t::TransformationTuple, x)
    ts = t.transformations
    sum(map((t,ix) -> logjac(t, x[ix]), ts, transformation_indexes(ts)))
end

inverse(t::TransformationTuple, y::Tuple) =
    vcat(map(inverse, t.transformations, y)...)

# wrapper

"""
$TYPEDEF

Wrap a transformation to achieve some specialized functionality.

Supports `length`, [`get_transformation`](@ref), and other methods depending on
the subtype.
"""
abstract type TransformationWrapper <: Function end

# loglikelihood wrapper

"""
    TransformLogLikelihood(ℓ, transformation::Union{Tuple, GroupedTransformation})

    TransformLogLikelihood(ℓ, transformations...)

Return a callable that

1. transforms its vector argument using a grouped transformation to a set of
   values,

2. calls `ℓ` (which should return a scalar) with this tuple.

3. returns the result above corrected by the log Jacobians.

Useful when `ℓ` is a log-likelihood function with a restricted domain, and
`transformations` is used to trasform to this domain from ``ℝ^n``.

See also [`get_transformation`](@ref), [`get_distribution`](@ref),
`Distributions.logpdf`, and [`logpdf_in_domain`](@ref).
"""
struct TransformLogLikelihood{L, T <: GroupedTransformation} <:
        TransformationWrapper
    loglikelihood::L
    transformation::T
end

TransformLogLikelihood(L, T::Tuple) =
    TransformLogLikelihood(L, TransformationTuple(T))

TransformLogLikelihood(L, ts::ContinuousTransformation...) =
    TransformLogLikelihood(L, TransformationTuple(ts))

"""
    $SIGNATURES

Return the log likelihood function.
"""
get_loglikelihood(t::TransformLogLikelihood) = t.loglikelihood

@forward TransformLogLikelihood.transformation length

(f::TransformLogLikelihood)(x) =
    f.loglikelihood(f.transformation(x)) + logjac(f.transformation, x)

function show(io::IO, f::TransformLogLikelihood)
    print(io, "TransformLogLikelihood of length $(length(f)), with ")
    print(io, f.transformation)
end

# transform distributions

"""
    TransformDistribution(distribution, transformation)

Given a `transformation` and a `distribution`, create a transformed distribution
object that has the distribution of `transformation(x)` with `x ∼ distribution`.

The transformation object is callable with the same syntax as
`transformation`. It also supports methods `rand`, `length`.

See also [`logpdf_in_domain`](@ref) and [`logpdf_in_image`](@ref).
"""
struct TransformDistribution{D <: Union{ContinuousUnivariateDistribution,
                                        ContinuousMultivariateDistribution},
                             T <: Union{UnivariateTransformation,
                                        GroupedTransformation}} <:
                                            TransformationWrapper
    distribution::D
    transformation::T
    function TransformDistribution(distribution::D, transformation::T) where
        {D <: ContinuousUnivariateDistribution, T <: UnivariateTransformation}
        new{D, T}(distribution, transformation)
    end
    function TransformDistribution(distribution::D, transformation::T) where
        {D <: ContinuousMultivariateDistribution, T <: GroupedTransformation}
        @argcheck(length(distribution) == length(transformation),
                  "Incompatible lengths.")
        new{D,T}(distribution, transformation)
    end
end

get_transformation(t::TransformDistribution) = t.transformation

"""
    $SIGNATURES

Return the wrapped distribution.
"""
get_distribution(t::TransformDistribution) = t.distribution

(t::TransformDistribution)(x) = t.transformation(x)

rand(t::TransformDistribution) = t(rand(t.distribution))

@forward TransformDistribution.distribution length

"""
    $SIGNATURES

The log pdf for a transformed distribution at `t(x)` in image, calculated in the
domain without performing the transformation.

The log pdf is adjusted with the log determinant of the Jacobian, ie the
following holds:

```julia logpdf_in_image(t, t(x)) == logpdf_in_domain(t, x) ```

See [`logpdf_in_image`](@ref).

!!! note

    Typical usage of this function would be drawing some random `x`s from the
    contained distribution (possibly also used for some other purpose), and
    obtaining the log pdfs at `t(y)` with the same values.
"""
function logpdf_in_domain(t::TransformDistribution, x)
    # NOTE adjusting by -logjac because of the derivative of the inverse rule
    logpdf(t.distribution, x) - logjac(t.transformation, x)
end

"""
    $SIGNATURES

The log pdf for a transformed distribution at `y` in image.

See also [`logpdf_in_domain`](@ref).
"""
function logpdf_in_image(t::TransformDistribution, y)
    x = inverse(t.transformation, y)
    logpdf_in_domain(t, x)
end

# misc utilities
# NOTE: may be experimental, removed/redesigned
# NOTE: code written for simplicity, not speed, benchmark & optimize later

"""
    $SIGNATURES

Make a container for ungrouped results.

The first argument, `Vector` or `Array`, determines the result of ungrouping.

`representative_elt` is a representative element, used to determine element
type and stucture, and `len` is the length of the result.
"""
function make_ungrouped_container(::Type{Vector},
                                  representative_elt::AbstractArray{T},
                                  len) where T
    [Vector{T}(undef, len) for _ in CartesianIndices(axes(representative_elt))]
end

function make_ungrouped_container(::Type{Array},
                                  representative_elt::AbstractArray{T},
                                  len) where T
    Array{T}(undef, len, size(representative_elt)...)
end

make_ungrouped_container(T, representative_elt::Tuple, len) =
    map(x -> make_ungrouped_container(T, x, len), representative_elt)

make_ungrouped_container(_, representative_elt::T, len) where {T <: Real} =
    Vector{T}(undef, len)

"""
    $SIGNATURES

Save `elt` in `container` at index `i`. `T` determines the layout of
`container`, conformably with [`make_ungrouped_container`](@ref).
"""
ungroup_elt!(T, container::NTuple{N, Any}, elt::NTuple{N, Any}, i) where N =
    foreach((c, e) -> ungroup_elt!(T, c, e, i), container, elt)

function ungroup_elt!(::Type{Vector}, container::AbstractArray,
                      elt::AbstractArray, i)
    @argcheck axes(container) == axes(elt)
    for j in 1:length(elt)
        ungroup_elt!(Vector, container[j], elt[j], i)
    end
end

ungroup_elt!(::Type{Array}, container::AbstractArray,
             elt::AbstractArray{T, N}, i) where {T, N} =
    container[i, fill(Colon(), N)...] .= elt

ungroup_elt!(_, container::AbstractVector, elt::Number, i) = container[i] = elt

"""
    $SIGNATURES

Map the vector `A` elementwise using `f`, and collect the results in arrays or
vectors organized as a single result from `f`. For example, if `f` returns
tuples, this function returns tuples of vectors.

The shape of the result is determined by the first argument, which can be
`Vector` or `Array`. The difference is that `Array` stacks arrays into arrays,
while `Vector` does not.

`Array` may be more useful for summaries along dimensions, eg `mean` and
similar.
"""
function ungrouping_map(T, f, A::AbstractVector)
    n = length(A)
    representative_elt = f(A[1])
    container = make_ungrouped_container(T, representative_elt, n)
    ungroup_elt!(T, container, representative_elt, 1)
    for i in 2:n
        ungroup_elt!(T, container, f(A[i]), i)
    end
    container
end

"""
    $SIGNATURES

Test if a numerical argument is made up of finite real numbers.
"""
all_finite(x::Real) = isfinite(x)

all_finite(A::AbstractArray{T}) where {T <: Real} = all(isfinite, A)

all_finite(t::Tuple) = all(all_finite, t)

all_finite(itr) = all(all_finite, itr)

end


text:
@scalar_rule(one(x), Zero())
@scalar_rule(zero(x), Zero())
@scalar_rule(abs2(x), Wirtinger(x', x))
@scalar_rule(log(x), inv(x))
@scalar_rule(log10(x), inv(x) / log(oftype(x, 10)))
@scalar_rule(log2(x), inv(x) / log(oftype(x, 2)))
@scalar_rule(log1p(x), inv(x + 1))
@scalar_rule(expm1(x), exp(x))

@scalar_rule(sin(x), cos(x))
@scalar_rule(cos(x), -sin(x))
@scalar_rule(sinpi(x), π * cospi(x))
@scalar_rule(cospi(x), -π * sinpi(x))
@scalar_rule(sind(x), (π / oftype(x, 180)) * cosd(x))
@scalar_rule(cosd(x), -(π / oftype(x, 180)) * sind(x))

@scalar_rule(asin(x), inv(sqrt(1 - x^2)))
@scalar_rule(acos(x), -inv(sqrt(1 - x^2)))
@scalar_rule(atan(x), inv(1 + x^2))
@scalar_rule(asec(x::Real), inv(abs(x) * sqrt(x^2 - 1)))
@scalar_rule(asec(x), inv(x^2 * sqrt(1 - x^-2)))
@scalar_rule(acsc(x::Real), -inv(abs(x) * sqrt(x^2 - 1)))
@scalar_rule(acsc(x), -inv(x^2 * sqrt(1 - x^-2)))
@scalar_rule(acot(x), -inv(1 + x^2))

@scalar_rule(asind(x), oftype(x, 180) / π / sqrt(1 - x^2))
@scalar_rule(acosd(x), -oftype(x, 180) / π / sqrt(1 - x^2))
@scalar_rule(atand(x), oftype(x, 180) / π / (1 + x^2))
@scalar_rule(asecd(x::Real), oftype(x, 180) / π / abs(x) / sqrt(x^2 - 1))
@scalar_rule(asecd(x), oftype(x, 180) / π / x^2 / sqrt(1 - x^-2))
@scalar_rule(acscd(x::Real), -oftype(x, 180) / π / abs(x) / sqrt(x^2 - 1))
@scalar_rule(acscd(x), -oftype(x, 180) / π / x^2 / sqrt(1 - x^-2))
@scalar_rule(acotd(x), -oftype(x, 180) / π / (1 + x^2))

@scalar_rule(sinh(x), cosh(x))
@scalar_rule(cosh(x), sinh(x))
@scalar_rule(tanh(x), sech(x)^2)
@scalar_rule(coth(x), -(csch(x)^2))

@scalar_rule(asinh(x), inv(sqrt(x^2 + 1)))
@scalar_rule(acosh(x), inv(sqrt(x^2 - 1)))
@scalar_rule(atanh(x), inv(1 - x^2))
@scalar_rule(asech(x), -inv(x * sqrt(1 - x^2)))
@scalar_rule(acsch(x::Real), -inv(abs(x) * sqrt(1 + x^2)))
@scalar_rule(acsch(x), -inv(x^2 * sqrt(1 + x^-2)))
@scalar_rule(acoth(x), inv(1 - x^2))

@scalar_rule(deg2rad(x), π / oftype(x, 180))
@scalar_rule(rad2deg(x), oftype(x, 180) / π)

@scalar_rule(conj(x), Wirtinger(Zero(), One()))
@scalar_rule(adjoint(x), Wirtinger(Zero(), One()))
@scalar_rule(transpose(x), One())

@scalar_rule(abs(x::Real), sign(x))
@scalar_rule(abs(x::Complex), Wirtinger(x' / 2Ω, x / 2Ω))
@scalar_rule(hypot(x::Real), sign(x))
@scalar_rule(hypot(x::Complex), Wirtinger(x' / 2Ω, x / 2Ω))
@scalar_rule(rem2pi(x, r::RoundingMode), (One(), DNE()))

@scalar_rule(+(x), One())
@scalar_rule(-(x), -1)
@scalar_rule(+(x, y), (One(), One()))
@scalar_rule(-(x, y), (One(), -1))
@scalar_rule(/(x, y), (inv(y), -(x / y / y)))
@scalar_rule(\(x, y), (-(y / x / x), inv(x)))
@scalar_rule(^(x, y), (y * x^(y - 1), Ω * log(x)))

@scalar_rule(inv(x), -Ω^2)
@scalar_rule(sqrt(x), inv(2 * Ω))
@scalar_rule(cbrt(x), inv(3 * Ω^2))
@scalar_rule(exp(x), Ω)
@scalar_rule(exp2(x), Ω * log(oftype(x, 2)))
@scalar_rule(exp10(x), Ω * log(oftype(x, 10)))

@scalar_rule(tan(x), 1 + Ω^2)
@scalar_rule(sec(x), Ω * tan(x))
@scalar_rule(csc(x), -Ω * cot(x))
@scalar_rule(cot(x), -(1 + Ω^2))
@scalar_rule(tand(x), (π / oftype(x, 180)) * (1 + Ω^2))
@scalar_rule(secd(x), (π / oftype(x, 180)) * Ω * tand(x))
@scalar_rule(cscd(x), -(π / oftype(x, 180)) * Ω * cotd(x))
@scalar_rule(cotd(x), -(π / oftype(x, 180)) * (1 + Ω^2))
@scalar_rule(sech(x), -tanh(x) * Ω)
@scalar_rule(csch(x), -coth(x) * Ω)

@scalar_rule(hypot(x::Real, y::Real), (x / Ω, y / Ω))
@scalar_rule(sincos(x), @setup((sinx, cosx) = Ω), cosx, -sinx)
@scalar_rule(atan(y, x), @setup(u = x^2 + y^2), (x / u, -y / u))

@scalar_rule(max(x, y), @setup(gt = x > y), (gt, !gt))
@scalar_rule(min(x, y), @setup(gt = x > y), (!gt, gt))
@scalar_rule(mod(x, y), @setup((u, nan) = promote(x / y, NaN16)),
             (ifelse(isint, nan, one(u)), ifelse(isint, nan, -floor(u))))
@scalar_rule(rem(x, y), @setup((u, nan) = promote(x / y, NaN16)),
             (ifelse(isint, nan, one(u)), ifelse(isint, nan, -trunc(u))))
@scalar_rule(angle(x::Complex), @setup(u = abs2(x)), Wirtinger(-im//2 * x' / u, im//2 * x / u))
@scalar_rule(angle(x::Real), Zero())
@scalar_rule(real(x::Complex), Wirtinger(1//2, 1//2))
@scalar_rule(real(x::Real), One())
@scalar_rule(imag(x::Complex), Wirtinger(-im//2, im//2))
@scalar_rule(imag(x::Real), Zero())

# product rule requires special care for arguments where `mul` is non-commutative

function frule(::typeof(*), x::Number, y::Number)
    function times_pushforward(_, Δx, Δy)
        return Δx * y + x * Δy
    end
    return x * y, times_pushforward
end

function rrule(::typeof(*), x::Number, y::Number)
    function times_pullback(ΔΩ)
        return (NO_FIELDS,  @thunk(ΔΩ * y'), @thunk(x' * ΔΩ))
    end
    return x * y, times_pullback
end

function frule(::typeof(identity), x)
    function identity_pushforward(_, ẏ)
        return ẏ
    end
    return x, identity_pushforward
end

function rrule(::typeof(identity), x)
    function identity_pullback(ȳ)
        return (NO_FIELDS, ȳ)
    end
    return x, identity_pullback
end


text:
export DGKS, ClassicalGramSchmidt, ModifiedGramSchmidt
export orthogonalize_and_normalize!

abstract type OrthogonalizationMethod end
struct DGKS <: OrthogonalizationMethod end
struct ClassicalGramSchmidt <: OrthogonalizationMethod end
struct ModifiedGramSchmidt <: OrthogonalizationMethod end

# Default to MGS, good enough for solving linear systems.
@inline orthogonalize_and_normalize!(V::StridedMatrix{T}, w::StridedVector{T}, h::StridedVector{T}) where {T} = orthogonalize_and_normalize!(V, w, h, ModifiedGramSchmidt)

function orthogonalize_and_normalize!(V::StridedMatrix{T}, w::StridedVector{T}, h::StridedVector{T}, ::Type{DGKS}) where {T}
    # Orthogonalize using BLAS-2 ops
    mul!(h, adjoint(V), w)
    BLAS.gemv!('N', -one(T), V, h, one(T), w)
    nrm = norm(w)

    # Constant used by ARPACK.
    η = one(real(T)) / √2

    projection_size = norm(h)

    # Repeat as long as the DGKS condition is satisfied
    # Typically this condition is true only once.
    while nrm < η * projection_size
        correction = adjoint(V)*w
        projection_size = norm(correction)
        # w = w - V * correction
        BLAS.gemv!('N', -one(T), V, correction, one(T), w)
        h .+= correction
        nrm = norm(w)
    end

    # Normalize; note that we already have norm(w).
    w .*= inv(nrm)

    nrm
end

function orthogonalize_and_normalize!(V::StridedMatrix{T}, w::StridedVector{T}, h::StridedVector{T}, ::Type{ClassicalGramSchmidt}) where {T}
    # Orthogonalize using BLAS-2 ops
    mul!(h, adjoint(V), w)
    BLAS.gemv!('N', -one(T), V, h, one(T), w)
    nrm = norm(w)

    # Normalize
    w .*= inv(nrm)

    nrm
end

function orthogonalize_and_normalize!(V::StridedVector{Vector{T}}, w::StridedVector{T}, h::StridedVector{T}, ::Type{ModifiedGramSchmidt}) where {T}
    # Orthogonalize using BLAS-1 ops
    for i = 1 : length(V)
        h[i] = dot(V[i], w)
        w .-= h[i] .* V[i]
    end

    # Normalize
    nrm = norm(w)
    w .*= inv(nrm)

    nrm
end

function orthogonalize_and_normalize!(V::StridedMatrix{T}, w::StridedVector{T}, h::StridedVector{T}, ::Type{ModifiedGramSchmidt}) where {T}
    # Orthogonalize using BLAS-1 ops and column views.
    for i = 1 : size(V, 2)
        column = view(V, :, i)
        h[i] = dot(column, w)
        w .-= h[i] .* column
    end

    nrm = norm(w)
    w .*= inv(nrm)

    nrm
end


text:
module VUMPS

include("cuda_patch.jl")
include("environment.jl")
include("fixedpoint.jl")
include("vumpsruntime.jl")
include("autodiff.jl")

end


text:
# Implements linear fitting (analagous to IDL's regress function)

# Dependent in the vector Y, which should have length N_data
# N_var independent variables are contained in the array X.
# X should have the size [N_var,N_data].

function regress(X,Y,sigY)
# Check that X has two dimensions, and that
# X and Y have are compatible (second dimension of X
# should equal the dimension of Y):
nx = ndims(X)
sx = size(X)
sy = size(Y)
if (nx == 2) & (sx[2] == sy[1]) then
# Set up array for computing coefficients:
  A = zeros(sx[1],sx[1])
  B = zeros(sx[1])
# Compute A = \sum_i \sum_j \sum_k X_{i,k} X_{j,k}/\sigY_k^2
  for i = 1:sx[1]
    for j = 1:sx[1]
      for k = 1:sx[2]
        A[i,j] += X[i,k]*X[j,k]/sigY[k]^2
      end
    end
# Compute B = \sum_i \sum_k X_{i,k} Y_k/\sigY_k^2
    for k = 1:sx[2]
      B[i] += X[i,k]*Y[k]/sigY[k]^2
    end
  end
# Covariance:
  covC = inv(A)
# Now, solve for the coefficients, A*C = B:
  C = \(A,B)
else
  println("Error")
end
return C,covC
end