text:
export substitute, SubstituteHandler

import Base: show, push!

mutable struct SubstituteActor{L, R, F, A} <: Actor{L}
    mapFn   :: F
    actor   :: A
    pending :: Union{Nothing, L}
    current :: Union{Nothing, R}
    handler
end

"""
    SubstituteHandler()

handler used to `release!` new values in `substitute` operator.

See also: [`substitute`](@ref)
"""
struct SubstituteHandler
    list :: List{SubstituteActor}

    SubstituteHandler() = new(List(SubstituteActor))
end

Base.push!(handler::SubstituteHandler, actor) = pushnode!(handler.list, actor)

release!(handler::SubstituteHandler) = foreach(release!, handler.list)

"""
    substitute(::Type{T}, mapFn::F, handler::SubstituteHandler) where { T, F <: Function }

This operator forces observable to substitute each emmited value with the latest computed value with the corresponding `handler` and `release!` function.
After `release!` on `handler` `substitute` operator computes new value with `mapFn` but does not emit it until next emission from source observable.
Always calls `mapFn` on first value from source observable.

# Producing

Stream of type `<: Subscribable{T}`

# Examples
```jldoctest
using Rocket

subject = Subject(Int)

handler = SubstituteHandler()
source  = subject |> substitute(String, i -> string("i = ", i), handler)

subscription = subscribe!(source, logger())

next!(subject, 1)
next!(subject, 2)
next!(subject, 3)

release!(handler)

next!(subject, 4)
next!(subject, 5)
next!(subject, 6)

unsubscribe!(subscription)
;

# output

[LogActor] Data: i = 1
[LogActor] Data: i = 1
[LogActor] Data: i = 1
[LogActor] Data: i = 3
[LogActor] Data: i = 3
[LogActor] Data: i = 3
```

See also: [`SubstituteHandler`](@ref)
"""
substitute(::Type{R}, mapFn::F, handler::SubstituteHandler) where { R, F <: Function } = SubstituteOperator{R, F}(mapFn, handler)

struct SubstituteOperator{R, F} <: RightTypedOperator{R}
    mapFn   :: F
    handler :: SubstituteHandler
end

function on_call!(::Type{L}, ::Type{R}, operator::SubstituteOperator{R, F}, source) where { L, R, F }
    return proxy(R, source, SubstituteProxy{L, R, F}(operator.mapFn, operator.handler))
end

operator_right(::SubstituteOperator{R}, ::Type{L}) where { L, R } = R

struct SubstituteProxy{L, R, F} <: ActorSourceProxy
    mapFn   :: F
    handler :: SubstituteHandler
end

actor_proxy!(::Type{R}, proxy::SubstituteProxy{L, R, F}, actor::A) where { L, R, F, A } = SubstituteActor{L, R, F, A}(proxy.mapFn, actor, nothing, nothing, nothing)

function release!(actor::SubstituteActor)
    if actor.pending !== nothing
        actor.current = actor.mapFn(actor.pending)
    end
    return nothing
end

function on_next!(actor::SubstituteActor{L}, data::L) where L
    actor.pending = data
    if actor.current === nothing
        release!(actor)
    end
    if actor.current !== nothing
        next!(actor.actor, actor.current)
    end
end

function on_error!(actor::SubstituteActor, err)
    remove(actor.handler)
    error!(actor.actor, err)
end

function on_complete!(actor::SubstituteActor)
    remove(actor.handler)
    complete!(actor.actor)
end

##

@subscribable struct SubstituteSource{L, S} <: Subscribable{L}
    source  :: S
    handler :: SubstituteHandler
end

source_proxy!(::Type{R}, proxy::SubstituteProxy{L, R}, source::S) where { L, R, S } = SubstituteSource{L, S}(source, proxy.handler)

function on_subscribe!(source::SubstituteSource, actor::SubstituteActor)
    handler       = push!(source.handler, actor)
    actor.handler = handler
    subscription  = subscribe!(source.source, actor)
    return SubstituteSubscription(handler, subscription)
end

struct SubstituteSubscription{H, S} <: Teardown
    handler      :: H
    subscription :: S
end

as_teardown(::Type{ <: SubstituteSubscription }) = UnsubscribableTeardownLogic()

function on_unsubscribe!(subscription::SubstituteSubscription)
    remove(subscription.handler)
    return unsubscribe!(subscription.subscription)
end

Base.show(io::IO, ::SubstituteOperator{R})   where {    R } = print(io, "SubstituteOperator($R)")
Base.show(io::IO, ::SubstituteProxy{L, R})   where { L, R } = print(io, "SubstituteProxy($L, $R)")
Base.show(io::IO, ::SubstituteActor{L, R})   where { L, R } = print(io, "SubstituteActor($L -> $R)")
Base.show(io::IO, ::SubstituteSource{R})     where R        = print(io, "SubstituteSource($R)")
Base.show(io::IO, ::SubstituteSubscription)                 = print(io, "SubstituteSubscription()")


text:
using Printf
using Turing
using CairoMakie
using ElectronDisplay

# Data generation

N = 100
x = range(0, 1, length=N)
y = @. 0.16x^2 + 4.5x + π
y = @. y * (1 + 0.1 * randn())

scatter(x, y)

# Linear model (y = mx + c)

@model function linear_model(x, y)
    m ~ Normal(0, 1)
    c ~ Normal(0, 1)
    σ ~ Exponential(1)

    N = length(x)
    for n in 1:N
        y[n] ~ Normal(m*x[n] + c, σ)
    end
end

chain_linear = sample(linear_model(x, y), NUTS(0.65), 10000, progress=true)

function y_linear(chain, n, x)
    m = chain[:m][n]
    c = chain[:c][n]
    return @. m * x + c
end

function animate_chain(chain, y_model, x, y; filepath, N_samples=1000, N_history=50, max_alpha=0.5)
    n = Observable(1)

    # Compute all predictions to plot.
    yₚ = []
    for i in 1:N_samples
        push!(yₚ, y_model(chain, i, x))
    end

    fig = Figure()

    title = @lift "Iteration $($n)"
    ax = Axis(fig[1, 1], xlabel="x", ylabel="y"; title)

    # Create N_history lines with different alpha that we can change the data on.
    line_plots = []
    for i in 1:N_history
        line_plot = lines!(x, yₚ[1], color=(:gray, max_alpha * i/N_history))
        push!(line_plots, line_plot)
    end

    # Plot the actual data points.
    scatter!(x, y)

    record(fig, filepath, 1:N_samples, framerate=30) do chain_iter
        @info "Animating $filepath: frame $chain_iter/$N_samples"
        n[] = chain_iter

        # Change the data of the line plots to show the last N_history predictions.
        for (i, j) in enumerate(chain_iter:-1:max(1, chain_iter-N_history+1))
            line_plots[i][2] = yₚ[j]
        end
    end

    return
end

animate_chain(chain_linear, y_linear, x, y, filepath="linear_model.mp4")

function plot_chains_and_densities(chain, vars; filepath)
    fig = Figure(resolution=(800, 200*length(vars)))

    for (n, var) in enumerate(vars)
        ax1 = Axis(fig[n, 1], ylabel="$var", xlabel="iteration")
        lines!(ax1, chain[var][:] |> skipmissing |> collect)

        ax2 = Axis(fig[n, 2], ylabel="P($var)", xlabel="$var", yticklabelsvisible=false)
        density!(ax2, chain[var][:] |> skipmissing |> collect)
    end

    save(filepath, fig, px_per_unit=2)

    return fig
end

plot_chains_and_densities(chain_linear, [:m, :c, :σ], filepath="linear_model_posteriors.png")

# Quadratic model (y = ax² + bx + c)

@model function quadratic_model(x, y)
    a ~ Normal(0, 1)
    b ~ Normal(0, 1)
    c ~ Normal(0, 1)
    σ ~ Exponential(1)

    N = length(x)
    for n in 1:N
        y[n] ~ Normal(a*x[n]^2 + b*x[n] + c, σ)
    end
end

chain_quadratic = sample(quadratic_model(x, y), NUTS(0.65), 10000, progress=true)

function y_quadratic(chain, n, x)
    a = chain[:a][n]
    b = chain[:b][n]
    c = chain[:c][n]
    return @. a*x^2 + b*x + c
end

animate_chain(chain_quadratic, y_quadratic, x, y, filepath="quadratic_model.mp4")

plot_chains_and_densities(chain_quadratic, [:a, :b, :c, :σ], filepath="quadratic_model_posteriors.png")

# Exponential model [y = α + β * exp(γ * x)]

@model function exponential_model(x, y)
    α ~ Normal(0, 1)
    β ~ Normal(0, 1)
    γ ~ Normal(0, 1)
    σ ~ Exponential(1)

    N = length(x)
    for n in 1:N
        y[n] ~ Normal(α + β * exp(γ*x[n]), σ)
    end
end

chain_exponential = sample(exponential_model(x, y), NUTS(0.65), 10000, progress=true)

function y_exponential(chain, n, x)
    α = chain[:α][n]
    β = chain[:β][n]
    γ = chain[:γ][n]
    return @. α + β * exp(γ * x)
end

animate_chain(chain_exponential, y_exponential, x, y, filepath="exponential_model.mp4")

plot_chains_and_densities(chain_exponential, [:α, :β, :γ, :σ], filepath="exponential_model_posteriors.png")

# Categorical model (picks between linear, quadratic, and exponential)

@model function categorical_model1(x, y)
    N = length(x)

    pl ~ Beta(1, 1)
    M ~ Categorical([pl, 1-pl])

    if M == 1
        m ~ Normal(0, 1)
        k ~ Normal(0, 1)
        σl ~ Exponential(1)

        for n in 1:N
            y[n] ~ Normal(m*x[n] + k, σl)
        end
    else
        a ~ Normal(0, 1)
        b ~ Normal(0, 1)
        c ~ Normal(0, 1)
        σq ~ Exponential(1)

        for n in 1:N
            y[n] ~ Normal(a*x[n]^2 + b*x[n] + c, σq)
        end
    end
end

prior_chain_categorical1 = sample(categorical_model1(x, y), Prior(), 10000)
plot_chains_and_densities(prior_chain_categorical1, [:pl, :M, :a, :b, :c, :m, :k], filepath="categorical_model1_priors.png")

# Errors since variables defined in if statements are missing during inference?
sampler = Gibbs(PG(100, :M), NUTS(100, 0.65, :pl, :m, :k, :σl, :a, :b, :c, :σq))
chain_categorical1 = sample(categorical_model1(x, y), sampler, 1000, progress=true)

@model function categorical_model2(x, y)
    N = length(x)

    pl ~ Beta(1, 1)
    M ~ Categorical([pl, 1-pl])

    m ~ Normal(0, 1)
    k ~ Normal(0, 1)
    σl ~ Exponential(1)

    a ~ Normal(0, 1)
    b ~ Normal(0, 1)
    c ~ Normal(0, 1)
    σq ~ Exponential(1)

    if M == 1
        for n in 1:N
            y[n] ~ Normal(m*x[n] + k, σl)
        end
    else
        for n in 1:N
            y[n] ~ Normal(a*x[n]^2 + b*x[n] + c, σq)
        end
    end
end

prior_chain_categorical2 = sample(categorical_model2(x, y), Prior(), 10000)
plot_chains_and_densities(prior_chain_categorical2, [:pl, :M, :a, :b, :c, :m, :k], filepath="categorical_model2_priors.png")

# Inference fails? since it always picks M=1 or M=2 and doesn't sample from both models.
sampler = Gibbs(PG(10, :M), NUTS(100, 0.65, :pl, :m, :k, :σl, :a, :b, :c, :σq))
chain_categorical2 = sample(categorical_model2(x, y), sampler, 1000, progress=true)
plot_chains_and_densities(chain_categorical2, [:pl, :M, :m, :k, :σl, :a, :b, :c, :σq], filepath="categorical_model2_posteriors.png")

@model function categorical_model3(x, y)
    N = length(x)

    m ~ Normal(0, 1)
    k ~ Normal(0, 1)
    σl ~ Exponential(1)

    a ~ Normal(0, 1)
    b ~ Normal(0, 1)
    c ~ Normal(0, 1)
    σq ~ Exponential(1)

    pl ~ Beta(1, 1)

    if rand() < pl
        for n in 1:N
            y[n] ~ Normal(m*x[n] + k, σl)
        end
    else
        for n in 1:N
            y[n] ~ Normal(a*x[n]^2 + b*x[n] + c, σq)
        end
    end
end

# Inference is bad and chain is very correlated.
chain_categorical3 = sample(categorical_model3(x, y), NUTS(), 1000, progress=true)
plot_chains_and_densities(chain_categorical3, [:pl, :m, :k, :σl, :a, :b, :c, :σq], filepath="categorical_model3_posteriors.png")

function plot_model_comparison(chain, x, y)
    fig = Figure()
    ax = Axis(fig[1, 1])

    pl = mode(chain[:pl])

    m = mode(chain[:m])
    k = mode(chain[:k])

    a = mode(chain[:a])
    b = mode(chain[:b])
    c = mode(chain[:c])

    y_linear = @. m*x + k
    y_quadratic = @. a*x^2 + b*x + c

    lines!(ax, x, y_linear, label=@sprintf("linear (%d%%)", 100pl))
    lines!(ax, x, y_quadratic, label=@sprintf("quadratic (%d%%)", 100*(1-pl)))

    scatter!(x, y)
    axislegend(ax)

    return fig
end

plot_model_comparison(chain_categorical3, x, y)

chain_time(chain) = chain.info.stop_time - chain.info.start_time

chains = (chain_linear, chain_quadratic, chain_exponential, chain_categorical2, chain_categorical3)
chain_names = ("linear", "quadratic", "exponential", "categorical1", "categorical2")

open("inference_wall_clock_times.txt","w") do io
    for (chain, name) in zip(chains, chain_names)
        line = @sprintf("%s: %.1f seconds for %d samples", name, chain_time(chain), size(chain, 1))
        @info line
        println(io, line)
    end
end

ext: jl


text:
for (elty, relty, ext) in ((:Float32, :Float32, :s),
                           (:Float64, :Float64, :d),
                           (:Complex64, :Float32, :c),
                           (:Complex128, :Float64, :z))

    # Distributed sparse gemv
    for (trans, elenum) in (("", :NORMAL), ("t", :TRANSPOSE), ("c", :ADJOINT))

        f = Symbol("A", trans, "_mul_B!")

        @eval begin
            function ($f)(α::$elty, A::DistSparseMatrix{$elty}, x::DistMultiVec{$elty}, β::$elty, y::DistMultiVec{$elty})
                err = ccall(($(string("ElMultiplyDist_", ext)), libEl), Cuint,
                    (Cint, $elty, Ptr{Void}, Ptr{Void}, $elty, Ptr{Void}),
                    $elenum, α, A.obj, x.obj, β, y.obj)
                err == 0 || throw(ElError(err))
                return y
            end
        end
    end
end


text:
# Create summary for Summary_value (default constructor by
# protobuf is broken).
SummaryCollection(;kwargs...) = Summary(value=Vector{Summary_Value}(); kwargs...)
SummaryCollection(summaries::Vector{Summary_Value}; kwargs...) = Summary(value=summaries; kwargs...)
SummaryCollection(summary::Summary_Value; kwargs...) = Summary(value=[summary]; kwargs...)

function make_event(logger::TBLogger, summary::Summary;
                    step::Int=TensorBoardLogger.step(logger))
    return Event(wall_time=time(), summary=summary, step=step)
end

function write_event(file::IOStream, event::Event)
    data = PipeBuffer();
    writeproto(data, event)

    #header
    header     = collect(reinterpret(UInt8, [data.size]))
    crc_header = reinterpret(UInt8, UInt32[masked_crc32c(header)])
    crc_data   = reinterpret(UInt8, UInt32[masked_crc32c(data.data)])

    write(file, header)
    write(file, crc_header)
    write(file, data.data)
    write(file, crc_data)
    flush(file)
end

write_event(logger::TBLogger, event::Event) = write_event(logger.file, event)


text:
"""
constructDormandPrince()

Constructs the tableau object for the Dormand-Prince Order 4/5 method.
"""
function constructDormandPrince(T::Type = Float64)
  A = [0 0 0 0 0 0 0
      1//5 0 0 0 0 0 0
      3//40 9//40 0 0 0 0 0
      44//45 -56//15 32//9 0 0 0 0
      19372//6561 -25360//2187 64448//6561 -212//729 0 0 0
      9017//3168 -355//33 46732//5247 49//176 -5103//18656 0 0
      35//384 0 500//1113 125//192 -2187//6784 11//84 0]
  c = [0;1//5;3//10;4//5;8//9;1;1]
  α = [35//384;0;500//1113;125//192;-2187//6784;11//84;0]
  αEEst = [5179//57600;0;7571//16695;393//640;-92097//339200;187//2100;1//40]
  A = map(T,A)
  α = map(T,α)
  αEEst = map(T,αEEst)
  c = map(T,c)
  return(ExplicitRKTableau(A,c,α,5,αEEst=αEEst,adaptiveorder=4,fsal=true))
end

"""
ODE_DEFAULT_TABLEAU

Sets the default tableau for the ODE solver. Currently Dormand-Prince 4/5.
"""
const ODE_DEFAULT_TABLEAU = constructDormandPrince()

@inline ODE_DEFAULT_NORM(u) = sqrt(sum(abs2,u) / length(u))
@inline ODE_DEFAULT_NORM(u::Number) = abs(u)
@inline ODE_DEFAULT_ISOUTOFDOMAIN(t,u) = false
@inline ODE_DEFAULT_PROG_MESSAGE(dt,t,u) = "dt="*string(dt)*"\nt="*string(t)*"\nmax u="*string(maximum(abs.(u)))
@inline ODE_DEFAULT_UNSTABLE_CHECK(dt,t,u) = any(isnan,u)


text:
export TRCFile

DatasetManager.srcname_default(::Type{Source{C3DFile}}) = "c3d"
DatasetManager.srcext(::Type{Source{C3DFile}}) = ".c3d"

function DatasetManager.readsource(s::Source{C3DFile}; kwargs...)
    return readc3d(sourcepath(s); kwargs...)
end

struct TRCFile; end

function DatasetManager.readsource(s::Source{TRCFile}; select=nothing, drop=["Frame#"])
    colnames = open(sourcepath(s)) do io
        readline(io); readline(io); readline(io)
        ln = readline(io)
        split(ln, r"[,\s]")
    end

    trc = CSV.read(sourcepath(s), DataFrame; header=false, datarow=7)
    while (length(colnames)-2) % 3 !== 0
        pop!(colnames)
    end
    mkrs = Dict("Frame#" => copy(trc[!, 1]),
        "Time" => copy(trc[!, 2]),
        ( string(colnames[i]) => Matrix(trc[!, i:i+2]) for i in 3:3:lastindex(colnames) )...)

    !isnothing(select) && delete!.(Ref(mkrs), setdiff(select, keys(mkrs)))
    !isnothing(drop) && delete!.(Ref(mkrs), drop)

    return mkrs
end

DatasetManager.dependencies(::Type{Source{TRCFile}}) = (Source{C3DFile},)
DatasetManager.srcname_default(::Type{Source{TRCFile}}) = "trc"
DatasetManager.srcext(::Type{Source{TRCFile}}) = ".trc"

function DatasetManager.generatesource(
    trial,
    src::Source{TRCFile},
    deps;
    f! = identity,
    kwargs...
)
    srcpath = sourcepath(src)
    mkpath(dirname(srcpath))

    c3dsrc = getsource(trial, only(filter(x-> Source{C3DFile} === x ||
        (x isa Pair && Source{C3DFile} ∈ x), deps)))
    c3dfl = readsource(c3dsrc)
    f!(c3dfl)
    trcpath = splitext(srcpath)[2] == ".trc" ? srcpath : srcpath*".trc"
    if c3dfl.groups[:POINT][String, :UNITS] == "mm"
        precision=2
    else
        precision=5
    end
    writetrc(trcpath, c3dfl; precision, kwargs...)

    return Source{TRCFile}(trcpath)
end


text:
using ROCKS
using Test

@testset "ROC, KS" begin
    @testset "ranks binning" begin
        x = [1, 1, 1, 2, 2, 2, 2, 2, 3, 3]
        xm = [1, 1, missing, 2, 2, 2, 3, 3, 3, 3]
        @test ranks(x) == [1, 1, 1, 5, 5, 5, 5, 5, 8, 8]
        @test ranks(x, rev=true) == [8, 8, 8, 4, 4, 4, 4, 4, 1, 1]
        @test ranks(x, rank=StatsBase.ordinalrank) == [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
        @test ranks(x, rank=StatsBase.competerank) == [0, 0, 0, 3, 3, 3, 3, 3, 8, 8]
        @test ranks(x, rank=StatsBase.denserank) == [0, 0, 0, 1, 1, 1, 1, 1, 2, 2]
        @test isequal(ranks(xm), [1, 1, missing, 4, 4, 4, 7, 7, 7, 7])
    end

    @testset "ks roc" begin
        t = [0, 0, 0, 1, 1, 1]
        x = [1, 1, 1, 2, 2, 2]

        @test kstest(t, x) ==
              (n = 6, n1 = 3, n0 = 3, baserate = 0.5, ks = 1.0, ksarg = 2, ksdep = 0.5)
        @test roc(t, x) == (conc = 9, tied = 0, disc = 0, auc = 1.0, gini = 1.0)
    end
end


text:
immutable Objective1 <: Cut{2}; groups; end
immutable Objective2 <: Cut{2}; groups; end
immutable Objective3 <: Cut{2}; groups; end
immutable Objective4 <: Cut{2}; groups; end
export Objective1, Objective2, Objective3, Objective4

function getcut(prob::Problem, vars::Variables, ::Objective1, i, j)
	Lˣ, Lʸ = prob.W, prob.H
	ℓˣ, ℓʸ = vars.ℓˣ, vars.ℓʸ
	zˣ, zʸ = vars.zˣ, vars.zʸ
	dˣ, dʸ = vars.dˣ, vars.dʸ
	@LinearConstraints(begin
		dˣ[i,j] ≥ 0.5*(ℓˣ[i]+ℓˣ[j]) - Lˣ*(1-zˣ[i,j]-zˣ[j,i])
		dʸ[i,j] ≥ 0.5*(ℓʸ[i]+ℓʸ[j]) - Lʸ*(1-zʸ[i,j]-zʸ[j,i])
	end)
end

function getcut(prob::Problem, vars::Variables, ::Objective2, i, j)
	Lˣ,  Lʸ =  prob.W,   prob.H
	lbˣ, lbʸ = prob.wlb, prob.hlb
	cˣ, cʸ = vars.cˣ, vars.cʸ
	ℓˣ, ℓʸ = vars.ℓˣ, vars.ℓʸ
	zˣ, zʸ = vars.zˣ, vars.zʸ
	dˣ, dʸ = vars.dˣ, vars.dʸ
	@LinearConstraints(begin
		dˣ[i,j] ≥ cˣ[i] - cˣ[j] + ℓˣ[i] + lbˣ[j]*(zˣ[i,j]+zˣ[j,i]) - Lˣ*(1-zˣ[i,j])
		dˣ[i,j] ≥ cˣ[i] - cˣ[j] + ℓˣ[j] + lbˣ[i]*(zˣ[i,j]+zˣ[j,i]) - Lˣ*(1-zˣ[i,j])
		dˣ[i,j] ≥ cˣ[j] - cˣ[i] + ℓˣ[i] + lbˣ[j]*(zˣ[i,j]+zˣ[j,i]) - Lˣ*(1-zˣ[j,i])
		dˣ[i,j] ≥ cˣ[j] - cˣ[i] + ℓˣ[j] + lbˣ[i]*(zˣ[i,j]+zˣ[j,i]) - Lˣ*(1-zˣ[j,i])
		dʸ[i,j] ≥ cʸ[i] - cʸ[j] + ℓʸ[i] + lbʸ[j]*(zʸ[i,j]+zʸ[j,i]) - Lʸ*(1-zʸ[i,j])
		dʸ[i,j] ≥ cʸ[i] - cʸ[j] + ℓʸ[j] + lbʸ[i]*(zʸ[i,j]+zʸ[j,i]) - Lʸ*(1-zʸ[i,j])
		dʸ[i,j] ≥ cʸ[j] - cʸ[i] + ℓʸ[i] + lbʸ[j]*(zʸ[i,j]+zʸ[j,i]) - Lʸ*(1-zʸ[j,i])
		dʸ[i,j] ≥ cʸ[j] - cʸ[i] + ℓʸ[j] + lbʸ[i]*(zʸ[i,j]+zʸ[j,i]) - Lʸ*(1-zʸ[j,i])
	end)
end

function getcut(prob::Problem, vars::Variables, ::Objective3, i, j)
	lbˣ, lbʸ = prob.wlb, prob.hlb
	cˣ, cʸ = vars.cˣ, vars.cʸ
	zˣ, zʸ = vars.zˣ, vars.zʸ
	dˣ, dʸ = vars.dˣ, vars.dʸ
	@LinearConstraints(begin
		dˣ[i,j] ≥ cˣ[i] - cˣ[j] + (lbˣ[i]+lbˣ[j])*zˣ[i,j]
		dˣ[i,j] ≥ cˣ[j] - cˣ[i] + (lbˣ[i]+lbˣ[j])*zˣ[j,i]
		dʸ[i,j] ≥ cʸ[i] - cʸ[j] + (lbʸ[i]+lbʸ[j])*zʸ[i,j]
		dʸ[i,j] ≥ cʸ[j] - cʸ[i] + (lbʸ[i]+lbʸ[j])*zʸ[j,i]
	end)
end

function getcut(prob::Problem, vars::Variables, ::Objective4, i, j)
	Lˣ,  Lʸ =  prob.W,   prob.H
	lbˣ, lbʸ = prob.wlb, prob.hlb
	ℓˣ, ℓʸ = vars.ℓˣ, vars.ℓʸ
	zˣ, zʸ = vars.zˣ, vars.zʸ
	dˣ, dʸ = vars.dˣ, vars.dʸ
	@LinearConstraints(begin
		2dˣ[i,j] ≥ ℓˣ[i] - Lˣ*(1-zˣ[i,j]-zˣ[j,i]) + lbˣ[j]*(zˣ[i,j]+zˣ[j,i])
		2dˣ[i,j] ≥ ℓˣ[j] - Lˣ*(1-zˣ[i,j]-zˣ[j,i]) + lbˣ[i]*(zˣ[i,j]+zˣ[j,i])
	end)
end

immutable ThreeBoxObjective1 <: Cut{3}; groups; end
immutable ThreeBoxObjective2 <: Cut{3}; groups; end
immutable ThreeBoxObjective3 <: Cut{3}; groups; end
immutable ThreeBoxObjective4 <: Cut{3}; groups; end
export ThreeBoxObjective1, ThreeBoxObjective2, ThreeBoxObjective3, ThreeBoxObjective4

function getcut(prob::Problem, vars::Variables, ::ThreeBoxObjective1, ii, jj, kk)
	Lˣ, Lʸ = prob.W, prob.H
	lbˣ, lbʸ = prob.wlb, prob.hlb
	ℓˣ, ℓʸ = vars.ℓˣ, vars.ℓʸ
	zˣ, zʸ = vars.zˣ, vars.zʸ
	dˣ, dʸ = vars.dˣ, vars.dʸ
	cuts = LinearConstraint[]
	for (_i,_j) in combinations([ii,jj,kk],2) # the two outer ones
		for (i,j) in permutations([_i,_j]) # get both orderings
			k = setdiff([ii,jj,kk],[i,j])[1]
			append!(cuts, @LinearConstraints(begin
				dˣ[i,j] ≥ 0.5*(ℓˣ[i]+ℓˣ[j]) - Lˣ*(1-zˣ[i,j]-zˣ[j,i]) + lbˣ[k]*(zˣ[i,k]+zˣ[k,j]-1)
				dʸ[i,j] ≥ 0.5*(ℓʸ[i]+ℓʸ[j]) - Lʸ*(1-zʸ[i,j]-zʸ[j,i]) + lbʸ[k]*(zʸ[i,k]+zʸ[k,j]-1)
			end))
		end
	end
	cuts
end

function getcut(prob::Problem, vars::Variables, ::ThreeBoxObjective2, ii, jj, kk)
	Lˣ,  Lʸ =  prob.W,   prob.H
	lbˣ, lbʸ = prob.wlb, prob.hlb
	cˣ, cʸ = vars.cˣ, vars.cʸ
	ℓˣ, ℓʸ = vars.ℓˣ, vars.ℓʸ
	zˣ, zʸ = vars.zˣ, vars.zʸ
	dˣ, dʸ = vars.dˣ, vars.dʸ
	cuts = LinearConstraint[]
	for (_i,_j) in combinations([ii,jj,kk],2) # the two outer ones
		for (i,j) in permutations([_i,_j]) # get both orderings
			k = setdiff([ii,jj,kk],[i,j])[1]
			append!(cuts, @LinearConstraints(begin
				dˣ[i,j] ≥ cˣ[i] - cˣ[j] + ℓˣ[i] + lbˣ[j]*(zˣ[i,j]+zˣ[j,i]) - Lˣ*(1-zˣ[i,j]) + lbˣ[k]*(zˣ[i,k]+zˣ[k,j]-1)
				dˣ[i,j] ≥ cˣ[i] - cˣ[j] + ℓˣ[j] + lbˣ[i]*(zˣ[i,j]+zˣ[j,i]) - Lˣ*(1-zˣ[i,j]) + lbˣ[k]*(zˣ[i,k]+zˣ[k,j]-1)
				dʸ[i,j] ≥ cʸ[i] - cʸ[j] + ℓʸ[i] + lbʸ[j]*(zʸ[i,j]+zʸ[j,i]) - Lʸ*(1-zʸ[i,j]) + lbʸ[k]*(zʸ[i,k]+zʸ[k,j]-1)
				dʸ[i,j] ≥ cʸ[i] - cʸ[j] + ℓʸ[j] + lbʸ[i]*(zʸ[i,j]+zʸ[j,i]) - Lʸ*(1-zʸ[i,j]) + lbʸ[k]*(zʸ[i,k]+zʸ[k,j]-1)
			end))
		end
	end
	cuts
end

function getcut(prob::Problem, vars::Variables, ::ThreeBoxObjective3, ii, jj, kk)
	lbˣ, lbʸ = prob.wlb, prob.hlb
	cˣ, cʸ = vars.cˣ, vars.cʸ
	zˣ, zʸ = vars.zˣ, vars.zʸ
	dˣ, dʸ = vars.dˣ, vars.dʸ
	cuts = LinearConstraint[]
	for (_i,_j) in combinations([ii,jj,kk],2) # the two outer ones
		for (i,j) in permutations([_i,_j]) # get both orderings
			k = setdiff([ii,jj,kk],[i,j])[1]
			append!(cuts, @LinearConstraints(begin
				dˣ[i,j] ≥ cˣ[i] - cˣ[j] + (lbˣ[i]+lbˣ[j])*zˣ[i,j] + lbˣ[k]*(zˣ[i,k]+zˣ[k,j]-1)
				dʸ[i,j] ≥ cʸ[i] - cʸ[j] + (lbʸ[i]+lbʸ[j])*zʸ[i,j] + lbʸ[k]*(zʸ[i,k]+zʸ[k,j]-1)
			end))
		end
	end
	cuts
end

function getcut(prob::Problem, vars::Variables, ::ThreeBoxObjective4, ii, jj, kk)
	Lˣ,  Lʸ =  prob.W,   prob.H
	lbˣ, lbʸ = prob.wlb, prob.hlb
	ℓˣ, ℓʸ = vars.ℓˣ, vars.ℓʸ
	zˣ, zʸ = vars.zˣ, vars.zʸ
	dˣ, dʸ = vars.dˣ, vars.dʸ
	cuts = LinearConstraint[]
	for (_i,_j) in combinations([ii,jj,kk],2) # the two outer ones
		for (i,j) in permutations([_i,_j]) # get both orderings
			k = setdiff([ii,jj,kk],[i,j])[1]
			append!(cuts, @LinearConstraints(begin
				2dˣ[i,j] ≥ ℓˣ[i] - Lˣ*(1-zˣ[i,j]-zˣ[j,i]) + lbˣ[j]*(zˣ[i,j]+zˣ[j,i]) + 2lbˣ[k]*(zˣ[i,k]+zˣ[k,j]-1)
				2dˣ[i,j] ≥ ℓˣ[j] - Lˣ*(1-zˣ[i,j]-zˣ[j,i]) + lbˣ[i]*(zˣ[i,j]+zˣ[j,i]) + 2lbʸ[k]*(zʸ[i,k]+zʸ[k,j]-1)
			end))
		end
	end
	cuts
end


text:
MOI.supports(::Optimizer, ::MOI.NLPBlock) = true

function walk_and_strip_variable_index!(expr::Expr)
    for i in 1:length(expr.args)
        if expr.args[i] isa MOI.VariableIndex
            expr.args[i] = expr.args[i].value
        end
        walk_and_strip_variable_index!(expr.args[i])
    end
    return
end

walk_and_strip_variable_index!(not_expr) = nothing

function MOI.set(model::Optimizer, ::MOI.NLPBlock, nlp_data::MOI.NLPBlockData)
    if model.nlp_block_data !== nothing
        error("Nonlinear block already set; cannot overwrite. Create a new model instead.")
    end
    model.nlp_block_data = nlp_data

    nlp_eval = nlp_data.evaluator

    MOI.initialize(nlp_eval, [:ExprGraph])

    if nlp_data.has_objective
        # according to test: test_nonlinear_objective_and_moi_objective_test
        # from MOI 0.10.9, linear objectives are just ignores if the noliena exists
        # if model.inner.objective_expr !== nothing
            # error("Two objectives set: One linear, one nonlinear.")
        # end
        obj = verify_support(MOI.objective_expr(nlp_eval))
        walk_and_strip_variable_index!(obj)

        model.inner.objective_expr = obj
        # if obj == :NaN
        #     model.inner.objective_expr = 0.0
        # else
        # end
    end

    for i in 1:length(nlp_data.constraint_bounds)
        expr = verify_support(MOI.constraint_expr(nlp_eval, i))
        lb = nlp_data.constraint_bounds[i].lower
        ub = nlp_data.constraint_bounds[i].upper
        @assert expr.head == :call
        if expr.args[1] == :(==)
            @assert lb == ub == expr.args[3]
        elseif expr.args[1] == :(<=)
            @assert lb == -Inf
            lb = nothing
            @assert ub == expr.args[3]
        elseif expr.args[1] == :(>=)
            @assert lb == expr.args[3]
            @assert ub == Inf
            ub = nothing
        else
            error("Unexpected expression $expr.")
        end
        expr = expr.args[2]
        walk_and_strip_variable_index!(expr)
        push!(model.inner.constraint_info, ConstraintInfo(expr, lb, ub, nothing))
    end
    return
end


text:
using Test, ColorSchemes, ColorSchemeTools, FileIO, Colors
using ImageMagick, QuartzImageIO

function run_all_tests()

    @testset "basic functions" begin
        # load existing scheme from ColorSchemes.jl
        hok = ColorSchemes.hokusai

        @test length(hok) == 32

        # image file located here in the test directory

        # create a colorscheme from image file, default is 10
        hokusai_test = ColorSchemeTools.extract(dirname(@__FILE__) * "/hokusai.jpg")

        @test length(hokusai_test) == 10

        # extract colors and weights
        c, w =  ColorSchemeTools.extract_weighted_colors(dirname(@__FILE__) * "/hokusai.jpg", 10, 10, 0.01; shrink = 4)

        @test length(c) == 10
        @test length(w) == 10

        # test that sampling schemes yield different values
        @test get(hokusai_test, 0.0) != get(hokusai_test, 0.5)

        # test sort
        @test ColorSchemeTools.sortcolorscheme(hokusai_test, rev=true) != ColorSchemeTools.sortcolorscheme(hokusai_test)

        # create weighted palette; there is some unpredictability here... :)
        csw = colorscheme_weighted(c, w, 37)
        @test 36 <= length(csw) <= 38

        # default is 50
        csw = colorscheme_weighted(c, w)
        @test length(csw) == 50

        # save as Julia text file
        colorscheme_to_text(ColorSchemes.hokusai, "hokusai_test_version", "hokusai_as_text.jl")

        # and read it in as text
        open("hokusai_as_text.jl") do f
            lines = readlines(f)
            @test occursin("loadcolorscheme(", lines[1])
            @test occursin("Colors.RGB{Float64}(0.1112499123425623", lines[4])
        end
    end

    @testset "getinverse tests" begin
        getinverse(ColorSchemes.colorschemes[:leonardo], RGB(1, 0, 0))
        getinverse(ColorScheme([Colors.RGB(0,0,0), Colors.RGB(1,1,1)]),  Colors.RGB(.5,.5,.5))
        cs = ColorScheme(range(Colors.RGB(0,0,0), stop=Colors.RGB(1,1,1), length=5))
        gi = getinverse(cs, cs[3])
        @test gi == 0.5
    end

    @testset "convert to scheme tests" begin
        # Add color to a grayscale image.
        # now using ColorScheme objects and .colors accessors
        red_cs = ColorScheme(range(RGB(0,0,0), stop=RGB(1,0,0), length=11))
        gray_cs = ColorScheme(range(RGB(0,0,0), stop=RGB(1,1,0), length=11))
        vs = [getinverse(gray_cs, p) for p in red_cs.colors]
        cs = ColorScheme([RGB(v, v, v) for v in vs])
        rcs = [get(red_cs, p) for p in vs]

        new_img = convert_to_scheme(red_cs, gray_cs.colors)
        # TODO
        # This is broken.. It should be way more specific. See next test.
        @test all(.≈(new_img, red_cs.colors, atol=0.5))

        # Should be able to uniquely match each increasing color with the next
        # increasing color in the new scale.
        red_cs = ColorScheme(range(RGB(0,0,0), stop=RGB(1,1,1)))
        blue_scale_img = range(RGB(0,0,0), stop=RGB(0,0,1))
        new_img = convert_to_scheme(red_cs, blue_scale_img)
        @test_broken unique(new_img) == new_img
    end

    @testset "make_colorscheme tests" begin
        cs = make_colorscheme([colorant"red", colorant"green", colorant"blue"], 20)
        @test cs[1] == RGB{Float64}(1.0,0.0,0.0)
        @test cs[end] == RGB{Float64}(0.0,0.0,1.0)
        colorscheme_to_text(cs, "rgb_scheme", "rgb_scheme.jl")
    end
end

function run_minimum_tests()
    @testset "basic minimum tests" begin
        # load scheme
        hok = ColorSchemes.hokusai

        @test length(hok) == 32

        # test sort
        @test ColorSchemeTools.sortcolorscheme(hok, rev=true) != ColorSchemeTools.sortcolorscheme(hok)

        # save as text
        ColorSchemeTools.colorscheme_to_text(hok, "hokusai_test_version", "hokusai_as_text.jl")

        @test filesize("hokusai_as_text.jl") > 2000

        open("hokusai_as_text.jl") do f
            lines = readlines(f)
            @test occursin("Colors.RGB{Float64}(0.1112499123425623", lines[4])
        end

        # convert an Array{T,2} to an RGB image
        tmp = get(ColorSchemes.leonardo, rand(10, 10))
        @test typeof(tmp) == Array{ColorTypes.RGB{Float64}, 2}

        # test conversion with default clamp
        x = [0.0 1.0 ; -1.0 2.0]
        y = get(ColorSchemes.leonardo, x)
        @test y[1,1] == y[2,1]
        @test y[1,2] == y[2,2]

        # test conversion with symbol clamp
        y2 = get(ColorSchemes.leonardo, x, :clamp)
        @test y2 == y

        # test conversion with symbol extrema
        y2 = get(ColorSchemes.leonardo, x, :extrema)
        @test y2[2,1] == y[1,1]   # Minimum now becomes one edge of ColorScheme
        @test y2[2,2] == y[1,2]   # Maximum now becomes other edge of ColorScheme
        @test y2[1,1] !== y2[2,1] # Inbetween values or now different

        # test conversion with manually supplied range
        y3 = get(ColorSchemes.leonardo, x, (-1.0, 2.0))
        @test y3 == y2

        # test with steplen (#17)
        r  = range(0, stop=5, length=10)
        y  = get(ColorSchemes.leonardo, r)
        y2 = get(ColorSchemes.leonardo, collect(r))
        @test y == y2

        # test for specific value
        val = 0.2
        y   = get(ColorSchemes.leonardo, [val])
        y2  = get(ColorSchemes.leonardo, val)
        @test y2 == y[1]
        end
end

if get(ENV, "ColorSchemeTools_KEEP_TEST_RESULTS", false) == "true"
    # they changed mktempdir in v1.3
    if VERSION <= v"1.2"
        cd(mktempdir())
    else
        cd(mktempdir(cleanup=false))
    end
    @info("...Keeping the results in: $(pwd())")
    @info("..running minimum tests")
    run_minimum_tests()
    @info("..running all tests")
    run_all_tests()
    @info("Test images saved in: $(pwd())")
else
    mktempdir() do tmpdir
    cd(tmpdir) do
        @info("..running tests in: $(pwd())")
        @info("..but not keeping the results")
        @info("..because you didn't do: ENV[\"ColorSchemeTools_KEEP_TEST_RESULTS\"] = \"true\"")
        @info("..running minimum tests")
        run_minimum_tests()
        @info("..running all tests")
        run_all_tests()
        @info("..Test images weren't saved. To see the test images, next time do this before running:")
        @info(" ENV[\"ColorSchemeTools_KEEP_TEST_RESULTS\"] = \"true\"")
        end
    end
end


text:
using .Dictionaries

# Monoid / Alternative
# ====================

TypeClasses.neutral(T::Type{<:AbstractDictionary}) = T()

# generic combine/⊕ for Dict: using ⊕ on the elements when needed
"""
```jldoctest
julia> using TypeClasses, Dictionaries

julia> dict1 = Dictionary([:a, :b, :c], ["1", "2", "3"])
3-element Dictionaries.Dictionary{Symbol, String}
 :a │ "1"
 :b │ "2"
 :c │ "3"

julia> dict2 = Dictionary([:b, :c, :d], ["4", "5", "6"])
3-element Dictionaries.Dictionary{Symbol, String}
 :b │ "4"
 :c │ "5"
 :d │ "6"

julia> dict1 ⊕ dict2
4-element Dictionaries.Dictionary{Symbol, String}
 :a │ "1"
 :b │ "24"
 :c │ "35"
 :d │ "6"
```
"""
function TypeClasses.combine(d1::AbstractDictionary, d2::AbstractDictionary)
    d1_combined = map(pairs(d1)) do (key, value)
        if key in keys(d2)
            value ⊕ d2[key]
        else
            value
        end
    end
    remaining_keys = setdiff(keys(d2), keys(d1))
    merge(d1_combined, getindices(d2, remaining_keys))
end

# generic orelse/⊘ for Dict
"""
    orelse(d1::Dict, d2::Dict) -> Dict

Following the orelse semantics on Option values, the first value is retained, and the second is dropped.
Hence this is the flipped version of `Base.merge`.
"""
TypeClasses.orelse(d1::AbstractDictionary, d2::AbstractDictionary) = merge(d2, d1)

# Functor / Applicative / Monad
# =============================

# adapted from the Scala Cats implementation https://github.com/typelevel/cats/blob/main/core/src/main/scala/cats/instances/map.scala

# map/foreach are already defined

# there is no TypeClasses.pure implementation because we cannot come up with an arbitrary key

"""
```jldoctest
julia> using TypeClasses, Dictionaries

julia> dict1 = Dictionary(["a", "b", "c"], [1, 2, 3])
3-element Dictionaries.Dictionary{String, Int64}
 "a" │ 1
 "b" │ 2
 "c" │ 3

julia> dict2 = Dictionary(["b", "c", "d"], [2, 3, 4])
3-element Dictionaries.Dictionary{String, Int64}
 "b" │ 2
 "c" │ 3
 "d" │ 4

mapn(+, dict1, dict2)
2-element Dictionaries.Dictionary{String, Int64}
 "b" │ 4
 "c" │ 6
```
"""
# TypeClasses.ap(f::AbstractDictionary, d::AbstractDictionary) = TypeClasses.default_ap_having_map_flatmap(f, d)

"""
```jldoctest
julia> using TypeClasses, Dictionaries

julia> dict = Dictionary(["a", "b", "c"], [1, 2, 3])
3-element Dictionaries.Dictionary{String, Int64}
 "a" │ 1
 "b" │ 2
 "c" │ 3

julia> f(x) = Dictionary(["b", "c", "d"], [10x, 20x, 30x])
f (generic function with 1 method)

julia> flatmap(f, dict)
2-element Dictionaries.Dictionary{String, Int64}
 "b" │ 20
 "c" │ 60
```
"""
function TypeClasses.flatmap(f, d::AbstractDictionary)
    returntype = Core.Compiler.return_type(f, Tuple{eltype(d)})
    result = similar(d, eltype(returntype))
    isused = fill(false, d)

    for (key, value) in pairs(d)
        subdict = f(value)
        # only add the key if key appears in both Dictionaries
        if key in keys(subdict)
            result[key] = subdict[key]
            isused[key] = true
        end
    end
    # finally we may have lost a couple of keys
    getindices(result, findall(isused))
end

# FlipTypes
# =========

"""
```jldoctest
julia> using TypeClasses, Dictionaries

julia> d1 = dictionary([:a => Option(1), :b => Option(2)])
2-element Dictionaries.Dictionary{Symbol, Identity{Int64}}
 :a │ Identity(1)
 :b │ Identity(2)

julia> flip_types(d1)
Identity({:a = 1, :b = 2})

julia> d2 = dictionary([:a => Option(1), :b => Option()])
2-element Dictionaries.Dictionary{Symbol, Option{Int64}}
 :a │ Identity(1)
 :b │ Const(nothing)

julia> flip_types(d2)
Const(nothing)
```
"""
function TypeClasses.flip_types(a::A) where A <: AbstractDictionary
    iter = pairs(a)
    first = iterate(iter)
    constructdict(key, c) = fill(c, getindices(a, Indices([key])))  # using fill to create the correct type automatically

    if first === nothing
      # only in this case we actually need `pure(eltype(A))` and `neutral(A)`
      # for non-empty sequences everything works for Types without both
      pure(eltype(A), neutral(Base.typename(A).wrapper))
    else
      (key, b), state = first
      start = map(c -> constructdict(key, c), b)  # we can only combine on ABC
      Base.foldl(Iterators.rest(iter, state); init = start) do acc, (key, b)
        mapn(acc, b) do acc′, c  # working in applicative context B
          acc′ ⊕ constructdict(key, c)  # combining on AC
        end
      end
    end
end


text:
struct CellQuadratures
    quads::Any
    celltoquad::Any
    ncells::Any
    function CellQuadratures(quads, celltoquad)
        nphase, ncells = size(celltoquad)
        @assert nphase == 2
        @assert all(celltoquad .>= 0)
        @assert all(celltoquad .<= length(quads))
        new(quads, celltoquad, ncells)
    end
end

function CellQuadratures(mesh, levelset, numuniformqp, numcutqp)
    polyorder = order(levelset)
    dim = dimension(mesh)
    interpgrad = InterpolatingPolynomial(dim, dim, polyorder)
    numcells = number_of_cells(mesh)
    tpq = tensor_product_quadrature(dim, numuniformqp)

    quads = [tpq]

    celltoquad = zeros(Int, 2, numcells)
    xL, xR = [-1.0, -1.0], [1.0, 1.0]

    for cellid = 1:numcells
        s = cell_sign(mesh, cellid)
        if s == +1
            celltoquad[1, cellid] = 1
        elseif s == -1
            celltoquad[2, cellid] = 1
        elseif s == 0
            nodeids = nodal_connectivity(background_mesh(mesh), cellid)
            load_coefficients!(levelset, cellid)
            update_interpolating_gradient!(interpgrad, interpolater(levelset))

            try
                pquad = area_quadrature(
                    interpolater(levelset),
                    interpgrad,
                    +1,
                    xL,
                    xR,
                    numcutqp,
                    numsplits = 2,
                )
                push!(quads, pquad)
                celltoquad[1, cellid] = length(quads)

                nquad = area_quadrature(
                    interpolater(levelset),
                    interpgrad,
                    -1,
                    xL,
                    xR,
                    numcutqp,
                    numsplits = 2,
                )
                push!(quads, nquad)
                celltoquad[2, cellid] = length(quads)
            catch e
                pquad = area_quadrature(
                    interpolater(levelset),
                    interpgrad,
                    +1,
                    xL,
                    xR,
                    numcutqp,
                    numsplits = 3,
                )
                push!(quads, pquad)
                celltoquad[1, cellid] = length(quads)

                nquad = area_quadrature(
                    interpolater(levelset),
                    interpgrad,
                    -1,
                    xL,
                    xR,
                    numcutqp,
                    numsplits = 3,
                )
                push!(quads, nquad)
                celltoquad[2, cellid] = length(quads)
            end
        else
            error("Expected cellsign ∈ {-1,0,+1}, got cellsign = $s")
        end
    end
    return CellQuadratures(quads, celltoquad)
end

function CellQuadratures(cutmesh, levelset, numqp)
    return CellQuadratures(cutmesh, levelset, numqp, numqp)
end

function Base.getindex(vquads::CellQuadratures, s, cellid)
    row = cell_sign_to_row(s)
    idx = vquads.celltoquad[row, cellid]
    idx > 0 ||
        error("Cell $cellid, cellsign $s, does not have a cell quadrature")
    return vquads.quads[vquads.celltoquad[row, cellid]]
end

function Base.show(io::IO, cellquads::CellQuadratures)
    ncells = cellquads.ncells
    nuniquequads = length(cellquads.quads)
    str = "CellQuadratures\n\tNum. Cells: $ncells\n\tNum. Unique Quadratures: $nuniquequads"
    print(io, str)
end

function uniform_cell_quadrature(vquads::CellQuadratures)
    return vquads.quads[1]
end

function update_cell_quadrature!(cellquads::CellQuadratures, s, cellid, quad)
    row = cell_sign_to_row(s)
    idx = cellquads.celltoquad[row, cellid]
    cellquads.quads[idx] = quad
end


text:
@static if VERSION < v"1.6.0-DEV"
    macro load_preference(name, default)
        esc(default)
    end
else
    using Preferences: Preferences, @load_preference
    set_preferences!(args...; kwargs...) = Preferences.set_preferences!(@__MODULE__, args...; kwargs...)
end

macro pref_log_destruction(ex)
    if @load_preference("LOG_DESTRUCTION", "false") == "true"
        ex = quote
            @ccall jl_safe_printf("Finalizing $handle"::Cstring)::Cvoid
            was_destroyed = $ex
            msg =
                was_destroyed ? ":\e[32m destroyed\e[m\n" :
                ":\e[33m nothing to do\e[m\n"
            @ccall jl_safe_printf(msg::Cstring)::Cvoid
        end
    end
    esc(ex)
end

macro pref_log_refcount(ex)
    if @load_preference("LOG_REFCOUNT", "false") == "true"
        ex = quote
            val = Int(handle.refcount.value)
            @ccall jl_safe_printf("Refcount: $handle: $val "::Cstring)::Cvoid
            $ex
            val_after = Int(handle.refcount.value)
            if val_after < val
                @ccall jl_safe_printf("\e[31m↓\e[m"::Cstring)::Cvoid
            else
                @ccall jl_safe_printf("\e[32m↑\e[m"::Cstring)::Cvoid
            end
            @ccall jl_safe_printf(" $val_after\n"::Cstring)::Cvoid
        end
    end
    esc(ex)
end


text:
# constants.jl  -  Physical constants

# 2016-07 - Updated accuracy from https://github.com/DataWookie/PhysicalConstants.jl
const h = 6.62606896e-34;                          # Js = kg m2 / s
const hbar = const ħ = 1.05457162825e-34;          # Js = kg m2 / s

const eV = const q = const ElectronVolt = 1.602176487e-19;    # kg m2 / s2
const MassElectron = 9.10938188e-31;                          # kg

const ε0 = const VacuumPermittivity = 8.854187817e-12;        # A2 s4 / kg m3
const VacuumPermeability = 1.25663706144e-6;                  # kg m / a2 s2
const Debye = 3.33564095198e-30;                              # A s2 / m2
const Gauss = 1e-4;                                           # kg / A s2

const Rydberg = 2.17987196968e-18;                            # J = kg m2 / s2
const Hartree = 2*Rydberg;                                    # J = kg m2 / s2
const Boltzmann = const kB =  1.3806504e-23;                  # J/K = kg m2 / K s2

# convenience for conversions
const hbar_in_eV = const ħeV = hbar/eV #eV energy units
const kB_in_eV = kB/eV # in units of eV
const Hartree_in_eV = Hartree/eV # eV
const Rydberg_in_eV = Rydberg/eV

# Dodgy unit conversions, done dirt cheap.
const Å=1E-10 # Angstrom


text:
module Analysis
    include("uncertainty.jl")
end


text:
## header to provide surface and context
using Cairo
c = CairoRGBSurface(256,256);
cr = CairoContext(c);

save(cr);
set_source_rgb(cr,0.8,0.8,0.8);    # light gray
rectangle(cr,0.0,0.0,256.0,256.0); # background
fill(cr);
restore(cr);

save(cr);
## original example, following here

image = read_from_png("data/mulberry.png");
w = image.width;
h = image.height;

translate(cr, 128.0, 128.0);
rotate(cr, 45* pi/180);
scale(cr, 256.0/w, 256.0/h);
translate(cr, -0.5*w, -0.5*h);

set_source_surface(cr, image, 0, 0);
paint(cr);
#surface_destroy (image);

## mark picture with current date
restore(cr);
move_to(cr,0.0,12.0);
set_source_rgb(cr, 0,0,0);
show_text(cr,Libc.strftime(time()));
write_to_png(c,"sample_image.png");


text:
function UpdatableFrame(arg0::Frame, arg1::Transform, arg2::JString)
    return UpdatableFrame((Frame, Transform, JString), arg0, arg1, arg2)
end

function UpdatableFrame(arg0::Frame, arg1::Transform, arg2::JString, arg3::jboolean)
    return UpdatableFrame((Frame, Transform, JString, jboolean), arg0, arg1, arg2, arg3)
end

function update_transform(obj::UpdatableFrame, arg0::Frame, arg1::Frame, arg2::Transform, arg3::AbsoluteDate)
    return jcall(obj, "updateTransform", void, (Frame, Frame, Transform, AbsoluteDate), arg0, arg1, arg2, arg3)
end


text:
@doc doc"""
    PowerManifold{TM<:Manifold, TSize<:Tuple} <: Manifold

Power manifold $M^{N_1 \times N_2 \times \dots \times N_n}$ with power geometry.
It is represented by an array-like structure with $n$ dimensions and sizes
$N_1, N_2, \dots, N_n$, along each dimension.
`TSize` statically defines the number of elements along each axis.
For example, a manifold-valued time series would be represented by a power
manifold with $n$ equal to 1 and $N_1$ equal to the number of samples.
A manifold-valued image (for example in diffusion tensor imaging) would
be represented by a two-axis power manifold ($n=2$) with $N_1$ and $N_2$
equal to width and height of the image.

While the size of the manifold is static, points on the power manifold
would not be represented by statically-sized arrays. Operations on small
power manifolds might be faster if they are represented as [`ProductManifold`](@ref).

# Constructor

    PowerManifold(M, N_1, N_2, ..., N_n)

generates the power manifold $M^{N_1 \times N_2 \times \dots \times N_n}$.
"""
struct PowerManifold{TM<:Manifold, TSize} <: Manifold
    manifold::TM
end

@generated function rep_size_to_colons(rep_size::Tuple)
    N = length(rep_size.parameters)
    return ntuple(i -> Colon(), N)
end

function PowerManifold(manifold::Manifold, size::Int...)
    return PowerManifold{typeof(manifold), Tuple{size...}}(manifold)
end

function get_iterator(M::PowerManifold{<:Manifold, Tuple{N}}) where N
    return 1:N
end

@generated function get_iterator(M::PowerManifold{<:Manifold, SizeTuple}) where SizeTuple
    return Base.product(map(Base.OneTo, tuple(SizeTuple.parameters...))...)
end

@inline function _read(rep_size::Tuple, x::AbstractArray, i::Int)
    return _read(rep_size, x, (i,))
end

@inline function _read(rep_size::Tuple, x::AbstractArray, i::Tuple)
    return view(x, rep_size_to_colons(rep_size)..., i...)
end

@inline function _read(rep_size::Tuple, x::HybridArray, i::Tuple)
    return x[rep_size_to_colons(rep_size)..., i...]
end

@inline function _write(rep_size::Tuple, x::AbstractArray, i::Int)
    return _write(rep_size, x, (i,))
end

@inline function _write(rep_size::Tuple, x::AbstractArray, i::Tuple)
    return view(x, rep_size_to_colons(rep_size)..., i...)
end

function isapprox(M::PowerManifold, x, y; kwargs...)
    result = true
    rep_size = representation_size(M.manifold)
    for i in get_iterator(M)
        result &= isapprox(M.manifold,
            _read(rep_size, x, i),
            _read(rep_size, y, i);
            kwargs...)
    end
    return result
end

function isapprox(M::PowerManifold, x, v, w; kwargs...)
    result = true
    rep_size = representation_size(M.manifold)
    for i in get_iterator(M)
        result &= isapprox(M.manifold,
            _read(rep_size, x, i),
            _read(rep_size, v, i),
            _read(rep_size, w, i);
            kwargs...)
    end
    return result
end

function representation_size(M::PowerManifold{<:Manifold, TSize}) where TSize
    return (representation_size(M.manifold)..., size_to_tuple(TSize)...)
end

function manifold_dimension(M::PowerManifold{<:Manifold, TSize}) where TSize
    return manifold_dimension(M.manifold)*prod(size_to_tuple(TSize))
end

struct PowerMetric <: Metric end

function det_local_metric(M::MetricManifold{PowerManifold, PowerMetric}, x::AbstractArray)
    result = one(eltype(x))
    rep_size = representation_size(M.manifold)
    for i in get_iterator(M)
        result *= det_local_metric(M.manifold,
            _read(rep_size, x, i))
    end
    return result
end

function inner(M::PowerManifold, x, v, w)
    result = zero(eltype(v))
    rep_size = representation_size(M.manifold)
    for i in get_iterator(M)
        result += inner(M.manifold,
            _read(rep_size, x, i),
            _read(rep_size, v, i),
            _read(rep_size, w, i))
    end
    return result
end

function norm(M::PowerManifold, x, v)
    sum_squares = zero(eltype(v))
    rep_size = representation_size(M.manifold)
    for i in get_iterator(M)
        sum_squares += norm(M.manifold,
            _read(rep_size, x, i),
            _read(rep_size, v, i))^2
    end
    return sqrt(sum_squares)
end

function distance(M::PowerManifold, x, y)
    sum_squares = zero(eltype(x))
    rep_size = representation_size(M.manifold)
    for i in get_iterator(M)
        sum_squares += distance(M.manifold,
            _read(rep_size, x, i),
            _read(rep_size, y, i))^2
    end
    return sqrt(sum_squares)
end

function exp!(M::PowerManifold, y, x, v)
    rep_size = representation_size(M.manifold)
    for i in get_iterator(M)
        exp!(M.manifold,
            _write(rep_size, y, i),
            _read(rep_size, x, i),
            _read(rep_size, v, i))
    end
    return y
end

function log!(M::PowerManifold, v, x, y)
    rep_size = representation_size(M.manifold)
    for i in get_iterator(M)
        log!(M.manifold,
            _write(rep_size, v, i),
            _read(rep_size, x, i),
            _read(rep_size, y, i))
    end
    return v
end

function injectivity_radius(M::PowerManifold, x)
    radius = 0.0
    initialized = false
    rep_size = representation_size(M.manifold)
    for i in get_iterator(M)
        cur_rad = injectivity_radius(M.manifold,
            _read(rep_size, x, i))
        if initialized
            radius = min(cur_rad, radius)
        else
            radius = cur_rad
            initialized = true
        end
    end
    return radius
end

function injectivity_radius(M::PowerManifold)
    radius = 0.0
    initialized = false
    for i in get_iterator(M)
        cur_rad = injectivity_radius(M.manifold)
        if initialized
            radius = min(cur_rad, radius)
        else
            radius = cur_rad
            initialized = true
        end
    end
    return radius
end

"""
    PowerRetraction(retraction::AbstractRetractionMethod)

Power retraction based on `retraction`. Works on [`PowerManifold`](@ref).
"""
struct PowerRetraction{TR<:AbstractRetractionMethod} <: AbstractRetractionMethod
    retraction::TR
end

function retract!(M::PowerManifold, y, x, v, method::PowerRetraction)
    rep_size = representation_size(M.manifold)
    for i in get_iterator(M)
        retract!(M.manifold,
            _write(rep_size, y, i),
            _read(rep_size, x, i),
            _read(rep_size, v, i),
            method.retraction)
    end
    return y
end

"""
    InversePowerRetraction(inverse_retractions::AbstractInverseRetractionMethod...)

Power inverse retraction of `inverse_retractions`.
Works on [`PowerManifold`](@ref).
"""
struct InversePowerRetraction{TR<:AbstractInverseRetractionMethod} <: AbstractInverseRetractionMethod
    inverse_retraction::TR
end

function inverse_retract!(M::PowerManifold, v, x, y, method::InversePowerRetraction)
    rep_size = representation_size(M.manifold)
    for i in get_iterator(M)
        inverse_retract!(M.manifold,
            _write(rep_size, v, i),
            _read(rep_size, x, i),
            _read(rep_size, y, i),
            method.inverse_retraction)
    end
    return v
end

function flat!(M::PowerManifold, v::FVector{CotangentSpaceType}, x, w::FVector{TangentSpaceType})
    rep_size = representation_size(M.manifold)
    for i in get_iterator(M)
        flat!(M.manifold,
            FVector(CotangentSpace, _write(rep_size, v.data, i)),
            _read(rep_size, x, i),
            FVector(TangentSpace, _read(rep_size, w.data, i)))
    end
    return v
end

function sharp!(M::PowerManifold, v::FVector{TangentSpaceType}, x, w::FVector{CotangentSpaceType})
    rep_size = representation_size(M.manifold)
    for i in get_iterator(M)
        sharp!(M.manifold,
            FVector(TangentSpace, _write(rep_size, v.data, i)),
            _read(rep_size, x, i),
            FVector(CotangentSpace, _read(rep_size, w.data, i)))
    end
    return v
end

"""
    check_manifold_point(M::ProductManifold, x; kwargs...)

Check whether `x` is a valid point on the [`ProductManifold`](@ref) `M`.

The tolerance for the last test can be set using the ´kwargs...`.
"""
function check_manifold_point(M::PowerManifold, x; kwargs...)
    rep_size = representation_size(M.manifold)
    for i in get_iterator(M)
        imp = check_manifold_point(M.manifold,
            _read(rep_size, x, i); kwargs...)
        imp === nothing || return imp
    end
    return nothing
end

"""
    check_tangent_vector(M::ProductManifold, x, v; kwargs... )

Check whether `v` is a tangent vector to `x` on the [`ProductManifold`](@ref)
`M`, i.e. atfer [`check_manifold_point`](@ref)`(M, x)`, and all projections to
base manifolds must be respective tangent vectors.

The tolerance for the last test can be set using the ´kwargs...`.
"""
function check_tangent_vector(M::PowerManifold, x, v; kwargs...)
    mpe = check_manifold_point(M, x)
    if mpe !== nothing
        return mpe
    end
    rep_size = representation_size(M.manifold)
    for i in get_iterator(M)
        imp = check_tangent_vector(M.manifold,
            _read(rep_size, x, i),
            _read(rep_size, v, i); kwargs...)
        imp === nothing || return imp
    end
    return nothing
end

"""
    PowerPointDistribution(M::PowerManifold, distribution)

Power distribution on manifold `M`, based on `distribution`.
"""
struct PowerPointDistribution{TM<:PowerManifold, TD<:MPointDistribution, TX} <: MPointDistribution{TM}
    manifold::TM
    distribution::TD
    x::TX
end

function support(d::PowerPointDistribution)
    return MPointSupport(d.manifold)
end

function rand(rng::AbstractRNG, d::PowerPointDistribution)
    x = similar(d.x)
    _rand!(rng, d, x)
    return x
end

function _rand!(rng::AbstractRNG, d::PowerPointDistribution, x::AbstractArray)
    rep_size = representation_size(d.manifold.manifold)
    for i in get_iterator(d.manifold)
        _rand!(rng,
            d.distribution,
            _write(rep_size, x, i))
    end
    return x
end

"""
    PowerFVectorDistribution([type::VectorBundleFibers], [x], distr)

Generates a random vector at point `x` from vector space (a fiber of a tangent
bundle) of type `type` using the power distribution of `distr`.

Vector space type and `x` can be automatically inferred from distribution `distr`.
"""
struct PowerFVectorDistribution{TSpace<:VectorBundleFibers{<:VectorSpaceType, <:PowerManifold}, TD<:FVectorDistribution, TX} <: FVectorDistribution{TSpace, TX}
    type::TSpace
    x::TX
    distribution::TD
end

function support(tvd::PowerFVectorDistribution)
    return FVectorSupport(tvd.type, tvd.x)
end

function rand(rng::AbstractRNG, d::PowerFVectorDistribution)
    fv = zero_vector(d.type, d.x)
    _rand!(rng, d, fv)
    return fv
end

function _rand!(rng::AbstractRNG, d::PowerFVectorDistribution, v::AbstractArray)
    rep_size = representation_size(d.type.M.manifold)
    for i in get_iterator(d.type.M)
        copyto!(d.distribution.x, _read(rep_size, d.x, i))
        _rand!(rng, d.distribution, _read(rep_size, v, i))
    end
    return v
end


text:
using MyDistanceSimilarityPKG

@testset "Teste 1" begin

end


text:
using ApproximateBayesianComputation

#Start Testing...
data = rand(Distributions.Normal(2.0, 0.5), 100)

summaryStatistics = mean(data)
nsteps = 10
nparticles = 200
ninitial = 1000

sample_prior, density_prior = make_model_prior(Distributions.Normal(0.0, 0.5))
sample_kernel, density_kernel = make_normal_kernel()
kernel_bandwidth = weightedstd2

forward_model(μ::Float64) = rand(Distributions.Normal(μ, 0.5), 100)
compute_distance(x::Array{Float64, 1}, y::Float64) = abs(mean(x) - y)

rank_distances = identity
shrink_threshold{A <: SingleMeasureAbcPmc}(abcpmc::A) = quantile_threshold(abcpmc, 0.3)

out = abc_pmc(summaryStatistics, nsteps, nparticles, ninitial, sample_prior,
              density_prior, sample_kernel, density_kernel, forward_model,
              compute_distance, rank_distances, kernel_bandwidth, shrink_threshold)


text:
################################################################################
#
#  Parent constructor
#
################################################################################

function SMatSpace(R::Ring, r::Int, c::Int; cached = true)
  T = elem_type(R)
  return SMatSpace{T}(R, r, c, cached)
end

################################################################################
#
#  Parent and base ring
#
################################################################################

base_ring(A::SMatSpace{T}) where {T} = A.base_ring::parent_type(T)

parent(A::SMat) = SMatSpace(base_ring(A), A.r, A.c)

base_ring(A::SMat{T}) where {T} = A.base_ring::parent_type(T)

function rows(A::SMat)
  @assert A.r == length(A.rows)
  return A.r
end

function cols(A::SMat)
  return A.c
end

size(A::SMat) = (rows(A), cols(A))

################################################################################
#
#  Hashing
#
################################################################################

function hash(A::SMat, h::UInt)
  return hash(A.r, hash(A.c, hash(A.rows, h)))
end

################################################################################
#
#  Comparison
#
################################################################################

function ==(x::SMat{T}, y::SMat{T}) where T
  parent(x) != parent(y) && error("Parents incompatible")
  return x.rows == y.rows
end

################################################################################
#
#  Sparsity
#
################################################################################

# this makes only sense for SMat{fmpz}
function sparsity(A::SMat{T}) where T
  return A.nnz/(A.r * A.c), nbits(maximum(abs, A))
end

################################################################################
#
#  String I/O
#
################################################################################

function show(io::IO, A::SMat{T}) where T
  print(io, "Sparse ", A.r, " x ", A.c, " matrix with ")
  print(io, A.nnz, " non-zero entries\n")
end

################################################################################
#
#  Empty sparse matrix (given base ring)
#
################################################################################

function SMat(R::T) where T <: Ring
  r = SMat{elem_type(R)}()
  r.base_ring = R
  return r
end

################################################################################
#
#  Copy
#
################################################################################

function copy(A::SMat{T}) where T
  return sub(A, 1:rows(A), 1:cols(A))
end

################################################################################
#
#  Entry/row access
#
################################################################################

function getindex(A::SMat{T}, i::Int, j::Int) where T
  if i in 1:A.r
    ra = A.rows[i]
    p = findfirst(x->x==j, ra.pos)
    if p != 0
      return ra.values[p]
    end
  end
  return zero(base_ring(A))
end

function getindex(A::SMat{T}, i::Int) where T
  if i in 1:A.r
    return A.rows[i]
  end
  return SRow{T}()
end

function setindex!(A::SMat{T}, b::SRow{T}, i::Int) where T
  A.rows[i] = b
end

################################################################################
#
#  Random row from sparse matrix
#
################################################################################

function randrow(A::SMat{T}) where T
  return rand(A.rows)
end

################################################################################
#
#  Sparse matrix from dense matrix
#
################################################################################

# The entries of the sparse matrix will be coerced into the ring R.
# It defaults to the base ring of the input matrix.
# If keepzrows is set, then zero rows will not be remove
function SMat(A::T; R::S = base_ring(A),
                     keepzrows::Bool = true) where {T <: MatElem, S <: Ring}

  m = SMat(R)
  m.c = cols(A)
  m.r = 0

  for i=1:rows(A)
    if iszero_row(A, i)
      if !keepzrows
        continue
      else
        r = SRow{elem_type(R)}()
        #push!(m.rows, SRow{elem_type(R)}())
      end
    else
      r = SRow{elem_type(R)}()
      for j = 1:cols(A)
        t = A[i, j]
        if t != 0
          m.nnz += 1
          push!(r.values, R(t))
          push!(r.pos, j)
        end
      end
    end
    push!(m.rows, r)
    m.r += 1
  end
  return m
end

doc"""
    SMat{T}(A::Array{T, 2}) -> SMat{T}

> Constructs the SMat (Hecke-sparse matrix) with coefficients of
> type T corresponding to A.
"""
function SMat(A::Array{T, 2}) where T <: RingElem
  length(A) == 0 && error("Cannot create sparse matrix from empty array")
  m = SMat(parent(A[1, 1]))
  m.c = Base.size(A, 2)
  m.r = 0
  for i=1:size(A, 1)
    if iszero_row(A, i)
      continue
    end
    r = SRow{T}()
    for j =1:size(A, 2)
      if A[i,j] != 0
        m.nnz += 1
        push!(r.values, T(A[i,j]))
        push!(r.pos, j)
      end
    end
    push!(m.rows, r)
    m.r += 1
  end
  return m
end

function (M::SMatSpace)(A::T; R::S = base_ring(A),
                                 keepzrows::Bool = true) where {T <: MatElem, S <: Ring}
  return SMat(A, R, keepzrows)
end

function (M::SMatSpace)(A::Array{T, 2}) where T <: MatElem
  return SMat(A)
end

# a faster version for nmod_mat -> SMat{T}
# it avoids the creation of elements in ResidueRing(FlintZZ, n)
doc"""
    SMat{S <: Ring}(A::nmod_mat; R::S = base_ring(A), keepzrows::Bool = false)

> "Lifts" the entries in $A$ to a sparse matrix over $R$.
"""
function SMat(A::nmod_mat; R::S = base_ring(A), keepzrows::Bool = false) where S <: Ring
  if false && R == base_ring(A)
    return _SMat(A, R = R)
  end

  m = SMat{elem_type(R)}()
  m.c = cols(A)
  m.r = 0

  for i=1:rows(A)
    if iszero_row(A, i)
      if !keepzrows
        continue
      else
        #push!(m.rows, SRow{elem_type(R)}())
        r = SRow(R)
      end
    else
      r = SRow(R)
      for j =1:cols(A)
        t = ccall((:nmod_mat_get_entry, :libflint), Base.GMP.Limb, (Ptr{nmod_mat}, Int, Int), &A, i - 1, j - 1)
        if t != 0
          m.nnz += 1
          push!(r.values, R(t))
          push!(r.pos, j)
        end
      end
    end
    push!(m.rows, r)
    m.r += 1
  end
  return m
end

################################################################################
#
#  Modular reduction
#
################################################################################

doc"""
***
    mod_sym!(A::SMat{fmpz}, n::fmpz)

> Inplace reduction of all entries of $A$ modulo $n$ to the symmetric residue
> system.
"""
function mod_sym!(A::SMat{fmpz}, b::fmpz)
  for i=A
    mod_sym!(i, b)
  end
end

################################################################################
#
#  Ring change
#
################################################################################

doc"""
    SMat(A::SMat{fmpz}, n::Int) -> SMat{nmod}

> Converts $A$ to be a sparse matrix (row) over $Z/nZ$
"""
function SMat(A::SMat{fmpz}, n::Int)
  R = ResidueRing(FlintZZ, n, cached=false)
  return SMat(A, R)
end

doc"""
    SMat(A::SMat{fmpz}, R::Ring) -> SMat{elem_type(R)}
    SRow(A::SMat{fmpz}, R::Ring) -> SRow{elem_type(R)}

> Convert the matrix (row) $A$ to be over $R$.
"""
function SMat(A::SMat{fmpz}, R::T) where T <: Ring
  z = SMat(R)
  z.r = A.r
  z.c = A.c
  for r in A
    rz = SRow(r, R)
    if length(rz.pos) != 0
      push!(z.rows, rz)
      z.nnz += length(rz.pos)
    end
  end
  return z
end

function SMat(A::SMat{fmpz}, n::fmpz)
  R = ResidueRing(FlintZZ, n, cached=false)
  return SMat(A, R)
end

################################################################################
#
#  Transpose
#
################################################################################

doc"""
***
    transpose(A::SMat) -> SMat

> Returns the transpose of $A$.
"""
function transpose(A::SMat)
  R = base_ring(A)
  B = SMat(base_ring(A))
  n = rows(A)
  m = cols(A)
  B.rows = Array{SRow{elem_type(R)}}(m)
  for i=1:m
    B.rows[i] = SRow(R)
  end
  for i = 1:length(A.rows)
    for j = 1:length(A.rows[i].pos)
      r = A.rows[i].pos[j]
      push!(B.rows[r].pos, i)
      push!(B.rows[r].values, A.rows[i].values[j])
    end
  end
  B.c = n
  B.r = m
  B.nnz = A.nnz

  return B
end

################################################################################
#
#  Make sparse matrices iteratable
#
################################################################################

function endof(A::SMat)
  return length(A.rows)
end

function start(A::SMat)
  return 1
end

function next(A::SMat, st::Int)
  return A.rows[st], st + 1
end

function done(A::SMat, st::Int)
  return st > rows(A)
end

function length(A::SMat)
  return rows(A)
end

################################################################################
#
#  Multiplication
#
################################################################################

# (dense Array{T, 1}) * SMat{T} as (dense Array{T, 1})
# inplace
function mul!(c::Array{T, 1}, A::SMat{T}, b::Array{T, 1}) where T
  assert(length(b) == cols(A))
  assert(length(c) == rows(A))
  for i = 1:length(A.rows)
    s = 0
    I = A.rows[i]
    for j=1:length(I.pos)
      s += I.values[j]*b[I.pos[j]]
    end
    c[i] = s
  end
  return c
end

# (dense Array{T, 1}) * SMat{T} as (dense Array{T, 1})
function mul(A::SMat{T}, b::Array{T, 1}) where T
  assert(length(b) == cols(A))
  c = Array{T}(rows(A))
  mul!(c, A, b)
  return c
end

# - (dense Array{S, 1}) * SMat{T} as (dense Array{S, 1}) modulo mod::S
# - Inplace
# - Reduction as the last step, no intermediate reductions.
function mul_mod!(c::Array{S, 1}, A::SMat{T}, b::Array{S, 1}, mod::S) where {S, T}
  assert( length(b) == cols(A))
  assert( length(c) == rows(A))
  for i = 1:length(A.rows)
    s = 0
    I = A.rows[i]
    for j=1:length(I.pos)
      s += S(I.values[j]) * b[I.pos[j]]
    end
    c[i] = s % mod
  end
  return c
end

# - (dense Array{S, 1}) * SMat{T} as (dense Array{S, 1}) modulo mod::S
# - Inplace
# - Intermediate reductions.
function mul_mod_big!(c::Array{S, 1}, A::SMat{T}, b::Array{S, 1}, mod::S) where {S, T}
  assert(length(b) == cols(A))
  assert(length(c) == rows(A))
  for i = 1:length(A.rows)
    s = 0
    I = A.rows[i]
    for j=1:length(I.pos)
      s += S(I.values[j] *b[I.pos[j]] % mod)
    end
    c[i] = s % mod
  end
  return c
end

# - SMat{T} * Array{T, 2} as Array{T, 2}
# - Inplace
function mul!(c::Array{T, 2}, A::SMat{T}, b::Array{T, 2}) where T
  sz = size(b)
  assert(sz[1] == cols(A))
  tz = size(c)
  assert(tz[1] == rows(A))
  assert(tz[2] == sz[2])
  for m = 1:sz[2]
    for i = 1:length(A.rows)
      s = 0
      for j=1:length(A.rows[i].pos)
        s += A.rows[i].values[j] * b[A.rows[i].pos[j], m]
      end
      c[i, m] = s
    end
  end
  return c
end

# - SMat{T} * Array{T, 2} as Array{T, 2}
function mul(A::SMat{T}, b::Array{T, 2}) where T
  sz = size(b)
  assert(sz[1] == cols(A))
  c = Array{T}(sz[1], sz[2])
  return mul!(c, A, b)
end

# - SMat{T} * fmpz_mat as fmpz_mat
# - Inplace
function mul!(c::fmpz_mat, A::SMat{T}, b::fmpz_mat) where T
  assert(rows(b) == cols(A))
  assert(rows(c) == rows(A))
  assert(cols(c) == cols(b))
  for m = 1:cols(b)
    for i = 1:length(A.rows)
      s = 0
      for j=1:length(A.rows[i].pos)
        s += A.rows[i].values[j] * b[A.rows[i].pos[j], m]
      end
      c[i, m] = s
    end
  end
  return c
end

# - SMat{T} * fmpz_mat as fmpz_mat
function mul(A::SMat{T}, b::fmpz_mat) where T
  assert(rows(b) == cols(A))
  c = zero_matrix(FlintZZ, rows(A), cols(b))
  return mul!(c, A, b)
end

# - SMat{T} * SMat{T} as MatElem{T}
function mul(A::SMat{T}, B::SMat{T}) where T
  @assert A.c == B.r
  C = zero_matrix(base_ring(A), A.r, B.c)
  for i=1:A.r
    for j=1:B.c
      C[i,j] = mul(A[i], B[j])
    end
  end
  return C
end

# - SMat{nmod} * SMat{nmod} as nmod_mat
function mul(A::SMat{nmod}, B::SMat{nmod})
  @assert A.c == B.r
  C = zero_matrix(base_ring(A), A.r, B.c)
  for i=1:A.r
    for j=1:B.c
      C[i,j] = mul(A[i], B[j])
    end
  end
  return C
end

# - SRow{T} * SMat{T} as SRow{T}
function mul(A::SRow{T}, B::SMat{T}) where T
  C = SRow{T}()
  for (p, v) = A
    C = add_scaled_row(B[p], C, v)
  end
  return C
end

################################################################################
#
#  Addition
#
################################################################################

function +(A::SMat{T}, B::SMat{T}) where T
  C = SMat(base_ring(A))
  m = min(rows(A), rows(B))
  for i=1:m
    push!(C, A[i]+B[i])
  end
  for i=m+1:rows(A)
    push!(C, A[i])
  end
  for i=m+1:rows(B)
    push!(C, B[i])
  end
  return C
end

function -(A::SMat{T}, B::SMat{T}) where T
  C = SMat{T}()
  m = min(rows(A), rows(B))
  for i=1:m
    push!(C, A[i]-B[i])
  end
  for i=m+1:rows(A)
    push!(C, A[i])
  end
  if m < rows(B)
    n = base_ring(B)(-1)
    for i=m+1:rows(B)
      push!(C, n*B[i])
    end
  end
  return C
end

#TODO: this is a problem in anther ops as well: a length zero row does not have a base_ring
function -(A::SRow{T}, B::SRow{T}) where T
  if length(A) == 0
    if length(B) == 0
      return A
    else
      return add_scaled_row(B, A, base_ring(B)(-1))
    end
  end
  return add_scaled_row(B, A, base_ring(A)(-1))
end

################################################################################
#
#  Scalar multiplication
#
################################################################################

# The following runs into an infinite recursion in case T = fmpz
#function *{T}(b::fmpz, A::SMat{T})
#  return base_ring(A)(b)*A
#end

function *(b::T, A::SMat{T}) where T
  B = SMat(base_ring(A))
  if iszero(b)
    return B
  end
  for a = A
    push!(B, b*a)
  end
  return B
end

function *(b::Integer, A::SMat{T}) where T
  return base_ring(A)(b)*A
end

################################################################################
#
#  Submatrix
#
################################################################################

function sub(A::SMat{T}, r::UnitRange, c::UnitRange) where T
  B = SMat(base_ring(A))
  B.nnz = 0
  for i=r
    rw = SRow{T}()
    ra = A.rows[i]
    for j=1:length(ra.values)
      if ra.pos[j] in c
        push!(rw.values, ra.values[j])
        push!(rw.pos, ra.pos[j]-c.start+1)
      end
    end
    if true || length(rw.pos)>0
      push!(B.rows, rw)
      B.nnz += length(rw.pos)
    end
  end
  B.r = length(r)
  B.c = length(c)
  return B
end

################################################################################
#
#  Valence
#
################################################################################

doc"""
  valence_mc{T}(A::SMat{T}; extra_prime = 2, trans = Array{SMatSLP_add_row{T}, 1}()) -> T

  Uses a Monte-Carlo alorithm to  compute the valence of A. The valence is the valence of the minimal polynomial f of A'*A, thus the last non-zero coefficient,
  typically f(0).

  The valence is computed modulo various primes until the computation
  stabilises for extra_prime many.

  trans, if given, is  a SLP (straight-line-program) in GL(n, Z). Then
  the valence of trans * A  is computed instead.
"""
function valence_mc(A::SMat{T}; extra_prime = 2, trans = Array{SMatSLP_add_row{T}, 1}()) where T
  # we work in At * A (or A * At) where we choose the smaller of the 2
  # matrices
  if false && cols(A) > rows(A)
    At = A
    A = transpose(A)
  else
    At = transpose(A)
  end
  if maximum(abs, A) > 2^20
    mm = mul_mod_big!
    println("mul big case")
  else
    mm = mul_mod!
    println("mul small case")
  end
  c1 = Array{Int}(cols(A))
  c2 = Array{Int}(rows(A))

  for i=1:cols(A)
    c1[i] = Int(rand(-10:10))
  end

  c = copy(c1) ## need the starting value for the other primes as well

  p = next_prime(2^30)
  k = FiniteField(p)
  d = 10
  v = Array{typeof(k(1)), 1}()
  push!(v, k(c1[1]))
  while true
    while length(v) <= d
      mm(c2, A, c1, p)
      if length(trans) >0
        apply_t!(c2, trans)
        apply!(c2, trans)
      end
      mm(c1, At, c2, p)
      push!(v, k(c1[1]))
    end
    f = berlekamp_massey(v)
    fl = true
    if !fl || degree(f) > div(d, 2)-1
      d += 10
      continue
    end
    df = degree(f)
    println("Poly degree is $df, dims $(rows(A)) x $(cols(A))")

    V = fmpz(leading_coefficient(f))
    pp = fmpz(p)

    v = Array{typeof(k(1))}(2*degree(f)+1)
    while true
      p = next_prime(p)
      println(p)
      k = FiniteField(p)
      copy!(c1, c)
      v[1] = k(c1[1])
      @time for i=1:2*degree(f)
        mm(c2, A, c1, p)
        if length(trans) >0
          apply_t!(c2, trans)
          apply!(c2, trans)
        end
        mm(c1, At, c2, p)
        v[i+1] = k(c1[1])
      end
      f = berlekamp_massey(v)
      fl = true
      if !fl || degree(f) > df
        println("initial guess was wrong...")
        d += 10
        return -1
      end

      Vn = crt(V, pp, fmpz(leading_coefficient(f)), fmpz(p))
      pp *= p
      if 2*Vn > pp
        Vn = Vn - pp
      end
      if Vn == V
        extra_prime -= 1
        if extra_prime == 0
          return V
        end
      end
      V = Vn
    end
  end
end

function valence_mc(A::SMat{T}, p::Int) where T
  # we work in At * A (or A * At) where we choose the smaller of the 2
  # matrices
  if false && cols(A) > rows(A)
    At = A
    A = transpose(A)
  else
    At = transpose(A)
  end
  if maximum(abs, A) > 2^20
    mm = mul_mod_big!
    println("mul big case")
  else
    mm = mul_mod!
    println("mul small case")
  end
  c1 = Array{Int}(cols(A))
  c2 = Array{Int}(rows(A))

  for i=1:cols(A)
    c1[i] = Int(rand(-10:10))
  end

  c = copy(c1) ## need the starting value for the other primes as well

  k = FiniteField(p)
  d = 10
  v = Array{typeof(k(1)), 1}()
  push!(v, k(c1[1]))
  while true
    while length(v) <= d
      mm(c2, A, c1, p)
      mm(c1, At, c2, p)
      push!(v, k(c1[1]))
    end
    f = berlekamp_massey(v)
    fl = true
    if !fl || degree(f) > div(d, 2)-1
      d += 10
      continue
    end
    df = degree(f)
    println("Poly degree is $df, dims $(rows(A)) x $(cols(A))")
    return f
  end
end

################################################################################
#
#  Vertical concatentation
#
################################################################################

doc"""
***
    vcat!{T}(A::SMat{T}, B::SMat{T})

> Vertically joins $A$ and $B$ inplace, that is, the rows of $B$ are
> appended to $A$.
"""
function vcat!(A::SMat{T}, B::SMat{T}) where T
  @assert length(A.rows) == A.r
  @assert length(B.rows) == B.r
  A.r += B.r
  A.c = max(A.c, B.c)
  A.nnz += B.nnz
  A.rows = vcat(A.rows, B.rows)
  @assert length(A.rows) == A.r
end

doc"""
***
    vcat{T}(A::SMat{T}, B::SMat{T})

> Vertically joins $A$ and $B$.
"""
function vcat(A::SMat{T}, B::SMat{T}) where T
  @assert length(A.rows) == A.r
  @assert length(B.rows) == B.r
  C = copy(A)
  vcat!(C, B)
  return C
end

################################################################################
#
#  Horizontal concatentation
#
################################################################################

doc"""
***
    hcat!{T}(A::SMat{T}, B::SMat{T})

> Horizontally concatenates $A$ and $B$, inplace, changing $A$.
"""
function hcat!(A::SMat{T}, B::SMat{T}) where T
  o = A.c
  A.c += B.c
  nnz = A.nnz
  for i=1:min(rows(A), rows(B))
    for (p,v) in B[i]
      push!(A[i].pos, p+o)
      push!(A[i].values, v)
    end
  end
  for i=min(rows(A), rows(B))+1:rows(B)
    sr = SRow{T}()
    for (p,v) in B[i]
      push!(sr.pos, p+o)
      push!(sr.values, v)
    end
    push!(A, sr)
  end
  A.nnz = nnz + B.nnz #A.nnz may have changed - if B is longer
end

doc"""
***
    hcat{T}(A::SMat{T}, B::SMat{T})

> Horizontally concatenates $A$ and $B$.
"""
function hcat(A::SMat{T}, B::SMat{T}) where T
  C = copy(A)
  hcat!(C, B)
  return C
end

################################################################################
#
#  Append sparse row to sparse matrix
#
################################################################################

function push!(A::SMat{T}, B::SRow{T}) where T
  push!(A.rows, B)
  A.r += 1
  @assert length(A.rows) == A.r
  A.nnz += length(B.pos)
  if length(B.pos) > 0
    A.c = max(A.c, B.pos[end])
  end
end

################################################################################
#
#  Conversion to fmpz_mat
#
################################################################################

doc"""
***
    fmpz_mat{T <: Integer}(A::SMat{T})

> The same matix $A$, but as an fmpz_mat.
> Requires a conversion from the base ring of $A$ to $\mathbf ZZ$.
"""
function fmpz_mat(A::SMat{T}) where T <: Integer
  B = zero_matrix(FlintZZ, A.r, A.c)
  for i = 1:length(A.rows)
    ra = A.rows[i]
    for j = 1:length(ra.pos)
      B[i, ra.pos[j]] = ra.values[j]
    end
  end
  return B
end

doc"""
***
    fmpz_mat(A::SMat{fmpz})

> The same matix $A$, but as an fmpz_mat.
"""
function fmpz_mat(A::SMat{fmpz})
  B = zero_matrix(FlintZZ, A.r, A.c)
  for i = 1:length(A.rows)
    ra = A.rows[i]
    for j = 1:length(ra.pos)
      B[i, ra.pos[j]] = ra.values[j]
    end
  end
  return B
end

################################################################################
#
#  Maximum, minimum and norms
#
################################################################################

doc"""
***
    norm2(A::SRow{fmpz})

> The square of the euclidean norm of $A$.
"""
function norm2(A::SRow{fmpz})
  return sum([x*x for x= A.values])
end

doc"""
***
    hadamard_bound2(A::SMat{fmpz})

> The square of the product of the norms of the rows of $A$.
"""
function hadamard_bound2(A::SMat{fmpz})
  return prod([norm2(x) for x=A])
end

doc"""
***
    maximum(abs, A::SMat{fmpz}) -> fmpz

  Finds the largest, in absolute value, entry of $A$.
"""
function maximum(::typeof(abs), A::SMat{fmpz})
  m = abs(A.rows[1].values[1])
  for i in A.rows
    for j in i.values
      if cmpabs(m, j) < 0
        m = j
      end
    end
  end
  return abs(m)
end

doc"""
***
    maximum(A::SMat{fmpz}) -> fmpz

> Finds the largest entry of $A$.
"""
function maximum(A::SMat{fmpz})
  m = A.rows[1].values[1]
  for i in A.rows
    for j in i.values
      if cmp(m, j) < 0
        m = j
      end
    end
  end
  return m
end

doc"""
***
    minimum(A::SMat{fmpz}) -> fmpz

> Finds the smallest entry of $A$.
"""
function minimum(A::SMat{fmpz})
  m = fmpz(0)
  for i in A.rows
    for j in i.values
      if cmp(m, j) > 0
        m = j
      end
    end
  end
  return m
end

################################################################################
#
#  Test upper triangular shape
#
################################################################################

doc"""
    isupper_triangular(A::SMat)

> Returns true iff $A$ is upper triangular.
"""
function isupper_triangular(A::SMat)
  for i=2:A.r
    if A[i-1].pos[1] >= A[i].pos[1]
      return false
    end
  end
  return true
end

################################################################################
#
#  Sparse to dense
#
################################################################################

function Nemo.matrix(A::SMat)
  R = base_ring(A)
  M = zero_matrix(R, rows(A), cols(A))
  for i=1:rows(A)
    for j=1:length(A.rows[i].pos)
      M[i, A.rows[i].pos[j]] = A.rows[i].values[j]
    end
  end
  return M
end

################################################################################
#
#  Identity matrices
#
################################################################################

doc"""
    id{T}(::Type{SMat{T}}, n::Int) -> SMat{T}

> The $n\times n$ identity matrix as a SMat of type T.
"""
function id{T}(::Type{SMat{T}}, n::Int)
  A = SMat{T}()
  for i=1:n
    push!(A, SRow{T}([(i, T(1))]))
  end
  return A
end

doc"""
    id{S}(::Type{SMat}, R::S, n::Int) -> SMat{elem_type(R)}

> The $n \times n$ identity over $R$ as a SMat.
> Necessary if $T(1)$ for the type $T$ does not work.
"""
function id(::Type{SMat}, R::Ring, n::Int)
  T = elem_type(R)
  A = SMat(R)
  for i=1:n
    push!(A, SRow{T}([(i, R(1))]))
  end
  return A
end

doc"""
   identity_matrix(::Type{SMat}, R::Ring, n::Int)
   identity_matrix(::Type{MatElem}, R::Ring, n::Int)
> Create a sparse (resp. dense) $n$ times $n$ identity matrix over $R$.
"""
function identity_matrix(::Type{SMat}, R::Ring, n::Int)
  return id(SMat, R, n)
end

function identity_matrix(::Type{MatElem}, R::Ring, n::Int)
  return identity_matrix(R, n)
end

doc"""
   zero_matrix(::Type{SMat}, R::Ring, n::Int)
   zero_matrix(::Type{SMat}, R::Ring, n::Int, m::Int)
   zero_matrix(::Type{MatElem}, R::Ring, n::Int)
   zero_matrix(::Type{MatElem}, R::Ring, n::Int, m::Int)
> Create a sparse (resp. dense) $n$ times $n$ (resp. $n$ times $m$) zero matrix over $R$.
"""
function zero_matrix(::Type{SMat}, R::Ring, n::Int)
  S = SMat(R)
  S.rows = [SRow(R)() for i=1:n]
  S.c = S.r = n
  return S
end

function zero_matrix(::Type{SMat}, R::Ring, n::Int, m::Int)
  S = SMat(R)
  S.rows = [SRow(R)() for i=1:n]
  S.r = n
  S.c = m
  return S
end

function zero_matrix(::Type{MatElem}, R::Ring, n::Int)
  return zero_matrix(R, n)
end

function zero_matrix(::Type{MatElem}, R::Ring, n::Int, m::Int)
  return zero_matrix(R, n, m)
end

function Base.cat(n::Int, A::SMat...)
  if n==1
    return vcat(A...)
  elseif n==2
    return hcat(A...)
  else
    error("dims must be 1 or 2")
  end
end

function Base.vcat(A::SMat...)
  B = copy(A[1])
  for i=2:length(A)
    for r = A[i].rows
      push!(B, copy(r))
    end
  end
  return A
end

#base case
function Base.hcat(A::SMat)
  return A
end

function Base.hcat(A::SMat, B::SMat, C::SMat...)
  return hcat(hcat(A, B), C...)
end

function Base.cat(dims::Tuple{Int, Int}, A::SMat...)
  if dims != (1,2)
    error("dims must be (1, 2)")
  end
  B = copy(A[1])
  c = B.c
  for i=2:length(A)
    for j=1:rows(A[i])
      R = SRow(base_ring(B))
      for (k,v) = A[i].rows[j]
        push!(R.pos, k+c)
        push!(R.values, v)
      end
      push!(B, R)
    end
    c += A[i].c
  end
  return B
end

doc"""
    isid{T}(A::SMat{T})

> Tests if $A$ is the $n \times n$ identity.
"""
function isid(A::SMat{T}) where T
  if A.c != A.r
    return false
  end
  for i = 1:A.r
    if length(A.rows[i].pos) != 1
      return false
    end
    if A.rows[i].pos[1] != i ||
       !isone(A.rows[i].values[1])
      return false
    end
  end
  return true
end

################################################################################
#
#  Serialization
#
################################################################################

doc"""
  toNemo(io::IOStream, A::SMat; name = "A")

  Prints the SMat as a julia-program into the file corresponding to io.
  The file can be included to get the matrix.
  `name` controls the variable name of the matrix.
"""
function toNemo(io::IOStream, A::SMat; name = "A")
  T = typeof(A.rows[1].values[1])
  println(io, name, " = SMat{$T}()")
  for i=A.rows
    print(io, "push!($(name), SRow{$T}(Tuple{Int, $T}[");
    for j=1:length(i.pos)-1
      print(io, "($(i.pos[j]), $(i.values[j])), ")
    end
    println(io, "($(i.pos[end]), $(i.values[end]))]))")
  end
end

doc"""
  toNemo(io::String, A::SMat; name = "A")

  Prints the SMat as a julia-program into the file named io.
  The file can be included to get the matrix.
  `name` controls the variable name of the matrix.
"""
function toNemo(f::String, A::SMat; name = "A")
  io = open(f, "w")
  toNemo(io, A, name=name)
  close(io)
end

################################################################################
#
#  Conversion to julia types
#
################################################################################

doc"""
***
    sparse{T}(A::SMat{T}) -> sparse{T}

> The same matrix, but as a sparse matrix of julia type.
"""
function sparse(A::SMat{T}) where T
  I = Array{Int}(A.nnz)
  J = Array{Int}(A.nnz)
  V = Array{T}(A.nnz)
  i = 1
  for r = 1:rows(A)
    for j=1:length(A.rows[r].pos)
      I[i] = r
      J[i] = A.rows[r].pos[j]
      V[i] = A.rows[r].values[j]
      i += 1
    end
  end
  return sparse(I, J, V)
end

doc"""
***
    Array{T}(A::SMat{T}) -> Array{T, 2}

> The same matrix, but as a two-dimensional julia array.
"""
function Array(A::SMat{T}) where T
  R = zero_matrix(base_ring(A), A.r, A.c)
           # otherwise, most entries will be #undef
           # at least if T is a flint-type
  for i=1:rows(A)
    for j=1:length(A.rows[i].pos)
      R[i,A.rows[i].pos[j]] = A.rows[i].values[j]
    end
  end
  return R
end


text:
# This file is part of Kpax3. License is MIT.

include("basic_functions.jl")
include("partition_functions.jl")
include("io.jl")


text:
# Test the different anderson accelerator types

using COSMO, SparseArrays, LinearAlgebra, Test

anderson_types = [
    AndersonAccelerator{Float64, Type1, RollingMemory, NoRegularizer};
    AndersonAccelerator{Float64, Type1, RestartedMemory, NoRegularizer};
    AndersonAccelerator{Float64, Type1, RollingMemory, TikonovRegularizer};
    AndersonAccelerator{Float64, Type1, RestartedMemory, TikonovRegularizer};
    AndersonAccelerator{Float64, Type1, RollingMemory, FrobeniusNormRegularizer};
    AndersonAccelerator{Float64, Type1, RestartedMemory, FrobeniusNormRegularizer};
    AndersonAccelerator{Float64, Type2{NormalEquations}, RollingMemory, NoRegularizer};
    AndersonAccelerator{Float64, Type2{NormalEquations}, RestartedMemory, NoRegularizer};
    AndersonAccelerator{Float64, Type2{NormalEquations}, RollingMemory, TikonovRegularizer};
    AndersonAccelerator{Float64, Type2{NormalEquations}, RestartedMemory, TikonovRegularizer};
    AndersonAccelerator{Float64, Type2{NormalEquations}, RollingMemory, FrobeniusNormRegularizer};
    AndersonAccelerator{Float64, Type2{NormalEquations}, RestartedMemory, FrobeniusNormRegularizer};
    AndersonAccelerator{Float64, Type2{QRDecomp}, RestartedMemory, NoRegularizer}
]

@testset "AndersonAccelerator" begin

    @testset "Types" begin
        q = [1; 1.];
        P = sparse([4. 1; 1 2]);
        A = [1. 1; 1 0; 0 1];
        l = [1.; 0; 0];
        u = [1; 0.7; 0.7];

        Aa = [-A; A]
        ba = [u; -l]
        constraint1 = COSMO.Constraint(Aa, ba, COSMO.Nonnegatives);

        for at in anderson_types
            settings = COSMO.Settings(accelerator = at);
            model = COSMO.Model();
            assemble!(model, P, q, constraint1, settings = settings);
            res = COSMO.optimize!(model);
            @test res.status == :Solved
        end
    end

    @testset "Types" begin
        q = [1; 1.];
        P = sparse([4. 1; 1 2]);
        A = [1. 1; 1 0; 0 1];
        l = [1.; 0; 0];
        u = [1; 0.7; 0.7];

        Aa = [-A; A]
        ba = [u; -l]
        constraint1 = COSMO.Constraint(Aa, ba, COSMO.Nonnegatives);

        for at in anderson_types
            settings = COSMO.Settings(accelerator = at);
            model = COSMO.Model();
            assemble!(model, P, q, constraint1, settings = settings);
            res = COSMO.optimize!(model);
            @test res.status == :Solved
        end
    end

end