text:
##
using Pkg
pkg"activate ."
pkg"add JuMP"
pkg"add Interact"
ENV["GUROBI_HOME"] = "/home/franklin/Documentos/gurobi/gurobi903/linux64"
pkg"add Gurobi"
pkg"add Plots"
pkg"add Random"
pkg"add DataFrames"
pkg"add Latexify"
##
using Plots; gr(size = (600, 400))
## função que calcula a distância Euclidiana
dist_euclid(x1, x2, y1, y2) = ((x1-y1)^2+(x2-y2)^2)^0.5

##
distMin = 10
npts = [2,2]

##
# Essa função gera pontos aleatórios com uma pequena perturbação em uma malha
# de [m n] pontos. O parâmetro distMin é a distância mínima entre os pontos
# horizontalmente e verticalmente.
using Random
function gerapontos(distMin, npts)
    Random.seed!(1)
    Coord = [0 0]
    for i = 1 : npts[1], j = 1:npts[2]
        if i!=1 || j!=1
            Coord = [Coord; [distMin*(i-1)+2*rand() distMin*(j-1)+2*rand()]]
        end
    end
    np = length(Coord[:,1])
    return Coord, np
end

##
# Essa função calcula as distâncias entre os pontos, atribui adjacências entre os pontos
# mais próximos, e Cria uma demanda aleatória nas arestas criadas.
function DistAdjaDemand(Coord, np)
    Distancias = [dist_euclid(Coord[i,1], Coord[i,2], Coord[j,1], Coord[j,2])
        for i in 1:np, j in 1:np];
    MenDist = fill(0.0, (np, 3))
    for i = 1:np
        vet_temp = sort(Distancias[i,:])
        MenDist[i,1:3] = vet_temp[2:4]
    end
    Adjacencias = fill(0, (np, np))
    Demanda = fill(0, (np, np))
    for i in 1:np
        for j in i:np
            if i != j
                if Distancias[i,j]>MenDist[i,3]
                    Distancias[i,j] = 1e4
                    Distancias[j,i] = 1e4
                else
                    Adjacencias[i,j] = 1
                    Adjacencias[j,i] = 1
                    Demanda[i,j] = rand(50:100)
                    Demanda[j,i] = Demanda[i,j]
                end
            end
        end
    end
    return Distancias, Adjacencias, Demanda
end
##
# Cria o Plot das adjacências criadas pela função anterior
function plotAdja(Adjacencias, Coord, np)
    pl1=plot()
    for i = 1:np, j = 1:np
            if Adjacencias[i,j] == 1
                plot!(pl1, [Coord[i,1], Coord[j,1]], [Coord[i,2], Coord[j,2]], c=:blue,
                      l=:arrow)
            end
    end
    scatter!(pl1, Coord[:,1], Coord[:,2], c="black", leg = false, ms = 3)
    return pl1
end
##
# Função que verifica a factibilidade do problema, comparando arestas e demandas com a capacidade
# e número de caminhões e quantidade de arestas possíveis de visitação.
function VerificaCenario(Distancias, Demanda, Adjacencias, W, nIJ, nT, nV)
    C = Distancias
    q = Demanda
    adja = Adjacencias

    if sum(adja) <= nT*nV
        println("Arestas - OK")
    else
        println("Mais arestas a serem percorridas que Veículos")
        return 1
    end

    if sum(q/2) <= W*nV # se True, Factível
        println("Capacidade - OK")
    else
        println("Demanda maior que capacidade dos caminhões")
        return 2
    end
    return 0
end
##
using JuMP, Gurobi
##
function ResolveModelo1(Distancias, Demanda, Adjacencias, W, nIJ, nT, nP, TLimit=180)
    C = Distancias
    q = Demanda
    adja = Adjacencias

    model = Model(optimizer_with_attributes(Gurobi.Optimizer, "TimeLimit" => TLimit))
    @variable(model, x[i=1:nIJ, j=1:nIJ, t=1:nT, p=1:nP] ≥ 0, Int) #Int
    @variable(model, xt[i=1:nIJ, j=1:nIJ, p=1:nP] ≥ 0)
    @variable(model, l[i=1:nIJ, j=1:nIJ, t=1:nT, p=1:nP] ≥ 0, Bin) #Bin
    @variable(model, lt[i=1:nIJ, j=1:nIJ, p=1:nP] ≥ 0)
    @variable(model, f[i=1:nIJ, j=1:nIJ, p=1:nP] ≥ 0) # Incluir t?

    num_var_total = length(all_variables(model))

    @objective(model, Min, sum(x[i,j,t,p] * (1+t/nT) * C[i,j] for i=1:nIJ, j=1:nIJ, t=1:nT, p=1:nP));
    @constraint(model, fluxo[i=1:nIJ,p=1:nP],
                sum(x[i,j,t,p] for j=1:nIJ, t=1:nT) ==
                sum(x[j,i,t,p] for j=1:nIJ, t=1:nT))
    @constraint(model, temp_x[i=1:nIJ,j=1:nIJ,p=1:nP],
                sum(x[i,j,t,p] for t=1:nT) == xt[i,j,p])
    @constraint(model, temp_x2[t=1:nT,p=1:nP],
                sum(x[i,j,t,p] for i=1:nIJ, j=1:nIJ) <=1)
    @constraint(model, arc_atend[i=1:nIJ,j=1:nIJ],
                sum(l[i,j,t,p] + l[j,i,t,p] for p=1:nP, t=1:nT) == adja[i,j])
    @constraint(model, temp_l[i=1:nIJ, j=1:nIJ, p=1:nP],
                sum(l[i,j,t,p] for t=1:nT) == lt[i,j,p])
    @constraint(model, pass_demanda[i=1:nIJ,j=1:nIJ, p=1:nP],
                sum(x[i,j,t,p] for t=1:nT) >=
                sum(l[i,j,t,p] for t=1:nT))
    @constraint(model, capac[p=1:nP],
                sum(l[i,j,t,p]*q[i,j] for i=1:nIJ,j=1:nIJ,t=1:nT) <= W)
    @constraint(model, contr_fluxo[i=2:nIJ,p=1:nP],
                sum(f[i,k,p] for k=1:nIJ) - sum(f[k,i,p] for k=1:nIJ) ==
                sum(lt[i,j,p] for j=1:nIJ))
    @constraint(model, contr_fluxo2[i=1:nIJ, j=1:nIJ ,p=1:nP],
                f[i,j,p] <= nIJ^2*xt[i,j,p]);
    optimize!(model)

    if termination_status(model) == MOI.OPTIMAL
        SolOt = round.(Int, value.(xt))
        SolStatus = "ótima"
        ObjOt = objective_value(model)
        SimplexIter = simplex_iterations(model)
        Gap = relative_gap(model)
        SolTime = solve_time(model)
    elseif termination_status(model) == MOI.TIME_LIMIT && has_values(model)
        SolOt = round.(Int, value.(xt))
        SolStatus = "Limitada pelo tempo"
        ObjOt = objective_value(model)
        SimplexIter = simplex_iterations(model)
        Gap = relative_gap(model)
        SolTime = solve_time(model)
    else
        error("The model was not solved correctly.")
    end

    return SolOt, num_var_total, SolOt, SolStatus, ObjOt, SimplexIter, Gap, SolTime
end
##
function ResolveModelo2(Distancias, Demanda, Adjacencias, W, nIJ, nT, nP, TLimit=300)
    C = Distancias
    q = Demanda
    adja = Adjacencias

    model = Model(optimizer_with_attributes(Gurobi.Optimizer, "TimeLimit" => TLimit))
    @variable(model, x[i=1:nIJ, j=1:nIJ, t=1:nT, p=1:nP] ≥ 0, Int) #Int
    @variable(model, xt[i=1:nIJ, j=1:nIJ, p=1:nP] ≥ 0)
    @variable(model, l[i=1:nIJ, j=1:nIJ, t=1:nT, p=1:nP] ≥ 0, Bin) #Bin
    @variable(model, lt[i=1:nIJ, j=1:nIJ, p=1:nP] ≥ 0)
    @variable(model, f[i=1:nIJ, j=1:nIJ, p=1:nP] ≥ 0) # Incluir t?

    num_var_total = length(all_variables(model))

    @expression(model, pco[i=1:nIJ, j=1:nIJ, t=1:nT, p=1:nP],
        sum(l[i,j,t-1,p] * q[i,j] for i=1:nIJ, j=1:nIJ, t=2:nT, p=1:nP))

    @objective(model, Min, sum(x[i,j,t,p] * C[i,j] *
        (1+ pco[i,j,t,p]/W) for i=1:nIJ, j=1:nIJ, t=1:nT, p=1:nP));

    @constraint(model, fluxo[i=1:nIJ,p=1:nP],
                sum(x[i,j,t,p] for j=1:nIJ, t=1:nT) ==
                sum(x[j,i,t,p] for j=1:nIJ, t=1:nT))
    @constraint(model, temp_x[i=1:nIJ,j=1:nIJ,p=1:nP],
                sum(x[i,j,t,p] for t=1:nT) == xt[i,j,p])
    @constraint(model, temp_x2[t=1:nT,p=1:nP],
                sum(x[i,j,t,p] for i=1:nIJ, j=1:nIJ) <=1)
    @constraint(model, arc_atend[i=1:nIJ,j=1:nIJ],
                sum(l[i,j,t,p] + l[j,i,t,p] for p=1:nP, t=1:nT) == adja[i,j])
    @constraint(model, temp_l[i=1:nIJ, j=1:nIJ, p=1:nP],
                sum(l[i,j,t,p] for t=1:nT) == lt[i,j,p])
    @constraint(model, pass_demanda[i=1:nIJ,j=1:nIJ, p=1:nP],
                sum(x[i,j,t,p] for t=1:nT) >=
                sum(l[i,j,t,p] for t=1:nT))
    @constraint(model, capac[p=1:nP],
                sum(l[i,j,t,p]*q[i,j] for i=1:nIJ,j=1:nIJ,t=1:nT) <= W)
    @constraint(model, contr_fluxo[i=2:nIJ,p=1:nP],
                sum(f[i,k,p] for k=1:nIJ) - sum(f[k,i,p] for k=1:nIJ) ==
                sum(lt[i,j,p] for j=1:nIJ))
    @constraint(model, contr_fluxo2[i=1:nIJ, j=1:nIJ ,p=1:nP],
                f[i,j,p] <= nIJ^2*xt[i,j,p]);
    optimize!(model)

    if termination_status(model) == MOI.OPTIMAL
        SolOt = round.(Int, value.(xt))
        SolStatus = "ótima"
        ObjOt = objective_value(model)
        SimplexIter = simplex_iterations(model)
        Gap = relative_gap(model)
        SolTime = solve_time(model)
    elseif termination_status(model) == MOI.TIME_LIMIT && has_values(model)
        SolOt = round.(Int, value.(xt))
        SolStatus = "Limitada pelo tempo"
        ObjOt = objective_value(model)
        SimplexIter = simplex_iterations(model)
        Gap = relative_gap(model)
        SolTime = solve_time(model)
    else
        error("The model was not solved correctly.")
    end

    return SolOt, num_var_total, SolOt, SolStatus, ObjOt, SimplexIter, Gap, SolTime
end
##
# Função que plota os resultados do modelo
function plotResult(x, Coord, nIJ, nP)
    Colors = [:red, :blue, :green, :yellow, :magenta]
    pl2=plot()
    for p=1:nP
        for i = 1:nIJ, j = 1:nIJ
            if (x[i,j,p] > 0)
                pl2=plot!([Coord[i,1]+0.2*p, Coord[j,1]+0.2*p],
                          [Coord[i,2]+0.2*p, Coord[j,2]+0.2*p],
                          ls =:dash, c=Colors[p], lw = 2, lab="")
            end
        end
    end
    pl2=scatter!(Coord[:,1], Coord[:,2], c="black", leg = false, ms = 3)
    pl2
end
## Cria o Data-Frame das soluções
using DataFrames
Tabela_Resultados = DataFrame(
    Cenários = String[],
    Adjac = Int64[],
    Demanda = Float64[],
    Capac_Viagens = Int64[],
    Periodos = Int64[],
    Veiculos = Int64[],
    Variaveis = Int64[],
    Simplex_iters = Float64[],
    Time_elapsed = Float64[],
    Solucao = String[],
    Gap = Float64[])
##

## Cenário 1
distMin = 10
npts = [3, 3]
Coord1, np1 = gerapontos(distMin, npts)
Distancias1, Adjacencias1, Demanda1 = DistAdjaDemand(Coord1, np1)
W = 500
nIJ = np1
nT = 16
nV = 2
VerificaCenario(Distancias1, Demanda1, Adjacencias1, W, nIJ, nT, nV)
##
sum(Adjacencias1), sum(Demanda1)/2
##
sol_x, num_var_total, SolOt, SolStatus, ObjOt, SimplexIter, Gap, SolTime =
    ResolveModelo2(Distancias1, Demanda1, Adjacencias1, W, nIJ, nT, nV);
##
push!(Tabela_Resultados, ["Cenário 1", sum(Adjacencias1), sum(Demanda1)/2,
                          W, nT, nV, num_var_total, SimplexIter, SolTime, SolStatus, Gap])
##
show(Tabela_Resultados, allcols=true)
##
plot(plotAdja(Adjacencias1, Coord1, np1), ratio =1 )
##
plot(plotResult(sol_x, Coord1, nIJ, nV), ratio = 1)
##

## Cenário 2
distMin = 10
npts = [5, 2]
Coord1, np1 = gerapontos(distMin, npts)
Distancias1, Adjacencias1, Demanda1 = DistAdjaDemand(Coord1, np1)
W = 600
nIJ = np1
nT = 16
nV = 2
VerificaCenario(Distancias1, Demanda1, Adjacencias1, W, nIJ, nT, nV)
##
sum(Adjacencias1), sum(Demanda1)/2
##
sol_x, num_var_total, SolOt, SolStatus, ObjOt, SimplexIter, Gap, SolTime =
    ResolveModelo2(Distancias1, Demanda1, Adjacencias1, W, nIJ, nT, nV);
##
push!(Tabela_Resultados, ["Cenário 2", sum(Adjacencias1), sum(Demanda1)/2,
                          W, nT, nV, num_var_total, SimplexIter, SolTime, SolStatus, Gap])
##
show(Tabela_Resultados, allcols=true)
##
plot(plotAdja(Adjacencias1, Coord1, np1), ratio = 1)
##
plot(plotResult(sol_x, Coord1, nIJ, nV), ratio = 1)
##

## Cenário 3
distMin = 10
npts = [3, 4]
Coord1, np1 = gerapontos(distMin, npts)
Distancias1, Adjacencias1, Demanda1 = DistAdjaDemand(Coord1, np1)
W = 700
nIJ = np1
nT = 19
nV = 2
VerificaCenario(Distancias1, Demanda1, Adjacencias1, W, nIJ, nT, nV)
##
sum(Adjacencias1), sum(Demanda1)/2
##
sol_x, num_var_total, SolOt, SolStatus, ObjOt, SimplexIter, Gap, SolTime =
    ResolveModelo2(Distancias1, Demanda1, Adjacencias1, W, nIJ, nT, nV);
##
push!(Tabela_Resultados, ["Cenário 3", sum(Adjacencias1), sum(Demanda1)/2,
                          W, nT, nV, num_var_total, SimplexIter, SolTime, SolStatus, Gap])
##
show(Tabela_Resultados, allcols=true)
##
plot(plotAdja(Adjacencias1, Coord1, np1), ratio = 1)
##
plot(plotResult(sol_x, Coord1, nIJ, nV), ratio = 1)
##

using Latexify
latexify(Tabela_Resultados[:, 1:7]; env=:table, latex=false)
##
df = select(Tabela_Resultados, Not(["Veiculos", "Adjac", "Demanda",
        "Capac_Viagens", "Periodos", "Variaveis"]))
latexify(df; env=:table, latex=false)

## Cenário 4
distMin = 10
npts = [3,5]
Coord1, np1 = gerapontos(distMin, npts)
Distancias1, Adjacencias1, Demanda1 = DistAdjaDemand(Coord1, np1)
W = 1600
nIJ = np1
nT = 44
nV = 1
VerificaCenario(Distancias1, Demanda1, Adjacencias1, W, nIJ, nT, nV)
##
sum(Adjacencias1), sum(Demanda1)/2
##
sol_x, num_var_total, SolOt, SolStatus, ObjOt, SimplexIter, Gap, SolTime =
    ResolveModelo2(Distancias1, Demanda1, Adjacencias1, W, nIJ, nT, nV);
##
push!(Tabela_Resultados, ["Cenário 4", sum(Adjacencias1), sum(Demanda1)/2,
                          W, nT, nV, num_var_total, SimplexIter, SolTime, SolStatus, Gap])
##
show(Tabela_Resultados, allcols=true)
##
plot(plotAdja(Adjacencias1, Coord1, np1), plotResult(sol_x, Coord1, nIJ, nV), ratio = 1)
##

## Cenário 5- Rodar só esse já dá pau no Kernel
distMin = 10
npts = [3,5]
Coord1, np1 = gerapontos(distMin, npts)
Distancias1, Adjacencias1, Demanda1 = DistAdjaDemand(Coord1, np1)
W = 900
nIJ = np1
nT = 23
nV = 2
VerificaCenario(Distancias1, Demanda1, Adjacencias1, W, nIJ, nT, nV)
##
sum(Adjacencias1), sum(Demanda1)/2
##
sol_x, num_var_total, SolOt, SolStatus, ObjOt, SimplexIter, Gap, SolTime =
    ResolveModelo2(Distancias1, Demanda1, Adjacencias1, W, nIJ, nT, nV);
##
push!(Tabela_Resultados, ["Cenário 5", sum(Adjacencias1), sum(Demanda1)/2,
                          W, nT, nV, num_var_total, SimplexIter, SolTime, SolStatus, Gap])
##
show(Tabela_Resultados, allcols=true)
##
plot(plotAdja(Adjacencias1, Coord1, np1), plotResult(sol_x, Coord1, nIJ, nV), ratio = 1)
##

max_stars_repo_path: Projeto 1 - Modelo2.jl

max_issues_repo_path: Projeto 1 - Modelo2.jl

max_forks_repo_path: Projeto 1 - Modelo2.jl


text:
module TDD
include("graph.jl")
include("interval.jl")

# Write your package code here.

end

ext: jl


text:
using Gadfly
using Colors
using Base.Test

# This should throw a MethodError, not a stackoverflow
@test_throws MethodError Scale.color_discrete_manual([RGBA(0,0,0,0.1)])

plot([sin, cos], 0, 25, Scale.color_discrete_manual("red", "purple"))


text:
import Cartesian

function testcounter(m, n, N)
    sloop = 0
    tbase = @elapsed begin
        for k = 1:N
            for j = 1:n
                for i = 1:m
                    sloop += 1
                end
            end
        end
    end
    sz = [m,n]
    scounter = 0
    t = @elapsed begin
        for k = 1:N
            @Cartesian.forcartesian c sz begin
                scounter += 1
            end
        end
    end
    @assert scounter == sloop
    tbase, t
end

testcounter(2,2,1)
tbase, t = testcounter(500,300,100)
# @assert t < 10tbase

include("testloops.jl")


text:
struct AnimationTrack{T}
    name::String
    jstype::String
    events::Vector{Pair{Int, T}} # frame => value
end

AnimationTrack{T}(name::String, jstype::String) where {T} = AnimationTrack(name, jstype, Pair{Int, T}[])

function Base.insert!(track::AnimationTrack, frame::Integer, value)
    i = searchsortedfirst(track.events, frame; by=first)
    if i <= length(track.events) && first(track.events[i]) == frame
        track.events[i] = frame => value
    else
        insert!(track.events, i, frame => value)
    end
    return track
end

function Base.merge!(a::AnimationTrack{T}, others::AnimationTrack{T}...) where T
    l = length(a.events)
    for other in others
        @assert other.name == a.name
        @assert other.jstype == a.jstype
        l += length(other.events)
    end
    events = similar(a.events, l)
    mergesorted!(events, a.events, Iterators.flatten((other.events for other in others)))
    combine!(events, by=first, combine=(a, b) -> b)
    resize!(a.events, length(events))
    copyto!(a.events, events)
    return a
end

@with_kw struct AnimationClip
    tracks::Dict{String, AnimationTrack} = Dict{String, AnimationTrack}()
    fps::Int = 30
    name::String = "default"
end

function Base.merge!(a::AnimationClip, others::AnimationClip...)
    for other in others
        @assert other.fps == a.fps
        merge!(merge!, a.tracks, other.tracks) # merge tracks recursively
    end
    return a
end

struct Animation
    clips::Dict{Path, AnimationClip}
    default_framerate::Int
end

"""
Create a new animation. See [`atframe`](@ref) to adjust object poses or properties
in that animation.

$(TYPEDSIGNATURES)
"""
Animation(fps::Int=30) = Animation(Dict{Path, AnimationClip}(), fps)

"""
Merge multiple animations, storing the result in `a`.

$(TYPEDSIGNATURES)
"""
function Base.merge!(a::Animation, others::Animation...)
    for other in others
        @assert a.default_framerate == other.default_framerate
        merge!(merge!, a.clips, other.clips) # merge clips recursively
    end
    return a
end

"""
Merge two or more animations, returning a new animation.

$(TYPEDSIGNATURES)

The animations may involve the same properties or different properties
(animations of the same property on the same path will have their events
interleaved). All animations must have the same framerate.
"""
Base.merge(a::Animation, others::Animation...) = merge!(Animation(a.default_framerate), a, others...)

"""
Convert the `.tar` file of still images produced by the meshcat "record" feature
into a video.

$(TYPEDSIGNATURES)

To record an animation which has been played in the meshcat visualizer, click
"Open Controls", then navigate to the "Animations" folder and click "record".
This will step through the animation frame-by-frame and then prompt you to
download a `.tar` file containing the resulting frames. You can then pass the
path to that `.tar` file to this function to produce a video.

This uses FFMPEG.jl internally, so you do *not* need to have installed ffmpeg
manually.
"""
function convert_frames_to_video(tar_file_path::AbstractString, output_path::AbstractString="output.mp4"; framerate=60, overwrite=false)
    output_path = abspath(output_path)
    if !isfile(tar_file_path)
        error("Could not find the input file $tar_file_path")
    end
    if isfile(output_path) && !overwrite
        error("The output path $output_path already exists. To overwrite that file, you can pass `overwrite=true` to this function")
    end

    mktempdir() do tmpdir
        run(unpack_cmd(tar_file_path, tmpdir, ".tar", nothing))
        cmd = ["-r", string(framerate), "-i", "%07d.png", "-vcodec", "libx264", "-preset", "slow", "-crf", "18"]
        if overwrite
            push!(cmd, "-y")
        end
        push!(cmd, output_path)

        cd(tmpdir) do
            FFMPEG.exe(cmd...)
        end
    end
    @info("Saved output as $output_path")
    return output_path
end


text:
export open_video_out, close_video_out!, get_codec_name

mutable struct VideoWriter{T<:GraphType}
    format_context::AVFormatContextPtr
    codec_context::AVCodecContextPtr
    frame_graph::T
    packet::AVPacketPtr
    stream_index0::Int
    scanline_major::Bool
    next_index::Int
end

graph_input_frame(r::VideoWriter) = graph_input_frame(r.frame_graph)

graph_output_frame(r::VideoWriter) = graph_output_frame(r.frame_graph)

isopen(w::VideoWriter) = check_ptr_valid(w.format_context, false)

function encode_mux!(writer::VideoWriter, flush = false)
    pkt = writer.packet
    av_init_packet(pkt)
    frame = graph_output_frame(writer)
    if flush
        fret = avcodec_send_frame(writer.codec_context, C_NULL)
    else
        fret = avcodec_send_frame(writer.codec_context, frame)
    end
    if fret < 0 && !in(fret, [-Libc.EAGAIN, VIO_AVERROR_EOF])
        error("Error $fret sending a frame for encoding")
    end

    pret = Cint(0)
    while pret >= 0
        pret = avcodec_receive_packet(writer.codec_context, pkt)
        if pret == -Libc.EAGAIN || pret == VIO_AVERROR_EOF
            break
        elseif pret < 0
            error("Error $pret during encoding")
        end
        if pkt.duration == 0
            codec_pts_duration = round(
                Int,
                1 / (
                    convert(Rational, writer.codec_context.framerate) *
                    convert(Rational, writer.codec_context.time_base)
                ),
            )
            pkt.duration = codec_pts_duration
        end
        stream = writer.format_context.streams[writer.stream_index0+1]
        av_packet_rescale_ts(pkt, writer.codec_context.time_base, stream.time_base)
        pkt.stream_index = writer.stream_index0
        ret = av_interleaved_write_frame(writer.format_context, pkt)
        # No packet_unref, av_interleaved_write_frame now owns packet.
        ret != 0 && error("Error muxing packet")
    end
    if !flush && fret == -Libc.EAGAIN && pret != VIO_AVERROR_EOF
        fret = avcodec_send_frame(writer.codec_context, frame)
        if fret < 0 && fret != VIO_AVERROR_EOF
            error("Error $fret sending a frame for encoding")
        end
    end
    return pret
end

function prepare_video_frame!(writer::VideoWriter, img, index)
    dstframe = graph_output_frame(writer)
    ret = av_frame_make_writable(dstframe)
    ret < 0 && error("av_frame_make_writable() error")
    dstframe.pts = index
    transfer_img_buf_to_frame!(graph_input_frame(writer), img, writer.scanline_major)
    return execute_graph!(writer)
end

execute_graph!(writer::VideoWriter) = exec!(writer.frame_graph)

"""
    write(writer::VideoWriter, img)
    write(writer::VideoWriter, img, index)

Prepare frame `img` for encoding, encode it, mux it, and either cache it or
write it to the file described by `writer`. `img` must be the same size and
element type as the size and element type that was used to create `writer`.
If `index` is provided, it must start at zero and increment monotonically.
"""
function write(writer::VideoWriter, img, index::Int)
    isopen(writer) || error("VideoWriter is closed for writing")
    prepare_video_frame!(writer, img, index)
    encode_mux!(writer)
    return writer.next_index = index + 1
end
write(writer::VideoWriter, img) = write(writer::VideoWriter, img, writer.next_index)

"""
    close_video_out!(writer::VideoWriter)

Write all frames cached in `writer` to the video container that it describes,
and then close the file. Once all frames in a video have been added to `writer`,
then it must be closed with this function to flush any cached frames to the file,
and then finally close the file and release resources associated with `writer`.
"""
function close_video_out!(writer::VideoWriter)
    if check_ptr_valid(writer.format_context, false)
        encode_mux!(writer, true) # flush
        format_context = writer.format_context
        ret = av_write_trailer(format_context)
        ret < 0 && error("Could not write trailer")
        if format_context.oformat.flags & AVFMT_NOFILE == 0
            pb_ptr = field_ptr(format_context, :pb)
            ret = @preserve writer avio_closep(pb_ptr)
            ret < 0 && error("Could not free AVIOContext")
        end
    end
    # Free allocated memory through finalizers
    writer.format_context = AVFormatContextPtr(C_NULL)
    writer.codec_context = AVCodecContextPtr(C_NULL)
    writer.frame_graph = null_graph(writer.frame_graph)
    writer.packet = AVPacketPtr(C_NULL)
    writer.stream_index0 = -1
    return writer
end

function get_array_from_avarray(ptr::Union{Ptr{T},Ref{T},NestedCStruct{T}}, term; make_copy = true) where {T}
    check_ptr_valid(ptr, false) || return Vector{T}()
    i = 1
    el = ptr[i]
    while el != term
        i += 1
        el = ptr[i]
    end
    len = i - 1
    if make_copy
        dst = Vector{T}(undef, len)
        unsafe_copyto!(dst, ptr, len)
    else
        dst = unsafe_wrap(Array, ptr, len)
    end
    return dst
end

function determine_best_encoding_format(target_pix_fmt, transfer_pix_fmt, codec, loss_flags = 0)
    if target_pix_fmt === nothing
        @preserve codec begin
            encoding_pix_fmt, losses =
                _vio_determine_best_pix_fmt(transfer_pix_fmt, codec.pix_fmts; loss_flags = loss_flags)
        end
    else
        @preserve codec begin
            codec_pix_fmts = get_array_from_avarray(codec.pix_fmts, AV_PIX_FMT_NONE; make_copy = false)
            codec_name = unsafe_string(codec.name)
            target_pix_fmt in codec_pix_fmts ||
                throw(ArgumentError("Pixel format $target_pix_fmt not compatible with codec $codec_name"))
        end
        encoding_pix_fmt = target_pix_fmt
    end
    return encoding_pix_fmt
end

function create_encoding_frame_graph(
    transfer_pix_fmt,
    encoding_pix_fmt,
    width,
    height,
    transfer_colorspace_details,
    dst_color_primaries,
    dst_color_trc,
    dst_colorspace,
    dst_color_range,
    use_vio_gray_transform,
    swscale_options;
    sws_color_options::OptionsT = (;),
)
    if use_vio_gray_transform
        frame_graph = GrayTransform()
        set_basic_frame_properties!(frame_graph.srcframe, width, height, transfer_pix_fmt)
        bit_depth, padded_bits_per_pixel = pix_fmt_to_bits_per_pixel(encoding_pix_fmt)
        frame_graph.src_depth = frame_graph.dst_depth = bit_depth
        frame_graph.srcframe.color_range = transfer_colorspace_details.color_range
    elseif transfer_pix_fmt == encoding_pix_fmt
        frame_graph = AVFramePtr()
    else
        frame_graph = SwsTransform(
            width,
            height,
            transfer_pix_fmt,
            transfer_colorspace_details.color_primaries,
            transfer_colorspace_details.color_range,
            encoding_pix_fmt,
            dst_color_primaries,
            dst_color_range,
            sws_color_options,
            swscale_options,
        )
        set_basic_frame_properties!(frame_graph.srcframe, width, height, transfer_pix_fmt)
    end
    dstframe = graph_output_frame(frame_graph)
    dstframe.color_range = dst_color_range
    dstframe.colorspace = dst_colorspace
    dstframe.color_trc = dst_color_trc
    dstframe.color_primaries = dst_color_primaries
    set_basic_frame_properties!(dstframe, width, height, encoding_pix_fmt)
    return frame_graph
end

function maybe_configure_codec_context_colorspace_details!(codec_context, colorspace_details)
    if codec_context.colorspace == AVCOL_SPC_UNSPECIFIED
        codec_context.colorspace = colorspace_details.colorspace
    end
    if codec_context.color_primaries == AVCOL_PRI_UNSPECIFIED
        codec_context.color_primaries = colorspace_details.color_primaries
    end
    if codec_context.color_trc == AVCOL_TRC_UNSPECIFIED
        codec_context.color_trc = colorspace_details.color_trc
    end
    if codec_context.color_range == AVCOL_RANGE_UNSPECIFIED
        codec_context.color_range = colorspace_details.color_range
    end
    return nothing
end

"""
    get_codec_name(writer) -> String

Get the name of the encoder used by `writer`.
"""
function get_codec_name(w::VideoWriter)
    if check_ptr_valid(w.codec_context, false) && check_ptr_valid(w.codec_context.codec, false)
        return unsafe_string(w.codec_context.codec.name)
    else
        return "None"
    end
end

function VideoWriter(
    filename::AbstractString,
    ::Type{T},
    sz::NTuple{2,Integer};
    codec_name::Union{AbstractString,Nothing} = nothing,
    framerate::Real = 24,
    scanline_major::Bool = false,
    container_options::OptionsT = (;),
    container_private_options::OptionsT = (;),
    encoder_options::OptionsT = (;),
    encoder_private_options::OptionsT = (;),
    swscale_options::OptionsT = (;),
    target_pix_fmt::Union{Nothing,Cint} = nothing,
    pix_fmt_loss_flags = 0,
    input_colorspace_details = nothing,
    allow_vio_gray_transform = true,
    sws_color_options::OptionsT = (;),
    thread_count::Union{Nothing,Int} = nothing,
) where {T}
    framerate > 0 || error("Framerate must be strictly positive")

    if haskey(encoder_options, :priv_data)
        throw(
            ArgumentError(
                """The field `priv_data` is no longer supported. Either reorganize as a flat NamedTuple or Dict,
e.g. `encoder_options=(color_range=2, crf=0, preset=\"medium\")` to rely on auto routing of generic and private
options, or pass the private options to `encoder_private_options` explicitly""",
            ),
        )
    end
    if !is_eltype_transfer_supported(T)
        throw(ArgumentError("Encoding arrays with eltype $T not yet supported"))
    end
    if scanline_major
        width, height = sz
    else
        height, width = sz
    end
    if isodd(width) || isodd(height)
        throw(ArgumentError("Encoding error: Frame dims must be a multiple of two"))
    end
    transfer_pix_fmt = get_transfer_pix_fmt(T)

    if input_colorspace_details === nothing
        transfer_colorspace_details =
            get(VIO_DEFAULT_TRANSFER_COLORSPACE_DETAILS, transfer_pix_fmt, VIO_DEFAULT_COLORSPACE_DETAILS)
    else
        transfer_colorspace_details = input_colorspace_details
    end

    format_context = output_AVFormatContextPtr(filename)
    default_codec = codec_name === nothing
    if default_codec
        codec_enum = format_context.oformat.video_codec
        codec_enum == AV_CODEC_ID_NONE && error("No default codec found")
        codec_p = avcodec_find_encoder(codec_enum)
        check_ptr_valid(codec_p, false) || error("Default codec not found")
    else
        codec_p = avcodec_find_encoder_by_name(codec_name)
        check_ptr_valid(codec_p, false) || error("Codec '$codec_name' not found")
    end
    codec = AVCodecPtr(codec_p)
    if !check_ptr_valid(codec.pix_fmts, false)
        error("Codec has no supported pixel formats")
    end
    encoding_pix_fmt = determine_best_encoding_format(target_pix_fmt, transfer_pix_fmt, codec, pix_fmt_loss_flags)
    if !default_codec
        ret = avformat_query_codec(format_context.oformat, codec.id, libffmpeg.FF_COMPLIANCE_NORMAL)
        ret == 1 || error("Format not compatible with codec $codec_name")
    end
    codec_context = AVCodecContextPtr(codec)
    codec_context.width = width
    codec_context.height = height
    framerate_rat = Rational(framerate)
    target_timebase = 1 / framerate_rat
    codec_context.time_base = target_timebase
    codec_context.framerate = framerate_rat
    codec_context.pix_fmt = encoding_pix_fmt
    codec_context.thread_count = thread_count === nothing ? codec_context.thread_count : thread_count

    if format_context.oformat.flags & AVFMT_GLOBALHEADER != 0
        codec_context.flags |= AV_CODEC_FLAG_GLOBAL_HEADER
    end

    set_class_options(format_context, container_options)
    if check_ptr_valid(format_context.oformat.priv_class, false)
        set_class_options(format_context.priv_data, container_private_options)
    elseif !isempty(container_private_options)
        @warn "This container format does not support private options, and will be ignored"
    end
    set_class_options(codec_context, encoder_options)
    set_class_options(codec_context.priv_data, encoder_private_options)

    sigatomic_begin()
    lock(VIO_LOCK)
    ret = avcodec_open2(codec_context, codec, C_NULL)
    unlock(VIO_LOCK)
    sigatomic_end()
    ret < 0 && error("Could not open codec: Return code $(ret)")

    stream_p = avformat_new_stream(format_context, C_NULL)
    check_ptr_valid(stream_p, false) || error("Could not allocate output stream")
    stream = AVStreamPtr(stream_p)
    stream_index0 = stream.index
    stream.time_base = codec_context.time_base
    stream.avg_frame_rate = 1 / convert(Rational, stream.time_base)
    stream.r_frame_rate = stream.avg_frame_rate
    ret = avcodec_parameters_from_context(stream.codecpar, codec_context)
    ret < 0 && error("Could not set parameters of stream")

    if format_context.oformat.flags & AVFMT_NOFILE == 0
        pb_ptr = field_ptr(format_context, :pb)
        ret = @preserve format_context avio_open(pb_ptr, filename, AVIO_FLAG_WRITE)
        ret < 0 && error("Could not open file $filename for writing")
    end
    avformat_write_header(format_context, C_NULL)
    ret < 0 && error("Could not write header")

    use_vio_gray_transform =
        transfer_pix_fmt == encoding_pix_fmt &&
        allow_vio_gray_transform &&
        transfer_pix_fmt in VIO_GRAY_SCALE_TYPES &&
        transfer_colorspace_details.color_range != codec_context.color_range
    if !use_vio_gray_transform && transfer_pix_fmt == encoding_pix_fmt
        maybe_configure_codec_context_colorspace_details!(codec_context, transfer_colorspace_details)
    end
    frame_graph = create_encoding_frame_graph(
        transfer_pix_fmt,
        encoding_pix_fmt,
        width,
        height,
        transfer_colorspace_details,
        codec_context.color_primaries,
        codec_context.color_trc,
        codec_context.colorspace,
        codec_context.color_range,
        use_vio_gray_transform,
        swscale_options;
        sws_color_options = sws_color_options,
    )
    packet = AVPacketPtr()

    return VideoWriter(format_context, codec_context, frame_graph, packet, Int(stream_index0), scanline_major, 0)
end

VideoWriter(filename, img::AbstractMatrix{T}; kwargs...) where {T} =
    VideoWriter(filename, img_params(img)...; kwargs...)

"""
    open_video_out(filename, ::Type{T}, sz::NTuple{2, Integer};
                   <keyword arguments>) -> writer
    open_video_out(filename, first_img::Matrix; ...)
    open_video_out(f, ...; ...)

Open file `filename` and prepare to encode a video stream into it, returning
object `writer` that can be used to encode frames. The size and element type of
the video can either be specified by passing the first frame of the movie
`first_img`, which will not be encoded, or instead the element type `T` and
2-tuple size `sz`. If the size is explicitly specified, the first element will
be the height, and the second width, unless keyword argument `scanline_major = true`, in which case the order is reversed. Both height and width must be even.
The element type `T` must be one of the supported element types, which is any
key of `VideoIO.VIO_DEF_ELTYPE_PIX_FMT_LU`, or instead the `Normed` or
`Unsigned` type for a corresponding `Gray` element type. The container type will
be inferred from `filename`.

Frames are encoded with[ `write`](@ref), which must use frames with
the same size, element type, and obey the same value of `scanline_major`. The
video must be closed once all frames are encoded with
[`close_video_out!`](@ref).

If called with a function as the first argument, `f`, then the function will be
called with the writer object `writer` as its only argument. This `writer` object
will be closed once the call is complete, regardless of whether or not an error
occurred.

# Keyword arguments

  - `codec_name::Union{AbstractString, Nothing} = nothing`: Name of the codec to
    use. Must be a name accepted by [FFmpeg](https://ffmpeg.org/), and
    compatible with the chosen container type, or `nothing`, in which case the
    codec will be automatically selected by FFmpeg based on the container.
  - `framerate::Real = 24`: Framerate of the resulting video.
  - `scanline_major::Bool = false`: If `false`, then Julia arrays are assumed to
    have frame height in the first dimension, and frame width on the second. If
    `true`, then pixels that adjacent to eachother in the same scanline (i.e.
    horizontal line of the video) are assumed to be adjacent to eachother in
    memory. `scanline_major = true` videos must be `StridedArray`s with unit
    stride in the first dimension. For normal arrays, this corresponds to a
    matrix where frame width is in the first dimension, and frame height is in
    the second.
  - `container_options::OptionsT = (;)`: A `NamedTuple` or `Dict{Symbol, Any}`
    of options for the container. Must correspond to option names and values
    accepted by [FFmpeg](https://ffmpeg.org/).
  - `container_private_options::OptionsT = (;)`: A `NamedTuple` or
    `Dict{Symbol, Any}` of private options for the container. Must correspond
    to private options names and values accepted by
    [FFmpeg](https://ffmpeg.org/) for the chosen container type.
  - `encoder_options::OptionsT = (;)`: A `NamedTuple` or `Dict{Symbol, Any}` of
    options for the encoder context. Must correspond to option names and values
    accepted by [FFmpeg](https://ffmpeg.org/).
  - `encoder_private_options::OptionsT = (;)`: A `NamedTuple` or
    `Dict{Symbol, Any}` of private options for the encoder context. Must
    correspond to private option names and values accepted by
    [FFmpeg](https://ffmpeg.org/) for the chosen codec specified by `codec_name`.
  - `swscale_options::OptionsT = (;)`: A `Namedtuple`, or `Dict{Symbol, Any}` of
    options for the swscale object used to perform color space scaling. Options
    must correspond with options for FFmpeg's
    [scaler](https://ffmpeg.org/ffmpeg-all.html#Scaler-Options) filter.
  - `target_pix_fmt::Union{Nothing, Cint} = nothing`: The pixel format that will
    be used to input data into the encoder. This can either by a
    `VideoIO.AV_PIX_FMT_*` value corresponding to a FFmpeg
    [`AVPixelFormat`]
    (https://ffmpeg.org/doxygen/4.1/pixfmt_8h.html#a9a8e335cf3be472042bc9f0cf80cd4c5),
    and must then be a format supported by the encoder, or instead `nothing`,
    in which case it will be chosen automatically by FFmpeg.
  - `pix_fmt_loss_flags = 0`: Loss flags to control how encoding pixel format is
    chosen. Only valid if `target_pix_fmt = nothing`. Flags must correspond to
    FFmpeg
    [loss flags]
    (http://ffmpeg.org/doxygen/2.5/pixdesc_8h.html#a445e6541dde2408332c216b8d0accb2d).
  - `input_colorspace_details = nothing`: Information about the color space
    of input Julia arrays. If `nothing`, then this will correspond to a
    best-effort interpretation of `Colors.jl` for the corresponding element type.
    To override these defaults, create a `VideoIO.VioColorspaceDetails` object
    using the appropriate `AVCOL_` definitions from FFmpeg, or use
    `VideoIO.VioColorspaceDetails()` to use the FFmpeg defaults. If data in the
    input Julia arrays is already in the mpeg color range, then set this to
    `VideoIO.VioColorspaceDetails()` to avoid additional scaling by `sws_scale`.
  - `allow_vio_gray_transform = true`: Instead of using `sws_scale` for gray data,
    use a more accurate color space transformation implemented in `VideoIO` if
    `allow_vio_gray_transform = true`. Otherwise, use `sws_scale`.
  - `sws_color_options::OptionsT = (;)`: Additional keyword arguments passed to
    [sws_setColorspaceDetails]
    (http://ffmpeg.org/doxygen/2.5/group__libsws.html#ga541bdffa8149f5f9203664f955faa040).
  - `thread_count::Union{Nothing, Int} = nothing`: The number of threads the codec is
    allowed to use, or `nothing` for default codec behavior. Defaults to `nothing`.

See also: [`write`](@ref), [`close_video_out!`](@ref)
"""
open_video_out(s::AbstractString, args...; kwargs...) = VideoWriter(s, args...; kwargs...)

function open_video_out(f, s::AbstractString, args...; kwargs...)
    writer = open_video_out(s, args...; kwargs...)
    try
        f(writer)
    finally
        close_video_out!(writer)
    end
end

img_params(img::AbstractMatrix{T}) where {T} = (T, size(img))

"""
    save(filename::String, imgstack; ...)

Create a video container `filename` and encode the set of frames `imgstack` into
it. `imgstack` must be an iterable of matrices and each frame must have the same
dimensions and element type.

Encoding options, restrictions on frame size and element type, and other
details are described in [`open_video_out`](@ref). All keyword arguments are
passed to `open_video_out`.

See also: [`open_video_out`](@ref), [`write`](@ref),
[`close_video_out!`](@ref)
"""
function save(filename::String, imgstack; kwargs...)
    open_video_out(filename, first(imgstack); kwargs...) do writer
        for img in imgstack
            write(writer, img)
        end
    end
    return nothing
end


text:
using ConstrainedDynamics

# Parameters
ex = [1.;0.;0.]

h = 1.
r = .05

vert11 = [0.;0.;h / 2]
vert12 = -vert11

# Initial orientation
phi = pi / 4
q1 = UnitQuaternion(RotX(phi))

# Links
origin = Origin{Float64}()
links = [Cylinder(r, h, h, color = RGBA(1., 0., 0.), name="body"*string(i)) for i = 1:2]

# Constraints
jointb1 = EqualityConstraint(Revolute(origin, links[1], ex; p2=vert11), name="jointb1")
joint12 = EqualityConstraint(Revolute(links[1], links[2], ex; p1=vert12, p2=vert11), name="joint12")
constraints = [jointb1;joint12]

mech = Mechanism(origin, links, constraints;Δt = 0.01)
setPosition!(origin, links[1], p2 = vert11,Δq = q1)
setPosition!(links[1], links[2], p1 = vert12, p2 = vert11)

deactivate!(mech,4)
activate!(mech,4)

function doublependulum_disconnection_control!(mechanism,k)
    if k==75
        deactivate!(mechanism,4)
    end
end


text:
include("../scripts/perturbation_bin_file.jl")

function generate_tao_perturbation()
    url_reference_without_correctcrust = "/work/05880/tg851791/stampede2/specfem/20190115/tao_h_files/DATABASES_MPI"
    url_tao_binfile_smooth = "/work/05880/tg851791/stampede2/model/FWEA18_smooth"
    url_output = "/scratch/05880/tg851791/binfile/perturbation/tao_perturbation_smooth"
    nproc = 336
    nspec = 4960
    generate_perturbation(url_tao_binfile_smooth, url_reference_without_correctcrust, url_output, nproc, nspec)
end

function generate_min_perturbation()
    url_reference_without_correctcrust = "/work/05880/tg851791/stampede2/specfem/20190115/min_h_files/DATABASES_MPI"
    url_min_binfile = "/work/05880/tg851791/stampede2/model/EARA2014"
    url_output = "/scratch/05880/tg851791/binfile/perturbation/min_perturbation"
    nproc = 144
    nspec = 7236
    generate_perturbation(url_min_binfile, url_reference_without_correctcrust, url_output, nproc, nspec)
end

function generate_my_perturbation()
    url_reference_with_correctcrust = "/work/05880/tg851791/stampede2/specfem/20190115/simulation_taoreg/DATABASES_MPI"
    url_background_binfile = "/work/05880/tg851791/stampede2/specfem/20190115/simulation_taoreg/DATABASES_MPI"
    url_output = "/scratch/05880/tg851791/binfile/perturbation/background_perturbation"
    nproc = 441
    nspec = 4480
    generate_perturbation(url_background_binfile, url_reference_with_correctcrust, url_output, nproc, nspec)
end

function recover_tao()
    url_reference_with_correctcrust = "/work/05880/tg851791/stampede2/specfem/20190115/simulation_taoreg/DATABASES_MPI"
    url_tao_perturbation = "/scratch/05880/tg851791/binfile/interp/s362ani_addtao-pert"
    url_output = "/scratch/05880/tg851791/binfile/gll/bg-tao"
    nproc = 441
    nspec = 4480
    generate_real(url_tao_perturbation, url_reference_with_correctcrust, url_output, nproc, nspec)
end

function recover_tao_min()
    url_reference_with_correctcrust = "/work/05880/tg851791/stampede2/specfem/20190115/simulation_taoreg/DATABASES_MPI"
    url_min_perturbation = "/scratch/05880/tg851791/binfile/interp/s362ani_addtao_addmin-pert"
    url_output = "/scratch/05880/tg851791/binfile/gll/bg-tao-min"
    nproc = 441
    nspec = 4480
    generate_real(url_min_perturbation, url_reference_with_correctcrust, url_output, nproc, nspec)
end

function recover_min()
    url_reference_with_correctcrust = "/work/05880/tg851791/stampede2/specfem/20190115/simulation_taoreg/DATABASES_MPI"
    url_min_perturbation = "/scratch/05880/tg851791/binfile/interp/s362ani_addmin-pert"
    url_output = "/scratch/05880/tg851791/binfile/gll/bg-min"
    nproc = 441
    nspec = 4480
    generate_real(url_min_perturbation, url_reference_with_correctcrust, url_output, nproc, nspec)
end

function recover_min_tao()
    url_reference_with_correctcrust = "/work/05880/tg851791/stampede2/specfem/20190115/simulation_taoreg/DATABASES_MPI"
    url_min_perturbation = "/scratch/05880/tg851791/binfile/interp/s362ani_addmin_addtao-pert"
    url_output = "/scratch/05880/tg851791/binfile/gll/bg-min-tao"
    nproc = 441
    nspec = 4480
    generate_real(url_min_perturbation, url_reference_with_correctcrust, url_output, nproc, nspec)
end

max_stars_repo_name: ziyixi/SeisScripts

max_issues_repo_name: ziyixi/SeisScripts

max_forks_repo_name: ziyixi/SeisScripts


text:
"""
    ConvergenceState

A struct that summarizes the convergence state of a solution. The fields in this struct are:
 - `Δx`: the infinity norm of the change in the solution `x`.
 - `Δf`: the change in the objective value `f`.
 - `relΔf`: the ratio of change in the objective value `f`.
 - `kkt_residual`: the Karush-Kuhn-Tucker (KKT) residual of the solution. If [`ScaledKKTCriteria`](@ref) is used instead of [`KKTCriteria`](@ref), the `kkt_residual` will be divided by a factor.
 - `ipopt_residual`: the modified KKT residual used in the IPOPT solver.
 - `infeas`: maximum infeasibility amount. This is 0 if the solution is feasible.
 - `x_converged`: true if `Δx` is less than the `x` tolerance in [`MMAOptions`](@ref).
 - `fabs_converged`: true if `Δf` is less than the `f` tolerance in [`MMAOptions`](@ref).
 - `frel_converged`: true if `relΔf` is less than the `f` tolerance in [`MMAOptions`](@ref).
 - `kkt_converged`: true if the `kkt_residual` is less than the KKT tolerance in [`MMAOptions`](@ref).
 - `ipopt_converged`: true if the `ipopt_residual` is less than the KKT tolerance in [`MMAOptions`](@ref).
 - `infeas_converged`: true if `infeas` is less than the infeasibility tolerance in [`MMAOptions`](@ref).
 - `f_increased`: true if the objective value of the current solution is higher than that of the previous solution.
 - `converged`: true if the solution satisfies the convergence criteria of choice. See [`ConvergenceCriteria`](@ref) for more on the different convergence criteria available.
"""
@with_kw mutable struct ConvergenceState{T}
    Δx::T = Inf
    Δf::T = Inf
    relΔf::T = Inf
    kkt_residual::T = Inf
    ipopt_residual::T = Inf
    infeas::T = Inf
    x_converged::Bool = false
    fabs_converged::Bool = false
    frel_converged::Bool = false
    kkt_converged::Bool = false
    ipopt_converged::Bool = false
    infeas_converged::Bool = false
    f_increased::Bool = true
    converged::Bool = false
end
ConvergenceState(::Type{T}) where {T} = ConvergenceState{T}()

"""
    hasconverged(s::ConvergenceState)

Returns true if `s.converged` is true.
"""
hasconverged(s::ConvergenceState) = s.converged

"""
    Solution

A struct that stores all the information about a solution. The following are the fields of `Solution`:
 - `x`: the current primal solution
 - `prevx`: the previous primal solution
 - `λ`: the current dual solution
 - `f`: the objective value of `x`
 - `prevf`: the obejctive value of `prevx`
 - `∇f`: the gradient of the objective at `x`
 - `g`: the constraint function value at `x`
 - `∇g`: the Jacobian of the constraint functions at `x`
 - `convstate`: the convergence state of the solution
"""
@params mutable struct Solution
    prevx
    x
    λ
    prevf
    f
    ∇f
    g
    ∇g
    convstate::ConvergenceState
end

"""
    Solution(dualmodel, λ)

Construct an empty solution for the dual model `dualmodel` given a sample dual solution `λ`.
"""
function Solution(dualmodel, λ)
    prevx = copy(getxk(dualmodel))
    x = copy(prevx)
    λ = copy(λ)
    prevf = Inf
    fg = getfk(dualmodel)
    ∇fg = get∇fk(dualmodel)
    f = fg[1]
    g = fg[2:end]
    ∇f = ∇fg[1, :]
    ∇g = ∇fg[2:end, :]
    convstate = ConvergenceState()
    return Solution(prevx, x, λ, prevf, f, ∇f, g, ∇g, convstate)
end

"""
    hasconverged(s::Solution)

Returns true if `hasconverged(s.convstate)` is true.
"""
hasconverged(s::Solution) = hasconverged(s.convstate)

"""
    Tolerance

A struct specifying the different tolerances used to assess the convergence of the algorithms. The following are the fields of `Tolerance`:
 - `x`: the tolerance for `Δx` in [`ConvergenceState`](@ref). `x_converged` will be true if `Δx` is less than the `x` tolerance in `Tolerance`. This is used to assess convergence when the [`GenericCriteria`](@ref) is used as the convergence criteria.
 - `fabs`: the tolerance for `Δf` in [`ConvergenceState`](@ref). `f_converged` will be true if `Δf` is less than the `fabs` tolerance in `Tolerance`. This is used to assess convergence when the [`GenericCriteria`](@ref) is used as the convergence criteria.
 - `frel`: the tolerance for `relΔf` in [`ConvergenceState`](@ref). `f_converged` will be true if `relΔf` is less than the `frel` tolerance in `Tolerance`. This is used to assess convergence when the [`GenericCriteria`](@ref) is used as the convergence criteria.
 - `kkt`: the KKT tolerance. `kkt_converged` in [`ConvergenceState`](@ref) will be true if the `kkt_residual` is less than the KKT tolerance. And `ipopt_converged` in [`ConvergenceState`](@ref) will be true if `ipopt_residual` is less than the KKT tolerance. This is used to assess convergence when the [`KKTCriteria`](@ref), the [`ScaledKKTCriteria`](@ref) or [`IpoptCriteria`](@ref) criteria is used as the convergence criteria.
 - `infeas`: the maximum infeasibility tolerance. `infeas_converged` in [`ConvergenceState`](@ref) will be true if the maximum infeasibility is less than the infeasibility tolerance. This is used to assess convergence regardless of the convergence criteria used.

For more on convergence criteria, see [`GenericCriteria`](@ref), [`KKTCriteria`](@ref), [`ScaledKKTCriteria`](@ref) and [`IpoptCriteria`](@ref).
"""
struct Tolerance{Tx, Tf, Tkkt, Tinfeas}
    x::Tx
    fabs::Tf
    frel::Tf
    kkt::Tkkt
    infeas::Tinfeas
end
function Tolerance(;
    x = 0.0, f = 1e-5, fabs = f, frel = f,
    kkt = 1e-5, infeas = 1e-5,
)
    Tolerance(x, fabs, frel, kkt, infeas)
end
function (tol::Tolerance{<:Function, <:Function, <:Function})(i)
    return Tolerance(tol.x(i), tol.fabs(i), tol.frel(i), tol.kkt(i), tol.infeas)
end

function Base.:*(t::Tolerance, m::Real)
    return Tolerance(t.x * m, t.fabs * m, t.frel * m, t.kkt * m, t.infeas * m)
end

abstract type ConvergenceCriteria end

"""
    GenericCriteria

This is a generic convergence criteria that uses:
1. The maximum change in the solution, `Δx`,
2. The change in the objective value, `Δf`, and
3. The change percentage in the objective value, `Δf`, and
4. The maximum infeasibility `infeas`.
to assess convergence. More details are given in [`assess_convergence!`](@ref).
"""
struct GenericCriteria <: ConvergenceCriteria end

"""
    KKTCriteria

This convergence criteria uses the Karush-Kuhn-Tucker residual and maximum infeasibility to assess convergence. More details are given in [`assess_convergence!`](@ref).
"""
struct KKTCriteria <: ConvergenceCriteria end

"""
    IpoptCriteria

This convergence criteria uses a scaled version of the Karush-Kuhn-Tucker (KKT) residual and maximum infeasibility to assess convergence. This scaled KKT residual is used in the IPOPT nonlinear programming solver as explained in [this paper](http://cepac.cheme.cmu.edu/pasilectures/biegler/ipopt.pdf). More details are given in [`assess_convergence!`](@ref).
"""
struct IpoptCriteria <: ConvergenceCriteria end

"""
    ScaledKKTCriteria

This convergence criteria uses another scaled version of the Karush-Kuhn-Tucker (KKT) residual and maximum infeasibility to assess convergence. In particular if the objective was scaled by a factor `m`, the KKT residual will be scaled down by a factor `max(m, 1/m)`. The objective scaling is explained in [`optimize!`](@ref). This scaling was found to make the convergence criteria less sensitive to scale compared to using the traditional KKT residual. More details are given in [`assess_convergence!`](@ref).
"""
struct ScaledKKTCriteria <: ConvergenceCriteria end

"""
```
assess_convergence!(
    solution::Solution,
    model::AbstractModel,
    tol::Tolerance,
    criteria::ConvergenceCriteria,
)
```

Evaluates the convergence state `solution.convstate` given the current solution, `solution`, the tolerance, `tol`, and the convergence criteria `criteria`. `solution.convstate.converged` is then updated.

If `criteria` is an instance of `GenericCriteria`, `converged = (x_converged || f_converged) && infeas_converged`. `x_converged`, `f_converged` and `infeas_converged` are explained in [`Tolerance`](@ref). If `criteria` is an instance of `KKTCriteria` or `ScaledKKTCriteria`, `converged = kkt_converged && infeas_converged`. `kkt_converged` and `infeas_converged` are explained in [`Tolerance`](@ref). If `criteria` is an instance of `IpoptCriteria`, `converged = ipopt_converged && infeas_converged`. `ipopt_converged` and `infeas_converged` are explained in [`Tolerance`](@ref).
"""
function assess_convergence!(
    solution::Solution,
    model::AbstractModel,
    tol::Tolerance,
    criteria::ConvergenceCriteria,
)
    xtol, fabstol, freltol, kkttol, infeastol = tol.x, tol.fabs, tol.frel, tol.kkt, tol.infeas
    Δx, Δf, infeas = getresiduals(solution, model, GenericCriteria())
    relΔf = Δf / (abs(solution.f) + freltol)
    kkt_residual, infeas = getresiduals(solution, model, KKTCriteria())
    ipopt_residual, infeas = getresiduals(solution, model, IpoptCriteria())
    if show_residuals[]
        @show kkt_residual, ipopt_residual, kkttol
    end

    x_converged = Δx < xtol
    fabs_converged = Δf < fabstol
    frel_converged = relΔf < freltol
    if criteria isa ScaledKKTCriteria
        if debugging[]
            #@show get_objective_multiple(model)
        end
        m = get_objective_multiple(model)
        kkt_residual = kkt_residual / max(m, 1/m)
    end
    kkt_converged = kkt_residual < kkttol
    ipopt_converged = ipopt_residual < kkttol
    infeas_converged = infeas <= infeastol
    f_increased = solution.f > solution.prevf

    if criteria isa GenericCriteria
        converged = (x_converged || fabs_converged || frel_converged) && infeas_converged
    elseif criteria isa KKTCriteria || criteria isa ScaledKKTCriteria
        converged = kkt_converged && infeas_converged
    elseif criteria isa IpoptCriteria
        converged = ipopt_converged && infeas_converged
    else
        throw("Unsupported convergence criteria for MMA.")
    end
    @pack! solution.convstate = x_converged,
                                fabs_converged,
                                frel_converged,
                                kkt_converged,
                                ipopt_converged,
                                infeas_converged,
                                Δx,
                                Δf,
                                relΔf,
                                kkt_residual,
                                ipopt_residual,
                                infeas,
                                f_increased,
                                converged
    return solution
end

function getresiduals(solution::Solution, ::AbstractModel, ::GenericCriteria)
    @unpack prevx, x, prevf, f, g = solution
    Δx = maximum(abs(x[j] - prevx[j]) for j in 1:length(x))
    Δf = abs(f - prevf)
    infeas = length(g) == 0 ? zero(eltype(g)) : max(0, maximum(g))
    return Δx, Δf, infeas
end
function getresiduals(solution::Solution, model::AbstractModel, ::KKTCriteria)
    @unpack ∇f, g, ∇g, λ, x = solution
    xmin, xmax = getmin(model), getmax(model)
    T = eltype(x)
    res = maximum(1:length(x)) do j
        @views temp = ∇f[j] + dot(∇g[:,j], λ)
        if xmin[j] >= x[j]
            return abs(min(0, temp))
        elseif x[j] >= xmax[j]
            return max(0, temp)
        else
            return abs(temp)
        end
    end
    if debugging[]
        @show λ, g
    end
    res = length(g) == 0 ? res : max(res, maximum(abs.(λ .* g)))
    if debugging[]
        @show maximum(abs, g)
        @show maximum(abs, λ)
        @show maximum(x)
    end
    infeas = length(g) == 0 ? zero(eltype(g)) : max(maximum(g), 0)
    return res, infeas
end
function getresiduals(solution::Solution, model::AbstractModel, ::IpoptCriteria)
    @unpack ∇f, g, ∇g, λ, x = solution
    xmin, xmax = getmin(model), getmax(model)
    T = eltype(x)
    n, m, s = length(x), length(λ), zero(T)
    res = maximum(1:n) do j
        @views temp = ∇f[j] + dot(∇g[:,j], λ)
        if xmin[j] >= x[j]
            dj = temp
            s += max(dj, 0)
            return abs(min(0, dj))
        elseif x[j] >= xmax[j]
            yj = -temp
            s += max(yj, 0)
            return abs(min(0, yj))
        else
            return abs(temp)
        end
    end
    sd = max(100, (sum(abs, λ) + s) / (n + m)) / 100
    res = length(g) == 0 ? res : max(res, maximum(abs.(λ .* g)))
    res = res / sd
    infeas = length(g) == 0 ? zero(eltype(g)) : max(maximum(g), 0)
    if debugging[]
        println("Agg infeas = ", infeas)
    end
    return res, infeas
end


text:
using Test
using OrdinaryDiffEq
using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems;
importodeproblems();
import DiffEqProblemLibrary.ODEProblemLibrary:
    prob_ode_linear, prob_ode_2Dlinear, prob_ode_lotkavoltera, prob_ode_fitzhughnagumo

@testset "Test the different diffusion models" begin
    prob = prob_ode_fitzhughnagumo
    true_sol = solve(prob, Tsit5(), abstol=1e-12, reltol=1e-12)

    @testset "Time-Varying Diffusion" begin
        sol = solve(
            prob,
            EK0(diffusionmodel=:dynamic, smooth=false),
            dense=false,
            adaptive=false,
            dt=1e-4,
        )
        @test sol.u[end] ≈ true_sol.(sol.t)[end]
    end

    @testset "Time-Fixed Diffusion" begin
        sol = solve(
            prob,
            EK0(diffusionmodel=:fixed, smooth=false),
            dense=false,
            adaptive=false,
            dt=1e-4,
        )
        @test sol.u[end] ≈ true_sol.(sol.t)[end]
    end

    @testset "Time-Varying Diagonal Diffusion" begin
        sol = solve(
            prob,
            EK0(diffusionmodel=:dynamicMV, smooth=false),
            dense=false,
            adaptive=false,
            dt=1e-4,
        )
        @test sol.u[end] ≈ true_sol.(sol.t)[end]
    end

    @testset "Time-Fixed Diagonal Diffusion" begin
        sol = solve(
            prob,
            EK0(diffusionmodel=:fixedMV, smooth=false),
            dense=false,
            adaptive=false,
            dt=1e-4,
        )
        @test sol.u[end] ≈ true_sol.(sol.t)[end]
    end

    # @testset "Time-Fixed Diffusion MAP" begin
    #     sol = solve(prob, EK0(diffusionmodel=:fixedMAP, smooth=false),
    #                 dense=false, adaptive=false, dt=1e-4)
    #     @test sol.u[end] ≈ true_sol.(sol.t)[end]
    # end

end


text:
import PyCall

# Change that to whatever packages you need.
const PACKAGES = ["pygraphviz", "PyPDF2", "wntr"]

try
    for package in PACKAGES
        PyCall.pyimport(package)
    end
catch
    @warn("Python `pip` will be used to install Python modules.")
    try
        # Import pip.
        PyCall.pyimport("pip")
    catch
        # If pip is not found, install it.
        @warn("Pip not found on your sytstem. Downloading it.")
        get_pip = joinpath(dirname(@__FILE__), "get-pip.py")
        download("https://bootstrap.pypa.io/get-pip.py", get_pip)
        run(`$(PyCall.python) $get_pip --user`)
    end

    # don't upgrade... let pip work with whatever exists or install new
    # run(`$(PyCall.python) -m pip install --user --upgrade pip setuptools`)
    run(`$(PyCall.python) -m pip install --user $(PACKAGES)`)
end


text:
module Plotting

using AbstractPlotting
using AbstractPlotting.MakieLayout
using ColorSchemes
using LinearAlgebra
using StaticArrays

using ..Funs
using ..Manifolds
using ..SparseOps

################################################################################

export plot_manifold
function plot_manifold(mfd::Manifold{D,C,S}, filename=nothing) where {D,C,S}
    D::Int
    C::Int
    @assert 0 ≤ D ≤ C

    xmin = SVector{D}(minimum(x -> x[d], get_coords(mfd)) for d in 1:D)
    xmax = SVector{D}(maximum(x -> x[d], get_coords(mfd)) for d in 1:D)
    dx = xmax - xmin
    sz = norm(dx) / 100

    visible(xs) = true
    # visible(xs) = xs[1] ≥ 0.5 && xs[2] ≤ 0.2 && xs[3] ≥ 0.3

    if C == 2
        scene, layout = layoutscene(; resolution=(1024, 1024))
        laxis = layout[1, 1] = LAxis(scene; aspect=DataAspect())
        canvas = laxis
    elseif C == 3
        scene = Scene(; resolution=(1024, 1024))
        canvas = scene
    else
        error("C ∉ (2, 3)")
    end

    # Edges
    edges = SVector{C,S}[]
    for i in axes(get_simplices(mfd, 1), 2)
        sj = sparse_column_rows(get_simplices(mfd, 1), i)
        @assert length(sj) == 2
        xs1 = get_coords(mfd)[sj[1]]
        xs2 = get_coords(mfd)[sj[2]]
        if visible(xs1) && visible(xs2)
            push!(edges, xs1)
            push!(edges, xs2)
        end
    end
    linesegments!(canvas, edges; color=:green, linestyle=:solid, linewidth=3)

    # Vertices
    vertices = SVector{C,S}[]
    for i in axes(get_simplices(mfd, 0), 2)
        sj = sparse_column_rows(get_simplices(mfd, 0), i)
        @assert length(sj) == 1
        xs1 = get_coords(mfd)[sj[1]]
        if visible(xs1)
            push!(vertices, xs1)
        end
    end
    scatter!(canvas, vertices; markersize=sz, strokecolor=:red, strokewidth=5)

    # Dual edges
    dualedges = SVector{C,S}[]
    for i in axes(get_lookup(mfd, D, D - 1), 2)
        sj = sparse_column_rows(get_lookup(mfd, D, D - 1), i)
        for j in sj
            xs1 = get_dualcoords(mfd, D - 1)[i]
            xs2 = get_dualcoords(mfd, D)[j]
            if visible(xs1) && visible(xs2)
                push!(dualedges, xs1)
                push!(dualedges, xs2)
            end
        end
    end
    linesegments!(canvas, dualedges; color=:red, linestyle=:dash, linewidth=3)

    # Dual vertices
    dualvertices = SVector{C,S}[]
    for i in axes(get_dualcoords(mfd, D), 1)
        xs1 = get_dualcoords(mfd, D)[i]
        if visible(xs1)
            push!(dualvertices, xs1)
        end
    end
    scatter!(canvas, dualvertices; markersize=sz, strokecolor=:blue,
             strokewidth=5)

    if C == 2
        xmin -= dx / 10
        xmax += dx / 10
        limits!(canvas, xmin[1], xmax[1], xmin[2], xmax[2])
    end

    if filename ≢ nothing
        save(filename, scene)
    end
    return scene
end

################################################################################

export plot_function
export plot_function1d

function plot_function(fun::Fun{D,P,R,1,S,T},
                       filename=nothing) where {D,P,R,S,T}
    D::Int
    C = 1
    @assert 0 ≤ D ≤ C

    mfd = fun.manifold

    visible(xs) = true
    # visible(xs) = xs[1] ≥ 0.5 && xs[2] ≤ 0.2 && xs[3] ≥ 0.3

    scene, layout = layoutscene(; resolution=(1024, 1024))
    laxis = layout[1, 1] = LAxis(scene)
    canvas = laxis

    @assert P == Pr && R == 0

    points = [(get_coords(mfd)[i][1], fun.values[i])
              for i in axes(get_simplices(mfd, D), 1)]
    sort!(points)
    lines!(canvas, points; color=(:black, 0.6), linewidth=4)
    scatter!(canvas, points; markersize=4, strokecolor=:red, strokewidth=4)

    # xmin = SVector{D}(minimum(x -> x[d], get_coords(mfd)) for d in 1:D)
    # xmax = SVector{D}(maximum(x -> x[d], get_coords(mfd)) for d in 1:D)
    # dx = xmax - xmin
    # xmin -= dx / 10
    # xmax += dx / 10
    # limits!(canvas, xmin[1], xmax[1], xmin[2], xmax[2])

    if filename ≢ nothing
        save(filename, scene)
    end

    return scene
end

function plot_function1d(fun::Fun{D,P,R,2,S,T},
                         filename=nothing) where {D,P,R,S,T}
    D::Int
    C = 2
    @assert 0 ≤ R ≤ D ≤ C

    mfd = fun.manifold

    # visible(xs) = true
    visible(xs) = 0.5 - 0.001 ≤ xs[2] ≤ 0.5 + 0.001

    scene, layout = layoutscene(; resolution=(1024, 1024))
    laxis = layout[1, 1] = LAxis(scene)
    canvas = laxis

    @assert P == Pr

    if R == 0
        points = [(get_coords(mfd, R)[i][1],
                   fun.values[i] / get_volumes(mfd, R)[i])
                  for i in axes(get_simplices(mfd, R), 2)
                  if visible(get_coords(mfd, R)[i])]
        sort!(points)
        lines!(canvas, points; color=(:black, 0.6), linewidth=4)
        scatter!(canvas, points; markersize=4, strokecolor=:red, strokewidth=4)

        points = [(get_coords(mfd, R)[i][1],
                   fun.values[i] / get_volumes(mfd, R)[i])
                  for i in axes(get_simplices(mfd, R), 2)
                  if visible(get_coords(mfd, R)[i] - 1 / 8)]
        sort!(points)
        lines!(canvas, points; color=(:black, 0.6), linewidth=4)
        scatter!(canvas, points; markersize=4, strokecolor=:red, strokewidth=4)

        points = [(get_coords(mfd, R)[i][1],
                   fun.values[i] / get_volumes(mfd, R)[i])
                  for i in axes(get_simplices(mfd, R), 2)
                  if visible(get_coords(mfd, R)[i] - 1 / 16)]
        sort!(points)
        lines!(canvas, points; color=(:blue, 0.6), linewidth=4)
        scatter!(canvas, points; markersize=4, strokecolor=:purple,
                 strokewidth=4)
    elseif R == 1
        points = [(get_coords(mfd, R)[i][1],
                   fun.values[i] *
                   (get_coords(mfd, 0)[sparse_column_rows(get_simplices(mfd, R),
                                                          i)[2]] -
                    get_coords(mfd, 0)[sparse_column_rows(get_simplices(mfd, R),
                                                          i)[1]]) ⋅ (1, 0))
                  for i in axes(get_simplices(mfd, R), 2)
                  if visible(get_coords(mfd, R)[i])]
        sort!(points)
        lines!(canvas, points; color=(:black, 0.6), linewidth=4)
        scatter!(canvas, points; markersize=4, strokecolor=:red, strokewidth=4)

        points = [(get_coords(mfd, R)[i][1],
                   fun.values[i] *
                   (get_coords(mfd, 0)[sparse_column_rows(get_simplices(mfd, R),
                                                          i)[2]] -
                    get_coords(mfd, 0)[sparse_column_rows(get_simplices(mfd, R),
                                                          i)[1]]) ⋅ (1, 0))
                  for i in axes(get_simplices(mfd, R), 2)
                  if visible(get_coords(mfd, R)[i] - 1 / 16)]
        sort!(points)
        lines!(canvas, points; color=(:black, 0.6), linewidth=4)
        scatter!(canvas, points; markersize=4, strokecolor=:red, strokewidth=4)

        points = [(get_coords(mfd, R)[i][1],
                   fun.values[i] *
                   (get_coords(mfd, 0)[sparse_column_rows(get_simplices(mfd, R),
                                                          i)[2]] -
                    get_coords(mfd, 0)[sparse_column_rows(get_simplices(mfd, R),
                                                          i)[1]]) ⋅ (1, 0))
                  for i in axes(get_simplices(mfd, R), 2)
                  if visible(get_coords(mfd, R)[i] - 1 / 32)]
        sort!(points)
        lines!(canvas, points; color=(:blue, 0.6), linewidth=4)
        scatter!(canvas, points; markersize=4, strokecolor=:purple,
                 strokewidth=4)
    end

    # xmin = SVector{D}(minimum(x -> x[d], get_coords(mfd)) for d in 1:D)
    # xmax = SVector{D}(maximum(x -> x[d], get_coords(mfd)) for d in 1:D)
    # dx = xmax - xmin
    # xmin -= dx / 10
    # xmax += dx / 10
    # limits!(canvas, xmin[1], xmax[1], xmin[2], xmax[2])

    if filename ≢ nothing
        save(filename, scene)
    end

    return scene
end

function plot_function(fun::Fun{D,P,R,2,S,T},
                       filename=nothing) where {D,P,R,S,T}
    D::Int
    C = 2
    @assert 0 ≤ D ≤ C

    mfd = fun.manifold

    visible(xs) = true
    # visible(xs) = xs[1] ≥ 0.5 && xs[2] ≤ 0.2 && xs[3] ≥ 0.3

    scene, layout = layoutscene(; resolution=(1024, 1024))
    laxis = layout[1, 1] = LAxis(scene; aspect=DataAspect())
    canvas = laxis

    @assert P == Pr && R == 0

    color = fun.values.vec
    colormap = ColorSchemes.plasma.colors
    colorrange = (minimum(color), maximum(color))

    vertices = [get_coords(mfd)[i][d]
                for i in axes(get_simplices(mfd, D), 1), d in 1:D]
    connectivity = [Int(sparse_column_rows(get_simplices(mfd, D), j)[n])
                    for j in axes(get_simplices(mfd, D), 2), n in 1:(D + 1)]
    poly!(canvas, vertices, connectivity; color=color, colormap=colormap,
          colorrange=colorrange, strokecolor=(:black, 0.6), strokewidth=4)
    colorlegend!(colormap, colorrange)

    xmin = SVector{D}(minimum(x -> x[d], get_coords(mfd)) for d in 1:D)
    xmax = SVector{D}(maximum(x -> x[d], get_coords(mfd)) for d in 1:D)
    dx = xmax - xmin
    xmin -= dx / 10
    xmax += dx / 10
    limits!(canvas, xmin[1], xmax[1], xmin[2], xmax[2])

    if filename ≢ nothing
        save(filename, scene)
    end

    return scene
end

function plot_function(fun::Fun{D,P,R,3,S,T},
                       filename=nothing) where {D,P,R,S,T}
    D::Int
    C = 3
    @assert 0 ≤ D ≤ C

    mfd = fun.manifold

    visible(xs) = true
    # visible(xs) = xs[1] ≥ 0.5 && xs[2] ≤ 0.2 && xs[3] ≥ 0.3

    xmin = SVector{D}(minimum(x -> x[d], get_coords(mfd)) for d in 1:D)
    xmax = SVector{D}(maximum(x -> x[d], get_coords(mfd)) for d in 1:D)
    dx = xmax - xmin
    sz = norm(dx) / 100

    scene = Scene(; resolution=(1024, 1024))
    canvas = scene

    @assert P == Pr && R == 0

    color = fun.values.vec

    # Edges
    edges = SVector{C,S}[]
    for i in axes(get_simplices(mfd, 1), 2)
        sj = sparse_column_rows(get_simplices(mfd, 1), i)
        @assert length(sj) == 2
        xs1 = get_coords(mfd)[sj[1]]
        xs2 = get_coords(mfd)[sj[2]]
        if visible(xs1) && visible(xs2)
            push!(edges, xs1)
            push!(edges, xs2)
        end
    end
    linesegments!(canvas, edges; color=:green, linestyle=:solid, linewidth=3)

    # Vertices
    vertices = SVector{C,S}[]
    for i in axes(get_simplices(mfd, 0), 2)
        sj = sparse_column_rows(get_simplices(mfd, 0), i)
        @assert length(sj) == 1
        xs1 = get_coords(mfd)[sj[1]]
        if visible(xs1)
            push!(vertices, xs1)
        end
    end
    scatter!(canvas, vertices; markersize=sz, color=color, strokecolor=:red,
             strokewidth=10)

    if filename ≢ nothing
        save(filename, scene)
    end

    return scene
end

end


text:
# Autogenerated wrapper script for mdbtools_jll for i686-linux-gnu
export libmdb, libmdbsql, mbd_ver, mdb_array, mdb_count, mdb_export, mdb_header, mdb_hexdump, mdb_parsecsv, mdb_prop, mdb_schema, mdb_sql, mdb_tables

using Glib_jll
## Global variables
PATH = ""
LIBPATH = ""
LIBPATH_env = "LD_LIBRARY_PATH"
LIBPATH_default = ""

# Relative path to `libmdb`
const libmdb_splitpath = ["lib", "libmdb.so"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libmdb_path = ""

# libmdb-specific global declaration
# This will be filled out by __init__()
libmdb_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libmdb = "libmdb.so.2"

# Relative path to `libmdbsql`
const libmdbsql_splitpath = ["lib", "libmdbsql.so"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libmdbsql_path = ""

# libmdbsql-specific global declaration
# This will be filled out by __init__()
libmdbsql_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libmdbsql = "libmdbsql.so.2"

# Relative path to `mbd_ver`
const mbd_ver_splitpath = ["bin", "mdb-ver"]

# This will be filled out by __init__() for all products, as it must be done at runtime
mbd_ver_path = ""

# mbd_ver-specific global declaration
function mbd_ver(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(mbd_ver_path)
    end
end

# Relative path to `mdb_array`
const mdb_array_splitpath = ["bin", "mdb-array"]

# This will be filled out by __init__() for all products, as it must be done at runtime
mdb_array_path = ""

# mdb_array-specific global declaration
function mdb_array(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(mdb_array_path)
    end
end

# Relative path to `mdb_count`
const mdb_count_splitpath = ["bin", "mdb-count"]

# This will be filled out by __init__() for all products, as it must be done at runtime
mdb_count_path = ""

# mdb_count-specific global declaration
function mdb_count(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(mdb_count_path)
    end
end

# Relative path to `mdb_export`
const mdb_export_splitpath = ["bin", "mdb-export"]

# This will be filled out by __init__() for all products, as it must be done at runtime
mdb_export_path = ""

# mdb_export-specific global declaration
function mdb_export(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(mdb_export_path)
    end
end

# Relative path to `mdb_header`
const mdb_header_splitpath = ["bin", "mdb-header"]

# This will be filled out by __init__() for all products, as it must be done at runtime
mdb_header_path = ""

# mdb_header-specific global declaration
function mdb_header(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(mdb_header_path)
    end
end

# Relative path to `mdb_hexdump`
const mdb_hexdump_splitpath = ["bin", "mdb-hexdump"]

# This will be filled out by __init__() for all products, as it must be done at runtime
mdb_hexdump_path = ""

# mdb_hexdump-specific global declaration
function mdb_hexdump(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(mdb_hexdump_path)
    end
end

# Relative path to `mdb_parsecsv`
const mdb_parsecsv_splitpath = ["bin", "mdb-parsecsv"]

# This will be filled out by __init__() for all products, as it must be done at runtime
mdb_parsecsv_path = ""

# mdb_parsecsv-specific global declaration
function mdb_parsecsv(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(mdb_parsecsv_path)
    end
end

# Relative path to `mdb_prop`
const mdb_prop_splitpath = ["bin", "mdb-prop"]

# This will be filled out by __init__() for all products, as it must be done at runtime
mdb_prop_path = ""

# mdb_prop-specific global declaration
function mdb_prop(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(mdb_prop_path)
    end
end

# Relative path to `mdb_schema`
const mdb_schema_splitpath = ["bin", "mdb-schema"]

# This will be filled out by __init__() for all products, as it must be done at runtime
mdb_schema_path = ""

# mdb_schema-specific global declaration
function mdb_schema(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(mdb_schema_path)
    end
end

# Relative path to `mdb_sql`
const mdb_sql_splitpath = ["bin", "mdb-sql"]

# This will be filled out by __init__() for all products, as it must be done at runtime
mdb_sql_path = ""

# mdb_sql-specific global declaration
function mdb_sql(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(mdb_sql_path)
    end
end

# Relative path to `mdb_tables`
const mdb_tables_splitpath = ["bin", "mdb-tables"]

# This will be filled out by __init__() for all products, as it must be done at runtime
mdb_tables_path = ""

# mdb_tables-specific global declaration
function mdb_tables(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(mdb_tables_path)
    end
end

"""
Open all libraries
"""
function __init__()
    global artifact_dir = abspath(artifact"mdbtools")

    # Initialize PATH and LIBPATH environment variable listings
    global PATH_list, LIBPATH_list
    # From the list of our dependencies, generate a tuple of all the PATH and LIBPATH lists,
    # then append them to our own.
    foreach(p -> append!(PATH_list, p), (Glib_jll.PATH_list,))
    foreach(p -> append!(LIBPATH_list, p), (Glib_jll.LIBPATH_list,))

    global libmdb_path = normpath(joinpath(artifact_dir, libmdb_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libmdb_handle = dlopen(libmdb_path)
    push!(LIBPATH_list, dirname(libmdb_path))

    global libmdbsql_path = normpath(joinpath(artifact_dir, libmdbsql_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libmdbsql_handle = dlopen(libmdbsql_path)
    push!(LIBPATH_list, dirname(libmdbsql_path))

    global mbd_ver_path = normpath(joinpath(artifact_dir, mbd_ver_splitpath...))

    push!(PATH_list, dirname(mbd_ver_path))
    global mdb_array_path = normpath(joinpath(artifact_dir, mdb_array_splitpath...))

    push!(PATH_list, dirname(mdb_array_path))
    global mdb_count_path = normpath(joinpath(artifact_dir, mdb_count_splitpath...))

    push!(PATH_list, dirname(mdb_count_path))
    global mdb_export_path = normpath(joinpath(artifact_dir, mdb_export_splitpath...))

    push!(PATH_list, dirname(mdb_export_path))
    global mdb_header_path = normpath(joinpath(artifact_dir, mdb_header_splitpath...))

    push!(PATH_list, dirname(mdb_header_path))
    global mdb_hexdump_path = normpath(joinpath(artifact_dir, mdb_hexdump_splitpath...))

    push!(PATH_list, dirname(mdb_hexdump_path))
    global mdb_parsecsv_path = normpath(joinpath(artifact_dir, mdb_parsecsv_splitpath...))

    push!(PATH_list, dirname(mdb_parsecsv_path))
    global mdb_prop_path = normpath(joinpath(artifact_dir, mdb_prop_splitpath...))

    push!(PATH_list, dirname(mdb_prop_path))
    global mdb_schema_path = normpath(joinpath(artifact_dir, mdb_schema_splitpath...))

    push!(PATH_list, dirname(mdb_schema_path))
    global mdb_sql_path = normpath(joinpath(artifact_dir, mdb_sql_splitpath...))

    push!(PATH_list, dirname(mdb_sql_path))
    global mdb_tables_path = normpath(joinpath(artifact_dir, mdb_tables_splitpath...))

    push!(PATH_list, dirname(mdb_tables_path))
    # Filter out duplicate and empty entries in our PATH and LIBPATH entries
    filter!(!isempty, unique!(PATH_list))
    filter!(!isempty, unique!(LIBPATH_list))
    global PATH = join(PATH_list, ':')
    global LIBPATH = join(vcat(LIBPATH_list, [joinpath(Sys.BINDIR, Base.LIBDIR, "julia"), joinpath(Sys.BINDIR, Base.LIBDIR)]), ':')

end  # __init__()


text:
function timearray(response::Requests.Response)
    #This function transform the Response object on a TimeArray

    # Split the data on every "\n"
    data = split(response.data, "\n")

    # Extract the head and body of the data
    head = data[1]
    body = data[2:end]

    # Parse body
    body[end] == "" ? pop!(body) : nothing # remove trailing empty string if it's there
    body      = [split(line, ",") for line in body] # split on comma

    ######### Timestamp
    # take the first row (assuming it's date)
    # TODO: regex query needed to catch edge cases
    dates     = [line[1] for line in body]
    timestamp = Date{ISOCalendar}[date(d) for d in dates] # parse dates

    ######### Values
    vals = [line[2:end] for line in body] # get rows 2 to the end
    vals = ss2float(vals) # SubString{ASCIIString}["88.94","89.48","88.8","89.23","30617089.0"]

    ######### Column names
    names = split(head, ",")[2:end] # Won't need the Date name (fist column) for TimeArray
    names = ASCIIString[name for name in names]

    return TimeArray(timestamp, vals, names)
end


text:
export Weighted, WeightedMatrix, ClampedWeighted, UnClampedWeighted, WeightOpt
export array, weights, lastlength, aname, wname, unclamp,
    normalise, normalise!, unnormalise, unnormalise!,
	totweight, wcopy, wcopy!, wglue, flatten, flatcopy, flatcopy!, svecs, CatView
export AbsVec, AbsMat, AbsArray

using Parameters

@with_kw_noshow struct WeightOpt
    norm::Bool = true
    clamp::Bool = false
    lo::Float64 = 0.0
    hi::Float64 = 1.0
    aname::String = "θ"
    wname::String = "p(θ)"
    like::Bool = false
end

mutable struct Weighted{T1,T2}
    array::T1
    weights::T2
    opt::WeightOpt

    function Weighted(array::T1, weights::T2, opt::WeightOpt) where {T1<:AbstractArray,T2<:AbstractVector}
        size(array,ndims(array)) == length(weights) || @error "length of weights much match size of last dimension of array"
        new{T1, T2}(array, weights, opt)
    end
end

const WeightedArray = Weighted
const WeightedVector = Weighted{<:AbstractVector}
const WeightedMatrix = Weighted{<:AbstractMatrix}

const MaybeWeightedArray = Union{AbstractArray, WeightedArray}
const MaybeWeightedVector = Union{AbstractVector, WeightedVector}
const MaybeWeightedMatrix = Union{AbstractMatrix, WeightedMatrix}

const AbsVec = AbstractVector
const AbsMat = AbstractMatrix
const AbsArray = AbstractArray

##### Constructors

using LinearAlgebra

"""
    x = Weighted(array, weights)
Created a Weighted Array, with given weights associated to the last dimension of array, e.g. to columns of a matrix.
Keyword `norm=true` specifies that `weights(x)` sum to 1 always; this is enforced on construction.
Calling `normalise!(x)` will re-enforce this, after scaling `2x` or `hcat()`-ing `x + y` or mutating.

    x = Weighted(array, weights, lo, hi)
Here every `θ ∈ array` is constrained `lo ≦ θ ≦ hi`, the same for all dimensions. (Weights are similarly `0 ≦ λ < ∞`.)
Calling `clamp!(x)` will re-enforce these box constraints, for instance after mutating `x[1,2] *= 3` or on `y = x .+ 0.4`.
"""
function Weighted(array::AbsArray, vector::AbsVec = ones(lastlength(array)); norm=true)
    weights = clamp.(vector,0,Inf)
    norm && rmul!(weights, 1/sum(weights))
    Weighted(array, weights, WeightOpt(norm=norm))
end

Weighted(array::AbsArray, lo::Real, hi::Real; norm=true) = Weighted(array, ones(lastlength(array)), lo, hi; norm=norm)

function Weighted(array::AbsArray, vector::AbsVec, lo::Real, hi::Real; norm=true)
    weights = clamp.(vector,0,Inf)
    norm && rmul!(weights, 1/sum(weights))
    Weighted(clamp.(array,lo,hi), weights, WeightOpt(norm=norm, clamp=true, lo=lo, hi=hi))
end

Weighted() = Weighted(Float64[])

WeightedMatrix(x::AbsMat, rest... ; kw...) = Weighted(x, rest...; kw...)
WeightedMatrix(x::AbsVec, rest... ; kw...) = Weighted(reshape(x,:,1), rest...; kw...)

"""
    hcat(x::WeightedMatrix, y, z...)
    x + y
These will `hcat` the arrays, and `vcat` the weights.
"""
function Base.hcat(x::Weighted, zz::Vararg{Union{Weighted, AbstractArray}})
    if length(x.array)==0
        if length(zz)==1 return zz[1]
        else return hcat(zz...)
        end
    end
    arr = hcat(x.array, [array(z) for z in zz]...)
    wei = vcat(x.weights, [sureweights(z) for z in zz]...)
    Weighted(arr, wei, x.opt)
end

sureweights(x::Weighted) = x.weights
sureweights(x) = begin n = lastlength(x); ones(n) ./ n end

Base.:+(x::Weighted, y::Union{AbsArray, Weighted}) = hcat(x,y)

# """
#     λ*Π
# For `Π::Weighted` this will scale `Π.weights` by `λ`.
# """
# Base.:*(λ::Number, x::Weighted) = Weighted(x.array, λ .* x.weights, x.opt)
# Base.:*(x::Weighted, λ::Number) = Weighted(x.array, λ .* x.weights, x.opt)
# """
#     λ*Π
# For `Π::Weighted` this will scale `Π.array` by `λ`. NB this is the OPPOSITE of what it did in 0.6, sorry.
# """
# Base.:*(x::Weighted, λ::Number) = x * λ
# Base.:*(x::Weighted, λ::Number) = rmul!(copy(x),λ)

using Setfield

wname(o::WeightOpt, s::Union{String,Symbol}) = set(o, @lens(_.wname), string(s))
aname(o::WeightOpt, s::Union{String,Symbol}) = set(o, @lens(_.aname), string(s))
addname(o::WeightOpt, s::Union{String,Symbol})  = set(o, @lens(_.aname), o.aname * string(s))
addlname(o::WeightOpt, s::Union{String,Symbol}) = set(o, @lens(_.aname), string(s) * o.aname)

function LinearAlgebra.rmul!(x::Weighted, λ::Number)
    rmul!(x.array, λ)
    if x.opt.clamp
        o = x.opt
        o = set(o, @lens(_.lo), λ * o.lo)
        o = set(o, @lens(_.hi), λ * o.hi)
        x.opt = o
    end
    x
end

##### Copying

Base.copy(x::Weighted) = Weighted(Base.copy(x.array), Base.copy(x.weights), x.opt)
Base.copyto!(x::Weighted, y::Weighted) = begin copyto!(x.array, y.array); copyto!(x.weights, y.weights); x end

wcopy!(x::Weighted, vec::AbsVec) = (x.weights = clamp.(vec,0,Inf); x)
wcopy(x::Weighted, vec::AbsVec) = Weighted(x.array, clamp.(vec,0,Inf), x.opt) ## only weights are copied!

##### Extractors

array(x::Weighted) = x.array
array(x) = x

# weights(x::NormWeighted) = x.weights .* (1/sum(x.weights))
# weights(x::UnNormWeighted) = x.weights
weights(x::Weighted) = x.opt.norm ? x.weights .* (1/sum(x.weights)) : x.weights ## sum(positive, x.weights) gives exotic problems

positive(x) = max(x,zero(x))

totweight(x::Weighted) = sum(positive, x.weights)
maxweight(x::Weighted) = maximum(x.weights)

Base.getindex(x::Weighted, rr::AbsVec{Int}, c::Colon) = Weighted(x.array[rr,c], x.weights,  x.opt) ## get some dimensions
Base.getindex(x::Weighted, r::Colon, cc::AbsVec{Int}) = Weighted(x.array[r,cc], x.weights[cc], x.opt) ## get some points

using Lazy: @forward

@forward Weighted.array  Base.size, Base.ndims, Base.eltype, Base.getindex, Base.setindex!,
    Base.iterate, Base.lastindex,
    Base.minimum, Base.maximum, Base.extrema

lastlength(x) = size(x, ndims(x))

using StaticArrays

svecs(x::Weighted{<:Matrix{T}}) where T = reinterpret(SVector{size(x,1),T}, vec(x.array)) ## not type-stable
svecs(x::Weighted{<:Matrix{T}}, vald::Val{D}) where {T,D} = (@assert D==size(x,1); reinterpret(SVector{D,T}, vec(x.array)))

##### Standardisers

clampdoc = """
    clamp!(x::Weighted)
    clamp(x::Weighted)
Always clamps weights to be positive, and if flag `clamp=true` is set in `x.opt`,
then clamps `x.array` using `lo,hi` from `x.opt`.

    clamp(x, lo, hi)
    clamp!(x, lo, hi)
These alter the flag `x.opt.clamp` & then proceed.
"""
@doc clampdoc
function Base.clamp!(x::Weighted)
    clamp!(x.weights, 0, Inf)
    x.opt.clamp && clamp!(x.array, x.opt.lo, x.opt.hi)
    x
end
@doc clampdoc
function Base.clamp(x::Weighted)
    weights = clamp.(x.weights, 0, Inf)
    array = x.opt.clamp ? clamp.(x.array, x.opt.lo, x.opt.hi) : x.array
    Weighted(array, weights, x.opt)
end

Base.clamp!(x::Weighted, lo::Real, hi::Real) = begin x.opt = clamp(x.opt, lo, hi); clamp!(x) end
Base.clamp(x::Weighted, lo::Real, hi::Real) = clamp!(copy(x), lo, hi)

function Base.clamp(o::WeightOpt, lo=0.0, hi=1.0)
    o = set(o, @lens(_.clamp), true)
    o = set(o, @lens(_.lo), lo)
    o = set(o, @lens(_.hi), hi)
    o
end
unclamp(o::WeightOpt) = set(o, @lens(_.clamp), false)

"""
    normalise(x) ## with an s, NB!
    normalise!(x)
Ensures weights are positive and sum to 1.
On `x::Weighted`... the mutating form checks whether `norm=true` in `x.opt`;
the copying (!) form sets this flag first.

	unnormalise(x)
	unnormalise!(x)
Just alters the flag `x.opt.norm`.
"""
normalise(x::AbsVec) = begin itot = 1/sum(positive,x); itot .* clamp.(x, 0,Inf) end
@doc @doc(normalise)
normalise!(x::AbsVec) = begin clamp!(x,0,Inf); rmul!(x, 1/sum(x)) end
# normalise!(x::NormWeighted) = begin normalise!(x.weights); x end
# normalise!(x::UnNormWeighted) = x
normalise!(x::Weighted) = begin x.opt.norm && normalise!(x.weights); x end

normalise(o::WeightOpt) = set(o, @lens(_.norm), true)
unnormalise(o::WeightOpt) = set(o, @lens(_.norm), false)

normalise(x::Weighted) = Weighted(copy(x.array), normalise(x.weights), normalise(x.opt))

unnormalise(x::Weighted) =  Weighted(copy(x.array), copy(x.weights), unnormalise(x.opt))
unnormalise!(x::Weighted) =  begin x.opt = unnormalise(x.opt); x end

##### Flatten

"""
    flatten(x::Weighted)
Makes a vector `[ reshaped(x.array,:) ; x.weights ]`.

    flatcopy!(x, v::Vector)
Copies numbers from `v = flatten(x)` back into `x`.

    flatcopy(x, v)
Uses only `size(x)` to know what shape `::Weighted` to make out of `v` (and `x.opt` for details).
"""
flatten(x::Weighted) = vcat(reshape(x.array,:), reshape(x.weights,:)) ## does copy

@doc @doc(flatten)
function flatcopy!(x::Weighted, flat::Vector{<:Number})
    x.array[:] .= flat[1:length(x.array)]
    x.weights[:] .= flat[end-length(x.weights)+1:end]
    x
end

@doc @doc(flatten)
function flatcopy(x::Weighted, flat::AbsVec) ## AbstractVector allows ReverseDiff.TrackedArray
    array = reshape(flat[1:length(x.array)], size(x))
    weights = flat[end-length(x.weights)+1:end]
    Weighted(array, weights, x.opt)
end

using CatViews
"""
    CatView(x::Weighted)
Gives a vector-like view of both `x.array` and `x.weights`, the last `d` components are weights.
"""
CatViews.CatView(x::Weighted) = CatView(x.array, x.weights)

##### Pretty

function Base.show(io::IO, x::Weighted) ## compact, re-digestable
    print(io, "Weighted( ")
    show(io, array(x))

    print(io, ", ")
    show(io, weights(x))

    x.opt.clamp && print(io, ", ", x.opt.lo, ", ", x.opt.hi )
    x.opt.norm  || print(io, "; norm=false")
    print(io, " )")
end

function Base.show(io::IO, o::WeightOpt) ## compact, re-digestable
	print(io, "WeightOpt(")
	o.norm || print(io, "norm=false, ")
	o.clamp && print(io, "clamp=true, lo=",o.lo,", hi=",o.hi,", ")
    o.aname !== "θ" && print(io, "aname=",o.aname,", ")
    o.wname != "p(θ)" && print(io, "wname=",o.wname,", ")
    o.like && print(io, "like=true")
	print(io, ")")
end

function Base.show(io::IO, m::MIME"text/plain", x::Weighted) ## full
    ioc = IOContext(io, :compact => true, :limit => true)

    print(io, "Weighted ", summary(x.array))
    if x.opt.like
        println(io, ", likelihood ", x.opt.aname,":")
    elseif x.opt.clamp
        println(io, ", clamped ", round(x.opt.lo, digits=3), " ≦ ", x.opt.aname, " ≦ ", round(x.opt.hi, digits=3), ":")
    else
        println(io, ", of unclamped ", x.opt.aname,":")
    end
    if ndims(x.array)==1 && eltype(x.array) <: Number
        Base.print_array(ioc, x.array')  ## print vector as a row
    else
        Base.print_array(ioc, x.array)
    end

    length(x.weights)==0 && return nothing ## don't remember why
    if x.opt.norm
        println(io, "\nwith normalised weights ",x.opt.wname,", ", summary(x.weights), ":")
    else
        println(io, "\nwith weights ",x.opt.wname,", ", summary(x.weights), ":")
    end
    Base.print_array(ioc, weights(x)')

    if x.opt.norm && !isapprox(totweight(x), 1)
        print(io, "\nlazily normalised by weights(), raw total = ", round(totweight(x); digits=3), ".")
    end
end

function Base.summary(io::IO, x::Weighted)
    print(io, "Weighted ", summary(x.array))
    x.opt.clamp && print(io, " clamped to ", round(x.opt.lo, digits=3), " ... ", round(x.opt.hi, digits=3), ",")
    if x.opt.norm
        print(io, " with normalised weights ", summary(x.weights))
    else
        print(io, " with weights ", summary(x.weights))
    end
end